Search

In their Commentary paper, Villaverde and Massonis (On testing structural identifiability by a simple scaling method: relying on scaling symmetries can be misleading) have commented on our paper in which we proposed a simple scaling method to test structural identifiability. Our scaling invariance method (SIM) tests for scaling symmetries only, and Villaverde and Massonis correctly show the SIM may fail to detect identifiability problems when a model has other types of symmetries. We agree with the limitations raised by these authors but, also, we emphasize that the method is still valuable for its applicability to a wide variety of models, its simplicity, and even as a tool to introduce the problem of identifiability to investigators with little training in mathematics.

This paper uses mathematical modeling to study the mechanisms of surround suppression in the primate visual cortex. We present a large-scale neural circuit model consisting of three interconnected components: LGN and two input layers (Layer 4Ca and Layer 6) of the primary visual cortex V1, covering several hundred hypercolumns. Anatomical structures are incorporated and physiological parameters from realistic modeling work are used. The remaining parameters are chosen to produce model outputs that emulate experimentally observed size-tuning curves. Our two main results are: (i) we discovered the character of the long-range connections in Layer 6 responsible for surround effects in the input layers; and (ii) we showed that a net-inhibitory feedback, i.e., feedback that excites I-cells more than E-cells, from Layer 6 to Layer 4 is conducive to producing surround properties consistent with experimental data. These results are obtained through parameter selection and model analysis. The effects of nonlinear recurrent excitation and inhibition are also discussed. A feature that distinguishes our model from previous modeling work on surround suppression is that we have tried to reproduce realistic lengthscales that are crucial for quantitative comparison with data. Due to its size and the large number of unknown parameters, the model is computationally challenging. We demonstrate a strategy that involves first locating baseline values for relevant parameters using a linear model, followed by the introduction of nonlinearities where needed. We find such a methodology effective, and propose it as a possibility in the modeling of complex biological systems., Author summary The visual cortex is a part of the cortex that processes visual signal. The visual signal from the retina is relayed through the lateral geniculate nucleus (LGN), and enters the primary visual cortex directly. The primary visual cortex consists of many layers. A phenomenon called surround suppression, which means a reduction of neuronal activities in response to stimuli of increasing size, is well-observed in all layers of the primate primary visual cortex. In this paper, we built a large scale mathematical model consisting of the LGN input and tens of thousands of neuron groups in Layer 4 and Layer 6 of the primary visual cortex to study the mechanism of surround suppression. This model is constrained by both anatomical structures and realistic physiological parameters. Remaining unknown parameters are reverse-engineered from experimental data. We showed how the character of the surround suppression depends on the long-range connections in the feedback layer, as well as the feedback current from Layer 6 to Layer 4.

Striatal oscillatory activity is associated with movement, reward, and decision-making, and observed in several interacting frequency bands. Local field potential recordings in rodent striatum show dopamine- and reward-dependent transitions between two states: a “spontaneous” state involving β (∼15-30 Hz) and low γ (∼40-60 Hz), and a state involving θ (∼4-8 Hz) and high γ (∼60-100 Hz) in response to dopaminergic agonism and reward. The mechanisms underlying these rhythmic dynamics, their interactions, and their functional consequences are not well understood. In this paper, we propose a biophysical model of striatal microcircuits that comprehensively describes the generation and interaction of these rhythms, as well as their modulation by dopamine. Building on previous modeling and experimental work suggesting that striatal projection neurons (SPNs) are capable of generating β oscillations, we show that networks of striatal fast-spiking interneurons (FSIs) are capable of generating δ/θ (ie, 2 to 6 Hz) and γ rhythms. Under simulated low dopaminergic tone our model FSI network produces low γ band oscillations, while under high dopaminergic tone the FSI network produces high γ band activity nested within a δ/θ oscillation. SPN networks produce β rhythms in both conditions, but under high dopaminergic tone, this β oscillation is interrupted by δ/θ-periodic bursts of γ-frequency FSI inhibition. Thus, in the high dopamine state, packets of FSI γ and SPN β alternate at a δ/θ timescale. In addition to a mechanistic explanation for previously observed rhythmic interactions and transitions, our model suggests a hypothesis as to how the relationship between dopamine and rhythmicity impacts motor function. We hypothesize that high dopamine-induced periodic FSI γ-rhythmic inhibition enables switching between β-rhythmic SPN cell assemblies representing the currently active motor program, and thus that dopamine facilitates movement in part by allowing for rapid, periodic shifts in motor program execution., Author summary Striatal oscillatory activity is associated with movement, reward, and decision-making, and observed in several interacting frequency bands. The mechanisms underlying these rhythmic dynamics, their interactions, and their functional consequences are not well understood. In this paper, we propose a biophysical model of striatal microcircuits that comprehensively describes the generation and interaction of striatal rhythms, as well as their modulation by dopamine. Our model suggests a hypothesis as to how the relationship between dopamine and rhythmicity impacts the function of the motor system, enabling rapid, periodic shifts in motor program execution.

During early human immunodeficiency virus (HIV) infection multiple CD8+ T cell responses are elicited almost simultaneously. These responses exert strong selective pressures on different parts of HIV’s genome, and select for mutations that escape recognition and are thus beneficial to the virus. Some studies reveal that the later these escape mutations emerge, the more slowly they go to fixation. This pattern of escape rate decrease(ERD) can arise by distinct mechanisms. In particular, in large populations with high beneficial mutation rates interference among different escape strains –an effect that can emerge in evolution with asexual reproduction and results in delayed fixation times of beneficial mutations compared to sexual reproduction– could significantly impact the escape rates of mutations. In this paper, we investigated how interference between these concurrent escape mutations affects their escape rates in systems with multiple epitopes, and whether it could be a source of the ERD pattern. To address these issues, we developed a multilocus Wright-Fisher model of HIV dynamics with selection, mutation and recombination, serving as a null-model for interference. We also derived an interference-free null model assuming initial neutral evolution before immune response elicitation. We found that interference between several equally selectively advantageous mutations can generate the observed ERD pattern. We also found that the number of loci, as well as recombination rates substantially affect ERD. These effects can be explained by the underexponential decline of escape rates over time. Lastly, we found that the observed ERD pattern in HIV infected individuals is consistent with both independent, interference-free mutations as well as interference effects. Our results confirm that interference effects should be considered when analyzing HIV escape mutations. The challenge in estimating escape rates and mutation-associated selective coefficients posed by interference effects cannot simply be overcome by improved sampling frequencies or sizes. This problem is a consequence of the fundamental shortcomings of current estimation techniques under interference regimes. Hence, accounting for the stochastic nature of competition between mutations demands novel estimation methodologies based on the analysis of HIV strains, rather than mutation frequencies., Author Summary Within-host evolution of human immunodeficiency virus (HIV) is shaped by strong immune responses mounted against the virus. Multiple CD8+ T cell populations, each recognizing a specific part of an HIV protein, simultaneously suppress HIV growth. Escape mutations that arise in HIV genome regions coding for these virion protein parts, impair CD8+ T cell recognition and are consequently strongly selected. The emergence and rise of these escape mutations exhibits an intriguing temporal pattern: the earlier an escape mutation arises, the faster it goes to fixation. This pattern is termed escape rate decrease (ERD). In this paper, we tested computationally whether interference, i.e. the coexistence of multiple genetically distinct HIV strains engaged in competitive interaction within the host, could be a possible source of ERD. As an alternative, we also mathematically derived the temporal pattern of escapes under interference-free conditions, and compared this with data. We found that interference between multiple beneficial mutations could generate ERD. However, ERD does not imply the presence of interference. Thus, more detailed data is required to unambiguously determine whether interference effects influence ERD generation. Nevertheless, interference should be considered when studying the within-host evolution of HIV. Ignoring its effects on population dynamics can severely underestimate the protective capacity of CD8+ T cells.

We describe Orbit Image Analysis, an open-source whole slide image analysis tool. The tool consists of a generic tile-processing engine which allows the execution of various image analysis algorithms provided by either Orbit itself or from other open-source platforms using a tile-based map-reduce execution framework. Orbit Image Analysis is capable of sophisticated whole slide imaging analyses due to several key features. First, Orbit has machine-learning capabilities. This deep learning segmentation can be integrated with complex object detection for analysis of intricate tissues. In addition, Orbit can run locally as standalone or connect to the open-source image server OMERO. Another important characteristic is its scale-out functionality, using the Apache Spark framework for distributed computing. In this paper, we describe the use of Orbit in three different real-world applications: quantification of idiopathic lung fibrosis, nerve fibre density quantification, and glomeruli detection in the kidney., Author summary Whole slide images (WSI) are digital scans of samples, e.g. tissue sections. It is very convenient to view samples in this digital form, and with the increasing computation power it can also be used for quantification. These images are often too large to be analysed with standard tools. To overcome this issue, we created on open-source tool called Orbit Image Analysis which divides the images into smaller parts and allows the analysis of it with either embedded algorithms or the integration of existing tools. It also provides mechanisms to process huge amounts of images in distributed computing environments, such as clusters or cloud infrastructures. In this paper we describe the Orbit system and demonstrate its application based on three real-word use-cases.

A biophysical model that captures molecular homeostatic control of ions at the perisynaptic cradle (PsC) is of fundamental importance for understanding the interplay between astroglial and neuronal compartments. In this paper, we develop a multi-compartmental mathematical model which proposes a novel mechanism whereby the flow of cations in thin processes is restricted due to negatively charged membrane lipids which result in the formation of deep potential wells near the dipole heads. These wells restrict the flow of cations to “hopping” between adjacent wells as they transverse the process, and this surface retention of cations will be shown to give rise to the formation of potassium (K+) and sodium (Na+) microdomains at the PsC. We further propose that a K+ microdomain formed at the PsC, provides the driving force for the return of K+ to the extracellular space for uptake by the neurone, thereby preventing K+ undershoot. A slow decay of Na+ was also observed in our simulation after a period of glutamate stimulation which is in strong agreement with experimental observations. The pathological implications of microdomain formation during neuronal excitation are also discussed., Author summary During periods of neuronal activity, ionic homeostasis in the surrounding extracellular space (ECS) is disturbed. To provide a healthy environment for continued neuronal function, excess ions such as potassium must be buffered away from the ECS; a vital supportive role provided by astrocyte cells. It has long been thought that astrocytes not only removed ions from the ECS but also transport them to other areas of the brain where their concentrations are lower. However, while our computational model simulations agree that astrocytes do remove these ions from the ECS they also show that these ions are mainly stored locally at the PsC to be returned to the ECS, thus restoring ionic homeostasis. Furthermore, we detail in this paper that this happens because of a previously overlooked biophysical phenomenon that is only dominant in thin astrocyte processes. The flow of these cations within thin processes is primarily by surface conduction where they experience the attraction of fixed negative charge at the membrane inner surface. This negative charge constrains cation movement along the surface and so their flow rate is restricted. Consequently, ions such as potassium that are released during neuronal excitation enter the PsC and are stored locally due to the low conductance pathway between the PsC and the astrocyte soma. Our simulations also show that this local build-up of K+ is returned to the ECS after the neuronal activity dies off which could potentially explain why K+ undershoot has not been observed; this result agrees with experimental observations. Moreover, the same mechanism can also explain the transient behaviour of Na+ ions whereby in thin processes a slow decay time constant is experimentally observed. These findings have important implications for the role of astrocytes in regulating neuronal excitability under physiological and pathological conditions, and therefore highlight the significance of the work presented in this paper.

