Search

Information flow within and between cells depends significantly on calcium (Ca2+) signaling dynamics. However, the biophysical mechanisms that govern emergent patterns of Ca2+ signaling dynamics at the organ level remain elusive. Recent experimental studies in developing Drosophila wing imaginal discs demonstrate the emergence of four distinct patterns of Ca2+ activity: Ca2+ spikes, intercellular Ca2+ transients, tissue-level Ca2+ waves, and a global “fluttering” state. Here, we used a combination of computational modeling and experimental approaches to identify two different populations of cells within tissues that are connected by gap junction proteins. We term these two subpopulations “initiator cells,” defined by elevated levels of Phospholipase C (PLC) activity, and “standby cells,” which exhibit baseline activity. We found that the type and strength of hormonal stimulation and extent of gap junctional communication jointly determine the predominate class of Ca2+ signaling activity. Further, single-cell Ca2+ spikes are stimulated by insulin, while intercellular Ca2+ waves depend on Gαq activity. Our computational model successfully reproduces how the dynamics of Ca2+ transients varies during organ growth. Phenotypic analysis of perturbations to Gαq and insulin signaling support an integrated model of cytoplasmic Ca2+ as a dynamic reporter of overall tissue growth. Further, we show that perturbations to Ca2+ signaling tune the final size of organs. This work provides a platform to further study how organ size regulation emerges from the crosstalk between biochemical growth signals and heterogeneous cell signaling states., Author summary Calcium (Ca2+) is a universal second messenger that regulates a myriad of cellular processes such as cell division, cell proliferation and apoptosis. Multiple patterns of Ca2+ signaling including single-cell spikes, multicellular Ca2+ transients, large-scale Ca2+ waves, and global “fluttering” have been observed in epithelial systems during organ development. Key molecular players and biophysical mechanisms involved in formation of these patterns during organ development are not well understood. In this work, we developed a generalized multicellular model of Ca2+ that captures all the key categories of Ca2+ activity as a function of key hormonal signals. Integration of model predictions and experiments reveals two subclasses of cell populations and demonstrates that Ca2+ signaling activity at the organ scale is defined by a general decrease in gap junction communication as an organ grows. Our experiments also reveal that a “goldilocks zone” of optimal Ca2+ activity is required to achieve optimal growth at the organ level.

The islets of Langerhans exist as multicellular networks that regulate blood glucose levels. The majority of cells in the islet are excitable, insulin-producing β-cells that are electrically coupled via gap junction channels. β-cells are known to display heterogeneous functionality. However, due to gap junction coupling, β-cells show coordinated [Ca2+] oscillations when stimulated with glucose, and global quiescence when unstimulated. Small subpopulations of highly functional β-cells have been suggested to control [Ca2+] dynamics across the islet. When these populations were targeted by optogenetic silencing or photoablation, [Ca2+] dynamics across the islet were largely disrupted. In this study, we investigated the theoretical basis of these experiments and how small populations can disproportionality control islet [Ca2+] dynamics. Using a multicellular islet model, we generated normal, skewed or bimodal distributions of β-cell heterogeneity. We examined how islet [Ca2+] dynamics were disrupted when cells were targeted via hyperpolarization or populations were removed; to mimic optogenetic silencing or photoablation, respectively. Targeted cell populations were chosen based on characteristics linked to functional subpopulation, including metabolic rate of glucose oxidation or [Ca2+] oscillation frequency. Islets were susceptible to marked suppression of [Ca2+] when ~10% of cells with high metabolic activity were hyperpolarized; where hyperpolarizing cells with normal metabolic activity had little effect. However, when highly metabolic cells were removed from the model, [Ca2+] oscillations remained. Similarly, when ~10% of cells with either the highest frequency or earliest elevations in [Ca2+] were removed from the islet, the [Ca2+] oscillation frequency remained largely unchanged. Overall, these results indicate small populations of β-cells with either increased metabolic activity or increased frequency are unable to disproportionately control islet-wide [Ca2+] via gap junction coupling. Therefore, we need to reconsider the physiological basis for such small β-cell populations or the mechanism by which they may be acting to control normal islet function., Author summary Many biological systems can be studied using network theory. How heterogeneous cell subpopulations come together to create complex multicellular behavior is of great value in understanding function and dysfunction in tissues. The pancreatic islet of Langerhans is a highly coupled structure that is important for maintaining blood glucose homeostasis. β-cell electrical activity is coordinated via gap junction communication. The function of the insulin-producing β-cell within the islet is disrupted in diabetes. As such, to understand the causes of islet dysfunction we need to understand how different cells within the islet contribute to its overall function via gap junction coupling. Using a computational model of β-cell electrophysiology, we investigated how small highly functional β-cell populations within the islet contribute to its function. We found that when small populations with greater functionality were introduced into the islet, they displayed signatures of this enhanced functionality. However, when these cells were removed, the islet, retained near-normal function. Thus, in a highly coupled system, such as an islet, the heterogeneity of cells allows small subpopulations to be dispensable, and thus their absence is unable to disrupt the larger cellular network. These findings can be applied to other electrical systems that have heterogeneous cell populations.

Mutations are known to cause perturbations in essential functional features of integral membrane proteins, including ion channels. Even restricted or point mutations can result in substantially changed properties of ion currents. The additive effect of these alterations for a specific ion channel can result in significantly changed properties of the action potential (AP). Both AP shortening and AP prolongation can result from known mutations, and the consequences can be life-threatening. Here, we present a computational method for identifying new drugs utilizing combinations of existing drugs. Based on the knowledge of theoretical effects of existing drugs on individual ion currents, our aim is to compute optimal combinations that can ‘repair’ the mutant AP waveforms so that the baseline AP-properties are restored. More specifically, we compute optimal, combined, drug concentrations such that the waveforms of the transmembrane potential and the cytosolic calcium concentration of the mutant cardiomyocytes (CMs) becomes as similar as possible to their wild type counterparts after the drug has been applied. In order to demonstrate the utility of this method, we address the question of computing an optimal drug for the short QT syndrome type 1 (SQT1). For the SQT1 mutation N588K, there are available data sets that describe the effect of various drugs on the mutated K+ channel. These published findings are the basis for our computational analysis which can identify optimal compounds in the sense that the AP of the mutant CMs resembles essential biomarkers of the wild type CMs. Using recently developed insights regarding electrophysiological properties among myocytes from different species, we compute optimal drug combinations for hiPSC-CMs, rabbit ventricular CMs and adult human ventricular CMs with the SQT1 mutation. Since the ‘composition’ of ion channels that form the AP is different for the three types of myocytes under consideration, so is the composition of the optimal drug., Author summary Poly-pharmacology (using multiple drugs to treat disease) has been proposed for improving cardiac anti-arrhythmic therapy for at least two decades. However, the specific arrhythmia contexts in which polytherapy is likely to be both safe and effective have remained elusive. Type 1 short QT syndrome (SQT1) is a rare form of cardiac arrhythmia that results from mutations to the human Ether-á-go-go Related Gene (hERG) potassium channel. Functionally, these mutations are remarkably consistent in that they permit the channel to open earlier during each heart beat. While hundreds of compounds are known to inhibit hERG channels, the specific effect of SQT1 mutations that allows for early channel opening also limits the ability of most of those compounds to correct SQT1 dysfunction. Here, we have applied a suite of ventricular cardiomyocyte computational models to ask whether polytherapy may offer a more effective therapeutic strategy in SQT1, and if so, what the likely characteristics of that strategy are. Our analyses suggest that simultaneous induction of late sodium current and partial hERG blockade offers a promising strategy. While no activators of late sodium current have been clinically approved, several experimental compounds are available and may provide a basis for interrogating this strategy. The method presented here can be used to compute optimal drug combinations provided that the effect of each drug on every relevant ion channel is known.

