MCell4 with BioNetGen: A Monte Carlo simulator of rule-based reaction-diffusion systems with Python interface.

Biochemical signaling pathways in living cells are often highly organized into spatially segregated volumes, membranes, scaffolds, subcellular compartments, and organelles comprising small numbers of interacting molecules. At this level of granularity stochastic behavior dominates, well-mixed continuum approximations based on concentrations break down and a particle-based approach is more accurate and more efficient. We describe and validate a new version of the open-source MCell simulation program (MCell4), which supports generalized 3D Monte Carlo modeling of diffusion and chemical reaction of discrete molecules and macromolecular complexes in solution, on surfaces representing membranes, and combinations thereof. The main improvements in MCell4 compared to the previous versions, MCell3 and MCell3-R, include a Python interface and native BioNetGen reaction language (BNGL) support. MCell4's Python interface opens up completely new possibilities for interfacing with external simulators to allow creation of sophisticated event-driven multiscale/multiphysics simulations. The native BNGL support, implemented through a new open-source library libBNG (also introduced in this paper), provides the capability to run a given BNGL model spatially resolved in MCell4 and, with appropriate simplifying assumptions, also in the BioNetGen simulation environment, greatly accelerating and simplifying model validation and comparison.

Synaptic Information Storage Capacity Measured With Information Theory.

Variation in the strength of synapses can be quantified by measuring the anatomical properties of synapses. Quantifying precision of synaptic plasticity is fundamental to understanding information storage and retrieval in neural circuits. Synapses from the same axon onto the same dendrite have a common history of coactivation, making them ideal candidates for determining the precision of synaptic plasticity based on the similarity of their physical dimensions. Here, the precision and amount of information stored in synapse dimensions were quantified with Shannon information theory, expanding prior analysis that used signal detection theory (Bartol et al., 2015). The two methods were compared using dendritic spine head volumes in the middle of the stratum radiatum of hippocampal area CA1 as well-defined measures of synaptic strength. Information theory delineated the number of distinguishable synaptic strengths based on nonoverlapping bins of dendritic spine head volumes. Shannon entropy was applied to measure synaptic information storage capacity (SISC) and resulted in a lower bound of 4.1 bits and upper bound of 4.59 bits of information based on 24 distinguishable sizes. We further compared the distribution of distinguishable sizes and a uniform distribution using Kullback-Leibler divergence and discovered that there was a nearly uniform distribution of spine head volumes across the sizes, suggesting optimal use of the distinguishable values. Thus, SISC provides a new analytical measure that can be generalized to probe synaptic strengths and capacity for plasticity in different brain regions of different species and among animals raised in different conditions or during learning. How brain diseases and disorders affect the precision of synaptic plasticity can also be probed.

Multiscale modeling of presynaptic dynamics from molecular to mesoscale.

Chemical synapses exhibit a diverse array of internal mechanisms that affect the dynamics of transmission efficacy. Many of these processes, such as release of neurotransmitter and vesicle recycling, depend strongly on activity-dependent influx and accumulation of Ca2+. To model how each of these processes may affect the processing of information in neural circuits, and how their dysfunction may lead to disease states, requires a computationally efficient modelling framework, capable of generating accurate phenomenology without incurring a heavy computational cost per synapse. Constructing a phenomenologically realistic model requires the precise characterization of the timing and probability of neurotransmitter release. Difficulties arise in that functional forms of instantaneous release rate can be difficult to extract from noisy data without running many thousands of trials, and in biophysical synapses, facilitation of per-vesicle release probability is confounded by depletion. To overcome this, we obtained traces of free Ca2+ concentration in response to various action potential stimulus trains from a molecular MCell model of a hippocampal Schaffer collateral axon. Ca2+ sensors were placed at varying distance from a voltage-dependent calcium channel (VDCC) cluster, and Ca2+ was buffered by calbindin. Then, using the calcium traces to drive deterministic state vector models of synaptotagmin 1 and 7 (Syt-1/7), which respectively mediate synchronous and asynchronous release in excitatory hippocampal synapses, we obtained high-resolution profiles of instantaneous release rate, to which we applied functional fits. Synchronous vesicle release occurred predominantly within half a micron of the source of spike-evoked Ca2+ influx, while asynchronous release occurred more consistently at all distances. Both fast and slow mechanisms exhibited multi-exponential release rate curves, whose magnitudes decayed exponentially with distance from the Ca2+ source. Profile parameters facilitate on different time scales according to a single, general facilitation function. These functional descriptions lay the groundwork for efficient mesoscale modelling of vesicular release dynamics.