Activities of transcription factors (TFs) are temporally modulated to regulate dynamic cellular processes, including development, homeostasis, and disease. Recent developments of bioinformatic tools have enabled the analysis of TF activities using transcriptome data. However, because these methods typically use exon-based target expression levels, the estimated TF activities have limited temporal accuracy. To address this, we proposed a TF activity measure based on intron-level information in time-series RNA-seq data, and implemented it to decode the temporal control of TF activities during dynamic processes. We showed that TF activities inferred from intronic reads can better recapitulate instantaneous TF activities compared to the exon-based measure. By analyzing public and our own time-series transcriptome data, we found that intron-based TF activities improve the characterization of temporal phasing of cycling TFs during circadian rhythm, and facilitate the discovery of two temporally opposing TF modules during T cell activation. Collectively, we anticipate that the proposed approach would be broadly applicable for decoding global transcriptional architecture during dynamic processes., Author summary Many health-related cellular processes, such as immune response and disease progression, involve dynamic changes of gene expression state, which are orchestrated by transcription factors. Dissecting the activities of transcription factors is thus important for understanding cellular processes and for interfering with dysregulated processes. Our ability to analyze transcription factor activities has been facilitated by genome-wide gene expression data from high-throughput assays such as RNA sequencing. Existing methods typically estimate transcription factor activities based on the expression levels of matured mRNAs of target genes. However, because the levels of matured mRNAs are affected by transcriptional and post-transcriptional regulatory activities, the estimated transcription factor activities may not faithfully recapitulate the regulatory activities of transcription factors. In this paper, we proposed and validated an alternative approach for analyzing transcription factor activities using the expression levels of unmatured mRNAs of target genes, allowing us to decode how transcription factor activities are temporally controlled during key biological processes. Our results provide insights into the temporal phasing of key circadian regulator activities in mouse liver, and uncover two temporally opposing modules of transcription factors that dictate the immune responses in T cells. Therefore, this approach can help understand the regulatory principles of dynamic cellular processes.

System identification techniques -- projection pursuit regression models (PPRs) and convolutional neural networks (CNNs) -- provide state-of-the-art performance in predicting visual cortical neurons' responses to arbitrary input stimuli. However, the constituent kernels recovered by these methods are often noisy and lack coherent structure, making it difficult to understand the underlying component features of a neuron's receptive field. In this paper, we show that using a dictionary of diverse kernels with complex shapes learned from natural scenes based on efficient coding theory, as the front-end for PPRs and CNNs can improve their performance in neuronal response prediction as well as algorithmic data efficiency and convergence speed. Extensive experimental results also indicate that these sparse-code kernels provide important information on the component features of a neuron's receptive field. In addition, we find that models with the complex-shaped sparse code front-end are significantly better than models with a standard orientation-selective Gabor filter front-end for modeling V1 neurons that have been found to exhibit complex pattern selectivity. We show that the relative performance difference due to these two front-ends can be used to produce a sensitive metric for detecting complex selectivity in V1 neurons., 22 pages

Working memory is a cognitive function involving the storage and manipulation of latent information over brief intervals of time, thus making it crucial for context-dependent computation. Here, we use a top-down modeling approach to examine network-level mechanisms of working memory, an enigmatic issue and central topic of study in neuroscience. We optimize thousands of recurrent rate-based neural networks on a working memory task and then perform dynamical systems analysis on the ensuing optimized networks, wherein we find that four distinct dynamical mechanisms can emerge. In particular, we show the prevalence of a mechanism in which memories are encoded along slow stable manifolds in the network state space, leading to a phasic neuronal activation profile during memory periods. In contrast to mechanisms in which memories are directly encoded at stable attractors, these networks naturally forget stimuli over time. Despite this seeming functional disadvantage, they are more efficient in terms of how they leverage their attractor landscape and paradoxically, are considerably more robust to noise. Our results provide new hypotheses regarding how working memory function may be encoded within the dynamics of neural circuits., Author summary The ability to remember information for brief periods of time before using it is a key human ability. For example, retaining a phone number for a few moments prior to entering it into a keypad. Such ability, known as working memory, enables many more complex functions such as planning and reasoning. In this paper, we use theory and computational modeling approaches to try and better understand how circuits and networks in the brain might be achieving working memory. Specifically, we construct hypothetical network models to perform a task that embodies essential aspects of working memory. We then dissect our model to reveal how its components—simulated neural units—interact with each other to represent and maintain information. It turns out that some models can achieve memory by maintaining a fixed representation over time, i.e., the units remain ‘still’ during memory. However, other models fluctuate their activity during memory in a particular way that is seemingly quite efficient and less susceptible to distraction. In total, our computational study provides new theory for how this form of memory might be implemented in brain networks.

Animals rely on sensory feedback to generate accurate, reliable movements. In many flying insects, strain-sensitive neurons on the wings provide rapid feedback that is critical for stable flight control. While the impacts of wing structure on aerodynamic performance have been widely studied, the impacts of wing structure on sensing are largely unexplored. In this paper, we show how the structural properties of the wing and encoding by mechanosensory neurons interact to jointly determine optimal sensing strategies and performance. Specifically, we examine how neural sensors can be placed effectively on a flapping wing to detect body rotation about different axes, using a computational wing model with varying flexural stiffness. A small set of mechanosensors, conveying strain information at key locations with a single action potential per wingbeat, enable accurate detection of body rotation. Optimal sensor locations are concentrated at either the wing base or the wing tip, and they transition sharply as a function of both wing stiffness and neural threshold. Moreover, the sensing strategy and performance is robust to both external disturbances and sensor loss. Typically, only five sensors are needed to achieve near-peak accuracy, with a single sensor often providing accuracy well above chance. Our results show that small-amplitude, dynamic signals can be extracted efficiently with spatially and temporally sparse sensors in the context of flight. The demonstrated interaction of wing structure and neural encoding properties points to the importance of understanding each in the context of their joint evolution., Author summary In addition to generating forces for flight, insect wings also serve an important role as sensory structures, providing rapid feedback about wing bending that is used to stabilize flight. While much is known about how wing structure affects aerodynamic performance, the effects of wing structure on sensing remain unexplored. Using a computational model of a flapping wing, we examine how sensing strategies depend on wing stiffness and sensor properties. We show that body rotations can be accurately detected with a small number of sensors on the wing across a wide range of conditions. Optimal sensor locations are clustered at either the wing base or wing tip, depending on a combination of wing stiffness and sensor properties. Moreover, sensing performance is robust to multiple kinds of perturbations. Our work provides a basis for understanding how wing structure impacts incoming sensory information during flight.

Voltage imaging enables monitoring neural activity at sub-millisecond and sub-cellular scale, unlocking the study of subthreshold activity, synchrony, and network dynamics with unprecedented spatio-temporal resolution. However, high data rates (>800MB/s) and low signal-to-noise ratios create bottlenecks for analyzing such datasets. Here we present VolPy, an automated and scalable pipeline to pre-process voltage imaging datasets. VolPy features motion correction, memory mapping, automated segmentation, denoising and spike extraction, all built on a highly parallelizable, modular, and extensible framework optimized for memory and speed. To aid automated segmentation, we introduce a corpus of 24 manually annotated datasets from different preparations, brain areas and voltage indicators. We benchmark VolPy against ground truth segmentation, simulations and electrophysiology recordings, and we compare its performance with existing algorithms in detecting spikes. Our results indicate that VolPy’s performance in spike extraction and scalability are state-of-the-art., Author summary Roughly 290 million electrical action potentials occur every second in the human brain, facilitating the propagation of signals among cells in the nervous system and driving most of our daily operations. New methods in brain imaging are emerging that have the speed and resolution to capture events in the brain at the pace at which neurons typically communicate. These methods measure voltage in neurons by using light, and therefore can access very detailed brain signaling patterns. However, the adoption of these methods by a larger community, and not a restricted set of experts, is limited by the lack of computational tools, thereby greatly hindering progress in this field. In this paper, we present VolPy, a software framework that greatly facilitates the preprocessing of this new type of imaging datasets. This pipeline incorporates efficient and optimized algorithms that can identify neurons and extract their activity with great accuracy. The presented software will make this new imaging modality accessible to a wide audience.

There are two distinct classes of cells in the primary visual cortex (V1): simple cells and complex cells. One defining feature of complex cells is their spatial phase invariance; they respond strongly to oriented grating stimuli with a preferred orientation but with a wide range of spatial phases. A classical model of complete spatial phase invariance in complex cells is the energy model, in which the responses are the sum of the squared outputs of two linear spatially phase-shifted filters. However, recent experimental studies have shown that complex cells have a diverse range of spatial phase invariance and only a subset can be characterized by the energy model. While several models have been proposed to explain how complex cells could learn to be selective to orientation but invariant to spatial phase, most existing models overlook many biologically important details. We propose a biologically plausible model for complex cells that learns to pool inputs from simple cells based on the presentation of natural scene stimuli. The model is a three-layer network with rate-based neurons that describes the activities of LGN cells (layer 1), V1 simple cells (layer 2), and V1 complex cells (layer 3). The first two layers implement a recently proposed simple cell model that is biologically plausible and accounts for many experimental phenomena. The neural dynamics of the complex cells is modeled as the integration of simple cells inputs along with response normalization. Connections between LGN and simple cells are learned using Hebbian and anti-Hebbian plasticity. Connections between simple and complex cells are learned using a modified version of the Bienenstock, Cooper, and Munro (BCM) rule. Our results demonstrate that the learning rule can describe a diversity of complex cells, similar to those observed experimentally., Author summary Many cortical functions originate from the learning ability of the brain. How the properties of cortical cells are learned is vital for understanding how the brain works. There are many models that explain how V1 simple cells can be learned. However, how V1 complex cells are learned still remains unclear. In this paper, we propose a model of learning in complex cells based on the Bienenstock, Cooper, and Munro (BCM) rule. We demonstrate that properties of receptive fields of complex cells can be learned using this biologically plausible learning rule. Quantitative comparisons between the model and experimental data are performed. Results show that model complex cells can account for the diversity of complex cells found in experimental studies. In summary, this study provides a plausible explanation for how complex cells can be learned using biologically plausible plasticity mechanisms. Our findings help us to better understand biological vision processing and provide us with insights into the general signal processing principles that the visual cortex employs to process visual information.

We propose a developmental model inspired by the cortico-basal system (CX-BG) for vocal learning in babies and for solving the correspondence mismatch problem they face when they hear unfamiliar voices, with different tones and pitches. This model is based on the neural architecture INFERNO standing for Iterative Free-Energy Optimization of Recurrent Neural Networks. Free-energy minimization is used for rapidly exploring, selecting and learning the optimal choices of actions to perform (eg sound production) in order to reproduce and control as accurately as possible the spike trains representing desired perceptions (eg sound categories). We detail in this paper the CX-BG system responsible for linking causally the sound and motor primitives at the order of a few milliseconds. Two experiments performed with a small and a large audio database show the capabilities of exploration, generalization and robustness to noise of our neural architecture in retrieving audio primitives during vocal learning and during acoustic matching with unheared voices (different genders and tones)., Author summary We designed a developmental architecture inspired by the cortico-basal system for early vocal learning. Our neural system explores, evaluates and strengthens the motor primitives that match the best the sound repertoire created also dynamically. After a babbling process in which the network tests and aligns pronounced sound and motor vocal tracks, it is used for listening to novel voices, solving the correspondence problem.