Biomineralization is the process by which organisms use minerals to harden their tissues and provide them with physical support. Biomineralizing cells concentrate the mineral in vesicles that they secret into a dedicated compartment where crystallization occurs. The dynamics of vesicle motion and the molecular mechanisms that control it, are not well understood. Sea urchin larval skeletogenesis provides an excellent platform for investigating the kinetics of mineral-bearing vesicles. Here we used lattice light-sheet microscopy to study the three-dimensional (3D) dynamics of calcium-bearing vesicles in the cells of normal sea urchin embryos and of embryos where skeletogenesis is blocked through the inhibition of Vascular Endothelial Growth Factor Receptor (VEGFR). We developed computational tools for displaying 3D-volumetric movies and for automatically quantifying vesicle dynamics. Our findings imply that calcium vesicles perform an active diffusion motion in both, calcifying (skeletogenic) and non-calcifying (ectodermal) cells of the embryo. The diffusion coefficient and vesicle speed are larger in the mesenchymal skeletogenic cells compared to the epithelial ectodermal cells. These differences are possibly due to the distinct mechanical properties of the two tissues, demonstrated by the enhanced f-actin accumulation and myosinII activity in the ectodermal cells compared to the skeletogenic cells. Vesicle motion is not directed toward the biomineralization compartment, but the vesicles slow down when they approach it, and probably bind for mineral deposition. VEGFR inhibition leads to an increase of vesicle volume but hardly changes vesicle kinetics and doesn’t affect f-actin accumulation and myosinII activity. Thus, calcium vesicles perform an active diffusion motion in the cells of the sea urchin embryo, with diffusion length and speed that inversely correlate with the strength of the actomyosin network. Overall, our studies provide an unprecedented view of calcium vesicle 3D-dynamics and point toward cytoskeleton remodeling as an important effector of the motion of mineral-bearing vesicles., Author summary Biomineralization is a widespread, fundamental process by which organisms use minerals to harden their tissues. Mineral-bearing vesicles were observed in biomineralizing cells and play an essential role in biomineralization, yet little is known about their three-dimensional (3D) dynamics. Here we quantify 3D-vesicle-dynamics during calcite skeleton formation in sea urchin larvae, using lattice-light-sheet microscopy. We discover that calcium vesicles perform a diffusive motion in both calcifying and non-calcifying cells of the embryo. The diffusion coefficient and vesicle speed are higher in the mesenchymal skeletogenic cells compared to the epithelial ectodermal cells. This difference is possibly due to the higher rigidity of the ectodermal cells as demonstrated by the enhanced signal of f-actin and myosinII activity in these cells compared to the skeletogenic cells. The motion of the vesicles in the skeletogenic cells, is not directed toward the biomineralization compartment but the vesicles slow down near it, possibly to deposit their content. Blocking skeletogenesis through the inhibition of Vascular Endothelial Growth Factor Receptor (VEGFR), increases vesicle volume but doesn’t change the diffusion mode and the cytoskeleton markers in the cells. Our studies reveal the active diffusive motion of mineral bearing vesicles that is apparently defined by the mechanical properties of the cells.

Departures of normal blood flow and metabolite distribution from the cerebral microvasculature into neuronal tissue have been implicated with age-related neurodegeneration. Mathematical models informed by spatially and temporally distributed neuroimage data are becoming instrumental for reconstructing a coherent picture of normal and pathological oxygen delivery throughout the brain. Unfortunately, current mathematical models of cerebral blood flow and oxygen exchange become excessively large in size. They further suffer from boundary effects due to incomplete or physiologically inaccurate computational domains, numerical instabilities due to enormous length scale differences, and convergence problems associated with condition number deterioration at fine mesh resolutions. Our proposed simple finite volume discretization scheme for blood and oxygen microperfusion simulations does not require expensive mesh generation leading to the critical benefit that it drastically reduces matrix size and bandwidth of the coupled oxygen transfer problem. The compact problem formulation yields rapid and stable convergence. Moreover, boundary effects can effectively be suppressed by generating very large replica of the cortical microcirculation in silico using an image-based cerebrovascular network synthesis algorithm, so that boundaries of the perfusion simulations are far removed from the regions of interest. Massive simulations over sizeable portions of the cortex with feature resolution down to the micron scale become tractable with even modest computer resources. The feasibility and accuracy of the novel method is demonstrated and validated with in vivo oxygen perfusion data in cohorts of young and aged mice. Our oxygen exchange simulations quantify steep gradients near penetrating blood vessels and point towards pathological changes that might cause neurodegeneration in aged brains. This research aims to explain mechanistic interactions between anatomical structures and how they might change in diseases or with age. Rigorous quantification of age-related changes is of significant interest because it might aide in the search for imaging biomarkers for dementia and Alzheimer’s disease., Author summary Brain function critically depends on the maintenance of physiological blood supply and metabolism in the cortex. Disturbances to adequate perfusion have been linked to age-related neurodegeneration. However, the precise correlation between age-related hemodynamic changes and the resulting decline in oxygen delivery is not well understood and has not been quantified. Therefore, we introduce a new compact, and therefore highly scalable, computational method for predicting the physiological relationship between hemodynamics and cortical oxygen perfusion for large sections of the cortical microcirculation. We demonstrate the novel mesh generation-free (MGF), multi-scale simulation approach through realistic in vivo case studies of cortical microperfusion in the mouse brain. We further validate mechanistic correlations and a quantitative relationship between blood flow and brain oxygenation using experimental data from cohorts of young, middle aged and old mouse brains. Our computational approach overcomes size and performance limitations of previous unstructured meshing techniques to enable the prediction of oxygen tension with a spatial resolution of least two orders of magnitude higher than previously possible. Our simulation results support the hypothesis that structural changes in the microvasculature induce hypoxic pockets in the aged brain that are absent in the healthy, young mouse.

Extracellular recording is an accessible technique used in animals and humans to study the brain physiology and pathology. As the number of recording channels and their density grows it is natural to ask how much improvement the additional channels bring in and how we can optimally use the new capabilities for monitoring the brain. Here we show that for any given distribution of electrodes we can establish exactly what information about current sources in the brain can be recovered and what information is strictly unobservable. We demonstrate this in the general setting of previously proposed kernel Current Source Density method and illustrate it with simplified examples as well as using evoked potentials from the barrel cortex obtained with a Neuropixels probe and with compatible model data. We show that with conceptual separation of the estimation space from experimental setup one can recover sources not accessible to standard methods., Author summary Every set of measurements is a window into reality rendering its incomplete or distorted picture. It is often difficult to relate the obtained representation of the world to underlying ground truth. Here we show, for brain electrophysiology, for arbitrary experimental setup (distribution of electrodes), and arbitrary analytical setup (function space of current source densities), that one can identify distributions of current sources which can be recovered precisely, and those which are invisible in the system. This shows what is and what is not observable in the studied system for a given setup, allows to improve the analysis results by modifying analytical setup, and facilitates interpretation of the measured sets of LFP, ECoG and EEG recordings. In particular we show that with conceptual separation of the estimation space from experimental setup one can recover source distributions not accessible to standard methods.