Efficient models of polymerization applied to FtsZ ring assembly in Escherichia coli.

High protein concentrations complicate modeling of polymer assembly kinetics by introducing structural complexity and a large variety of protein forms. We present a modeling approach that achieves orders of magnitude speed-up by replacing distributions of lengths and widths with their average counterparts and by introducing a hierarchical classification of species and reactions into sets. We have used this model to study FtsZ ring assembly in Escherichia coli The model's prediction of key features of the ring formation, such as time to reach the steady state, total concentration of FtsZ species in the ring, total concentration of monomers, and average dimensions of filaments and bundles, are all in agreement with the experimentally observed values. Besides validating our model against the in vivo observations, this study fills some knowledge gaps by proposing a specific structure of the ring, describing the influence of the total concentration in short and long kinetics processes, determining some characteristic mechanisms in polymer assembly regulation, and providing insights about the role of ZapA proteins, critical components for both positioning and stability of the ring.

Long-term potentiation expands information content of hippocampal dentate gyrus synapses.

An approach combining signal detection theory and precise 3D reconstructions from serial section electron microscopy (3DEM) was used to investigate synaptic plasticity and information storage capacity at medial perforant path synapses in adult hippocampal dentate gyrus in vivo. Induction of long-term potentiation (LTP) markedly increased the frequencies of both small and large spines measured 30 minutes later. This bidirectional expansion resulted in heterosynaptic counterbalancing of total synaptic area per unit length of granule cell dendrite. Control hemispheres exhibited 6.5 distinct spine sizes for 2.7 bits of storage capacity while LTP resulted in 12.9 distinct spine sizes (3.7 bits). In contrast, control hippocampal CA1 synapses exhibited 4.7 bits with much greater synaptic precision than either control or potentiated dentate gyrus synapses. Thus, synaptic plasticity altered total capacity, yet hippocampal subregions differed dramatically in their synaptic information storage capacity, reflecting their diverse functions and activation histories.

A multi-state model of the CaMKII dodecamer suggests a role for calmodulin in maintenance of autophosphorylation.

Ca2+/calmodulin-dependent protein kinase II (CaMKII) accounts for up to 2 percent of all brain protein and is essential to memory function. CaMKII activity is known to regulate dynamic shifts in the size and signaling strength of neuronal connections, a process known as synaptic plasticity. Increasingly, computational models are used to explore synaptic plasticity and the mechanisms regulating CaMKII activity. Conventional modeling approaches may exclude biophysical detail due to the impractical number of state combinations that arise when explicitly monitoring the conformational changes, ligand binding, and phosphorylation events that occur on each of the CaMKII holoenzyme's subunits. To manage the combinatorial explosion without necessitating bias or loss in biological accuracy, we use a specialized syntax in the software MCell to create a rule-based model of a twelve-subunit CaMKII holoenzyme. Here we validate the rule-based model against previous experimental measures of CaMKII activity and investigate molecular mechanisms of CaMKII regulation. Specifically, we explore how Ca2+/CaM-binding may both stabilize CaMKII subunit activation and regulate maintenance of CaMKII autophosphorylation. Noting that Ca2+/CaM and protein phosphatases bind CaMKII at nearby or overlapping sites, we compare model scenarios in which Ca2+/CaM and protein phosphatase do or do not structurally exclude each other's binding to CaMKII. Our results suggest a functional mechanism for the so-called "CaM trapping" phenomenon, wherein Ca2+/CaM may structurally exclude phosphatase binding and thereby prolong CaMKII autophosphorylation. We conclude that structural protection of autophosphorylated CaMKII by Ca2+/CaM may be an important mechanism for regulation of synaptic plasticity.