The partial success of tumor immunotherapy induced by checkpoint blockade, which is not antigen-specific, suggests that the immune system of some patients contain antigen receptors able to specifically identify tumor cells. Here we focused on T-cell receptor (TCR) repertoires associated with spontaneous breast cancer. We studied the alpha and beta chain CDR3 domains of TCR repertoires of CD4 T cells using deep sequencing of cell populations in mice and applied the results to published TCR sequence data obtained from human patients. We screened peripheral blood T cells obtained monthly from individual mice spontaneously developing breast tumors by 5 months. We then looked at identical TCR sequences in published human studies; we used TCGA data from tumors and healthy tissues of 1,256 breast cancer resections and from 4 focused studies including sequences from tumors, lymph nodes, blood and healthy tissues, and from single cell dataset of 3 breast cancer subjects. We now report that mice spontaneously developing breast cancer manifest shared, Public CDR3 regions in both their alpha and beta and that a significant number of women with early breast cancer manifest identical CDR3 sequences. These findings suggest that the development of breast cancer is associated, across species, with biomarker, exclusive TCR repertoires., Author summary More than 1,100 cancer treatments are currently in clinical testing. Of these, 295 are immune-oncology medicines and vaccines. For Checkpoint Immunotherapy alone, 52,539 patients were enrolled in 2017. Of these treatments, 108 are for breast cancer, the leading cancer diagnosed in women. Although checkpoint immunotherapy has yet to be proven beneficial in most types of breast cancer, current observations of which specific T cell go into the tumor, and how these cells are connected to the response to treatment, have lead to a new understating of the role of the immune system. Yet, and in spite of this intense research into the role of T cells in tumors, very little is known about how clones of T cells behave during tumor formation. The paper we present here studies this clonal behavior over time. We study it during mammary tumor development in mice. We see that a set of T cell clones, which comes from a diversity of different nucleotide sequences, stands out in its behavior compared with other T cell clones. Importantly, we show that these mouse T cell clones have a strong connection to human breast cancer. The specific identification of a subset of T cell clones, together with the behavior of these clones over time, is critical to our understanding of the effect of modifications to the immune system, such as the modifications made by immunotherapy treatments, and to our understanding of how to better control the behavior of the immune system.

Learning the reliability of different sources of rewards is critical for making optimal choices. However, despite the existence of detailed theory describing how the expected reward is learned in the basal ganglia, it is not known how reward uncertainty is estimated in these circuits. This paper presents a class of models that encode both the mean reward and the spread of the rewards, the former in the difference between the synaptic weights of D1 and D2 neurons, and the latter in their sum. In the models, the tendency to seek (or avoid) options with variable reward can be controlled by increasing (or decreasing) the tonic level of dopamine. The models are consistent with the physiology of and synaptic plasticity in the basal ganglia, they explain the effects of dopaminergic manipulations on choices involving risks, and they make multiple experimental predictions., Author Summary To maximize their chances for survival, animals need to base their decisions not only on the average consequences of chosen actions, but also on the variability of the rewards resulting from these actions. For example, when an animal’s food reserves are depleted, it should prefer to forage in an area where food is guaranteed over an area where the amount of food is higher on average but variable, thus avoiding the risk of starvation. To implement such policies, animals need to be able to learn about variability of rewards resulting from taking different actions. This paper proposes how such learning may be implemented in a circuit of subcortical nuclei called the basal ganglia. It also suggests how the information about reward uncertainty can be used during decision making, so that animals can make choices that not only maximize expected rewards but also minimize risks.

Since the seminal 1961 paper of Monod and Jacob, mathematical models of biomolecular circuits have guided our understanding of cell regulation. Model-based exploration of the functional capabilities of any given circuit requires systematic mapping of multidimensional spaces of model parameters. Despite significant advances in computational dynamical systems approaches, this analysis remains a nontrivial task. Here, we use a nonlinear system of ordinary differential equations to model oocyte selection in Drosophila, a robust symmetry-breaking event that relies on autoregulatory localization of oocyte-specification factors. By applying an algorithmic approach that implements symbolic computation and topological methods, we enumerate all phase portraits of stable steady states in the limit when nonlinear regulatory interactions become discrete switches. Leveraging this initial exact partitioning and further using numerical exploration, we locate parameter regions that are dense in purely asymmetric steady states when the nonlinearities are not infinitely sharp, enabling systematic identification of parameter regions that correspond to robust oocyte selection. This framework can be generalized to map the full parameter spaces in a broad class of models involving biological switches., Author summary Identification of qualitatively different regimes in models of biomolecular switches is essential for understanding dynamics of complex biological processes, including symmetry breaking in cells and cell networks. We demonstrate how topological methods, symbolic computation, and numerical simulations can be combined for systematic mapping of symmetry-broken states in a mathematical model of oocyte specification in Drosophila, a leading experimental system of animal oogenesis. Our algorithmic framework reveals global connectedness of parameter domains corresponding to robust oocyte specification and enables systematic navigation through multidimensional parameter spaces in a large class of biomolecular switches.

Gene expression atlases have transformed our understanding of the development, composition and function of human tissues. New technologies promise improved cellular or molecular resolution, and have led to the identification of new cell types, or better defined cell states. But as new technologies emerge, information derived on old platforms becomes obsolete. We demonstrate that it is possible to combine a large number of different profiling experiments summarised from dozens of laboratories and representing hundreds of donors, to create an integrated molecular map of human tissue. As an example, we combine 850 samples from 38 platforms to build an integrated atlas of human blood cells. We achieve robust and unbiased cell type clustering using a variance partitioning method, selecting genes with low platform bias relative to biological variation. Other than an initial rescaling, no other transformation to the primary data is applied through batch correction or renormalisation. Additional data, including single-cell datasets, can be projected for comparison, classification and annotation. The resulting atlas provides a multi-scaled approach to visualise and analyse the relationships between sets of genes and blood cell lineages, including the maturation and activation of leukocytes in vivo and in vitro. In allowing for data integration across hundreds of studies, we address a key reproduciblity challenge which is faced by any new technology. This allows us to draw on the deep phenotypes and functional annotations that accompany traditional profiling methods, and provide important context to the high cellular resolution of single cell profiling. Here, we have implemented the blood atlas in the open access Stemformatics.org platform, drawing on its extensive collection of curated transcriptome data. The method is simple, scalable and amenable for rapid deployment in other biological systems or computational workflows., Author summary Combining data from many different studies is an attractive way of capturing new aspects of the biology being studied. Biological variance attributable to cell type, cellular niche, origin, disease status or environmental stimuli is the basis of most small-n transcriptome studies. In aggregation, these promise to capture emergent dimensions of a biology that is not possible to view from any individual study. However biological signal is easily swamped by technical artifact, especially when data is generated on platforms with profoundly different data structures. This is the case when comparing microarray data to RNAseq, or RNAseq to single cell profiling. Consequently, transcriptome atlases are generally comprised from a small number of donors/conditions surveyed using one technology platform. In this paper we present a simple and scalable data integration method that is platform agnostic. We provide a proof-of-principle by constructing an atlas of blood cells that combines many data sets measured on different platforms, and that in combination, recapitulates the known blood hierarchy. The atlas provides a reference to compare external samples to, allowing users to benchmark new derivation or isolation methods. It also provides a reference point for new data types, such as the classification of single cells. The approach allows for FAIR data reuse and robust identification of molecular signatures across multiple studies and experimental conditions.

The brain consists of many interconnected networks with time-varying, partially autonomous activity. There are multiple sources of noise and variation yet activity has to eventually converge to a stable, reproducible state (or sequence of states) for its computations to make sense. We approached this problem from a control-theory perspective by applying contraction analysis to recurrent neural networks. This allowed us to find mechanisms for achieving stability in multiple connected networks with biologically realistic dynamics, including synaptic plasticity and time-varying inputs. These mechanisms included inhibitory Hebbian plasticity, excitatory anti-Hebbian plasticity, synaptic sparsity and excitatory-inhibitory balance. Our findings shed light on how stable computations might be achieved despite biological complexity. Crucially, our analysis is not limited to analyzing the stability of fixed geometric objects in state space (e.g points, lines, planes), but rather the stability of state trajectories which may be complex and time-varying., Author summary Stability is essential for any complex system including, and perhaps especially, the brain. The brain’s neural networks are highly dynamic and noisy. Activity fluctuates from moment to moment and can be highly variable. Yet it is critical that these networks reach a consistent state (or sequence of states) for their computations to make sense. Failures in stability have consequences ranging from mild (e.g incorrect decisions) to severe (disease states). In this paper we use tools from control theory and dynamical systems theory to find mechanisms which produce stability in recurrent neural networks (RNNs). We show that a kind of “unlearning” (inhibitory Hebbian and excitatory anti-Hebbian plasticity), balance of excitation and inhibition, and sparse anatomical connectivity all lead to stability. Crucially, we focus on the stability of neural trajectories. This is different from traditional studies of stability of fixed points or planes. We do not assess what trajectories our networks will follow but, rather, when these trajectories will all converge towards each other to achieve stability.

Electrical stimulation of neural systems is a key tool for understanding neural dynamics and ultimately for developing clinical treatments. Many applications of electrical stimulation affect large populations of neurons. However, computational models of large networks of spiking neurons are inherently hard to simulate and analyze. We evaluate a reduced mean-field model of excitatory and inhibitory adaptive exponential integrate-and-fire (AdEx) neurons which can be used to efficiently study the effects of electrical stimulation on large neural populations. The rich dynamical properties of this basic cortical model are described in detail and validated using large network simulations. Bifurcation diagrams reflecting the network’s state reveal asynchronous up- and down-states, bistable regimes, and oscillatory regions corresponding to fast excitation-inhibition and slow excitation-adaptation feedback loops. The biophysical parameters of the AdEx neuron can be coupled to an electric field with realistic field strengths which then can be propagated up to the population description. We show how on the edge of bifurcation, direct electrical inputs cause network state transitions, such as turning on and off oscillations of the population rate. Oscillatory input can frequency-entrain and phase-lock endogenous oscillations. Relatively weak electric field strengths on the order of 1 V/m are able to produce these effects, indicating that field effects are strongly amplified in the network. The effects of time-varying external stimulation are well-predicted by the mean-field model, further underpinning the utility of low-dimensional neural mass models., Author summary Weak electrical inputs to the brain in vivo using transcranial electrical stimulation or in isolated cortex in vitro can affect the dynamics of the underlying neural populations. However, it is poorly understood what the exact mechanisms are that modulate the activity of neural populations as a whole and why the responses are so diverse in stimulation experiments. Despite this, electrical stimulation techniques are being developed for the treatment of neurological diseases in humans. To better understand these interactions, it is often necessary to simulate and analyze very large networks of neurons, which can be computationally demanding. In this theoretical paper, we present a reduced model of coupled neural populations that represents a piece of cortical tissue. This efficient model retains the dynamical properties of the large network of neurons it is based on while being several orders of magnitude faster to simulate. Due to the biophysical properties of the neuron model, an electric field can be coupled to the population. We show that weak electric fields often used in stimulation experiments can lead to entrainment of neural oscillations on the population level, and argue that the responses critically depend on the dynamical state of the neural system.

Over the past two decades, substantial amount of work has been conducted to characterize different odorant receptors, neuroanatomy and odorant response properties of the early olfactory system of Drosophila melanogaster. Yet many odorant receptors remain only partially characterized, and the odorant transduction process and the axon hillock spiking mechanism of the olfactory sensory neurons (OSNs) have yet to be fully determined. Identity and concentration, two key characteristics of the space of odorants, are encoded by the odorant transduction process. Detailed molecular models of the odorant transduction process are, however, scarce for fruit flies. To address these challenges we advance a comprehensive model of fruit fly OSNs as a cascade consisting of an odorant transduction process (OTP) and a biophysical spike generator (BSG). We model odorant identity and concentration using an odorant-receptor binding rate tensor, modulated by the odorant concentration profile, and an odorant-receptor dissociation rate tensor, and quantitatively describe the mechanics of the molecular ligand binding/dissociation of the OTP. We model the BSG as a Connor-Stevens point neuron. The resulting spatio-temporal encoding model of the Drosophila antenna provides a theoretical foundation for understanding the neural code of both odorant identity and odorant concentration and advances the state-of-the-art in a number of ways. First, it quantifies on the molecular level the spatio-temporal level of complexity of the transformation taking place in the antennae. The concentration-dependent spatio-temporal code at the output of the antenna circuits determines the level of complexity of olfactory processing in the downstream neuropils, such as odorant recognition and olfactory associative learning. Second, the model is biologically validated using multiple electrophysiological recordings. Third, the model demonstrates that the currently available data for odorant-receptor responses only enable the estimation of the affinity of the odorant-receptor pairs. The odorant-dissociation rate is only available for a few odorant-receptor pairs. Finally, our model calls for new experiments for massively identifying the odorant-receptor dissociation rates of relevance to flies., Author summary Identity and concentration, intrinsically embedded in the odorant space, are two key characteristics of olfactory coding that define the level of complexity of neural processing throughout the olfactory system in the fruit fly. In this paper we advance a theoretical foundation for understanding these two characteristics by quantifying mathematically the odorant space and devising a biophysical model of the olfactory sensory neurons (OSNs). To validate our modeling approach, we propose and apply an algorithm to estimate the affinity value and the dissociation rate, the two characteristics that define odorant identity, of multiple odorant-receptor pairs. We then evaluate our model with a multitude of odorant waveforms and demonstrate that the model output reproduces the temporal responses of OSNs obtained from in vivo electrophysiology recordings. Furthermore, we evaluate the model at the OSN population level and quantify on the molecular level the spatio-temporal level of complexity of the transformation taking place between the odorant space and the OSNs. The resulting concentration-dependent spatio-temporal code determines the level of complexity of the input space driving olfactory processing in the downstream neuropils. Lastly, our model demonstrates that the currently available data for OSN responses only enables estimation of affinity value. This calls for new experiments for massively identifying the odorant-receptor dissociation rates of relevance to flies.