Lung cancer is still a leading cause of death worldwide. In recent years, knowledge has been obtained of the mechanisms modulating ion channel kinetics and thus of cell bioelectric properties, which is promising for oncological biomarkers and targets. The complex interplay of channel expression and its consequences on malignant processes, however, is still insufficiently understood. We here introduce the first approach of an in-silico whole-cell ion current model of a cancer cell, in particular of the A549 human lung adenocarcinoma, including the main functionally expressed ion channels in the plasma membrane as so far known. This hidden Markov-based model represents the electrophysiology behind proliferation of the A549 cell, describing its rhythmic oscillation of the membrane potential able to trigger the transition between cell cycle phases, and it predicts membrane potential changes over the cell cycle provoked by targeted ion channel modulation. This first A549 in-silico cell model opens up a deeper insight and understanding of possible ion channel interactions in tumor development and progression, and is a valuable tool for simulating altered ion channel function in lung cancer electrophysiology., Author summary Advances in the understanding of functional alterations at genetic, epigenetic or protein expression and the expanding knowledge in mechanisms modulating ion channel kinetics and thus the cells’ bioelectric properties have arisen as promising cancer biomarkers and oncological targets. Our hidden Markov-based in-silico cell model represents the electrophysiology behind proliferation of the A549 cell line, explaining the cell’s rhythmic oscillation from hyperpolarized to depolarized states of the membrane potential, able to trigger the transition between cell cycle phases. The model enables the prediction of membrane potential changes over the cell cycle provoked by targeted modulation of specific ion channels, leading to cell cycle promotion or interruption. We are encouraged that the availability of this first cancer cell model will provide profound insight into possible roles and interactions of ion channels in tumor development and progression, and may aid in the testing of research hypotheses in lung cancer electrophysiology.

Heart failure (HF) is associated with an increased propensity for atrial fibrillation (AF), causing higher mortality than AF or HF alone. It is hypothesized that HF-induced remodelling of atrial cellular and tissue properties promotes the genesis of atrial action potential (AP) alternans and conduction alternans that perpetuate AF. However, the mechanism underlying the increased susceptibility to atrial alternans in HF remains incompletely elucidated. In this study, we investigated the effects of how HF-induced atrial cellular electrophysiological (with prolonged AP duration) and tissue structural (reduced cell-to-cell coupling caused by atrial fibrosis) remodelling can have an effect on the generation of atrial AP alternans and their conduction at the cellular and one-dimensional (1D) tissue levels. Simulation results showed that HF-induced atrial electrical remodelling prolonged AP duration, which was accompanied by an increased sarcoplasmic reticulum (SR) Ca2+ content and Ca2+ transient amplitude. Further analysis demonstrated that HF-induced atrial electrical remodelling increased susceptibility to atrial alternans mainly due to the increased sarcoplasmic reticulum Ca2+-ATPase (SERCA) Ca2+ reuptake, modulated by increased phospholamban (PLB) phosphorylation, and the decreased transient outward K+ current (Ito). The underlying mechanism has been suggested that the increased SR Ca2+ content and prolonged AP did not fully recover to their previous levels at the end of diastole, resulting in a smaller SR Ca2+ release and AP in the next beat. These produced Ca2+ transient alternans and AP alternans, and further caused AP alternans and Ca2+ transient alternans through Ca2+→AP coupling and AP→Ca2+ coupling, respectively. Simulation of a 1D tissue model showed that the combined action of HF-induced ion channel remodelling and a decrease in cell-to-cell coupling due to fibrosis increased the heart tissue’s susceptibility to the formation of spatially discordant alternans, resulting in an increased functional AP propagation dispersion, which is pro-arrhythmic. These findings provide insights into how HF promotes atrial arrhythmia in association with atrial alternans., Author summary Atrial Fibrillation (AF) is the most common arrhythmia in adults, especially in the elderly, with the increased incidence of stroke being a major complication that increases morbidity and mortality. The occurrence of AF is often accompanied by heart failure (HF). AF and HF are also known to have the bidirectional relationship that AF worsens HF and HF promotes AF. HF can induce atrial remodelling, including electrical remodelling, atrial fibrosis, stretch and dilatation, and oxidative stress, in which many factors are associated with arrhythmogenic atrial alternans. HF-induced atrial remodelling varies during various stages and complications of HF, but possible mechanisms underlying their pro-susceptibility to alternans have not been completely elucidated. In this study, we investigated the effects of HF-induced atrial remodelling with prolonged action potential duration (APD) and decreased cell-to-cell coupling on susceptibility to atrial alternans. Simulation results showed that HF-induced an increase in sarcoplasmic reticulum Ca2+-ATPase (SERCA) Ca2+ reuptake caused by increased phospholamban phosphorylation and a decrease in transient outward K+ current played significant roles in the genesis of Ca2+ transient alternans and action potential alternans at the single-cell level. The HF-induced decline of cell-to-cell coupling and APD prolongation promoted the genesis of spatially discordant alternans in atrial tissue. This provides insights into how HF facilitates atrial arrhythmia in relation to atrial alternans.

The renal outer medullary potassium (ROMK) channel is essential for potassium transport in the kidney, and its dysfunction is associated with a salt-wasting disorder known as Bartter syndrome. Despite its physiological significance, we lack a mechanistic understanding of the molecular defects in ROMK underlying most Bartter syndrome-associated mutations. To this end, we employed a ROMK-dependent yeast growth assay and tested single amino acid variants selected by a series of computational tools representative of different approaches to predict each variants’ pathogenicity. In one approach, we used in silico saturation mutagenesis, i.e. the scanning of all possible single amino acid substitutions at all sequence positions to estimate their impact on function, and then employed a new machine learning classifier known as Rhapsody. We also used two additional tools, EVmutation and Polyphen-2, which permitted us to make consensus predictions on the pathogenicity of single amino acid variants in ROMK. Experimental tests performed for selected mutants in different classes validated the vast majority of our predictions and provided insights into variants implicated in ROMK dysfunction. On a broader scope, our analysis suggests that consolidation of data from complementary computational approaches provides an improved and facile method to predict the severity of an amino acid substitution and may help accelerate the identification of disease-causing mutations in any protein., Author summary As the number of sequenced human genomes rises, a major challenge is to identify which single amino acid variations in a protein affect function and predispose individuals to disease. While predictive algorithms are available for this purpose, a comparative analysis of recently developed algorithms has not been adequately performed, nor is it clear whether combining algorithms would improve predictive power. To this end, we compared the efficacy of three publicly available algorithms and applied the results to Bartter syndrome, a human disease for which numerous poorly-characterized single amino acid variants have been identified and for which there is no cure. In silico saturation mutagenesis, i.e., the computational prediction of pathogenesis for every possible amino acid substitution, allowed us to experimentally test predictions by measuring the activity of an ion channel linked to Bartter syndrome. Based on data from blinded experiments, we discovered that Rhapsody and EVmutation successfully predicted deleterious mutations. Moreover, Rhapsody—which takes into account evolutionary as well as structural and dynamic considerations—predicted that >90% of known Bartter syndrome mutations are deleterious. Overall, our data will aid investigators who wish to test single amino acid variants in any protein and aid biomedical researchers who wish to develop hypotheses on the potential severity of genetic variants uncovered from genome databases.