Efficient Multiscale Models of Polymer Assembly.

Protein polymerization and bundling play a central role in cell physiology. Predictive modeling of these processes remains an open challenge, especially when the proteins involved become large and their concentrations high. We present an effective kinetics model of filament formation, bundling, and depolymerization after GTP hydrolysis, which involves a relatively small number of species and reactions, and remains robust over a wide range of concentrations and timescales. We apply this general model to study assembly of FtsZ protein, a basic element in the division process of prokaryotic cells such as Escherichia coli, Bacillus subtilis, or Caulobacter crescentus. This analysis demonstrates that our model outperforms its counterparts in terms of both accuracy and computational efficiency. Because our model comprises only 17 ordinary differential equations, its computational cost is orders-of-magnitude smaller than the current alternatives consisting of up to 1000 ordinary differential equations. It also provides, to our knowledge, a new insight into the characteristics and functioning of FtsZ proteins at high concentrations. The simplicity and versatility of our model render it a powerful computational tool, which can be used either as a standalone descriptor of other biopolymers' assembly or as a component in more complete kinetic models.

Multi-state modeling of biomolecules.

Multi-state modeling of biomolecules refers to a series of techniques used to represent and compute the behavior of biological molecules or complexes that can adopt a large number of possible functional states. Biological signaling systems often rely on complexes of biological macromolecules that can undergo several functionally significant modifications that are mutually compatible. Thus, they can exist in a very large number of functionally different states. Modeling such multi-state systems poses two problems: the problem of how to describe and specify a multi-state system (the "specification problem") and the problem of how to use a computer to simulate the progress of the system over time (the "computation problem"). To address the specification problem, modelers have in recent years moved away from explicit specification of all possible states and towards rule-based formalisms that allow for implicit model specification, including the κ-calculus, BioNetGen, the Allosteric Network Compiler, and others. To tackle the computation problem, they have turned to particle-based methods that have in many cases proved more computationally efficient than population-based methods based on ordinary differential equations, partial differential equations, or the Gillespie stochastic simulation algorithm. Given current computing technology, particle-based methods are sometimes the only possible option. Particle-based simulators fall into two further categories: nonspatial simulators, such as StochSim, DYNSTOC, RuleMonkey, and the Network-Free Stochastic Simulator (NFSim), and spatial simulators, including Meredys, SRSim, and MCell. Modelers can thus choose from a variety of tools, the best choice depending on the particular problem. Development of faster and more powerful methods is ongoing, promising the ability to simulate ever more complex signaling processes in the future.

Short-term plasticity constrains spatial organization of a hippocampal presynaptic terminal.

Although the CA3-CA1 synapse is critically important for learning and memory, experimental limitations have to date prevented direct determination of the structural features that determine the response plasticity. Specifically, the local calcium influx responsible for vesicular release and short-term synaptic facilitation strongly depends on the distance between the voltage-dependent calcium channels (VDCCs) and the presynaptic active zone. Estimates for this distance range over two orders of magnitude. Here, we use a biophysically detailed computational model of the presynaptic bouton and demonstrate that available experimental data provide sufficient constraints to uniquely reconstruct the presynaptic architecture. We predict that for a typical CA3-CA1 synapse, there are ~70 VDCCs located 300 nm from the active zone. This result is surprising, because structural studies on other synapses in the hippocampus report much tighter spatial coupling. We demonstrate that the unusual structure of this synapse reflects its functional role in short-term plasticity (STP).

Computer simulations predict the impact of neuronal atrophy on the calcium dynamics in Huntington's disease.