Antigen recognition by T-cells is guided by the T-cell receptor (TCR) heterodimer formed by α and β chains. A huge diversity of TCR sequences should be maintained by the immune system in order to be able to mount an effective response towards foreign pathogens, so, due to cooperative binding of α and β chains to the pathogen, any constraints on chain pairing can have a profound effect on immune repertoire structure, diversity and antigen specificity. By integrating available structural data and paired chain sequencing results we were able to show that there are almost no constraints on pairing in TCRαβ complexes, allowing naive T-cell repertoire to reach the highest possible diversity. Additional analysis reveals that the specific choice of contacting amino acids can still have a profound effect on complex conformation. Moreover, antigen-driven selection can distort the uniform landscape of chain pairing, while small, yet significant, differences in the pairing can be attributed to various specialized T-cell subsets such as MAIT and iNKT T-cells, as well as other TCR sets specific to certain antigens., Author summary In the present paper we study chain pairing preferences in the T-cell receptor (TCR) heterodimer complex. The TCR molecule is formed by α and β chains and binding of both of these chains to an antigen presented by the major histocompatibility complex (MHC) molecule is required in order to trigger an immune response against foreign pathogens and neoantigens. We show that chain pairing in the TCR complex is nearly random ensuring a highly diverse set of TCRs required for recognition of a vast set of antigens. Our results also show that chain pairing preferences can nevertheless influence TCR complex geometry and biases in TCR chain pairing can be used to identify antigen-driven selection or selection towards specialized subsets of T-cells such as mucosal-associated and natural killer invariant T-cells.

Cells contain a thiol redox regulatory network to coordinate metabolic and developmental activities with exogenous and endogenous cues. This network controls the redox state and activity of many target proteins. Electrons are fed into the network from metabolism and reach the target proteins via redox transmitters such as thioredoxin (TRX) and NADPH-dependent thioredoxin reductases (NTR). Electrons are drained from the network by reactive oxygen species (ROS) through thiol peroxidases, e.g., peroxiredoxins (PRX). Mathematical modeling promises access to quantitative understanding of the network function and was implemented by using published kinetic parameters combined with fitting to known biochemical data. Two networks were assembled, namely the ferredoxin (FDX), FDX-dependent TRX reductase (FTR), TRX, fructose-1,6-bisphosphatase (FBPase) pathway with 2-cysteine PRX/ROS as oxidant, and separately the FDX, FDX-dependent NADP reductase (FNR), NADPH, NTRC-pathway for 2-CysPRX reduction. Combining both modules allowed drawing several important conclusions of network performance. The resting H2O2 concentration was estimated to be about 30 nM in the chloroplast stroma. The electron flow to metabolism exceeds that into thiol regulation of FBPase more than 7000-fold under physiological conditions. The electron flow from NTRC to 2-CysPRX is about 5.32-times more efficient than that from TRX-f1 to 2-CysPRX. Under severe stress (30 μM H2O2) the ratio of electron flow to the thiol network relative to metabolism sinks to 1:251 whereas the ratio of e- flow from NTRC to 2-CysPRX and TRX-f1 to 2-CysPRX rises up to 1:67. Thus, the simulation provides clues on experimentally inaccessible parameters and describes the functional state of the chloroplast thiol regulatory network., Author summary The state of the thiol redox regulatory network is a fundamental feature of all cells and determines metabolic and developmental processes. However, only some parameters are quantifiable in experiments. This paper establishes simplified mathematical models which enable simulation of electron flows through the regulatory system. This in turn allows for estimating rates and states of components of the network and to tentatively address previously unknown parameters such as the resting hydrogen peroxide levels or the expenditure of reductive power for regulation relative to metabolism. The establishment of such models for simulating the performance and dynamics of the redox regulatory network is of significance not only for photosynthesis but also, e.g., in bacterial and animal cells exposed to environmental stress or pathological disorders.

Identification of induced pluripotent stem (iPS) progenitor cells, the iPS forming cells in early stage of reprogramming, could provide valuable information for studying the origin and underlying mechanism of iPS cells. However, it is very difficult to identify experimentally since there are no biomarkers known for early progenitor cells, and only about 6 days after reprogramming initiation, iPS cells can be experimentally determined via fluorescent probes. What is more, the ratio of progenitor cells during early reprograming period is below 5%, which is too low to capture experimentally in the early stage. In this paper, we propose a novel computational approach for the identification of iPS progenitor cells based on machine learning and microscopic image analysis. Firstly, we record the reprogramming process using a live cell imaging system after 48 hours of infection with retroviruses expressing Oct4, Sox2 and Klf4, later iPS progenitor cells and normal murine embryonic fibroblasts (MEFs) within 3 to 5 days after infection are labeled by retrospectively tracing the time-lapse microscopic image. We then calculate 11 types of cell morphological and motion features such as area, speed, etc., and select best time windows for modeling and perform feature selection. Finally, a prediction model using XGBoost is built based on the selected six types of features and best time windows. Our model allows several missing values/frames in the sample datasets, thus it is applicable to a wide range of scenarios. Cross-validation, holdout validation and independent test experiments show that the minimum precision is above 52%, that is, the ratio of predicted progenitor cells within 3 to 5 days after viral infection is above 52%. The results also confirm that the morphology and motion pattern of iPS progenitor cells is different from that of normal MEFs, which helps with the machine learning methods for iPS progenitor cell identification., Author summary Identification of induced pluripotent stem (iPS) progenitor cells could provide valuable information for studying the origin and underlying mechanism of iPS cells. However, it is very difficult to identify experimentally since there are no biomarkers known for early progenitor cells, and only after about 6 days of induction, iPS cells can be experimentally determined via fluorescent probes. What is more, the percentage of the progenitor cells during the early induction period is below 5%, too low to capture experimentally in early stage. In this work, we proposed an approach for the identification of iPS progenitor cells, the iPS forming cells, based on machine learning and microscopic image analysis. The aim is to help biologists to enrich iPS progenitor cells during the early stage of induction, which allows experimentalists to select iPS progenitor cells with much higher probability, and furthermore to study the biomarkers which trigger the reprogramming process.

Cortical spreading depression (SD) is a spreading disruption of ionic homeostasis in the brain during which neurons experience complete and prolonged depolarizations. SD is the basis of migraine aura and is increasingly associated with many other brain pathologies. Here, we study the role of glutamate and NMDA receptor dynamics in the context of an ionic electrodiffusion model. We perform simulations in one (1D) and two (2D) spatial dimension. Our 1D simulations reproduce the “inverted saddle” shape of the extracellular voltage signal for the first time. Our simulations suggest that SD propagation depends on two overlapping mechanisms; one dependent on extracellular glutamate diffusion and NMDA receptors and the other dependent on extracellular potassium diffusion and persistent sodium channel conductance. In 2D simulations, we study the dynamics of spiral waves. We study the properties of the spiral waves in relation to the planar 1D wave, and also compute the energy expenditure associated with the recurrent SD spirals., Author summary Cortical spreading depression is a wave of neuronal silencing and ion concentration changes that sweeps slowly through the brain. It is the basis of migraine aura, in which a migraine patient sees a defect move through ones visual field 30 minutes to an hour prior to the headache attack. In this paper, we study the mechanisms by which cortical spreading depression travels as a wave through the brain. We construct a detailed mathematical model based on the physics of ion movement as well as what is known about the molecular players of cortical spreading depression. We build a simulation program that successfully solves the resulting set of highly coupled equations. We find that cortical spreading depression can propagate as a wave by distinct but overlapping mechanisms. We also simulate spiral cortical spreading depression waves and study their properties.

Pluripotent stem cells retain the developmental timing of their species of origin in vitro, an observation that suggests the existence of a cell-intrinsic developmental clock, yet the nature and machinery of the clock remain a mystery. We hypothesize that one possible component may lie in species-specific differences in the kinetics of transcriptional responses to differentiation signals. Using a liquid-handling robot, mouse and human pluripotent stem cells were exposed to identical neural differentiation conditions and sampled for RNA-sequencing at high frequency, every 4 or 10 minutes, for the first 10 hours of differentiation to test for differences in transcriptomic response rates. The majority of initial transcriptional responses occurred within a rapid window in the first minutes of differentiation for both human and mouse stem cells. Despite similarly early onsets of gene expression changes, we observed shortened and condensed gene expression patterns in mouse pluripotent stem cells compared to protracted trends in human pluripotent stem cells. Moreover, the speed at which individual genes were upregulated, as measured by the slopes of gene expression changes over time, was significantly faster in mouse compared to human cells. These results suggest that downstream transcriptomic response kinetics to signaling cues are faster in mouse versus human cells, and may offer a partial account for the vast differences in developmental rates across species., Author summary For regenerative medicine, a slow developmental clock operating in human pluripotent stem cells imposes severe obstacles to many pluripotent stem cell therapies. Producing the desired cell types with the necessary level of physiological maturity in clinically-relevant time frames will be essential if the therapy is aimed at being patient-specific. Thus, understanding what controls species-specific rates of differentiation and identifying ways to accelerate it would have a broad impact on improving the therapeutic value of human pluripotent stem cells. In this paper, we measured the immediate-early gene expression responses to neural differentiation signals in mouse and human pluripotent stem cells. We quantitatively summarized the gene expression patterns using a flexible model that allowed us to identify shared trends and dynamic changes across species. We found that responses occurred equally early across species, but expression patterns were prolonged in human cells and the rate of upregulation in gene expression was significantly faster in mouse cells. Our results represent a significant contribution in furthering the knowledge of early developmental transcriptional programs.

Respiratory syncytial virus (RSV) infection results in millions of hospitalizations and thousands of deaths each year. Variations in the adaptive and innate immune response appear to be associated with RSV severity. To investigate the host response to RSV infection in infants, we performed a systems-level study of RSV pathophysiology, incorporating high-throughput measurements of the peripheral innate and adaptive immune systems and the airway epithelium and microbiota. We implemented a novel multi-omic data integration method based on multilayered principal component analysis, penalized regression, and feature weight back-propagation, which enabled us to identify cellular pathways associated with RSV severity. In both airway and immune cells, we found an association between RSV severity and activation of pathways controlling Th17 and acute phase response signaling, as well as inhibition of B cell receptor signaling. Dysregulation of both the humoral and mucosal response to RSV may play a critical role in determining illness severity., Author summary This paper presents a novel approach to understanding the localized molecular responses to respiratory syncytial virus (RSV) and the system-level correlates of clinical outcomes. To do this, we developed a novel statistical method able to integrate high dimensional molecular data characterizing the host airway microbota and immune and nasal gene expression. We show that this integrative approach facilitates superior performance in estimating clinical outcome as opposed to any single data type. Using this approach, we identified both cell type-specific and shared biomarkers and regulatory pathways associated with RSV severity. Specifically, we identified an association between RSV severity, activation of pathways controlling Th17, and inhibition of B cell receptor signaling, which were present in both the site of infection airway and in peripheral immune cells. These results can guide future efforts to identify biomarkers for identifying or predicting illness severity following infant RSV infection. They may also be useful as biomarkers to inform the efficacy of future interventions (e.g., therapies) or preventative measures to suppress the rate of severe disease (e.g., vaccines).