This work reports a dynamical Markov state model of CLC-2 “fast” (pore) gating, based on 600 microseconds of molecular dynamics (MD) simulation. In the starting conformation of our CLC-2 model, both outer and inner channel gates are closed. The first conformational change in our dataset involves rotation of the inner-gate backbone along residues S168-G169-I170. This change is strikingly similar to that observed in the cryo-EM structure of the bovine CLC-K channel, though the volume of the intracellular (inner) region of the ion conduction pathway is further expanded in our model. From this state (inner gate open and outer gate closed), two additional states are observed, each involving a unique rotameric flip of the outer-gate residue GLUex. Both additional states involve conformational changes that orient GLUex away from the extracellular (outer) region of the ion conduction pathway. In the first additional state, the rotameric flip of GLUex results in an open, or near-open, channel pore. The equilibrium population of this state is low (∼1%), consistent with the low open probability of CLC-2 observed experimentally in the absence of a membrane potential stimulus (0 mV). In the second additional state, GLUex rotates to occlude the channel pore. This state, which has a low equilibrium population (∼1%), is only accessible when GLUex is protonated. Together, these pathways model the opening of both an inner and outer gate within the CLC-2 selectivity filter, as a function of GLUex protonation. Collectively, our findings are consistent with published experimental analyses of CLC-2 gating and provide a high-resolution structural model to guide future investigations., Author summary In contrast to sodium-, potassium-, and calcium-selective ion channels, the roles of chloride-selective ion channels in mammalian physiology have been studied much less and are not sufficiently understood, despite known associations of chloride-channel defects with a variety of pathological conditions. CLC-2 is a voltage-activated chloride channel (one of 9 human CLC homologs) with broad tissue and organ distribution. In this work, we use simulations to model the conformational dynamics of the CLC-2 chloride ion channel selectivity filter (SF), which is the part of the protein that controls whether the channel is in an ion-conducting or non-conducting state. Our analysis identifies four primary conformational states and a specific progression through these states. Our results are consistent with structural and functional data in the literature and provide a high-resolution model for guiding further studies of CLC-2. These results will inform our understanding of how CLC-2 regulates electrical activity and ion homeostasis in the many tissues where it is expressed.

Computational modeling of dopamine transmission is challenged by complex underlying mechanisms. Here we present a new computational model that (I) simultaneously regards release, diffusion and uptake of dopamine, (II) considers multiple terminal release events and (III) comprises both synaptic and volume transmission by incorporating the geometry of the synaptic cleft. We were able to validate our model in that it simulates concentration values comparable to physiological values observed in empirical studies. Further, although synaptic dopamine diffuses into extra-synaptic space, our model reflects a very localized signal occurring on the synaptic level, i.e. synaptic dopamine release is negligibly recognized by neighboring synapses. Moreover, increasing evidence suggests that cognitive performance can be predicted by signal variability of neuroimaging data (e.g. BOLD). Signal variability in target areas of dopaminergic neurons (striatum, cortex) may arise from dopamine concentration variability. On that account we compared spatio-temporal variability in a simulation mimicking normal dopamine transmission in striatum to scenarios of enhanced dopamine release and dopamine uptake inhibition. We found different variability characteristics between the three settings, which may in part account for differences in empirical observations. From a clinical perspective, differences in striatal dopaminergic signaling contribute to differential learning and reward processing, with relevant implications for addictive- and compulsive-like behavior. Specifically, dopaminergic tone is assumed to impact on phasic dopamine and hence on the integration of reward-related signals. However, in humans DA tone is classically assessed using PET, which is an indirect measure of endogenous DA availability and suffers from temporal and spatial resolution issues. We discuss how this can lead to discrepancies with observations from other methods such as microdialysis and show how computational modeling can help to refine our understanding of DA transmission., Author summary The dopaminergic system of the brain is very complex and affects various cognitive domains like memory, learning and motor control. Alterations have been observed e.g. in Parkinson’s or Huntington’s Disease, ADHD, addiction and compulsive disorders, such as pathological gambling and also in obesity. We present a new computational model that allows to simulate the process of dopamine transmission from dopaminergic neurons originated in source brain regions like the VTA to target areas such as the striatum on a synaptic and on a larger, volume-spanning level. The model can further be used for simulations of dopamine related diseases or pharmacological interventions. In general, computational modeling helps to extend our understanding, gained from empirical research, to situations where in vivo measurements are not feasible.

The ventral striatum (VS) is a central node within a distributed network that controls appetitive behavior, and neuromodulation of the VS has demonstrated therapeutic potential for appetitive disorders. Local field potential (LFP) oscillations recorded from deep brain stimulation (DBS) electrodes within the VS are a pragmatic source of neural systems-level information about appetitive behavior that could be used in responsive neuromodulation systems. Here, we recorded LFPs from the bilateral nucleus accumbens core and shell (subregions of the VS) during limited access to palatable food across varying conditions of hunger and food palatability in male rats. We used standard statistical methods (logistic regression) as well as the machine learning algorithm lasso to predict aspects of feeding behavior using VS LFPs. We were able to predict the amount of food eaten, the increase in consumption following food deprivation, and the type of food eaten. Further, we were able to predict whether the initiation of feeding was imminent up to 42.5 seconds before feeding began and classify current behavior as either feeding or not-feeding. In classifying feeding behavior, we found an optimal balance between model complexity and performance with models using 3 LFP features primarily from the alpha and high gamma frequencies. As shown here, unbiased methods can identify systems-level neural activity linked to domains of mental illness with potential application to the development and personalization of novel treatments., Author summary As neuropsychiatry begins to leverage the power of computational methods to understand disease states and to develop better therapies, it is vital that we acknowledge the trade-offs between model complexity and performance. We show that computational methods can elucidate a neural signature of feeding behavior and we show how these methods could be used to discover neural patterns related to other behaviors and reveal new potential therapeutic targets. Further, our results help to contextualize both the limitations and potential of applying computational methods to neuropsychiatry by showing how changing the data being used to train predictive models (e.g., population vs. individual data) can have a large impact on how model performance generalizes across time, internal states, and individuals.

Electrocorticography (ECoG) is becoming more prevalent due to improvements in fabrication and recording technology as well as its ease of implantation compared to intracortical electrophysiology, larger cortical coverage, and potential advantages for use in long term chronic implantation. Given the flexibility in the design of ECoG grids, which is only increasing, it remains an open question what geometry of the electrodes is optimal for an application. Conductive polymer, PEDOT:PSS, coated microelectrodes have an advantage that they can be made very small without losing low impedance. This makes them suitable for evaluating the required granularity of ECoG recording in humans and experimental animals. We used two-dimensional (2D) micro-ECoG grids to record intra-operatively in humans and during acute implantations in mouse with separation distance between neighboring electrodes (i.e., pitch) of 0.4 mm and 0.2/0.25 mm respectively. To assess the spatial properties of the signals, we used the average correlation between electrodes as a function of the pitch. In agreement with prior studies, we find a strong frequency dependence in the spatial scale of correlation. By applying independent component analysis (ICA), we find that the spatial pattern of correlation is largely due to contributions from multiple spatially extended, time-locked sources present at any given time. Our analysis indicates the presence of spatially structured activity down to the sub-millimeter spatial scale in ECoG despite the effects of volume conduction, justifying the use of dense micro-ECoG grids., Author summary Electrocorticography (ECoG) is a type of electrophysiological monitoring that uses electrodes placed directly on the exposed surface of the brain. ECoG is a promising technique for studying the brain, and EcoG signals can be used to control brain-computer interfaces. Advances have made it possible to record simultaneously with an increasing number of smaller, and more closely spaced electrodes. However, a property of electrical recording from outside the brain is that common signals appear on different electrodes at different locations, and this affects decisions about how to best distribute a limited number of electrodes to maximize the information that can be gathered. Large spacing of electrodes around one centimeter apart on the brain’s surface has proven useful for clinical and research use, but how much benefit there is to recording from more locations in a smaller area remains to be answered. We found that we can explain the commonality between the different locations as the combination of different patterns of brain activity that are present at multiple electrode locations, and that signals recorded from very closely spaced electrodes, around a millimeter or less apart, are able to identify patterns that are at this small scale.