One of the early hallmarks of Huntington's disease (HD) is neuronal cell atrophy, especially in the striatum, underlying motor dysfunction in HD. Here using a computer model, we have predicted the impact of cell shrinkage on calcium dynamics at the cellular level. Our model indicates that as cytosolic volume decreases, the amplitude of calcium transients increases and the endoplasmic reticulum (ER) becomes more leaky due to calcium-induced calcium release and a "toxic" positive feedback mechanism mediated by ryanodine receptors that greatly increases calcium release into the cytosol. The excessive calcium release from ER saturates the calcium buffering capacity of calbindin and forces further accumulation of free calcium in the cytosol and cellular compartments including mitochondria. This leads to imbalance of calcium in both cytosol and ER regions. Excessive calcium accumulation in the cytosol can damage the mitochondria resulting in metabolic dysfunction in the cell consistent with the pathology of HD. Our computational model points toward potential drug targets and can accelerate and greatly help the experimental studies of HD paving the way for treatments of patients suffering from HD.

A spatial model of autophosphorylation of CaMKII in a glutamatergic spine suggests a network-driven kinetic mechanism for bistable changes in synaptic strength.

Activation of N-methyl-D-aspartate-type glutamate receptors (NMDARs) at synapses in the CNS triggers changes in synaptic strength that underlie memory formation in response to strong synaptic stimuli. The primary target of Ca2+ flowing through NMDARs is Ca2+/calmodulin-dependent protein kinase II (CaMKII) which forms dodecameric holoenzymes that are highly concentrated at the postsynaptic site. Activation of CaMKII is necessary to trigger long-term potentiation of synaptic strength (LTP), and is prolonged by autophosphorylation of subunits within the holoenzyme. Here we use MCell4, an agent-based, stochastic, modeling platform to model CaMKII holoenzymes placed within a realistic spine geometry. We show how two mechanisms of regulation of CaMKII, 'Ca2+-calmodulin-trapping (CaM-trapping)' and dephosphorylation by protein phosphatase-1 (PP1) shape the autophosphorylation response during a repeated high-frequency stimulus. Our simulation results suggest that autophosphorylation of CaMKII does not constitute a bistable switch. Instead, prolonged but temporary, autophosphorylation of CaMKII may contribute to a biochemical-network-based 'kinetic proof-reading" mechanism that controls induction of synaptic plasticity.

Long-Term Potentiation Produces a Sustained Expansion of Synaptic Information Storage Capacity in Adult Rat Hippocampus.

Long-term potentiation (LTP) has become a standard model for investigating synaptic mechanisms of learning and memory. Increasingly, it is of interest to understand how LTP affects the synaptic information storage capacity of the targeted population of synapses. Here, structural synaptic plasticity during LTP was explored using three-dimensional reconstruction from serial section electron microscopy. Storage capacity was assessed by applying a new analytical approach, Shannon information theory, to delineate the number of functionally distinguishable synaptic strengths. LTP was induced by delta-burst stimulation of perforant pathway inputs to the middle molecular layer of hippocampal dentate granule cells in adult rats. Spine head volumes were measured as predictors of synaptic strength and compared between LTP and control hemispheres at 30 min and 2 hr after the induction of LTP. Synapses from the same axon onto the same dendrite were used to determine the precision of synaptic plasticity based on the similarity of their physical dimensions. Shannon entropy was measured by exploiting the frequency of spine heads in functionally distinguishable sizes to assess the degree to which LTP altered the number of bits of information storage. Outcomes from these analyses reveal that LTP expanded storage capacity; the distribution of spine head volumes was increased from 2 bits in controls to 3 bits at 30 min and 2.7 bits at 2 hr after the induction of LTP. Furthermore, the distribution of spine head volumes was more uniform across the increased number of functionally distinguishable sizes following LTP, thus achieving more efficient use of coding space across the population of synapses.

Anomalous diffusion of single particles in cytoplasm.