Regulation of cytosolic calcium (Ca2+) dynamics is fundamental to microglial function. Temporal and spatial Ca2+ fluxes are induced from a complicated signal transduction pathway linked to brain ionic homeostasis. In this paper, we develop a novel biophysical model of Ca2+ and sodium (Na+) dynamics in human microglia and evaluate the contribution of purinergic receptors (P2XRs) to both intracellular Ca2+ and Na+ levels in response to agonist/ATP binding. This is the first comprehensive model that integrates P2XRs to predict intricate Ca2+ and Na+ transient responses in microglia. Specifically, a novel compact biophysical model is proposed for the capture of whole-cell patch-clamp currents associated with P2X4 and P2X7 receptors, which is composed of only four state variables. The entire model shows that intricate intracellular ion dynamics arise from the coupled interaction between P2X4 and P2X7 receptors, the Na+/Ca2+ exchanger (NCX), Ca2+ extrusion by the plasma membrane Ca2+ ATPase (PMCA), and Ca2+ and Na+ leak channels. Both P2XRs are modelled as two separate adenosine triphosphate (ATP) gated Ca2+ and Na+ conductance channels, where the stoichiometry is the removal of one Ca2+ for the hydrolysis of one ATP molecule. Two unique sets of model parameters were determined using an evolutionary algorithm to optimise fitting to experimental data for each of the receptors. This allows the proposed model to capture both human P2X7 and P2X4 data (hP2X7 and hP2X4). The model architecture enables a high degree of simplicity, accuracy and predictability of Ca2+ and Na+ dynamics thus providing quantitative insights into different behaviours of intracellular Na+ and Ca2+ which will guide future experimental research. Understanding the interactions between these receptors and other membrane-bound transporters provides a step forward in resolving the qualitative link between purinergic receptors and microglial physiology and their contribution to brain pathology., Author summary Mathematical modelling and computer simulation are powerful tools by which we can analyse complex biological systems, particularly, neural phenomena involved in brain dysfunction. In this research, we develop a theoretical foundation for studying P2X-mediated calcium and sodium signalling in human microglial cells. Microglia, which are brain-resident macrophages, restructure their intracellular actin cytoskeleton to enable motility; this restructuring requires a complex molecular cascade involving a set of ionic channels, membrane-coupled receptors and cytosolic components. Recent studies highlight the importance for increasing our understanding of microglia physiology, since their functions play critical roles in both normal physiological and pathological dynamics of the brain. There is a need to develop reliable human cellular models to investigate the biology of microglia aimed at understanding the influence of purinergic signalling in brain dysfunction to provide novel drug discovery targets. In this work, a detailed mathematical model is built for the dynamics of human P2XRs in microglia. Subsequently, experimental whole-cell currents are used to derive P2X-mediated electrophysiology of human microglia (i.e. sodium and calcium dynamics, and membrane potential). Our predictions reveal new quantitative insights into P2XRs on how they regulate ionic concentrations in terms of physiological interactions and transient responses.

Information processing can leave distinct footprints on the statistics of neural spiking. For example, efficient coding minimizes the statistical dependencies on the spiking history, while temporal integration of information may require the maintenance of information over different timescales. To investigate these footprints, we developed a novel approach to quantify history dependence within the spiking of a single neuron, using the mutual information between the entire past and current spiking. This measure captures how much past information is necessary to predict current spiking. In contrast, classical time-lagged measures of temporal dependence like the autocorrelation capture how long—potentially redundant—past information can still be read out. Strikingly, we find for model neurons that our method disentangles the strength and timescale of history dependence, whereas the two are mixed in classical approaches. When applying the method to experimental data, which are necessarily of limited size, a reliable estimation of mutual information is only possible for a coarse temporal binning of past spiking, a so-called past embedding. To still account for the vastly different spiking statistics and potentially long history dependence of living neurons, we developed an embedding-optimization approach that does not only vary the number and size, but also an exponential stretching of past bins. For extra-cellular spike recordings, we found that the strength and timescale of history dependence indeed can vary independently across experimental preparations. While hippocampus indicated strong and long history dependence, in visual cortex it was weak and short, while in vitro the history dependence was strong but short. This work enables an information-theoretic characterization of history dependence in recorded spike trains, which captures a footprint of information processing that is beyond time-lagged measures of temporal dependence. To facilitate the application of the method, we provide practical guidelines and a toolbox., Author summary Even with exciting advances in recording techniques of neural spiking activity, experiments only provide a comparably short glimpse into the activity of only a tiny subset of all neurons. How can we learn from these experiments about the organization of information processing in the brain? To that end, we exploit that different properties of information processing leave distinct footprints on the firing statistics of individual spiking neurons. In our work, we focus on a particular statistical footprint: We quantify how much a single neuron’s spiking depends on its own preceding activity, which we call history dependence. By quantifying history dependence in neural spike recordings, one can, in turn, infer some of the properties of information processing. Because recording lengths are limited in practice, a direct estimation of history dependence from experiments is challenging. The embedding optimization approach that we present in this paper aims at extracting a maximum of history dependence within the limits set by a reliable estimation. The approach is highly adaptive and thereby enables a meaningful comparison of history dependence between neurons with vastly different spiking statistics, which we exemplify on a diversity of spike recordings. In conjunction with recent, highly parallel spike recording techniques, the approach could yield valuable insights on how hierarchical processing is organized in the brain.

[Abstract] Brain networks exhibit very variable and dynamical functional connectivity and flexible configurations of information exchange despite their overall fixed structure. Brain oscillations are hypothesized to underlie time-dependent functional connectivity by periodically changing the excitability of neural populations. In this paper, we investigate the role of the connection delay and the detuning between the natural frequencies of neural populations in the transmission of signals. Based on numerical simulations and analytical arguments, we show that the amount of information transfer between two oscillating neural populations could be determined by their connection delay and the mismatch in their oscillation frequencies. Our results highlight the role of the collective phase response curve of the oscillating neural populations for the efficacy of signal transmission and the quality of the information transfer in brain networks., [Author summary] Collective dynamics in brain networks are characterized by a coordinated activity of their constituent neurons that lead to brain oscillations. Many evidences highlight the role that brain oscillations play in signal transmission, the control of the effective communication between brain areas, and the integration of information processed by different specialized regions. Oscillations periodically modulate the excitability of neurons and determine the response of those areas receiving the signals. Based on the communication through coherence (CTC) theory, the adjustment of the phase difference between local oscillations of connected areas can specify the timing of exchanged signals and therefore, the efficacy of the communication channels. In this respect, an important factor is the delay in the transmission of signals from one region to another that affects the phase difference and timing, and consequently the impact of the signals. Despite this delay plays an essential role in CTC theory, its role has been mostly overlooked in previous studies. In this manuscript, we concentrate on the role that the connection delay and the oscillation frequency of the populations play in the signal transmission, and consequently in the effective connectivity, between two brain areas. Through extensive numerical simulations, as well as analytical results with reduced models, we show that these parameters have two essential impacts on the effective connectivity of neural networks: First, that the populations advancing in phase to others do not necessarily play the role of the information source; and second, that the amount and direction of information transfer dependents on the oscillation frequency of the populations., The work of AP, IF and CM was partially supported by the Spanish State Research Agency, through the Severo Ochoa and Maria de Maeztu Program for Centers and Units of Excellence in R&D (MDM-2017-0711) and the MINECO (Spain) through project TEC2016-80063-C3 (AEI/FEDER, UE). IF and CM acknowledge support from the Ministerio de Ciencia e Innovación through projects PID2019-111537GB-C21/AEI/10.13039/501100011033 and PID2019-111537GB-C22/AEI/10.13039/501100011033, respectively.

The olfactory system faces the difficult task of identifying an enormous variety of odors independent of their intensity. Primacy coding, where the odor identity is encoded by the receptor types that respond earliest, might provide a compact and informative representation that can be interpreted efficiently by the brain. In this paper, we analyze the information transmitted by a simple model of primacy coding using numerical simulations and statistical descriptions. We show that the encoded information depends strongly on the number of receptor types included in the primacy representation, but only weakly on the size of the receptor repertoire. The representation is independent of the odor intensity and the transmitted information is useful to perform typical olfactory tasks with close to experimentally measured performance. Interestingly, we find situations in which a smaller receptor repertoire is advantageous for discriminating odors. The model also suggests that overly sensitive receptor types could dominate the entire response and make the whole array useless, which allows us to predict how receptor arrays need to adapt to stay useful during environmental changes. Taken together, we show that primacy coding is more useful than simple binary and normalized coding, essentially because the sparsity of the odor representations is independent of the odor statistics, in contrast to the alternatives. Primacy coding thus provides an efficient odor representation that is independent of the odor intensity and might thus help to identify odors in the olfactory cortex., Author summary Humans can identify odors independent of their intensity. Experimental data suggest that this is accomplished by representing the odor identity by the earliest responding receptor types. Using theoretical modeling, we here show that such a primacy code outperforms alternative encodings and allows discriminating odors with close to experimentally measured performance. This performance depends strongly on the number of receptors considered in the primacy code, but the receptor repertoire size is less important. The model also suggests a strong evolutionary pressure on the receptor sensitivities, which could explain observed receptor copy number adaptations. By predicting psycho-physical experiments, the model will thus contribute to our understanding of the olfactory system.

Neuroscience models come in a wide range of scales and specificity, from mean-field rate models to large-scale networks of spiking neurons. There are potential trade-offs between simplicity and realism, versatility and computational speed. This paper is about large-scale cortical network models, and the question we address is one of scalability: would scaling down cell density impact a network’s ability to reproduce cortical dynamics and function? We investigated this problem using a previously constructed realistic model of the monkey visual cortex that is true to size. Reducing cell density gradually up to 50-fold, we studied changes in model behavior. Size reduction without parameter adjustment was catastrophic. Surprisingly, relatively minor compensation in synaptic weights guided by a theoretical algorithm restored mean firing rates and basic function such as orientation selectivity to models 10-20 times smaller than the real cortex. Not all was normal in the reduced model cortices: intracellular dynamics acquired a character different from that of real neurons, and while the ability to relay feedforward inputs remained intact, reduced models showed signs of deficiency in functions that required dynamical interaction among cortical neurons. These findings are not confined to models of the visual cortex, and modelers should be aware of potential issues that accompany size reduction. Broader implications of this study include the importance of homeostatic maintenance of firing rates, and the functional consequences of feedforward versus recurrent dynamics, ideas that may shed light on other species and on systems suffering cell loss., Author summary With the vast numbers of neurons in the cerebral cortex, models in neuroscience are, for practical reasons, often orders of magnitude smaller than the actual network. We examine in this article the scalability of cortical networks. We find that function and dynamics in a network depend on network size. For illustration, we use a previously constructed realistic model of monkey visual cortex. Reducing the number of cells in the model, we find that small changes in synaptic weights can help maintain firing rates. However, model characteristics change fundamentally in the reduced models. Neurons have abnormal intracellular dynamics. The model becomes dominated by feedforward inputs and is less capable of functions requiring neuronal interaction. Modelers need to be aware of the potential issues with reduced cortical network models.