In this paper, we present a novel Bayesian adaptive dual controller (ADC) for autonomously programming deep brain stimulation devices. We evaluated the Bayesian ADC’s performance in the context of reducing beta power in a computational model of Parkinson’s disease, in which it was tasked with finding the set of stimulation parameters which optimally reduced beta power as fast as possible. Here, the Bayesian ADC has dual goals: (a) to minimize beta power by exploiting the best parameters found so far, and (b) to explore the space to find better parameters, thus allowing for better control in the future. The Bayesian ADC is composed of two parts: an inner parameterized feedback stimulator and an outer parameter adjustment loop. The inner loop operates on a short time scale, delivering stimulus based upon the phase and power of the beta oscillation. The outer loop operates on a long time scale, observing the effects of the stimulation parameters and using Bayesian optimization to intelligently select new parameters to minimize the beta power. We show that the Bayesian ADC can efficiently optimize stimulation parameters, and is superior to other optimization algorithms. The Bayesian ADC provides a robust and general framework for tuning stimulation parameters, can be adapted to use any feedback signal, and is applicable across diseases and stimulator designs., Author summary Deep brain stimulation (DBS) is an effective therapy for treating motor symptoms of Parkinson’s disease. However, the clinical success of DBS relies on selecting stimulation parameters that both relieve symptoms while avoiding side effects. Currently, DBS devices are programmed using a laborious trial-and-error process, requiring multiple clinic visits over the course of months. As DBS leads and algorithms become more complex, it will become impossible to select optimal DBS parameters manually. There is a clear need for an intelligent, automated approach to parameter tuning. We present a novel Bayesian adaptive dual controller (ADC), which can autonomously tune stimulation parameters. It uses a feedback signal measured from the patient to quantify the efficacy of a set of stimulation parameters, and uses this information to intelligently find the parameters which work best for each individual patient. The Bayesian ADC has the potential to improve DBS efficacy and reduce clinic visits by efficiently finding the best stimulation parameters.

The propagation of epileptic seizure activity in the brain is a widespread pathophysiology that, in principle, should yield to intervention techniques guided by mathematical models of neuronal ensemble dynamics. During a seizure, neural activity will deviate from its current dynamical regime to one in which there are significant signal fluctuations. In silico treatments of neural activity are an important tool for the understanding of how the healthy brain can maintain stability, as well as of how pathology can lead to seizures. The hope is that, contained within the mathematical foundations of such treatments, there lie potential strategies for mitigating instabilities, e.g. via external stimulation. Here, we demonstrate that the dynamic causal modelling neuronal state equation generalises to a Fokker-Planck formalism if one extends the framework to model the ways in which activity propagates along the structural connections of neural systems. Using the Jacobian of this generalised state equation, we show that an initially unstable system can be rendered stable via a reduction in diffusivity–i.e., by lowering the rate at which neuronal fluctuations disperse to neighbouring regions. We show, for neural systems prone to epileptic seizures, that such a reduction in diffusivity can be achieved via external stimulation. Specifically, we show that this stimulation should be applied in such a way as to temporarily mirror the activity profile of a pathological region in its functionally connected areas. This counter-intuitive method is intended to be used pre-emptively–i.e., in order to mitigate the effects of the seizure, or ideally even prevent it from occurring in the first place. We offer proof of principle using simulations based on functional neuroimaging data collected from patients with idiopathic generalised epilepsy, in which we successfully suppress pathological activity in a distinct sub-network prior to seizure onset. Our hope is that this technique can form the basis for future real-time monitoring and intervention devices that are capable of treating epilepsy in a non-invasive manner., Author summary Epilepsy is a disease that affects over 50 million people worldwide. Current treatments include dangerous surgical procedures in which brain connections are severed, or even in which entire problem brain regions are removed. Pharmaceutical options are available, but only about one third of patients are responsive. However, even in these cases the drugs can cause such severe side effects that the patients sometimes choose to suffer seizures. We are proposing an innovative treatment of epilepsy that could be achieved by using non-invasive electrical stimulation. Specifically, we show that stimulation should be applied in such a way as to mirror the activity in a problem brain region, by targeting its neighbouring areas. This counterintuitive approach is based on a mathematical model in which this mirroring strategy is applied pre-emptively, i.e. long before the seizure has a chance to set in. The hope is that future clinical trials will be able to use this model to lessen the effect of seizures, or even prevent them from occurring in the first place.

Neurons in various regions of the brain generate spike bursts. While the number of spikes within a burst has been shown to carry information, information coding by interspike intervals (ISIs) is less well understood. In particular, a burst with k spikes has k−1 intraburst ISIs, and these k−1 ISIs could theoretically encode k−1 independent values. In this study, we demonstrate that such combinatorial coding occurs for retinal bursts. By recording ganglion cell spikes from isolated salamander retinae, we found that intraburst ISIs encode oscillatory light sequences that are much faster than the light intensity modulation encoded by the number of spikes. When a burst has three spikes, the two intraburst ISIs combinatorially encode the amplitude and phase of the oscillatory sequence. Analysis of trial-to-trial variability suggested that intraburst ISIs are regulated by two independent mechanisms responding to orthogonal oscillatory components, one of which is common to bursts with a different number of spikes. Therefore, the retina encodes multiple stimulus features by exploiting all degrees of freedom of burst spike patterns, i.e., the spike number and multiple intraburst ISIs., Author summary Neurons in various regions of the brain generate spike bursts. Bursts are typically composed of a few spikes generated within dozens of milliseconds, and individual bursts are separated by much longer periods of silence (~hundreds of milliseconds). Recent evidence indicates that the number of spikes in a burst, the interspike intervals (ISIs), and the overall duration of a burst, as well as the timing of burst onset, encode information. However, it remains unknown whether multiple ISIs within a single burst encode multiple input features. Here we demonstrate that such combinatorial ISI coding occurs for spike bursts in the retina. We recorded ganglion cell spikes from isolated salamander retinae stimulated with computer-generated movies. Visual response analyses indicated that multiple ISIs within a single burst combinatorially encode the phase and amplitude of oscillatory light sequences, which are different from the stimulus feature encoded by the spike number. The result demonstrates that the retina encodes multiple stimulus features by exploiting all degrees of freedom of burst spike patterns, i.e., the spike number and multiple intraburst ISIs. Because synaptic transmission in the visual system is highly sensitive to ISIs, the combinatorial ISI coding must have a major impact on visual information processing.