The crowded intracellular environment poses a formidable challenge to experimental and theoretical analyses of intracellular transport mechanisms. Our measurements of single-particle trajectories in cytoplasm and their random-walk interpretations elucidate two of these mechanisms: molecular diffusion in crowded environments and cytoskeletal transport along microtubules. We employed acousto-optic deflector microscopy to map out the three-dimensional trajectories of microspheres migrating in the cytosolic fraction of a cellular extract. Classical Brownian motion (BM), continuous time random walk, and fractional BM were alternatively used to represent these trajectories. The comparison of the experimental and numerical data demonstrates that cytoskeletal transport along microtubules and diffusion in the cytosolic fraction exhibit anomalous (nonFickian) behavior and posses statistically distinct signatures. Among the three random-walk models used, continuous time random walk provides the best representation of diffusion, whereas microtubular transport is accurately modeled with fractional BM.

Imbalance of ionic conductances contributes to diverse symptoms of demyelination.

Fast axonal conduction of action potentials in mammals relies on myelin insulation. Demyelination can cause slowed, blocked, desynchronized, or paradoxically excessive spiking that underlies the symptoms observed in demyelination diseases. The diversity and timing of such symptoms are poorly understood, often intermittent, and uncorrelated with disease progress. We modeled the effects of demyelination (and secondary remodeling) on intrinsic axonal excitability using Hodgkin-Huxley and reduced Morris-Lecar models. Simulations and analysis suggested a simple explanation for the breadth of symptoms and revealed that the ratio of sodium to leak conductance, g(Na)/g(L), acted as a four-way switch controlling excitability patterns that included spike failure, single spike transmission, afterdischarge, and spontaneous spiking. Failure occurred when this ratio fell below a threshold value. Afterdischarge occurred at g(Na)/g(L) just below the threshold for spontaneous spiking and required a slow inward current that allowed for two stable attractor states, one corresponding to quiescence and the other to repetitive spiking. A neuron prone to afterdischarge could function normally unless it was switched to its "pathological" attractor state; thus, although the underlying pathology may develop slowly by continuous changes in membrane conductances, a discontinuous change in axonal excitability can occur and lead to paroxysmal symptoms. We conclude that tonic and paroxysmal positive symptoms as well as negative symptoms may be a consequence of varying degrees of imbalance between g(Na) and g(L) after demyelination. The KCNK family of g(L) potassium channels may be an important target for new drugs to treat the symptoms of demyelination.

Mitochondrial morphology governs ATP production rate.

Life is based on energy conversion. In particular, in the nervous system, significant amounts of energy are needed to maintain synaptic transmission and homeostasis. To a large extent, neurons depend on oxidative phosphorylation in mitochondria to meet their high energy demand. For a comprehensive understanding of the metabolic demands in neuronal signaling, accurate models of ATP production in mitochondria are required. Here, we present a thermodynamically consistent model of ATP production in mitochondria based on previous work. The significant improvement of the model is that the reaction rate constants are set such that detailed balance is satisfied. Moreover, using thermodynamic considerations, the dependence of the reaction rate constants on membrane potential, pH, and substrate concentrations are explicitly provided. These constraints assure that the model is physically plausible. Furthermore, we explore different parameter regimes to understand in which conditions ATP production or its export are the limiting steps in making ATP available in the cytosol. The outcomes reveal that, under the conditions used in our simulations, ATP production is the limiting step and not its export. Finally, we performed spatial simulations with nine 3-D realistic mitochondrial reconstructions and linked the ATP production rate in the cytosol with morphological features of the organelles.

Morphological principles of neuronal mitochondria.

In the highly dynamic metabolic landscape of a neuron, mitochondrial membrane architectures can provide critical insight into the unique energy balance of the cell. Current theoretical calculations of functional outputs like adenosine triphosphate and heat often represent mitochondria as idealized geometries, and therefore, can miscalculate the metabolic fluxes. To analyze mitochondrial morphology in neurons of mouse cerebellum neuropil, 3D tracings of complete synaptic and axonal mitochondria were constructed using a database of serial transmission electron microscopy (TEM) tomography images and converted to watertight meshes with minimal distortion of the original microscopy volumes with a granularity of 1.64 nanometer isotropic voxels. The resulting in-silico representations were subsequently quantified by differential geometry methods in terms of the mean and Gaussian curvatures, surface areas, volumes, and membrane motifs, all of which can alter the metabolic output of the organelle. Finally, we identify structural motifs present across this population of mitochondria, which may contribute to future modeling studies of mitochondrial physiology and metabolism in neurons.