Circadian rhythms are generated by interlocked transcriptional-translational negative feedback loops (TTFLs), the molecular process implemented within a cell. The contributions, weighting and balancing between the multiple feedback loops remain debated. Dissociated, free-running dynamics in the expression of distinct clock genes has been described in recent experimental studies that applied various perturbations such as slice preparations, light pulses, jet-lag, and culture medium exchange. In this paper, we provide evidence that this “presumably transient” dissociation of circadian gene expression oscillations may occur at the single-cell level. Conceptual and detailed mechanistic mathematical modeling suggests that such dissociation is due to a weak interaction between multiple feedback loops present within a single cell. The dissociable loops provide insights into underlying mechanisms and general design principles of the molecular circadian clock., Author summary Circadian clocks are endogenous pacemakers that generate gene expression oscillations with a period of approximately 24h. They enable an organism to anticipate daily changes in light and temperature and allow to align physiological properties to the most beneficial time around the solar day. The suprachiasmatic nucleus (SCN) of the hypothalamus is the master circadian pacemaker in mammals that coordinates peripheral clocks throughout the body and even encodes seasons. Gene expression oscillations of circadian clock genes in this master pacemaker have been shown to dissociate after perturbations of the system such as light pulses and jet-lag. The underlying mechanism remains unknown. We show that this dissociation may occur even within a single cell of the pacemaker. Data-driven mathematical modeling suggests that the dissociation relies upon a weak interaction between interlocked gene-regulatory feedback loops. Differential responses of these feedback loops to light perturbations is consistent with the concept of morning and evening oscillators.

Dendritic cells (DCs) are the most effective professional antigen-presenting cell. They ferry antigen from the extremities to T cells and are essential for the initiation of an adaptive immune response. Despite interest in how DCs respond to chemical stimuli, there have been few attempts to model DC migration. In this paper, we simulate the motility of DCs by modeling the generation of forces by filopodia and a force balance on the cell. The direction of fliopodial extension is coupled to differential occupancy of cognate chemokine receptors across the cell. Our model simulates chemokinesis and chemotaxis in a variety of chemical and mechanical environments. Simulated DCs undergoing chemokinesis were measured to have a speed of 5.1 ± 0.07 μm·min-1 and a persistence time of 3.2 ± 0.46 min, consistent with experiment. Cells undergoing chemotaxis exhibited a stronger chemotactic response when exposed to lower average chemokine concentrations, also consistent with experiment. We predicted that when placed in two opposing gradients, cells will cluster in a line, which we call the “line of equistimulation;” this clustering has also been observed. We calculated the effect of varying gradient steepness on the line of equistimulation, with steeper gradients resulting in tighter clustering. Moreover, gradients are found to be most potent when cells are in a gradient of chemokine whose mean concentration is close to the binding of the Kd to the receptor, and least potent when the mean concentration is 0.1Kd. Comparing our simulations to experiment, we can give a quantitative measure of the strength of certain chemokines relative to others. Assigning the signal of CCL19 binding CCR7 a baseline strength of 1, we found CCL21 binding CCR7 had a strength of 0.28, and CXCL12 binding CXCR4 had a strength of 0.30. These differences emerge despite both chemokines having virtually the same Kd, suggesting a mechanism of signal amplification in DCs requiring further study., Author summary Dendritic cells harvest and display antigen to other immune cells, and motility is essential to their function. Dendritic cells use filopodia to pull themselves forward, and orient their filopodia based on signals received from chemokines. We developed a model of dendritic cell motion based on a force balance, in which pulling from filopodia (calculated from a well-established model of filopodial adhesion-clutch dynamics) is counterbalanced by adhesive friction, and where the angular orientation of filopodia is based on the strength of chemotactic signal. The model can explain directional motion in chemotactic gradients, as well as the effect of superposition of coincident gradients of chemokine on chemotaxis. The model also identifies conditions in which dendritic cells can be focused at a single location when confronted with countergradients of chemokine. The model makes specific prediction regarding dendritic cell motion and defines conditions for the design of extracorporeal devices for focusing the position of dendritic cells.

We use stochastic neural field theory to analyze the stimulus-dependent tuning of neural variability in ring attractor networks. We apply perturbation methods to show how the neural field equations can be reduced to a pair of stochastic nonlinear phase equations describing the stochastic wandering of spontaneously formed tuning curves or bump solutions. These equations are analyzed using a modified version of the bivariate von Mises distribution, which is well-known in the theory of circular statistics. We first consider a single ring network and derive a simple mathematical expression that accounts for the experimentally observed bimodal (or M-shaped) tuning of neural variability. We then explore the effects of inter-network coupling on stimulus-dependent variability in a pair of ring networks. These could represent populations of cells in two different layers of a cortical hypercolumn linked via vertical synaptic connections, or two different cortical hypercolumns linked by horizontal patchy connections within the same layer. We find that neural variability can be suppressed or facilitated, depending on whether the inter-network coupling is excitatory or inhibitory, and on the relative strengths and biases of the external stimuli to the two networks. These results are consistent with the general observation that increasing the mean firing rate via external stimuli or modulating drives tends to reduce neural variability., Author summary A topic of considerable current interest concerns the neural mechanisms underlying the suppression of cortical variability following the onset of a stimulus. Since trial-by-trial variability and noise correlations are known to affect the information capacity of neurons, such suppression could improve the accuracy of population codes. One of the main candidate mechanisms is the suppression of noise-induced transitions between multiple attractors, as exemplified by ring attractor networks. The latter have been used to model experimentally measured stochastic tuning curves of directionally selective middle temporal (MT) neurons. In this paper we show how the stimulus-dependent tuning of neural variability in ring attractor networks can be analyzed in terms of the stochastic wandering of spontaneously formed tuning curves or bumps in a continuum neural field model. The advantage of neural fields is that one can derive explicit mathematical expressions for the second-order statistics of neural activity, and explore how this depends on important model parameters, such as the level of noise, the strength of recurrent connections, and the input contrast.

Across the mammalian nervous system, neurotrophins control synaptic plasticity, neuromodulation, and neuronal growth. The neurotrophin Brain-Derived Neurotrophic Factor (BDNF) is known to promote structural and functional synaptic plasticity in the hippocampus, the cerebral cortex, and many other brain areas. In recent years, a wealth of data has been accumulated revealing the paramount importance of BDNF for neuronal function. BDNF signaling gives rise to multiple complex signaling pathways that mediate neuronal survival and differentiation during development, and formation of new memories. These different roles of BDNF for neuronal function have essential consequences if BDNF signaling in the brain is reduced. Thus, BDNF knock-out mice or mice that are deficient in BDNF receptor signaling via TrkB and p75 receptors show deficits in neuronal development, synaptic plasticity, and memory formation. Accordingly, BDNF signaling dysfunctions are associated with many neurological and neurodegenerative conditions including Alzheimer’s and Huntington’s disease. However, despite the widespread implications of BDNF-dependent signaling in synaptic plasticity in healthy and pathological conditions, the interplay of the involved different biochemical pathways at the synaptic level remained mostly unknown. In this paper, we investigated the role of BDNF/TrkB signaling in spike-timing dependent plasticity (STDP) in rodent hippocampus CA1 pyramidal cells, by implementing the first subcellular model of BDNF regulated, spike timing-dependent long-term potentiation (t-LTP). The model is based on previously published experimental findings on STDP and accounts for the observed magnitude, time course, stimulation pattern and BDNF-dependence of t-LTP. It allows interpreting the main experimental findings concerning specific biomolecular processes, and it can be expanded to take into account more detailed biochemical reactions. The results point out a few predictions on how to enhance LTP induction in such a way to rescue or improve cognitive functions under pathological conditions., PLoS Computational Biology, 15 (4), ISSN:1553-734X, ISSN:1553-7358