Calmodulin-dependent kinase II (CaMKII) has long been known to play an important role in learning and memory as well as long term potentiation (LTP). More recently it has been suggested that it might be involved in the time averaging of synaptic signals, which can then lead to the high precision of information stored at a single synapse. However, the role of the scaffolding molecule, neurogranin (Ng), in governing the dynamics of CaMKII is not yet fully understood. In this work, we adopt a rule-based modeling approach through the Monte Carlo method to study the effect of Ca2+ signals on the dynamics of CaMKII phosphorylation in the postsynaptic density (PSD). Calcium surges are observed in synaptic spines during an EPSP and back-propagating action potential due to the opening of NMDA receptors and voltage dependent calcium channels. Using agent-based models, we computationally investigate the dynamics of phosphorylation of CaMKII monomers and dodecameric holoenzymes. The scaffolding molecule, Ng, when present in significant concentration, limits the availability of free calmodulin (CaM), the protein which activates CaMKII in the presence of calcium. We show that Ng plays an important modulatory role in CaMKII phosphorylation following a surge of high calcium concentration. We find a non-intuitive dependence of this effect on CaM concentration that results from the different affinities of CaM for CaMKII depending on the number of calcium ions bound to the former. It has been shown previously that in the absence of phosphatase, CaMKII monomers integrate over Ca2+ signals of certain frequencies through autophosphorylation (Pepke et al, Plos Comp. Bio., 2010). We also study the effect of multiple calcium spikes on CaMKII holoenzyme autophosphorylation, and show that in the presence of phosphatase, CaMKII behaves as a leaky integrator of calcium signals, a result that has been recently observed in vivo. Our models predict that the parameters of this leaky integrator are finely tuned through the interactions of Ng, CaM, CaMKII, and PP1, providing a mechanism to precisely control the sensitivity of synapses to calcium signals. Author Summary not valid for PLOS ONE submissions., Author summary Neurons communicate with each other through synapses. The strength of a particular synapse is effectively the level of sensitivity of the postsynaptic neuron in response to firing of the presynaptic neuron. The process of changing synaptic strength is dubbed synaptic plasticity, a foundational aspect of learning and memory. In this work, we create a computational model of a calcium signaling pathway that sets off a chain reaction in CaMKII phosphorylation, eventually leading to synaptic plasticity. Computational modeling provides a unique way to tease apart and understand the non-intuitive results of interactions between the molecules involved. Our model successfully predicts the experimentally observed activation dynamics of this crucially important enzyme which is necessary for learning. These dynamics, along with other pathways, regulate the size of the synapse, which is known to be highly correlated with synaptic strength. In this work, we reveal quantitative characteristics of CaMKII activation for various stimuli, leading to important insights regarding the potential role of Neurogranin, a scaffolding protein in this pathway.

Intracellular calcium ([Ca2+]i) is a basic and ubiquitous cellular signal controlling a wide variety of biological processes. A remarkable example is the steering of sea urchin spermatozoa towards the conspecific egg by a spatially and temporally orchestrated series of [Ca2+]i spikes. Although this process has been an experimental paradigm for reproduction and sperm chemotaxis studies, the composition and regulation of the signalling network underlying the cytosolic calcium fluctuations are hitherto not fully understood. Here, we used a differential equations model of the signalling network to assess which set of channels can explain the characteristic envelope and temporal organisation of the [Ca2+]i-spike trains. The signalling network comprises an initial membrane hyperpolarisation produced by an Upstream module triggered by the egg-released chemoattractant peptide, via receptor activation, cGMP synthesis and decay. Followed by downstream modules leading to intraflagellar pH (pHi), voltage and [Ca2+]i fluctuations. The Upstream module outputs were fitted to kinetic data on cGMP activity and early membrane potential changes measured in bulk cell populations. Two candidate modules featuring voltage-dependent Ca2+-channels link these outputs to the downstream dynamics and can independently explain the typical decaying envelope and the progressive spacing of the spikes. In the first module, [Ca2+]i-spike trains require the concerted action of a classical CaV-like channel and a potassium channel, BK (Slo1), whereas the second module relies on pHi-dependent CatSper dynamics articulated with voltage-dependent neutral sodium-proton exchanger (NHE). We analysed the dynamics of these two modules alone and in mixed scenarios. We show that the [Ca2+]i dynamics observed experimentally after sustained alkalinisation can be reproduced by a model featuring the CatSper and NHE module but not by those including the pH-independent CaV and BK module or proportionate mixed scenarios. We conclude in favour of the module containing CatSper and NHE and highlight experimentally testable predictions that would corroborate this conclusion., Author summary Fertilisation in marine invertebrates, such as the sea urchin, occurs during broadcast spawning events in which males and females of co-localised species ejaculate sperm and spawn eggs synchronously. During these events, spermatozoa have to find and fertilise conspecific eggs in the midst of all the other ones, which is a remarkable navigation and mating choice achievement. Sperm cells do this by navigating towards the source of species-specific peptides released by the egg, steered by spatial and temporally orchestrated peaks in intracellular calcium concentration that trigger sudden reorientations. How these calcium spikes are regulated and timed remains elusive. Different calcium channels have been implicated by indirect experimental evidence giving rise to a complex network of putative interacting components. We gained insight into the structure and function of this network by modelling it as a set of candidate modules that could be studied separately. By using this ‘divide and conquer’ approach to the complexity of the network, we could characterise the potential dynamics of each module and confront these dynamics with specific quantitative data. Our results indicate that the channel mediating calcium signals in sea urchin sperm is likely CatSper, a calcium channel necessary for human male fertility.

The cardiac muscarinic receptor (M2R) regulates heart rate, in part, by modulating the acetylcholine (ACh) activated K+ current IK,ACh through dissociation of G-proteins, that in turn activate KACh channels. Recently, M2Rs were noted to exhibit intrinsic voltage sensitivity, i.e. their affinity for ligands varies in a voltage dependent manner. The voltage sensitivity of M2R implies that the affinity for ACh (and thus the ACh effect) varies throughout the time course of a cardiac electrical cycle. The aim of this study was to investigate the contribution of M2R voltage sensitivity to the rate and shape of the human sinus node action potentials in physiological and pathophysiological conditions. We developed a Markovian model of the IK,ACh modulation by voltage and integrated it into a computational model of human sinus node. We performed simulations with the integrated model varying ACh concentration and voltage sensitivity. Low ACh exerted a larger effect on IK,ACh at hyperpolarized versus depolarized membrane voltages. This led to a slowing of the pacemaker rate due to an attenuated slope of phase 4 depolarization with only marginal effect on action potential duration and amplitude. We also simulated the theoretical effects of genetic variants that alter the voltage sensitivity of M2R. Modest negative shifts in voltage sensitivity, predicted to increase the affinity of the receptor for ACh, slowed the rate of phase 4 depolarization and slowed heart rate, while modest positive shifts increased heart rate. These simulations support our hypothesis that altered M2R voltage sensitivity contributes to disease and provide a novel mechanistic foundation to study clinical disorders such as atrial fibrillation and inappropriate sinus tachycardia., Author summary Heart rate regulation is dependent upon a delicate interplay between parasympathetic and sympathetic nerve activity at the level of the sinus node. Acetylcholine slows the heart rate by activating the M2 muscarinic receptor (M2R) that, in turn, opens the acetylcholine-activated potassium channel (IK,ACh) to slow the firing of the sinus node. Surprisingly, the M2R is sensitive to membrane potential and undergoes conformational changes throughout the cardiac action potential that alter the affinity for acetylcholine, with secondary consequences for IK,ACh activity. Here, we investigated the contribution of M2R voltage sensitivity to the rate and shape of the human sinus node action potential in physiological and pathophysiological conditions, using a Markovian model of the IK,ACh channel integrated into a computational model of human sinus node. The computational model allowed us to assess the effects of potential genetic variants that alter specific properties of voltage sensitivity. Our results indicate that alterations in the M2R voltage sensitivity play a significant role in the physiology and pathophysiology of the human sinus node and atria. Our computational model is relevant for future studies aimed at the design and development of anti-arrhythmic agents that specifically target the unique voltage-sensitive properties of M2R.

