Figures
Abstract
Cognitive disorders, including Down syndrome (DS), present significant morphological alterations in neuron architectural complexity. However, the relationship between neuromorphological alterations and impaired brain function is not fully understood. To address this gap, we propose a novel computational model that accounts for the observed cell deformations in DS. The model consists of a cross-sectional layer of the mouse motor cortex, composed of 3000 neurons. The network connectivity is obtained by accounting explicitly for two single-neuron morphological parameters: the mean dendritic tree radius and the spine density in excitatory pyramidal cells. We obtained these values by fitting reconstructed neuron data corresponding to three mouse models: wild-type (WT), transgenic (TgDyrk1A), and trisomic (Ts65Dn). Our findings reveal a dynamic interplay between pyramidal and fast-spiking interneurons leading to the emergence of gamma activity (∼40 Hz). In the DS models this gamma activity is diminished, corroborating experimental observations and validating our computational methodology. We further explore the impact of disrupted excitation-inhibition balance by mimicking the reduction recurrent inhibition present in DS. In this case, gamma power exhibits variable responses as a function of the external input to the network. Finally, we perform a numerical exploration of the morphological parameter space, unveiling the direct influence of each structural parameter on gamma frequency and power. Our research demonstrates a clear link between changes in morphology and the disruption of gamma oscillations in DS. This work underscores the potential of computational modeling to elucidate the relationship between neuron architecture and brain function, and ultimately improve our understanding of cognitive disorders.
Author summary
The structural integrity of individual brain neurons and the intricate networks they form are fundamental to all brain functions, with structural anomalies directly linked to neurological disorders. Deciphering these links is a leading question in developmental disorders such as Down syndrome. Our work sheds light on the pivotal role that the structure of complex neural systems plays in shaping emergent network activity. In particular, our anatomically informed network modeling enables determining the extent to which specific deficits in neuronal architecture and connectivity perturb oscillatory patterns of activity.
Citation: Clusella P, Manubens-Gil L, Garcia-Ojalvo J, Dierssen M (2024) Modeling the impact of neuromorphological alterations in Down syndrome on fast neural oscillations. PLoS Comput Biol 20(7): e1012259. https://doi.org/10.1371/journal.pcbi.1012259
Editor: Daniel Bush, University College London, UNITED KINGDOM OF GREAT BRITAIN AND NORTHERN IRELAND
Received: January 13, 2024; Accepted: June 18, 2024; Published: July 5, 2024
Copyright: © 2024 Clusella et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All code and data in support of this publication are publicly available at https://github.com/pclus/neuromorphology.
Funding: PC and JGO have received funding from the Future and Emerging Technologies Programme (FET) of the European Union’s Horizon 2020 research and innovation programme (project NEUROTWIN, grant agreement No 101017716). JGO was also financially supported by the Spanish Ministry of Science and Innovation (project reference PID2021-127311NB-I00), the FEDER funds of the Spanish State Research Agency (project reference MICIN/AEI/10.13039/501100011033/FEDER), the Maria de Maeztu Programme for Units of Excellence in R\&D of the Spanish State Research Agency (project reference CEX2018-000792-M), and by the ICREA Academia program. LMG was supported by the National Natural Science Foundation of China (RFIS-I grant, reference 32350410413). MD acknowledges financial support from the Spanish Ministry of Science and Innovation (projects CPP2022-009659 and RTC2019-007329-1), the Spanish State Research Agency (project INTO-DS, reference PID2022-141900OB-I00), the Marato-TV3 Foundation (202212-30-31-32), and the European Comission Horizon 2020 programme (projects H2020-899986 and GO-DS21-848077). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Down Syndrome (DS), caused by the trisomy of chromosome 21, is associated with a wide spectrum of cognitive deficits [1], making it the most prevalent form of intellectual disability. Notably, abnormalities in the nervous system of individuals with DS manifest already at the single-neuron level. Mouse models of DS exhibit a significant reduction of dendritic tree branching and spine density when compared to control groups [2–4], features that have also been found in human postmortem tissue [5]. These morphological alterations are believed to play a significant role in the disruption of neural circuitry, ultimately contributing to the cognitive impairments associated with DS. Nevertheless, the precise mechanisms underlying the relationship between microscopic morphological alterations and mesoscopic brain dysfunction remain unknown.
Electrophysiological studies also show abnormal neural synchronicity in DS [6]. In particular, gamma rhythms (∼ 40 Hz) appear to be significantly reduced in both awake and anesthetized DS mouse models [7]. Alterations of these fast neural rhythms are not exclusive to DS, and have also been observed in other neuropathologies, such as Alzheimer’s disease [8, 9]. Gamma oscillations emerge from the collective activity of neural networks in both the hippocampus and cortex, and have been consistently associated with various cognitive processes, such as decision-making and memory tasks, across different species [10–17].
Thus we hypothesized that there is a direct link between microscopic circuitry abnormalities and functional deficits. However, this relation cannot be tested experimentally due to the presence of confounding factors in animal models. For instance, Ruiz-Mejias et al. [7] also reported a significant reduction of inhibitory connections targeting parvalbumin-positive interneurons in DS. This weakened recurrent inhibition was hinted as a potential cause for the reduction of gamma oscillations using a computational model, but the role of neuromorphological alterations was not explored in that study.
Here we propose a data-driven computational model of a simplified local neural network that incorporates some of the observed neuromorphological changes present in DS mouse models. We selected two different mouse models. The first model is trisomic for about two-thirds of the genes orthologous to human chromosome 21, (Ts(17(16))65Dn) [18], and is a well-characterized model for studying DS. We refer to this genotype as Ts65Dn for brevity. The second model (TgDyrk1A) [19], overexpresses only the dual-specificity tyrosine phosphorylation-regulated kinase 1a (Dyrk1a), a gene whose overexpression recapitulates the main neuronal architecture defects and cognitive impairments of the trisomy [20]. By integrating empirical data on dendritic complexity and spine density obtained from wild-type (WT), transgenic (TgDyrk1A), and trisomic mice (Ts65Dn), we construct a simplified cortical-layer structural model representing the synaptic connectivity of a neural network composed of point neurons, including pyramidal and fast-spiking interneurons. Simulations of these neural networks using Izhikevich dynamics [21] provide the functional differences between the three genotypes. Specifically, the model reproduces the deficit in gamma oscillations observed in DS animal models and allows us to test the role of reduced recurrent inhibition observed in [7]. Moreover, the scalable nature of the modeled morphologies enables us to explore the morphological space. This exploration includes values that do not correspond to specific animal models, allowing us to assess the impact of fabricated topologies on gamma rhythm generation.
Altogether, our study offers new insights into the complex interplay between neural morphology and network-level dynamics, thus contributing to our understanding of neurodevelopmental disorders.
Results
Our model consists of an in silico representation of a cross-section of layers II/III of the mouse motor cortex. In this neural network, we considered a total of N = 3037 neurons randomly distributed in a 2-dimensional square, with each side measuring 1500 μm. Details on the model construction, relevant parameters, and their grounding to the literature are outlined in the Methods section, but we briefly summarize the main aspects here.
The main challenge in our computational approach is to represent neural connectivity in a way that integrates single-cell morphology. Following the ideas of previous modeling studies [22, 23], we approach this problem by assuming simplified neuronal shapes paired with synaptic contact probability clouds based on experimental data. Each synthetic neuron is composed of a soma, a dendritic tree with variable size, and an axon. Axons are generated following a biased random walk starting from each neuron’s soma (see Fig 1). Whenever there is an intersection between a dendritic tree and an axon, a synapse is established according to a synaptic contact probability (SCP). The SCP function quantifies the likelihood of encountering a branch with a spine at a certain radial distance from the soma, denoted as r. Consequently, the SCP is influenced by the unique morphological attributes of each genotype (WT, Ts65Dn, and TgDyrk1A).
Small and large black circles represent the neurons’ soma and dendritic tree, respectively. The color gradient of the dendritic tree corresponds to the synaptic contact probability of each of the two neurons. Blue curves depict the axons, grown according to a biased random walk. Since the axon of neuron i overlaps with the dendritic tree of neuron j, a synapse might be established. The strength of the synapse depends on the length of the overlap and the value of the SCP along the coincident sites.
Synaptic contact probability from morphological data
To model the impact of neuromorphological alterations on the neural network topology in healthy and DS conditions we analyzed 18 individual pyramidal neurons from WT, Ts65Dn, and TgDyrk1A animal models (6 neurons for each genotype, see Methods) to calculate SCP based on single-cell morphology data.
Reconstructions of the analyzed neurons are displayed in Fig 2(a)–2(c). For each cell, we obtained the dendritic branching using a Sholl analysis [24]. The Sholl intersection profile is obtained by counting the number of dendritic branches at a given distance from the soma and is a key measure of dendritic complexity. Fig 2(d) shows the average number of intersections as a function of the distance from the soma r for each genotype. Neurons corresponding to the pathological conditions display significantly less branch density and shorter dendritic trees than the control condition (WT). Nonetheless, the shape of the branch density distribution, e.g. how branches are distributed or clustered in a particular area, exhibits a consistent similarity among the three cases, with variations primarily attributable to scaling factors. Using the WT case as the reference, we performed a non-linear fitting of the averaged Sholl intersection profile using the function
(1)
The choice of this function allows for an appropriate fit of the data with only two free parameters. The black curve in Fig 2(d) depicts the outcome of the fitting, with the resulting function parameters detailed in Table 1.
(a-c) Reconstructed neurons for each of the three animal models. (d) Sholl intersectional profile (branch pattern complexity) for the three different genotypes (WT, Ts65Dn, and TgDyrk1A). Each dataset corresponds to the average reconstruction of 6 different neurons. Black curve corresponds to fitting Eq (1) to the WT data (see Table 1). (e) Spine density for each different genotype as published in previous literature. The spine numbers are per 10 μm, thus the distributions are divided by 10 in the fitting procedure. Black curve corresponds to fitting Eq (2) to the WT data (see Table 1). (f) Resulting synaptic contact probability function. Circles, squares, and triangles obtained from the data presented in panels (d) and (e), continuous curves obtained from the nonlinear fitting of Eq (4) using the rescaling parameters α and β (see Table 2).
Next, we consider the dependency of spine density on the distance from the soma. Here we use data previously published in [2, 3] (see Fig 2(e)). Again, the maximal spine density remains comparable across trisomic, transgenic, and WT genotypes, but is influenced by the shorter dendritic trees in pathological conditions. We fit the WT spine density distribution using a 6th degree polynomial
(2) Table 1 contains the parameter values obtained from the fitting, and the black curve in Fig 2(e) shows the resulting function.
The product of the BD and SD functions provides the average spine density at a certain distance from the soma. These functions assume straight dendrites, whereas these are actually irregularly shaped in nature. Indeed, while the largest distance between a spine and the soma in the WT case is 218 μm, using a convex polygon fitting of reconstructed WT neurons, we found an average mean dendritic tree radius of (see Methods). To account for the actual shape of the dendritic tree, we rescale the radius in the BD and SD functions by a factor γ = 218/156.30. Finally, we divide by 2πr to account for the circular shape of the dendritic tree.
Altogether, for a typical WT neuron, the probability of finding a spine at a distance r from the soma is given by
(3)
This function is depicted in Fig 2(f) (see red continuous curve), together with the morphological data (see red circles).
In order to obtain a SCP distribution for Ts65Dn and TgDyrk1A, we exploit the fact that their spine and branch densities in Fig 2(d) and 2(e) follow a similar shape to the WT case up to scaling factors. Therefore, we consider the following generalization of the SCP:
(4)
Here, the parameter α determines an overall scaling of the SCP in comparison to the WT, whereas
provides the ratio of the mean dendritic tree radius in comparison to WT. We fit the expression in Eq (4) to the data corresponding to Ts65dn and TgDyrk1A (see blue squares and green triangles in Fig 2(c)) with α and β as free parameters. Table 2 contains the resulting values of the neuromorphological parameters, and continuous curves in Fig 2(f) display the resulting shape of the SCP.
The good agreement between the SCP model (Eq (4)) and the morphological data for the three genotypes substantiates the characterization of single neurons’ morphology alterations in DS with only two parameters, α and β.
Morphology changes impact network connectivity
We generate network topologies corresponding to each of the studied genotypes with the synaptic contact algorithm described in Methods and illustrated in Fig 1. Topology generation parameters are identical for all genotypes (WT, Ts65Dn, and TgDyrk1A) except for α and β as described in the previous section (see Table 2).
In order to characterize the main circuitry changes between the three genotypes, we analyze and compare the topologies in terms of network density, degree centrality, and synaptic strength (see Methods for definitions). To account for the variability on different network generations, we create and analyze 10 different networks for each mouse model. Table 3 contains the average quantities of the computed measures for each case, and Fig 3 displays the corresponding degree and weight distributions.
(a,b) In-degree (panel (a)) and out-degree (panel (b)) distributions corresponding to WT (red), Ts65Dn (blue), and TgDyrk1A (green). (c) Distribution of synaptic strengths, i.e., number of connections between the same two neurons. In all panels, lines correspond to the average of 10 network generations, with shaded regions indicating standard deviation. Standard deviation in panel (c) is too small to be visibly appreciated.
Each value indicates the average measure over a population of 10 networks generated with the same genotype parameters α and β. Values in parenthesis correspond to the sample standard deviation. Precise definitions for each measure are provided in the Methods section. The source code to generate network topologies is openly available at github.com/pclus/neuromorphology.
First, in terms of network density, the three genotypes provide sparse topologies, with less than 7% of all possible links present. Moreover, pathological genotypes show a significant decrease of network density with respect to control. The in-degree and out-degree distribution of the networks (see Fig 3(a) and 3(b)) also reflect this reduction of connectivity in the Ts65dn and TgDyrk1A models, with the average degree being almost half of that of the WT in both cases. In spite of the reduced density of TgDyrk1A compared to Ts65Dn, the degree distributions do not show significant changes between the two pathological cases.
Since in our model a presynaptic neuron can establish several contacts with the same postsynaptic neuron, we also take into account how such synaptic strength varies across genotypes (see Fig 3(c)). In this case, the differences between DS models and WT are less prominent, with only TgDyrk1A showing a decrease of average synaptic strength with respect to WT.
Overall, this analysis reflects that morphology changes incur a direct impact on the connectivities generated by our model. These differences are mainly reflected by an important decrease of number of pre and post-synaptic connections established by each neuron in the DS models, with the strength of such connections displaying only mild changes between genotypes. This scenario is consistent across network generations, as indicated by the small standard deviations shown in parenthesis in Table 3 and in the shaded regions of Fig 3.
Simulations recapitulate disrupted gamma activity in pathological conditions
Next, our aim is to investigate how the topological changes induced by morphology alterations affect the network at a functional level. We simulate the dynamics of each neuron using the Izhikevich model with parameters set for regular spiking (pyramidal) and fast-spiking (interneurons) [21] (see Methods). Importantly, each neuron in the network receives an independent train of external excitatory inputs following a Poisson shot process with frequency λ. The code to simulate the network is openly available at github.com/pclus/neuromorphology.
Fig 4(a)–4(f) show raster plots and mean-firing rate time series for each genotype obtained from network simulations with λ = 9 kHz. Single unit spike trends are rather irregular due to the three different stochastic sources in the model: randomness of the network generation, external inputs λ, and finite-size effects. Nonetheless, noisy collective oscillations emerge due to the interplay between the fast-spiking interneurons, the excitatory neurons, and the external excitatory input. This onset of gamma rhythmic activity in a noisy environment corresponds to the paradigmatic pyramidal-interneuron network gamma (PING) mechanism [13, 25–27]. S1, S2 and S3 Movies show such dynamic activity at both single-cell and collective levels.
(a-c) Raster plots showing the spike times of excitatory (black) and inhibitory (red) neurons for WT (panel (a)), Ts65Dn (panel (b)), and TgDyrk1A (panel (c)). (d-f) Mean firing rate of excitatory (black) and inhibitory (red) neurons corresponding to the raster plots shown in panels (a-c). (g) Power spectrum of LFP signals for an input frequency of λ = 9 spikes/ms. Each curve corresponds to the average of 100 spectra corresponding to 10 independent realizations of the noise for 10 different topologies. Shaded regions indicate the standard deviation among the samples. Red, blue, and green correspond to the morphological parameters of the WT, Ts65Dn, and TgDyrk1A cases, respectively (see Table 2). Dashed lines correspond to simulations with recurrent inhibitory synapses reduced to 0.3 of the original value. (h,i) Location of the peak power in the averaged LFP power spectra (h) and corresponding power (i) obtained for different values of the external firing rate λ. Circles and continuous lines indicate results using the default network parameters. Triangles and dashed lines correspond to reduced recurrent inhibition, as in panel (g). Error bars indicate the standard deviation over the 100 simulations pool.
Visual inspection of individual simulations in Fig 4(a)–4(f) already indicate a lack of synchronicity in the Ts65Dn and TgDyrk1A models. In order to properly compare the dynamics of the three cases, we capture the collective activity in each simulation by computing the local field potential (LFP) as the average network firing rate (see Methods). Furthermore, to test the consistency of the results against statistical fluctuations of the topology generation and external input simulation, we use 10 networks for each of the three neuronal genotypes, and each of them is simulated independently 10 times, resulting in a total pool of 100 time series for each parameter set.
Continuous lines in Fig 4(g) show the average power spectra of the LFP corresponding to each animal model for an input rate λ = 9 kHz. In all three genotypes, the spectrum shows a clear peak around 40 Hz, with very small deviation across different noise realizations, confirming the robustness of the gamma activity in the model. Nonetheless, the TgDyrk1A and Ts65Dn models display a clear reduction of gamma power compared to the WT case. This decrease in power is paired with a slight decrease in the peak frequency. These results align with empirical observations of reduced gamma activity in prefrontal cortex of TgDyrk1A mice with compared to WT [7].
Further assessment of the dynamical differences between the three genotypes is provided by the analysis of the interspike interval (ISI) distribution. S1(a) and S1(b) Fig show the average interspike interval distribution of each simulation for excitatory and inhibitory neurons respectively. Remarkably, the ISI distributions of pyramidal neurons show a broad profile, with a coefficient of variation around 1.1. Moreover, the three genotypes display a peak around 100 ms, indicating thus a preferred spiking period four times larger than the period of the collective gamma oscillation. On the other hand, inhibitory neurons do show a prominence of spikes around 25 ms, revealing thus the significance of interneuron activity for the emerge of gamma oscillations. The WT ISI distribution of inhibitory neurons also presents a larger peak at shorter periods, corresponding to repetitive firing of interneurons when the excitatory feedback is strong enough.
Next, we investigate the effects of the external firing rate λ on the network dynamics. Fig 4(h) and 4(i) show the peak frequency and power, respectively, for the three genotypes upon varying λ (solid lines in the two panels). Additionally S3 and S4 Figs show raster plots and mean firing rate activity for individual simulations with λ = 6 and 12 kHz respectively. For all explored values, the three models produce robust oscillatory activity with a peak frequency generally within the gamma range (30–80 Hz). For low values of λ, the network frequency shows similar behavior for all three mouse models, displaying an increase with external input up to λ ≈ 12 kHz (see Fig 4(h)). However, for larger values of external input, the monotonic dependence of the frequency on λ breaks down for the DS models. Such decline in frequency observed in TgDyrk1A and Ts65Dn for large λ values is concomitant with an elevation in variability across network realizations, suggesting a lack of uniformity in generating gamma oscillations in the pathological models.
Regarding the power of the neural oscillations (Fig 4(i)), the WT model demonstrates a notably higher gamma amplitude as compared to the DS models for λ > 7 kHz. These differences in power between genotypes remain mostly unchanged upon increasing λ, although all three models show a reduction of oscillatory coherence as the external input increases. This is consistent with the increase of single neuron spike irregularity, as captured by the average CV displayed in S1(c) and S1(d) Fig. Overall, the scenario remains similar to that displayed in Fig 4(g), and reproduces experimental findings observed in electrophysiology studies of the TgDyrk1A model [7].
Reduced recurrent inhibition might increase or reduce oscillatory activity
Histological analysis of TgDyrk1A and WT mice cortex shows a significant reduction of the inhibitory synapses acting upon interneurons in the DS model [7]. We expect this to influence the network dynamics, since parvalbumin-positive interneurons are known to modulate gamma activity in the cortex [13, 14, 27]. Moreover, recurrent inhibition, one of the key factors involved in fast collective activity [28], was probably a leading cause of the gamma impairment in the DS model [7] given the reduction of GABAergic contacts among interneurons.
In this section, we test the effect of reduced recurrent inhibition in our model. While in vivo, this perturbation of the network balance only occurs for DS animals, we deliberately test the effects of disrupted inhibition in all three genotypes. This allows us to compare the effects of morphology alterations and reduced inhibition separately, both of which coexist in DS animal models.
Dashed lines in Fig 4(g) show the average power spectral density (PSD) obtained from simulations in networks with inhibitory-to-inhibitory synaptic strength reduced to 30% for λ = 9 kHz. While the three genotypes still exhibit a prominent peak around 40 Hz, there is a substantial decrease in activity across all frequency bands when contrasted with the unperturbed models. Significantly, the distinctions in power between WT and DS models vanish, resulting in nearly identical spectra for all three. Raster plots and mean firing rates of individual simulations in S2, S3 and S4 Figs indicate that this reduction of power corresponds to an increased inhibitory activity, which in turns lowers the activity of excitatory neurons.
Once more, we test the soundness of this scenario upon changing λ. Dashed lines and triangles in Fig 4(h) show the peak frequency for the models with weakened recurrent GABAergic contacts. For low values of the external input, the main frequency of oscillation remains close to the unperturbed models. Moreover, disrupted inhibition rescues the drop in gamma frequency of the DS models reported in the previous section. Indeed, with higher values of λ all three genotypes exhibit faster oscillatory activity compared to their respective unperturbed models, with no discernible distinctions between the three genotypes.
The stimulating effect of decreased recurrent inhibition on the oscillatory frequency within the DS models differs from its effects on peak power. Dashed lines and triangles in Fig 4(i) show two different scenarios depending on whether the external input is smaller or larger than λ ≈ 7 kHz. With lower external activity levels, disrupted recurrent inhibition enhances oscillatory power, with the WT model consistently exhibiting greater power than the DS models. Conversely, for higher external activity levels, the power of the three modified models rapidly declines below the levels of the default models, with all three genotypes reaching a plateau for λ > 10 kHz. In this range, no substantial differences in power are discernible among the three genotypes. Furthermore, within this range, the impact of disrupted excitation-inhibition balance on gamma power in the WT appears to be twofold compared to the drop of gamma between the unperturbed WT and DS models. Since frequencies above 40 Hz require λ > 9 kHZ (Fig 4(h)), these results suggest that reducing recurrent inhibition notably impedes gamma synchronicity, in agreement with [7].
The abrupt reduction in peak power in the disrupted networks corresponds to a transition from a highly synchronized state to a regime in which neurons fire irregularly, while still exhibiting some degree of collective synchronous behavior, mainly through the interneuron population. These two forms of gamma activity have been identified in previous computational studies and are usually referred to as strong and weak gamma, respectively [29, 30]. The sensitivity of the network rhythmicity on the external input λ in the perturbed models highlights the importance of recurrent inhibition to obtain robust gamma rhythms in cortical neural networks.
Parameter exploration
The unified SCP (Eq (4)) function derived from the morphological data for the three animal models enables the investigation of the network activity generated by hypothetical neuron morphologies. In this context, we explore the influence of the spine density (α) and dendritic tree size (β) on the power and frequency of the gamma rhythm.
Fig 5 shows the outcome of the LFP signal obtained from numerical simulations in network topologies generated with specific values of the scaling parameters for the mean dendritic tree β and synaptic contact probability α. The peak of the gamma activity is observed at higher frequencies for networks with low values of α and high values of β. Conversely, the power of such gamma activity becomes larger when both morphological parameters are high. This dual relationship highlights that there is no specific region within the morphological space where both frequency and power can be maximized concurrently. Instead, the emergence of these fast oscillations is driven by the topological features of the networks in a nonlinear manner.
Results from numerical simulations of networks for external input λ = 10 kHz. Each pixel of the heatmap corresponds to the average of results obtained with 10 different network topologies, each simulated with 10 different realizations of the noise. Parameters corresponding to the studied genotypes are marked with black circles for reference.
When using the exact values of the animal models in this exploration, it becomes evident that the most prominent difference between the DS and the WT is given by a substantial reduction of the parameter β. This implies that the size of the dendritic tree exerts a pivotal influence on the gamma abnormalities when compared with the reduction of the SCP parameter α. Consequently, our computational results indicate that the gamma impairment in DS models is primarily attributed to the loss of synaptic connections among neurons rather than a decline in the overall strength of these connections. Nonetheless, the parameter exploration suggests that changes on SCP stronger than those observed in Ts65Dn and TgDyrk1A could also significantly alter gamma oscillations (see also S5 Fig).
Discussion
In this study, we aimed to explore the interplay between neuromorphological alterations and network dynamics in DS, using a data-driven computational model. Our findings provide compelling evidence that the incorporation of empirical data on dendritic complexity and spine density into the model enables the faithful replication of the reduction in gamma oscillations documented in DS animal models [7]. These results strongly indicate that the neuromorphological changes observed in DS are pivotal contributors to the disruption of gamma activity, providing valuable insights into the underlying mechanisms of network dynamics associated with DS.
An important challenge of this work has been to integrate single-cell morphology data into the generation of the neural network architecture. Few studies include morphology in neural network topology, and the options vary based on the desired biological realism: from sophisticated cloning of reconstructed morphologies paired with touch-detection rules [31–33], to generating random networks based mostly on axo-dendritic overlap [22, 34]. A mid-way approach consists of defining a probability region based on different morphological parameters [23, 35, 36].
Here we took this intermediate route by incorporating a synaptic contact probability (SCP) function into the model developed in [22]. Through an analysis of single-cell morphology data, we devised a quantitative model for the SCP that integrates dendritic complexity and spine density of different mouse models. Remarkably, our findings demonstrate that the three genotypes examined in this study (WT, trisomic, and single-gene transgenic TgDyrk1A) can be accurately represented by a single SCP function, with scaling factors accounting for their distinct characteristics. This result allowed us to explore the broader morphological landscape beyond the confines of specific animal models.
The simulations conducted in the study demonstrate the influence of diminished recurrent inhibition on network dynamics. The reduction of inhibitory connections, particularly those targeting parvalbumin-positive interneurons, previously reported in DS animal models [7], has been proposed as a potential cause for the observed reduction in gamma oscillation. The model’s implementation of reduced recurrent inhibition successfully reproduces the effects on gamma activity, underscoring its pivotal role in modulating network dynamics. Notably, perturbing excitation-inhibition balance rescues the drop in frequency in DS models but concurrently diminishes oscillatory power. These results align with previous empirical and modeling studies identifying recurrent inhibition as a key ingredient for the emergence of gamma oscillations [14, 27, 28, 37–39].
Furthermore, the study explores the morphological parameter space and demonstrates that there is no single regime in which both frequency and power of gamma oscillations can be maximized. Moreover, the reduction in dendritic tree size (β) appears to be a major factor contributing to the gamma abnormalities in DS, while the reduction in synaptic contact probability (α) plays a secondary role.
In summary, the data-driven computational model presented in this study successfully integrates neuromorphological alterations observed in DS animal models and reproduces the reduction in gamma oscillations. The findings support the notion that microscopic circuitry abnormalities contribute to the disruption of network dynamics in DS. The model provides a controlled framework to explore the impact of morphology alterations and elucidate their role in network synchronicity.
Model limitations
Our data-driven computational approach relies on several modeling assumptions that should be taken into account when drawing biological conclusions from the results. One major simplification pertains to the morphological model for network topologies, which simplifies neuronal shapes and may not fully capture the diversity and intricacy of real morphologies. Furthermore, the axon growth model doesn’t consider the precise axon terminal locations.
Another aspect is the dynamical model employed for the simulations. Here we considered idealized point neurons, which are simplified models that disregard the influence of detailed structures, such as dendrites and axons, on the dynamics of the cell. While these simplifications facilitate easier analysis and interpretation, they may not fully capture the effects that morphology can have on the system dynamics. To this end, compartmental models could allow the incorporation of specific relations such as the effects of smaller dendritic trees on the synaptic delays.
Despite these limitations, our data-driven computational model is an important milestone in understanding the importance of neuromorphological alterations in neurodevelopmental diseases. Further work should allow us to improve and further validate our assumptions based on empirical data, as well as focus on other brain areas which are known to be also affected in DS, such as the hippocampus.
Methods
Data gathering
We used previously published single-neuron 2D tracings of cortical layer II/III pyramidal neurons. Those were traced as part of previous studies of our team, following experimental procedures detailed in [40–42]. Briefly, cells were injected with Lucifer Yellow, immuno-stained with a biotinylated secondary antibody and biotin–horseradish peroxidase complex. Tissue sections were imaged with brightfield microscopy and traced using camera lucida microscope attachment. Specifically, the details for the tracing of WT and Ts65Dn neurons can be found at [2], and those for TgDyrk1A neurons can be found at [3]. We did not find differences in neuronal morphology between the WT strains of trisomic and transgenic mice, and thus we consider WT parameters to be the same for the two DS mouse models, thus allowing comparisons among the three genotypes. All data used in this study is openly available at github.com/pclus/neuromorphology.
Topology generation
The generation of the network topology has been largely inspired by the model developed by Orlandi et. al. [22]. We included substantial modifications to their proposal to account for specific morphological variants of single neurons. A C code of the resulting algorithm is openly available at github.com/pclus/neuromorphology.
Using a full 3D representation of the cortical layer II/III with realistic neuronal densities is not possible due to computational limitations [22, 43]. Thus we considered a thin 2-dimensional cross-section of the cortical layer with a thickness of a single cell soma, 16 μm.
The neuronal density of the synthetic circuit is set to 1350 neurons/mm2, obtained by multiplying the neuronal density in the 3D layer by the width of the thin layer modeled. Based on this assumption, we placed randomly N = 3037 neurons in a 2-dimensional square with a side length of 1.5 mm. To mitigate the effects of imposing a too-small spatial domain, we provided the circuit with periodic boundary conditions. The spatial resolution of the in silico layer is 1 μm.
In the model, each neuron has three components: the soma, the dendritic tree, and the axon (see Fig 1 for a schematic representation). All neurons have identical soma, which are modeled as circles of radius RS = 16 μm. The center of each neuron’s soma is randomly distributed on the 2-dimensional layer, and overlapping somas are not allowed. The dendritic tree of each neuron is also modeled as a circle, but in this case, we considered variability among neurons. The radius of each neuron’s dendritic tree is a random number drawn from a Gaussian distribution with mean and standard deviation σ = 40. The mean dendritic tree radius
is one of the main control parameters of the study. As explained in the Results section, we use the WT mean dendritic tree radius
as a reference. We measured this quantity by calculating the 2D convex hull of the reconstructed WT dendritic trees. Specifically, we used the “boundary” function in MATLAB with shrink factor s = 0. To obtain the mean radius we assumed that the area obtained with the convex hull forms a circle for each tree.
The axon is modeled as a biased random walk starting from the center of the soma with a random direction. After it has grown 10 μm, it modifies direction by θ degrees, where θ is a random variable chosen from a Gaussian distribution with zero mean and standard deviation σ = π/30 radians. The total length of each axon is obtained from a Rayleigh distribution with mean . This mean axon length has been chosen to simulate only local horizontal connections (given by the local axonal tree in layers II/III), following experimental observations in the mouse M2 cortical layer II/III [43], and disregarding horizontal patchy connections, connections between cortical layers, and interhemispheric projections.
If the axon of neuron i overlaps the dendritic tree (but not the soma) of neuron j, then a synapse is established with a probability p that depends on the distance r from the overlap point to the neuron’s j soma. Such dependency is given by the function p = SCP(r;α, β) where α and β are parameters that depend on the morphological variables. In particular, α is the ratio between the synaptic contact probability of each mouse model and the wild-type (WT) value, and β is the ratio between the dendritic tree radius of each mouse model and WT, i.e., where
is the mean dendritic tree radius of a WT neuron. The derivation of the synaptic contact probability function SCP(r) and the values of α and β for the different animal models are detailed in the Results section.
Each grid square presenting an overlap between an axon and a dendritic tree can generate a synapse, thus each pair of neurons might have multiple, in some cases several, synapses. In other words, the resulting architecture of connections among neurons is a directed weighted network given by the weight matrix W = (wij). Autapses are not allowed.
Topology measures
We analyze the network topologies generated by our model using standard measures of complex network analysis:
- Network density: Proportion of links present in the network with respect to the total number of possible links:
(5) where A = (aij) is the network adjacency matrix, i.e.,
(6)
- In and out-degree: For each node i, the in-degree
is the number of pre-synaptic neurons with a contact on i. The out-degree
corresponds to the number of post-synaptic neurons connected by i. They are computed as
(7)
- Synaptic strength: Number of established connections between the same two neurons, wij.
Neuronal dynamics
Several models for spiking neuron dynamics exist in the literature. Here we use the model proposed in [21], due to its apt trade-off between dynamical richness and computational efficiency. The C code used to simulate the network dynamics is openly available at github.com/pclus/neuromorphology. We consider a network composed of pyramidal and fast-spiking interneurons only. Since the topology generation algorithm does not differentiate between the two types, each neuron is set as either excitatory or inhibitiory randomly at the beginning of each simulation with a 80%-20% proportion [26].
The dynamics of the i-th neuron in the network is ruled by the following ordinary differential equations,
(8)
(9)
where vi is the membrane voltage potential, and ui is a recovery variable accounting for the dynamics of the different ion channels. When the voltage of neuron i reaches a threshold of v(thr) = 30, the neuron emits a spike and the voltage is reset to vi ← c, whereas ui ← ui + d. Upon tuning the system parameters a, b, and d, one can obtain different dynamical behaviors for each neuron. Here we focus on Regular Spiking (RS) for pyramidal neurons and Fast Spiking (FS) for interneurons [21]. Table 4 indicates the numerical value for the parameters corresponding to each neuron type.
The voltage of each neuron is influenced by the synaptic inputs coming from the excitatory neurons of the network, , inhibitory neurons
, and glutamatergic inputs coming from outside the network, Iext. In all cases, these inputs have the form
(10)
with
being a reversal potential and
is a time-dependent conductance, and τ0 = 1 ms is a fixed synaptic delay. Upon receiving a spike, the neurotransmitter-activated ion channels open and close following an exponential decay. For the recurrent AMPA and GABA connectivities this reads
(11)
where W = (wij) is the weight matrix of the network,
is the time at which neuron m emitted its k-th spike, H is the Heaviside step function, τsyn is the decay time, gsyn is the strength of the synapse. The incoming signals from outside of the modeled layer consist of excitatory exponential pulses only
These inputs account for post-synaptic potentials from other cortical layers, as well as other brain regions projecting to layer 2/3. The spiking times t(k) are drawn from a Poissonian shot process with frequency λ. We use λ as the main control parameter to test the oscillatory response of the network. Values for gsyn, τsyn, and
for the three synaptic types are given in Table 5.
In order to capture the network activity, we model the local field potential (LFP) as the network average firing rate (with time bins of 0.1 ms):
(12)
We opted to use the firing rate as a proxy to capture the actual network activity in our model. Other options, which are based on the AMPA and GABA currents of the network, could be susceptible to parameter changes such as when we model disrupted networks by reducing gGABA targeting inhibitory neurons. Power Spectral Density of LFP time series have been computed using the GNU Scientific Library [44] by first applying a fast Fourier transform algorithm and then reducing the noise in the spectra through a Gaussian filter.
The firing activity of individual neurons and their regularity can be assessed by the interspike interval (ISI) distribution and the corresponding coefficient of variation (CV). For each neuron i we compute the time between consecutive spikes as . For periodically firing neurons, the ISI distribution approaches a delta function located at the period. On the other hand, a broad ISI distribution reflects more complex, irregular, spiking patterns. To characterize the regularity of the firing of each population type, we compute the ISI coefficient of variation as CV = σ/μ, where μ and σ correspond to the mean and standard deviation of
respectively. We perform this analysis for pyramidal and inhibitory neurons separately.
Supporting information
S1 Fig. Interspike interval and firing regularity for different network morphologies.
(a,b) Interspike interval (ISI) distribution corresponding to pyramidal neurons (a) and inhibitiory neurons (b) for λ = 9 kHz. Each curve corresponds to the average of 100 histograms corresponding to 10 independent realizations of the noise for 10 different topologies. Shaded regions indicate the standard deviation among the samples. Red, blue, and green correspond to the morphological parameters of the WT, Ts65Dn, and TgDyrk1A cases, respectively. Dashed lines correspond to simulations with recurrent inhibitory synapses reduced to 0.3 of the original value. (c-f) Coefficient of variation (CV) of the ISI distribution of pyramidal neurons (c) and inhibitory neurons (d) for different values of external input λ. Symbols correspond to the average CV of 100 ISI distributions, and errorbars indicate the respective standard deviation. Panels (c) and (d) correspond to the default network parameters, whereas panels (e) and (f) correspond to simulations with recurrent inhibition reduced to 0.3 of the original value.
https://doi.org/10.1371/journal.pcbi.1012259.s001
(PDF)
S2 Fig. Gamma activity dynamics for an external input λ = 9 kHz.
Raster plots and mean firing rate of the pyramidal (black) and inhibitory interneurons (red) for the WT (panels (a) and (d)), Ts65Dn (panels (b) and (e))), and TgDyrk1A (panels (c) and (f)) morphologies. Panels (a-c) correspond to simulations with unperturbed recurrent inhibition (same as in Fig 4(a)–4(f)), and panels (d-f) correspond to recurrent inhibitory synapses reduced to 0.3 of the original value.
https://doi.org/10.1371/journal.pcbi.1012259.s002
(PDF)
S3 Fig. Gamma activity dynamics for an external input λ = 6 kHz.
Raster plots and mean firing rate of the pyramidal (black) and inhibitory interneurons (red) for the WT (panels (a) and (d)), Ts65Dn (panels (b) and (e))), and TgDyrk1A (panels (c) and (f)) morphologies. Panels (a-c) correspond to simulations with unperturbed recurrent inhibition, and panels (d-f) correspond to recurrent inhibitory synapses reduced to 0.3 of the original value.
https://doi.org/10.1371/journal.pcbi.1012259.s003
(PDF)
S4 Fig. Gamma activity dynamics for an external input λ = 12 kHz.
Raster plots and mean firing rate of the pyramidal (black) and inhibitory interneurons (red) for the WT (panels (a) and (d)), Ts65Dn (panels (b) and (e))), and TgDyrk1A (panels (c) and (f)) morphologies. Panels (a-c) correspond to simulations with unperturbed recurrent inhibition, and panels (d-f) correspond to recurrent inhibitory synapses reduced to 0.3 of the original value.
https://doi.org/10.1371/journal.pcbi.1012259.s004
(PDF)
S5 Fig. Gamma activity dynamics corresponding to selected fabricated morphologies.
(a-d) Raster plots and mean firing rate of the pyramidal (black) and inhibitory interneurons (red) for network topologies generated with different values of the SCP scaling parameter α and the scaling of the mean dendritic tree size with respect to WT β. Rest of the parameters as in Fig 5. (e,f) Average power spectra of the LFP signal corresponding to Fig 5 for specific values of α and β.
https://doi.org/10.1371/journal.pcbi.1012259.s005
(PDF)
S1 Movie. Simulation of the model corresponding to WT morphology.
Top: spatial representation of the cortex slice model. Each circle represents a neuron, whose spikes are marked with red (excitatory neurons) and blue (inhibitory neurons) for a short time interval. Bottom: Time series of the average firing rate of excitatory (red) and inhibitory (blue) neurons. Simulation parameters as in Fig 4(a) (λ = 9 kHz).
https://doi.org/10.1371/journal.pcbi.1012259.s006
(MP4)
S2 Movie. Simulation of the model corresponding to Ts65Dn morphology.
Top: spatial representation of the cortex slice model. Each circle represents a neuron, whose spikes are marked with red (excitatory neurons) and blue (inhibitory neurons) for a short time interval. Bottom: Time series of the average firing rate of excitatory (red) and inhibitory (blue) neurons. Simulation parameters as in Fig 4(a) (λ = 9 kHz).
https://doi.org/10.1371/journal.pcbi.1012259.s007
(MP4)
S3 Movie. Simulation of the model corresponding to TgDyrk1A morphology.
Top: spatial representation of the cortex slice model. Each circle represents a neuron, whose spikes are marked with red (excitatory neurons) and blue (inhibitory neurons) for a short time interval. Bottom: Time series of the average firing rate of excitatory (red) and inhibitory (blue) neurons. Simulation parameters as in Fig 4(a) (λ = 9 kHz).
https://doi.org/10.1371/journal.pcbi.1012259.s008
(MP4)
References
- 1. Lott IT, Dierssen M. Cognitive Deficits and Associated Neurological Complications in Individuals with Down’s Syndrome. The Lancet Neurology. 2010;9(6):623–633. pmid:20494326
- 2. Dierssen M. Alterations of Neocortical Pyramidal Cell Phenotype in the Ts65Dn Mouse Model of Down Syndrome: Effects of Environmental Enrichment. Cerebral Cortex. 2003;13(7):758–764. pmid:12816891
- 3. Martinez de Lagran M, Benavides-Piccione R, Ballesteros-Yanez I, Calvo M, Morales M, Fillat C, et al. Dyrk1A Influences Neuronal Morphogenesis Through Regulation of Cytoskeletal Dynamics in Mammalian Cortical Neurons. Cerebral Cortex. 2012;22(12):2867–2877. pmid:22215728
- 4. Dierssen M, Sierra C, Sabariego M, Fernández-Blanco A, Cruciani S, Zamora-Moratalla A, et al. The lncRNA Snhg11, a New Candidate Contributing to Neurogenesis, Plasticity and Memory Deficits in Down Syndrome. In Review; 2023. pmid:37841843
- 5. Marin-Padilla M. Pyramidal Cell Abnormalities in the Motor Cortex of a Child with Down’s Syndrome. A Golgi Study. Journal of Comparative Neurology. 1976;167(1):63–81. pmid:131810
- 6. Ruiz-Mejias M. Outer Brain Oscillations in Down Syndrome. Frontiers in Systems Neuroscience. 2019;13:17. pmid:31139056
- 7. Ruiz-Mejias M, Martinez de Lagran M, Mattia M, Castano-Prat P, Perez-Mendez L, Ciria-Suarez L, et al. Overexpression of Dyrk1A, a Down Syndrome Candidate, Decreases Excitability and Impairs Gamma Oscillations in the Prefrontal Cortex. The Journal of Neuroscience. 2016;36(13):3648–3659. pmid:27030752
- 8. Verret L, Mann EO, Hang GB, Barth AMI, Cobos I, Ho K, et al. Inhibitory Interneuron Deficit Links Altered Network Activity and Cognitive Dysfunction in Alzheimer Model. Cell. 2012;149(3):708–721. pmid:22541439
- 9. Mably AJ, Colgin LL. Gamma oscillations in cognitive disorders. Current opinion in neurobiology. 2018;52:182–187. pmid:30121451
- 10. Ward LM. Synchronous neural oscillations and cognitive processes. Trends in cognitive sciences. 2003;7(12):553–559. pmid:14643372
- 11.
Buzsáki G. Rhythms of the Brain. Oxford University Press; 2006.
- 12. Hanks TD, Summerfield C. Perceptual decision making in rodents, monkeys, and humans. Neuron. 2017;93(1):15–31. pmid:28056343
- 13. Whittington MA, Cunningham MO, LeBeau FEN, Racca C, Traub RD. Multiple Origins of the Cortical Gamma Rhythm. Developmental Neurobiology. 2011;71(1):92–106. pmid:21154913
- 14. Buzsáki G, Wang XJ. Mechanisms of Gamma Oscillations. Annual Review of Neuroscience. 2012;35(1):203–225. pmid:22443509
- 15. Başar E. A Review of Gamma Oscillations in Healthy Subjects and in Cognitive Impairment. International Journal of Psychophysiology. 2013;90(2):99–117. pmid:23892065
- 16. Hansen BJ, Dragoi V. Adaptation-Induced Synchronization in Laminar Cortical Circuits. Proceedings of the National Academy of Sciences. 2011;108(26):10720–10725. pmid:21659632
- 17. Brunet NM, Bosman CA, Vinck M, Roberts M, Oostenveld R, Desimone R, et al. Stimulus Repetition Modulates Gamma-Band Synchronization in Primate Visual Cortex. Proceedings of the National Academy of Sciences. 2014;111(9):3626–3631. pmid:24554080
- 18.
Davisson Muriel T A EC Schmidt Cecilia. Segmental trisomy for murine chromosome 16: a new system for studying Down syndrome. In: Patterson David E CJ, editor. Molecular genetics of chromosome 21 and Down syndrome. Wiley-Liss; 1990. p. 263–280.
- 19. Altafaj X, Dierssen M, Baamonde C, Martí E, Visa J, Guimerà J, et al. Neurodevelopmental delay, motor abnormalities and cognitive deficits in transgenic mice overexpressing Dyrk1A (minibrain), a murine model of Down’s syndrome. Human Molecular Genetics. 2001;10(18):1915–1923. pmid:11555628
- 20. Martínez de Lagrán M, Altafaj X, Gallego X, Martí E, Estivill X, Sahún I, et al. Motor phenotypic alterations in TgDyrk1a transgenic mice implicate DYRK1A in Down syndrome motor dysfunction. Neurobiology of Disease. 2004;15(1):132–142. pmid:14751778
- 21. Izhikevich EM. Simple model of spiking neurons. IEEE Transactions on Neural Networks. 2003;14(6):1569–1572. pmid:18244602
- 22. Orlandi JG, Soriano J, Alvarez-Lacalle E, Teller S, Casademunt J. Noise focusing and the emergence of coherent activity in neuronal cultures. Nature Physics. 2013;9(9):582–590.
- 23. Gandolfi D, Mapelli J, Solinas S, De Schepper R, Geminiani A, Casellato C, et al. A Realistic Morpho-Anatomical Connection Strategy for Modelling Full-Scale Point-Neuron Microcircuits. Scientific Reports. 2022;12(1):13864. pmid:35974119
- 24. Sholl DA. Dendritic Organization in the Neurons of the Visual and Motor Cortices of the Cat. Journal of Anatomy. 1953;87(4):387–406. pmid:13117757
- 25. Brunel N. Dynamics of Sparsely Connected Networks of Excitatory and Inhibitory Spiking Neurons. Journal of Computational Neuroscience. 2000;8(3):183–208. pmid:10809012
- 26. Brunel N, Wang XJ. What Determines the Frequency of Fast Network Oscillations With Irregular Neural Discharges? I. Synaptic Dynamics and Excitation-Inhibition Balance. Journal of Neurophysiology. 2003;90(1):415–430. pmid:12611969
- 27. Wang XJ. Neurophysiological and Computational Principles of Cortical Rhythms in Cognition. Physiological Reviews. 2010;90(3):1195–1268. pmid:20664082
- 28. Bartos M, Vida I, Jonas P. Synaptic Mechanisms of Synchronized Gamma Oscillations in Inhibitory Interneuron Networks. Nature Reviews Neuroscience. 2007;8(1):45–56. pmid:17180162
- 29. Börgers C, Epstein S, Kopell NJ. Background Gamma Rhythmicity and Attention in Cortical Local Circuits: A Computational Study. Proceedings of the National Academy of Sciences. 2005;102(19):7002–7007. pmid:15870189
- 30. Lee S, Jones SR. Distinguishing Mechanisms of Gamma Frequency Oscillations in Human Current Source Signals Using a Computational Model of a Laminar Neocortical Network. Frontiers in Human Neuroscience. 2013;7. pmid:24385958
- 31. Reimann MW, King JG, Muller EB, Ramaswamy S, Markram H. An Algorithm to Predict the Connectome of Neural Microcircuits. Frontiers in Computational Neuroscience. 2015;9. pmid:26500529
- 32. Markram H, Muller E, Ramaswamy S, Reimann MW, Abdellah M, Sanchez CA, et al. Reconstruction and Simulation of Neocortical Microcircuitry. Cell. 2015;163(2):456–492. pmid:26451489
- 33. Udvary D, Harth P, Macke JH, Hege HC, De Kock CPJ, Sakmann B, et al. The Impact of Neuron Morphology on Cortical Network Architecture. Cell Reports. 2022;39(2):110677. pmid:35417720
- 34. Hjorth JJJ, Kozlov A, Carannante I, Frost Nylén J, Lindroos R, Johansson Y, et al. The Microcircuits of Striatum in Silico. Proceedings of the National Academy of Sciences. 2020;117(17):9554–9565. pmid:32321828
- 35. Pyka M, Klatt S, Cheng S. Parametric Anatomical Modeling: A Method for Modeling the Anatomical Layout of Neurons and Their Projections. Frontiers in Neuroanatomy. 2014;8. pmid:25309338
- 36. Casali S, Marenzi E, Medini C, Casellato C, D’Angelo E. Reconstruction and Simulation of a Scaffold Model of the Cerebellar Network. Frontiers in Neuroinformatics. 2019;13:37. pmid:31156416
- 37. Whittington MA, Traub RD, Jefferys JGR. Synchronized oscillations in interneuron networks driven by metabotropic glutamate receptor activation. Nature. 1995;373(6515):612–615. pmid:7854418
- 38. Brunel N, Hakim V. Fast Global Oscillations in Networks of Integrate-and-Fire Neurons with Low Firing Rates. Neural Computation. 1999;11(7):1621–1671. pmid:10490941
- 39. Whittington MA, Traub RD, Kopell N, Ermentrout B, Buhl EH. Inhibition-based rhythms: experimental and mathematical observations on network dynamics. International Journal of Psychophysiology. 2000;38(3):315–336. pmid:11102670
- 40. Buhl EH, Schlote W. Intracellular Lucifer Yellow Staining and Electron Microscopy of Neurones in Slices of Fixed Epitumourous Human Cortical Tissue. Acta Neuropathologica. 1987;75(2):140–146. pmid:3434221
- 41. Einstein G. Intracellular Injection of Lucifer Yellow into Cortical Neurons in Lightly Fixed Sections and Its Application to Human Autopsy Material. Journal of Neuroscience Methods. 1988;26(2):95–103. pmid:3216684
- 42. Elston GN, Pow DV, Calford MB. Neuronal Composition and Morphology in Layer IV of Two Vibrissal Barrel Subfields of Rat Cortex. Cerebral Cortex. 1997;7(5):422–431. pmid:9261572
- 43. Voges N, Schüz A, Aertsen A, Rotter S. A Modeler’s View on the Spatial Structure of Intrinsic Horizontal Connectivity in the Neocortex. Progress in Neurobiology. 2010;92(3):277–292. pmid:20685378
- 44. Galassi M, Gough B. GNU Scientific Library: Reference Manual. Network Theory; 2003.