Deconvolution of heterogeneous bulk tumor samples into distinct cellular populations is an important yet challenging problem, particularly when only partial references are available. A common approach to dealing with this problem is to deconvolve the mixed signals using available references and leverage the remaining signal as a new cell component. However, as indicated in our simulation, such an approach tends to over-estimate the proportions of known cell types and fails to detect novel cell types. Here, we propose PREDE, a partial reference-based deconvolution method using an iterative non-negative matrix factorization algorithm. Our method is verified to be effective in estimating cell proportions and expression profiles of unknown cell types based on simulated datasets at a variety of parameter settings. Applying our method to TCGA tumor samples, we found that proportions of pure cancer cells better indicate different subtypes of tumor samples. We also detected several cell types for each cancer type whose proportions successfully predicted patient survival. Our method makes a significant contribution to deconvolution of heterogeneous tumor samples and could be widely applied to varieties of high throughput bulk data. PREDE is implemented in R and is freely available from GitHub (https://xiaoqizheng.github.io/PREDE)., Author summary Tumor tissues are mixtures of different cell types. Identification and quantification of constitutional cell types within tumor tissues are important tasks in cancer research. The problem can be readily solved using regression-based methods if reference signals are available. But in most clinical applications, only partial references are available, which significantly reduces the deconvolution accuracy of the existing regression-based methods. In this paper, we propose a partial-reference based deconvolution model, PREDE, integrating the non-negative matrix factorization framework with an iterative optimization strategy. We conducted comprehensive evaluations for PREDE using both simulation and real data analyses, demonstrating better performance of our method than other existing methods.

Calcium imaging has been widely used for measuring spiking activities of neurons. When using calcium imaging, we need to extract summarized information from the raw movie beforehand. Recent studies have used matrix deconvolution for this preprocessing. However, such an approach can neither directly estimate the generative mechanism of spike trains nor use stimulus information that has a strong influence on neural activities. Here, we propose a new deconvolution method for calcium imaging using marked point processes. We consider that the observed movie is generated from a probabilistic model with marked point processes as hidden variables, and we calculate the posterior of these variables using a variational inference approach. Our method can simultaneously estimate various kinds of information, such as cell shape, spike occurrence time, and tuning curve. We apply our method to simulated and experimental data to verify its performance., Author summary Calcium imaging is a promising technique that enables the observation of the activities of large neural populations as a movie. Since the measured movie is a large-scale dataset containing a considerable amount of information, we need to apply a preprocessing procedure to extract crucial information from the raw movie for the analysis that follows. Recent studies have adopted matrix decomposition to decompose the observed movie into the product of two matrices: one consisting of cell shapes and the other consisting of calcium florescent time series. This approach can estimate cell locations and activities simultaneously; however, it cannot express some aspects of neural population codes. For instance, this approach cannot incorporate other covariates that may affect the neural population activities. In this paper, we propose a new statistical model for calcium imaging movies and an estimation procedure for this model. To express the random occurrence of spikes occurring in the movie, our model adopts a marked point process, which is used to express sequences of events to which certain characteristic values are attached. Our model can estimate cell shapes, spikes, and tuning curves of cells directly without any additional preprocessing procedure, and it also improves the estimation accuracy compared to the conventional approach.

Stochastic simulation has been a powerful tool for studying the dynamics of gene regulatory networks, particularly in terms of understanding how cell-phenotype stability and fate-transitions are impacted by noisy gene expression. However, gene networks often have dynamics characterized by multiple attractors. Stochastic simulation is often inefficient for such systems, because most of the simulation time is spent waiting for rare, barrier-crossing events to occur. We present a rare-event simulation-based method for computing epigenetic landscapes and phenotype-transitions in metastable gene networks. Our computational pipeline was inspired by studies of metastability and barrier-crossing in protein folding, and provides an automated means of computing and visualizing essential stationary and dynamic information that is generally inaccessible to conventional simulation. Applied to a network model of pluripotency in Embryonic Stem Cells, our simulations revealed rare phenotypes and approximately Markovian transitions among phenotype-states, occurring with a broad range of timescales. The relative probabilities of phenotypes and the transition paths linking pluripotency and differentiation are sensitive to global kinetic parameters governing transcription factor-DNA binding kinetics. Our approach significantly expands the capability of stochastic simulation to investigate gene regulatory network dynamics, which may help guide rational cell reprogramming strategies. Our approach is also generalizable to other types of molecular networks and stochastic dynamics frameworks., Author summary Cell phenotypes are controlled by complex interactions between genes, proteins, and other molecules within a cell, along with signals from the cell’s environment. Gene regulatory networks (GRNs) describe these interactions mathematically. In principle, a GRN model can produce a map of possible cell phenotypes and phenotype-transitions, potentially informing experimental strategies for controlling cell phenotypes. Such a map could have a profound impact on many medical fields, ranging from stem cell therapies to wound healing. However, analytical solution of GRN models is virtually impossible, except for the smallest networks. Instead, time course trajectories of GRN dynamics can be simulated using specialized algorithms. However, these methods suffer from the difficulty of studying rare events, such as the spontaneous transitions between cell phenotypes that can occur in Embryonic Stem Cells or cancer cells. In this paper, we present a method to expand current stochastic simulation algorithms for the sampling of rare phenotypes and phenotype-transitions. The output of the computational pipeline is a simplified network of a few stable phenotypes, linked by potential transitions with quantified probabilities. This simplified network gives an intuitive representation of cell phenotype-transition dynamics, which could be useful for understanding how molecular processes impact cellular responses and aid interpretation of experimental data.

Visual neurons respond to static images with specific dynamics: neuronal responses sum sub-additively over time, reduce in amplitude with repeated or sustained stimuli (neuronal adaptation), and are slower at low stimulus contrast. Here, we propose a simple model that predicts these seemingly disparate response patterns observed in a diverse set of measurements–intracranial electrodes in patients, fMRI, and macaque single unit spiking. The model takes a time-varying contrast time course of a stimulus as input, and produces predicted neuronal dynamics as output. Model computation consists of linear filtering, expansive exponentiation, and a divisive gain control. The gain control signal relates to but is slower than the linear signal, and this delay is critical in giving rise to predictions matched to the observed dynamics. Our model is simpler than previously proposed related models, and fitting the model to intracranial EEG data uncovers two regularities across human visual field maps: estimated linear filters (temporal receptive fields) systematically differ across and within visual field maps, and later areas exhibit more rapid and substantial gain control. The model is further generalizable to account for dynamics of contrast-dependent spike rates in macaque V1, and amplitudes of fMRI BOLD in human V1., Author summary This paper contributes to modeling and understanding the neuronal dynamics of visual cortex in four ways. First, we proposed a model that describes stimulus-driven neuronal dynamics in a simple and intuitive way. Second, we applied the model to intracranial EEG data and found regularities of response dynamics across and within human visual field maps. Third, the model was generalizable across different ways of measuring brain activity, allowing us to potentially link the sources underlying diverse measurements. Fourth, we comprehensively summarized existing models of neuronal dynamics, and identified effective components that give rise to accurate prediction.

Seizures occur in a recurrent manner with intermittent states of interictal and ictal discharges (IIDs and IDs). The transitions to and from IDs are determined by a set of processes, including synaptic interaction and ionic dynamics. Although mathematical models of separate types of epileptic discharges have been developed, modeling the transitions between states remains a challenge. A simple generic mathematical model of seizure dynamics (Epileptor) has recently been proposed by Jirsa et al. (2014); however, it is formulated in terms of abstract variables. In this paper, a minimal population-type model of IIDs and IDs is proposed that is as simple to use as the Epileptor, but the suggested model attributes physical meaning to the variables. The model is expressed in ordinary differential equations for extracellular potassium and intracellular sodium concentrations, membrane potential, and short-term synaptic depression variables. A quadratic integrate-and-fire model driven by the population input current is used to reproduce spike trains in a representative neuron. In simulations, potassium accumulation governs the transition from the silent state to the state of an ID. Each ID is composed of clustered IID-like events. The sodium accumulates during discharge and activates the sodium-potassium pump, which terminates the ID by restoring the potassium gradient and thus polarizing the neuronal membranes. The whole-cell and cell-attached recordings of a 4-AP-based in vitro model of epilepsy confirmed the primary model assumptions and predictions. The mathematical analysis revealed that the IID-like events are large-amplitude stochastic oscillations, which in the case of ID generation are controlled by slow oscillations of ionic concentrations. The IDs originate in the conditions of elevated potassium concentrations in a bath solution via a saddle-node-on-invariant-circle-like bifurcation for a non-smooth dynamical system. By providing a minimal biophysical description of ionic dynamics and network interactions, the model may serve as a hierarchical base from a simple to more complex modeling of seizures., Author summary In pathological conditions of epilepsy, the functioning of the neural network crucially depends on the ionic concentrations inside and outside neurons. A number of factors that affect neuronal activity is large. That is why the development of a minimal model that reproduces typical seizures could structure further experimental and analytical studies of the pathological mechanisms. Here, on a base of known biophysical models, we present a simple population-type model that includes only four principal variables, the extracellular potassium concentration, the intracellular sodium concentration, the membrane potential and the synaptic resource diminishing due to short-term synaptic depression. A simple modeled neuron is used as an observer of the population activity. We validate the model assumptions with in vitro experiments. Our model reproduces ictal and interictal events, where the latter result in bursts of spikes in single neurons, and the former represent the cluster of spike bursts. Mathematical analysis reveals that the bursts are spontaneous large-amplitude oscillations, which may cluster after a saddle-node on invariant circle bifurcation in the pro-epileptic conditions. Our consideration has significant bearing in understanding pathological neuronal network dynamics.

Neuronal information processing is regulated by fast and localized fluctuations of brain states. Brain states reliably switch between distinct spatiotemporal signatures at a network scale even though they are composed of heterogeneous and variable rhythms at a cellular scale. We investigated the mechanisms of this network control in a conductance-based population model that reliably switches between active and oscillatory mean-fields. Robust control of the mean-field properties relies critically on a switchable negative intrinsic conductance at the cellular level. This conductance endows circuits with a shared cellular positive feedback that can switch population rhythms on and off at a cellular resolution. The switch is largely independent from other intrinsic neuronal properties, network size and synaptic connectivity. It is therefore compatible with the temporal variability and spatial heterogeneity induced by slower regulatory functions such as neuromodulation, synaptic plasticity and homeostasis. Strikingly, the required cellular mechanism is available in all cell types that possess T-type calcium channels but unavailable in computational models that neglect the slow kinetics of their activation., Author summary Brain information processing involves electrophysiological signals at multiple temporal and spatial timescales, from the single neuron level to whole brain areas. A fast and local control of these signals by neurochemicals called neuromodulators is essential in complex tasks such as movement initiation and attentional focus. The neuromodulators act at the cellular scale to control signals that propagate at potentially much larger scales. The present paper highlights the critical role of a cellular switch of excitability for the fast and localized control of cellular and network states. By turning ON and OFF the cellular switch, neuromodulators can robustly switch large populations between distinct network states. We stress the importance of controlling the switch at a cellular level and independently of the connectivity to allow for tunable spatiotemporal signatures of the network states.

Single-cell RNA sequencing (scRNA-seq) has been widely applied to discover new cell types by detecting sub-populations in a heterogeneous group of cells. Since scRNA-seq experiments have lower read coverage/tag counts and introduce more technical biases compared to bulk RNA-seq experiments, the limited number of sampled cells combined with the experimental biases and other dataset specific variations presents a challenge to cross-dataset analysis and discovery of relevant biological variations across multiple cell populations. In this paper, we introduce a method of variance-driven multitask clustering of single-cell RNA-seq data (scVDMC) that utilizes multiple single-cell populations from biological replicates or different samples. scVDMC clusters single cells in multiple scRNA-seq experiments of similar cell types and markers but varying expression patterns such that the scRNA-seq data are better integrated than typical pooled analyses which only increase the sample size. By controlling the variance among the cell clusters within each dataset and across all the datasets, scVDMC detects cell sub-populations in each individual experiment with shared cell-type markers but varying cluster centers among all the experiments. Applied to two real scRNA-seq datasets with several replicates and one large-scale droplet-based dataset on three patient samples, scVDMC more accurately detected cell populations and known cell markers than pooled clustering and other recently proposed scRNA-seq clustering methods. In the case study applied to in-house Recessive Dystrophic Epidermolysis Bullosa (RDEB) scRNA-seq data, scVDMC revealed several new cell types and unknown markers validated by flow cytometry. MATLAB/Octave code available at https://github.com/kuanglab/scVDMC., Author summary scRNA-seq enables detailed profiling of heterogeneous cell populations and can be used to reveal lineage relationships or discover new cell types. In the literature, there has been little effort directed towards developing computational methods for cross-population transcriptome analysis of multiple single-cell populations. The cross-cell-population clustering problem is different from the traditional clustering problem because single-cell populations can be collected from different patients, different samples of a tissue, or different experimental replicates. The accompanying biological and technical variation tends to dominate the signals for clustering the pooled single cells from the multiple populations. In this work, we have developed a multitask clustering method to address the cross-population clustering problem. The method simultaneously clusters each individual cell population and controls variance among the cell-type cluster centers within each cell population and across the cell populations. We demonstrate that our multitask clustering method significantly improves clustering accuracy and marker discovery in three public scRNA-seq datasets and also apply the method to an in-house Recessive Dystrophic Epidermolysis Bullosa (RDEB) dataset. Our results make it evident that multitask clustering is a promising new approach for cross-population analysis of scRNA-seq data.

Many diverse biological systems are described by randomly moving particles that can be captured by traps in their environment. Examples include neurotransmitters diffusing in the synaptic cleft before binding to receptors and prey roaming an environment before capture by predators. In most cases, the traps cannot capture particles continuously. Rather, each trap must wait a transitory “recharge” time after capturing a particle before additional captures. This recharge time is often overlooked. In the case of instant recharge, the average number of particles captured before they escape grows linearly in the total number of particles. In stark contrast, we prove that for any nonzero recharge time, the average number of captured particles grows at most logarithmically in the total particle number. This is a fundamental effect of recharge, as it holds under very general assumptions on particle motion and spatial domain. Furthermore, we characterize the parameter regime in which a given recharge time will dramatically affect a system, allowing researchers to easily verify if they need to account for recharge in their specific system. Finally, we consider a few examples, including a neural system in which recharge reduces neurotransmitter bindings by several orders of magnitude., Author summary Consider particles that are released into an environment (think diffusing molecules or plankton), and suppose that there are traps in the environment. How many particles will be captured by the traps before they escape? In a standard model, the number of captured particles is proportional to the initial number released. In this paper, we show that for a more realistic model of a trap (one in which traps must recharge after every capture), the number of captures is proportional to the logarithm of the initial number released. That means that if 106 particles are released, only about 6 will be captured. We prove this result mathematically, and then consider a number of applications, including neuronal synapses and ambush predators.

Pairwise maximum-entropy models have been used in neuroscience to predict the activity of neuronal populations, given only the time-averaged correlations of the neuron activities. This paper provides evidence that the pairwise model, applied to experimental recordings, would produce a bimodal distribution for the population-averaged activity, and for some population sizes the second mode would peak at high activities, that experimentally would be equivalent to 90% of the neuron population active within time-windows of few milliseconds. Several problems are connected with this bimodality: 1. The presence of the high-activity mode is unrealistic in view of observed neuronal activity and on neurobiological grounds. 2. Boltzmann learning becomes non-ergodic, hence the pairwise maximum-entropy distribution cannot be found: in fact, Boltzmann learning would produce an incorrect distribution; similarly, common variants of mean-field approximations also produce an incorrect distribution. 3. The Glauber dynamics associated with the model is unrealistically bistable and cannot be used to generate realistic surrogate data. This bimodality problem is first demonstrated for an experimental dataset from 159 neurons in the motor cortex of macaque monkey. Evidence is then provided that this problem affects typical neural recordings of population sizes of a couple of hundreds or more neurons. The cause of the bimodality problem is identified as the inability of standard maximum-entropy distributions with a uniform reference measure to model neuronal inhibition. To eliminate this problem a modified maximum-entropy model is presented, which reflects a basic effect of inhibition in the form of a simple but non-uniform reference measure. This model does not lead to unrealistic bimodalities, can be found with Boltzmann learning, and has an associated Glauber dynamics which incorporates a minimal asymmetric inhibition., Author summary Networks of interacting units are ubiquitous in various fields of biology; e.g. gene regulatory networks, neuronal networks, social structures. If a limited set of observables is accessible, maximum-entropy models provide a way to construct a statistical model for such networks, under particular assumptions. The pairwise maximum-entropy model only uses the first two moments among those observables, and can be interpreted as a network with only pairwise interactions. If correlations are on average positive, we here show that the maximum entropy distribution tends to become bimodal. In the application to neuronal activity this is a problem, because the bimodality is an artefact of the statistical model and not observed in real data. This problem could also affect other fields in biology. We here explain under which conditions bimodality arises and present a solution to the problem by introducing a collective negative feedback, corresponding to a modified maximum-entropy model. This result may point to the existence of a homeostatic mechanism active in the system that is not part of our set of observable units.

Simultaneous changes in ion concentrations, glutamate, and cell volume together with exchange of matter between cell network and vasculature are ubiquitous in numerous brain pathologies. A complete understanding of pathological conditions as well as normal brain function, therefore, hinges on elucidating the molecular and cellular pathways involved in these mostly interdependent variations. In this paper, we develop the first computational framework that combines the Hodgkin–Huxley type spiking dynamics, dynamic ion concentrations and glutamate homeostasis, neuronal and astroglial volume changes, and ion exchange with vasculature into a comprehensive model to elucidate the role of glutamate uptake in the dynamics of spreading depolarization (SD)—the electrophysiological event underlying numerous pathologies including migraine, ischemic stroke, aneurysmal subarachnoid hemorrhage, intracerebral hematoma, and trauma. We are particularly interested in investigating the role of glutamate in the duration and termination of SD caused by K+ perfusion and oxygen-glucose deprivation. Our results demonstrate that glutamate signaling plays a key role in the dynamics of SD, and that impaired glutamate uptake leads to recovery failure of neurons from SD. We confirm predictions from our model experimentally by showing that inhibiting astrocytic glutamate uptake using TFB-TBOA nearly quadruples the duration of SD in layers 2-3 of visual cortical slices from juvenile rats. The model equations are either derived purely from first physical principles of electroneutrality, osmosis, and conservation of particles or a combination of these principles and known physiological facts. Accordingly, we claim that our approach can be used as a future guide to investigate the role of glutamate, ion concentrations, and dynamics cell volume in other brain pathologies and normal brain function., Author summary Pathological conditions such as seizure, migraine, traumatic brain injury, and stroke are associated with extreme changes in ion concentrations and glutamate, cell swelling, and heavy exchange of matter between neurons, glia, and vasculature. However, current experimental tools are capable of measuring only a few of these variables, which necessitates the development of biophysically relevant models. This study provides a comprehensive computational framework derived from first physical principles and physiological facts that enables us to investigate a wide range of key variables in SD caused by K+ perfusion and oxygen-glucose deprivation. While we use SD as an example, our approach can be extended to other pathological conditions and normal brain function.

Neural circuit development requires that synapses be formed between appropriate neurons. In addition, for a hierarchical network, successful development involves a sequencing of developmental events. It has been suggested that one mechanism that helps speed up development of proper connections is an early overproduction of synapses. Using a computational model of synapse development, such as adaptive synaptogenesis, it is possible to study such overproduction and its role in speeding up development; it is also possible to study other outcomes of synapse overproduction that are seemingly new to the literature. With a fixed number of neurons, adaptive synaptogenesis can control the speed of synaptic development in two ways: by altering the rate constants of the adaptive processes or by altering the initial number of rapidly but non-selectively accrued synapses. Using either mechanism, the simulations reveal that synapse overproduction appears as an unavoidable concomitant of rapid adaptive synaptogenesis. However, the shortest development times, which always produces the greatest amount of synapse overproduction, reduce adult performance by three measures: energy use, discrimination error rates, and proportional neuron allocation. Thus, the results here lead to the hypothesis that the observed speed of neural network development represents a particular inter-generational compromise: quick development benefits parental fecundity while slow development benefits offspring fecundity., Author summary The brain processes information by building and adjusting neural connections, and as an organism learns about the environment, it incorporates associations between features of the environment into these adjustable connections. Organisms, however, are not born with these connections; synapses must develop from a connectionless set of neurons. An important aspect of normal brain development is the timely development of the brain and of hierarchically arranged brain regions. In this paper, we show that an unsupervised algorithm, called adaptive synaptogenesis, builds neurons with stable connections and allows these neurons to appropriately discriminate patterns. Here, we show that larger rates of synaptogenesis and synaptic modification speed up development of stable connections. However, such a speed-up produces larger, steady-state energy-costs. We conjecture that these synaptic modification rates have evolved as a compromise via natural selection. In this regard, the evolution-theoretic idea of parent-offspring conflict arises as a developmental consideration; that is, time to stable connectivity must balance with subsequent neural performance. Importantly, we illustrate that the widely-observed phenomenon of synapse overproduction during development can be understood as a compromise between speedy and efficient development.

Much of the research on cannabinoids (CBs) has focused on their effects at the molecular and synaptic level. However, the effects of CBs on the dynamics of neural circuits remains poorly understood. This study aims to disentangle the effects of CBs on the functional dynamics of the hippocampal Schaffer collateral synapse by using data-driven nonparametric modeling. Multi-unit activity was recorded from rats doing an working memory task in control sessions and under the influence of exogenously administered tetrahydrocannabinol (THC), the primary CB found in marijuana. It was found that THC left firing rate unaltered and only slightly reduced theta oscillations. Multivariate autoregressive models, estimated from spontaneous spiking activity, were then used to describe the dynamical transformation from CA3 to CA1. They revealed that THC served to functionally isolate CA1 from CA3 by reducing feedforward excitation and theta information flow. The functional isolation was compensated by increased feedback excitation within CA1, thus leading to unaltered firing rates. Finally, both of these effects were shown to be correlated with memory impairments in the working memory task. By elucidating the circuit mechanisms of CBs, these results help close the gap in knowledge between the cellular and behavioral effects of CBs., Author summary Research into cannabinoids (CBs) over the last several decades has found that they induce a large variety of oftentimes opposing effects on various neuronal receptors and processes. Due to this plethora of effects, disentangling how CBs influence neuronal circuits has proven challenging. This paper contributes to our understanding of the circuit level effects of CBs by using data driven modeling to examine how THC affects the input-output relationship in the Schaffer collateral synapse in the hippocampus. It was found that THC functionally isolated CA1 from CA3 by reducing feedforward excitation and theta information flow while simultaneously increasing feedback excitation within CA1. By elucidating the circuit mechanisms of CBs, these results help close the gap in knowledge between the cellular and behavioral effects of CBs.

Cells use surface receptors to estimate concentrations of external ligands. Limits on the accuracy of such estimations have been well studied for pairs of ligand and receptor species. However, the environment typically contains many ligands, which can bind to the same receptors with different affinities, resulting in cross-talk. In traditional rate models, such cross-talk prevents accurate inference of concentrations of individual ligands. In contrast, here we show that knowing the precise timing sequence of stochastic binding and unbinding events allows one receptor to provide information about multiple ligands simultaneously and with a high accuracy. We show that such high-accuracy estimation of multiple concentrations can be realized with simple structural modifications of the familiar kinetic proofreading biochemical network diagram. We give two specific examples of such modifications. We argue that structural and functional features of real cellular biochemical sensory networks in immune cells, such as feedforward and feedback loops or ligand antagonism, sometimes can be understood as solutions to the accurate multi-ligand estimation problem., Author summary Cells live in chemically complex environments with many different chemical ligands around them. Can cells estimate concentrations of more ligands than they have receptor types? In this paper, we show that, surprisingly, the answer is “yes”, and the estimation can be implemented with simple biochemical components already present in many cells. Therefore, cells may “know” a lot more about their environment and thus may be able to implement more complex and accurate response strategies than was previously thought.

The Notch pathway plays a vital role in determining whether cells in the intestinal epithelium adopt a secretory or an absorptive phenotype. Cell fate specification is coordinated via Notch’s interaction with the canonical Wnt pathway. Here, we propose a new mathematical model of the Notch and Wnt pathways, in which the Hes1 promoter acts as a hub for pathway crosstalk. Computational simulations of the model can assist in understanding how healthy intestinal tissue is maintained, and predict the likely consequences of biochemical knockouts upon cell fate selection processes. Chemical reaction network theory (CRNT) is a powerful, generalised framework which assesses the capacity of our model for monostability or multistability, by analysing properties of the underlying network structure without recourse to specific parameter values or functional forms for reaction rates. CRNT highlights the role of β-catenin in stabilising the Notch pathway and damping oscillations, demonstrating that Wnt-mediated actions on the Hes1 promoter can induce dynamic transitions in the Notch system, from multistability to monostability. Time-dependent model simulations of cell pairs reveal the stabilising influence of Wnt upon the Notch pathway, in which β-catenin- and Dsh-mediated action on the Hes1 promoter are key in shaping the subcellular dynamics. Where Notch-mediated transcription of Hes1 dominates, there is Notch oscillation and maintenance of fate flexibility; Wnt-mediated transcription of Hes1 favours bistability akin to cell fate selection. Cells could therefore regulate the proportion of Wnt- and Notch-mediated control of the Hes1 promoter to coordinate the timing of cell fate selection as they migrate through the intestinal epithelium and are subject to reduced Wnt stimuli. Furthermore, mutant cells characterised by hyperstimulation of the Wnt pathway may, through coupling with Notch, invert cell fate in neighbouring healthy cells, enabling an aberrant cell to maintain its neighbours in mitotically active states., Author summary Epithelial cells which line the intestine form finger-shaped structures called crypts; these undergo a process of renewal at the base, causing cells to migrate upwards until they die and are sloughed off into the gut. Much of our understanding of how crypts function rests upon two processes: proliferation, in which cells divide to produce ‘daughter cells’; and differentiation, in which cells become progressively more specialised as they migrate along the crypt axis and mature. Coordinated proliferation and differentiation enable each crypt to renew itself and to produce a range of specialised cell types essential to its healthy functioning. In this paper we build a mathematical model for two reaction pathways which regulate proliferation and differentiation. We use this model to explore how crosstalk between these pathways in cell pairs influences the generation of distinct cell fates in intestinal tissues. By modifying our model to represent abnormal, ‘mutant’ cells, we investigate abnormalities typical of early colorectal cancer. Computational simulation of our model identifies an important region of crosstalk in our reaction network which determines whether cells adopt the same fate as one another, or different fates. Our model may prove useful for realistic simulations of whole crypts in the future.

Across the nervous system, certain population spiking patterns are observed far more frequently than others. A hypothesis about this structure is that these collective activity patterns function as population codewords–collective modes–carrying information distinct from that of any single cell. We investigate this phenomenon in recordings of ∼150 retinal ganglion cells, the retina’s output. We develop a novel statistical model that decomposes the population response into modes; it predicts the distribution of spiking activity in the ganglion cell population with high accuracy. We found that the modes represent localized features of the visual stimulus that are distinct from the features represented by single neurons. Modes form clusters of activity states that are readily discriminated from one another. When we repeated the same visual stimulus, we found that the same mode was robustly elicited. These results suggest that retinal ganglion cells’ collective signaling is endowed with a form of error-correcting code–a principle that may hold in brain areas beyond retina., Author Summary Neurons in most parts of the nervous system represent and process information in a collective fashion, yet the nature of this collective code is poorly understood. An important constraint placed on any such collective processing comes from the fact that individual neurons’ signaling is prone to corruption by noise. The information theory and engineering literatures have studied error-correcting codes that allow individual noise-prone coding units to “check” each other, forming an overall representation that is robust to errors. In this paper, we have analyzed the population code of one of the best-studied neural systems, the retina, and found that it is structured in a manner analogous to error-correcting schemes. Indeed, we found that the complex activity patterns over ~150 retinal ganglion cells, the output neurons of the retina, could be mapped onto collective code words, and that these code words represented precise visual information while suppressing noise. In order to analyze this coding scheme, we introduced a novel quantitative model of the retinal output that predicted neural activity patterns more accurately than existing state-of-the-art approaches.