Bacterial and human voltage-gated sodium channels (Navs) exhibit similar cation selectivity, despite their distinct EEEE and DEKA selectivity filter signature sequences. Recent high-resolution structures for bacterial Navs have allowed us to learn about ion conduction mechanisms in these simpler homo-tetrameric channels, but our understanding of the function of their mammalian counterparts remains limited. To probe these conduction mechanisms, a model of the human Nav1.2 channel has been constructed by grafting residues of its selectivity filter and external vestibular region onto the bacterial NavRh channel with atomic-resolution structure. Multi-μs fully atomistic simulations capture long time-scale ion and protein movements associated with the permeation of Na+ and K+ ions, and their differences. We observe a Na+ ion knock-on conduction mechanism facilitated by low energy multi-carboxylate/multi-Na+ complexes, akin to the bacterial channels. These complexes involve both the DEKA and vestibular EEDD rings, acting to draw multiple Na+ into the selectivity filter and promote permeation. When the DEKA ring lysine is protonated, we observe that its ammonium group is actively participating in Na+ permeation, presuming the role of another ion. It participates in the formation of a stable complex involving carboxylates that collectively bind both Na+ and the Lys ammonium group in a high-field strength site, permitting pass-by translocation of Na+. In contrast, multiple K+ ion complexes with the DEKA and EEDD rings are disfavored by up to 8.3 kcal/mol, with the K+-lysine-carboxylate complex non-existent. As a result, lysine acts as an electrostatic plug that partially blocks the flow of K+ ions, which must instead wait for isomerization of lysine downward to clear the path for K+ passage. These distinct mechanisms give us insight into the nature of ion conduction and selectivity in human Nav channels, while uncovering high field strength carboxylate binding complexes that define the more general phenomenon of Na+-selective ion transport in nature., Author summary Ion channels can rapidly and selectively conduct an ionic species, essential for the firing of neurons, where sodium and potassium channels respond to changes in membrane potential to release stores of sodium and potassium ions in succession. The ability of a protein pore to discriminate between these two nearly identical ions has remained an intriguing problem for decades. In particular, the origins of sodium selectivity have been obscured by diverse protein chemistries that exhibit sodium-selective conduction in prokaryotes and eukaryotes. Here we use multi-microsecond atomistic simulations to observe and contrast the permeation mechanisms of sodium and potassium ions in model human and bacterial sodium channels. These channels exhibit shared features of conduction, centered on the involvement of charged protein groups that form complexes with the smaller sodium ion. In the human channel model, we observe a special role for the signature lysine that allows sodium to pass, but partially blocks potassium. As sodium channels are vital to heartbeat, sensation, muscle contraction and brain activity, these insights could assist developments in improved therapeutics for disorders such as epilepsy and chronic pain, with mechanisms relevant to a range of ion transport processes in biology and materials.

The combination of Transcranial Magnetic Stimulation (TMS) with Electroencephalography (EEG) exposes the brain’s global response to localized and abrupt stimulations. However, large electric artifacts are induced in the EEG by the TMS, obscuring crucial stages of the brain’s response. Artifact removal is commonly performed by data processing techniques. However, an experimentally verified physical model for the origin and structure of the TMS-induced discharge artifacts, by which these methods can be justified or evaluated, is still lacking. We re-examine the known contribution of the skin in creating the artifacts, and outline a detailed model for the relaxation of the charge accumulated at the electrode-gel-skin interface due to the TMS pulse. We then experimentally validate implications set forth by the model. We find that the artifacts decay like a power law in time rather than the commonly assumed exponential. In fact, the skin creates a power-law decay of order 1 at each electrode, which is turned into a power law of order 2 by the reference electrode. We suggest an artifact removal method based on the model which can be applied from times after the pulse as short as 2 milliseconds onwards to expose the full EEG from the brain. The method can separate the capacitive discharge artifacts from those resulting from cranial muscle activation, demonstrating that the capacitive effect dominates at short times. Overall, our insight into the physical process allows us to accurately access TMS-evoked EEG responses that directly follow the TMS pulse, possibly opening new opportunities in TMS-EEG research., Author summary The combined use of transcranial magnetic stimulation (TMS) with electroencephalography (EEG) is rapidly becoming a useful diagnostic and research tool for measuring the brain’s global response to localized and abrupt stimulations. However, large electric artifacts that are induced in the EEG by the TMS must be removed from the EEG signal, since they are unrelated to brain activity and obscure crucial stages of the brain’s response. Here we show that a simple physical model of skin impedance explains the observed artifact decay, and enables the efficient removal of the artifact signal in an automatable manner. We thus discern the EEG response precisely from about 2 milliseconds after the TMS pulse, showing that the artifact is dominated by capacitive rather than muscular effects.

Short-term synaptic depression, caused by depletion of releasable neurotransmitter, modulates the strength of neuronal connections in a history-dependent manner. Quantifying the statistics of synaptic transmission requires stochastic models that link probabilistic neurotransmitter release with presynaptic spike-train statistics. Common approaches are to model the presynaptic spike train as either regular or a memory-less Poisson process: few analytical results are available that describe depressing synapses when the afferent spike train has more complex, temporally correlated statistics such as bursts. Here we present a series of analytical results—from vesicle release-site occupancy statistics, via neurotransmitter release, to the post-synaptic voltage mean and variance—for depressing synapses driven by correlated presynaptic spike trains. The class of presynaptic drive considered is that fully characterised by the inter-spike-interval distribution and encompasses a broad range of models used for neuronal circuit and network analyses, such as integrate-and-fire models with a complete post-spike reset and receiving sufficiently short-time correlated drive. We further demonstrate that the derived post-synaptic voltage mean and variance allow for a simple and accurate approximation of the firing rate of the post-synaptic neuron, using the exponential integrate-and-fire model as an example. These results extend the level of biological detail included in models of synaptic transmission and will allow for the incorporation of more complex and physiologically relevant firing patterns into future studies of neuronal networks., Author summary Synapses between neurons transmit signals with strengths that vary with the history of their activity, over scales from milliseconds to decades. Short-term changes in synaptic strength modulate and sculpt ongoing neuronal activity, whereas long-term changes underpin memory formation. Here we focus on changes of strength over timescales of less than a second caused by transitory depletion of the neurotransmitters that carry signals across the synapse. Neurotransmitters are stored in small vesicles that release their contents, with a certain probability, when the presynaptic neuron is active. Once a vesicle has been used it is replenished after a variable delay. There is therefore a complex interaction between the pattern of incoming signals to the synapse and the probablistic release and restock of packaged neurotransmitter. Here we extend existing models to examine how correlated synaptic activity is transmitted through synapses and affects the voltage fluctuations and firing rate of the target neuron. Our results provide a framework that will allow for the inclusion of biophysically realistic synaptic behaviour in studies of neuronal circuits.