Modelling vesicular release at hippocampal synapses.

We study local calcium dynamics leading to a vesicle fusion in a stochastic, and spatially explicit, biophysical model of the CA3-CA1 presynaptic bouton. The kinetic model for vesicle release has two calcium sensors, a sensor for fast synchronous release that lasts a few tens of milliseconds and a separate sensor for slow asynchronous release that lasts a few hundred milliseconds. A wide range of data can be accounted for consistently only when a refractory period lasting a few milliseconds between releases is included. The inclusion of a second sensor for asynchronous release with a slow unbinding site, and thereby a long memory, affects short-term plasticity by facilitating release. Our simulations also reveal a third time scale of vesicle release that is correlated with the stimulus and is distinct from the fast and the slow releases. In these detailed Monte Carlo simulations all three time scales of vesicle release are insensitive to the spatial details of the synaptic ultrastructure. Furthermore, our simulations allow us to identify features of synaptic transmission that are universal and those that are modulated by structure.

Markov Chain Abstractions of Electrochemical Reaction-Diffusion in Synaptic Transmission for Neuromorphic Computing.

Progress in computational neuroscience toward understanding brain function is challenged both by the complexity of molecular-scale electrochemical interactions at the level of individual neurons and synapses and the dimensionality of network dynamics across the brain covering a vast range of spatial and temporal scales. Our work abstracts an existing highly detailed, biophysically realistic 3D reaction-diffusion model of a chemical synapse to a compact internal state space representation that maps onto parallel neuromorphic hardware for efficient emulation at a very large scale and offers near-equivalence in input-output dynamics while preserving biologically interpretable tunable parameters.

NMDAR-dependent long-term depression is associated with increased short term plasticity through autophagy mediated loss of PSD-95.

Long-term depression (LTD) of synaptic strength can take multiple forms and contribute to circuit remodeling, memory encoding or erasure. The generic term LTD encompasses various induction pathways, including activation of NMDA, mGlu or P2X receptors. However, the associated specific molecular mechanisms and effects on synaptic physiology are still unclear. We here compare how NMDAR- or P2XR-dependent LTD affect synaptic nanoscale organization and function in rodents. While both LTDs are associated with a loss and reorganization of synaptic AMPARs, only NMDAR-dependent LTD induction triggers a profound reorganization of PSD-95. This modification, which requires the autophagy machinery to remove the T19-phosphorylated form of PSD-95 from synapses, leads to an increase in AMPAR surface mobility. We demonstrate that these post-synaptic changes that occur specifically during NMDAR-dependent LTD result in an increased short-term plasticity improving neuronal responsiveness of depressed synapses. Our results establish that P2XR- and NMDAR-mediated LTD are associated to functionally distinct forms of LTD.

A comparison of deterministic and stochastic simulations of neuronal vesicle release models.

We study the calcium-induced vesicle release into the synaptic cleft using a deterministic algorithm and MCell, a Monte Carlo algorithm that tracks individual molecules. We compare the average vesicle release probability obtained using both algorithms and investigate the effect of the three main sources of noise: diffusion, sensor kinetics and fluctuations from the voltage-dependent calcium channels (VDCCs). We find that the stochastic opening kinetics of the VDCCs are the main contributors to differences in the release probability. Our results show that the deterministic calculations lead to reliable results, with an error of less than 20%, when the sensor is located at least 50 nm from the VDCCs, corresponding to microdomain signaling. For smaller distances, i.e. nanodomain signaling, the error becomes larger and a stochastic algorithm is necessary.