Calmodulin (CaM) is a calcium sensing protein that regulates the function of a large number of proteins, thus playing a crucial part in many cell signaling pathways. CaM has the ability to bind more than 300 different target peptides in a Ca2+-dependent manner, mainly through the exposure of hydrophobic residues. How CaM can bind a large number of targets while retaining some selectivity is a fascinating open question. Here, we explore the mechanism of CaM selective promiscuity for selected target proteins. Analyzing enhanced sampling molecular dynamics simulations of Ca2+-bound and Ca2+-free CaM via spectral clustering has allowed us to identify distinct conformational states, characterized by interhelical angles, secondary structure determinants and the solvent exposure of specific residues. We searched for indicators of conformational selection by mapping solvent exposure of residues in these conformational states to contacts in structures of CaM/target peptide complexes. We thereby identified CaM states involved in various binding classes arranged along a depth binding gradient. Binding Ca2+ modifies the accessible hydrophobic surface of the two lobes and allows for deeper binding. Apo CaM indeed shows shallow binding involving predominantly polar and charged residues. Furthermore, binding to the C-terminal lobe of CaM appears selective and involves specific conformational states that can facilitate deep binding to target proteins, while binding to the N-terminal lobe appears to happen through a more flexible mechanism. Thus the long-ranged electrostatic interactions of the charged residues of the N-terminal lobe of CaM may initiate binding, while the short-ranged interactions of hydrophobic residues in the C-terminal lobe of CaM may account for selectivity. This work furthers our understanding of the mechanism of CaM binding and selectivity to different target proteins and paves the way towards a comprehensive model of CaM selectivity., Author summary Calmodulin is a protein involved in the regulation of a variety of cell signaling pathways. It acts by making usually calcium-insensitive proteins sensitive to changes in the calcium concentration inside the cell. Its two lobes bind calcium and allow the energetically unfavorable exposure of hydrophobic residues to the aqueous environment which can then bind target proteins. The mechanisms behind the simultaneous specificity and variation of target protein binding is yet unknown but will aid understanding of the calcium-signaling and regulation that occur in many of our cellular processes. Here, we used molecular dynamics simulations and data analysis techniques to investigate what effect calcium has on the binding modes of calmodulin. The simulations and analyses allow us to observe and differentiate specific states. One domain of calmodulin is shown to be selective with binding involving short-distance interactions between hydrophobic residues, while the other binds target proteins through a more flexible mechanism involving long-distance electrostatic interactions.

Recorded potentials in the extracellular space (ECS) of the brain is a standard measure of population activity in neural tissue. Computational models that simulate the relationship between the ECS potential and its underlying neurophysiological processes are commonly used in the interpretation of such measurements. Standard methods, such as volume-conductor theory and current-source density theory, assume that diffusion has a negligible effect on the ECS potential, at least in the range of frequencies picked up by most recording systems. This assumption remains to be verified. We here present a hybrid simulation framework that accounts for diffusive effects on the ECS potential. The framework uses (1) the NEURON simulator to compute the activity and ionic output currents from multicompartmental neuron models, and (2) the electrodiffusive Kirchhoff-Nernst-Planck framework to simulate the resulting dynamics of the potential and ion concentrations in the ECS, accounting for the effect of electrical migration as well as diffusion. Using this framework, we explore the effect that ECS diffusion has on the electrical potential surrounding a small population of 10 pyramidal neurons. The neural model was tuned so that simulations over ∼100 seconds of biological time led to shifts in ECS concentrations by a few millimolars, similar to what has been seen in experiments. By comparing simulations where ECS diffusion was absent with simulations where ECS diffusion was included, we made the following key findings: (i) ECS diffusion shifted the local potential by up to ∼0.2 mV. (ii) The power spectral density (PSD) of the diffusion-evoked potential shifts followed a 1/f2 power law. (iii) Diffusion effects dominated the PSD of the ECS potential for frequencies up to several hertz. In scenarios with large, but physiologically realistic ECS concentration gradients, diffusion was thus found to affect the ECS potential well within the frequency range picked up in experimental recordings., Author Summary When electrical potentials are measured in the extracellular space (ECS) of the brain, they are interpreted as a signature of neural signalling. The relationship between the ECS potentials and the underlying neuronal processes is often studied with the aid of computer models. The ECS potential is typically assumed not to be affected by diffusive currents in the ECS, and existing models therefore neglect diffusion. However, there may be scenarios where this assumption does not hold. Here, we present a new computational model which explicitly models ion-concentration dynamics in the ECS surrounding a neural population, and which allows us to quantify the effect that diffusive currents have on the ECS potential. Using this model, we simulate a scenario where a population of pyramidal neurons is active over a long time, and produces large, but realistic concentration gradients in the ECS. In this scenario, diffusive currents are found to influence the ECS potential at frequency components as high as ten hertz. Unlike previously believed, we thus predict that there are scenarios where recorded local field potentials (LFPs) are likely to contain signatures not only of neural activity, but also of ECS diffusion.

Hyperpolarization-activated cyclic nucleotide-modulated (HCN) channels control electrical rhythmicity in specialized brain and heart cells. We quantitatively analysed voltage-dependent activation of homotetrameric HCN2 channels and its modulation by the second messenger cAMP using global fits of hidden Markovian models to complex experimental data. We show that voltage-dependent activation is essentially governed by two separable voltage-dependent steps followed by voltage-independent opening of the pore. According to this model analysis, the binding of cAMP to the channels exerts multiple effects on the voltage-dependent gating: It stabilizes the open pore, reduces the total gating charge from ~8 to ~5, makes an additional closed state outside the activation pathway accessible and strongly accelerates the ON-gating but not the OFF-gating. Furthermore, the open channel has a much slower computed OFF-gating current than the closed channel, in both the absence and presence of cAMP. Together, these results provide detailed new insight into the voltage- and cAMP-induced activation gating of HCN channels., Author summary Hyperpolarization-activated cyclic nucleotide-modulated (HCN) channels are tetrameric voltage-controlled ion channels in the cell membrane of specialized nerve and heart cells. Their main function is to generate a so-called pacemaker current which plays a key role in the generation of electrical rhythmicity. A special messenger molecule, cAMP, synthesized within these cells at sympathetic stimulation, can bind to these channels, thereby enhancing channel opening evoked by voltage. The mechanism of this dual activation is still poorly understood. Here we quantified this duality of activation for HCN2 channels by globally fitting hidden Markovian state models to extensive sets of data. We propose that activation of this tetrameric channel requires for a full description only two voltage-dependent steps that are followed by a voltage-independent opening step of the channel pore. According to this model analysis cAMP exerts multiple effects on channel activation: It notably accelerates the charge movement of the voltage-dependent steps and reduces the number of the involved electrical charges. Furthermore, it introduces an additional closed state and stabilizes the open pore. Together, our results provide new insight into the duality of voltage- and cAMP-induced activation of HCN channels.