Deterministic Limit of Intracellular Calcium Spikes

Spatial modeling algorithms for reactions and transport in biological cells

Biological cells rely on precise spatiotemporal coordination of biochemical reactions to control their functions. Such cell signaling networks have been a common focus for mathematical models, but they remain challenging to simulate, particularly in realistic cell geometries. Here we present Spatial Modeling Algorithms for Reactions and Transport (SMART), a software package that takes in high-level user specifications about cell signaling networks and then assembles and solves the associated mathematical systems. SMART uses state-of-the-art finite element analysis, via the FEniCS Project software, to efficiently and accurately resolve cell signaling events over discretized cellular and subcellular geometries. We demonstrate its application to several different biological systems, including yes-associated protein (YAP)/PDZ-binding motif (TAZ) mechanotransduction, calcium signaling in neurons and cardiomyocytes, and ATP generation in mitochondria. Throughout, we utilize experimentally derived realistic cellular geometries represented by well-conditioned tetrahedral meshes. These scenarios demonstrate the applicability, flexibility, accuracy and efficiency of SMART across a range of temporal and spatial scales.

Ca2+ puffs underlie adhesion-triggered Ca2+ microdomains in T cells

Ca2+ signalling is pivotal in T cell activation, an essential process in adaptive immune responses. Key to this activation are Ca2+ microdomains, which are transient increases in cytosolic Ca2+ concentration occurring within narrow regions between the endoplasmic reticulum (ER) and the plasma membrane (PM), lasting a few tens of milliseconds. Adhesion Dependent Ca2+ Microdomains (ADCM) rely on store-operated Ca2+ entry (SOCE) via the ORAI/STIM system. The nanometric scale at which these microdomains form poses challenges for direct experimental observation. Following the previous work of Gil et al. [1], which introduced a three-dimensional model of the ER-PM junction, this study combines a detailed description of the Ca2+ fluxes at the junction with stochastic dynamics of a cluster of D-myo-inositol 1,4,5 trisphosphate receptors (IP3R) located in the ER surrounding the junction. Because the consideration of Ca2+ release through the IP3R calls for the simulation of a portion of the cytoplasm considerably larger than the junction, our study also investigates the spatial distribution of PMCAs, revealing their likely localization outside the ER-PM junction. Simulations indicate that Ca2+ puffs implying the opening of 2-6 IP3Rs create ADCMs by provoking local depletions of ER Ca2+ stimulating Ca2+ entry through the ORAI1 channels. Such conditions allow the reproduction of the amplitude, duration and spatial extent of the observed ADCMs. By integrating advanced computational techniques with insights from experimental studies, our approach provides valuable information on the mechanisms governing early Ca2+ signalling in T cell activation, paving the way for a deeper understanding of immune responses.

Biophysical Modeling of Synaptic Plasticity

Dendritic spines are small, bulbous compartments that function as postsynaptic sites and undergo intense biochemical and biophysical activity. The role of the myriad signaling pathways that are implicated in synaptic plasticity is well studied. A recent abundance of quantitative experimental data has made the events associated with synaptic plasticity amenable to quantitative biophysical modeling. Spines are also fascinating biophysical computational units because spine geometry, signal transduction, and mechanics work in a complex feedback loop to tune synaptic plasticity. In this sense, ideas from modeling cell motility can inspire us to develop multiscale approaches for predictive modeling of synaptic plasticity. In this article, we review the key steps in postsynaptic plasticity with a specific focus on the impact of spine geometry on signaling, cytoskeleton rearrangement, and membrane mechanics. We summarize the main experimental observations and highlight how theory and computation can aid our understanding of these complex processes.

The clock in growing hyphae and their synchronization in Neurospora crassa

Utilizing a microfluidic chip with serpentine channels, we inoculated the chip with an agar plug with Neurospora crassa mycelium and successfully captured individual hyphae in channels. For the first time, we report the presence of an autonomous clock in hyphae. Fluorescence of a mCherry reporter gene driven by a clock-controlled gene-2 promoter (ccg-2p) was measured simultaneously along hyphae every half an hour for at least 6 days. We entrained single hyphae to light over a wide range of day lengths, including 6,12, 24, and 36 h days. Hyphae tracked in individual serpentine channels were highly synchronized (K = 0.60-0.78). Furthermore, hyphae also displayed temperature compensation properties, where the oscillation period was stable over a physiological range of temperatures from 24 °C to 30 °C (Q10 = 1.00-1.10). A Clock Tube Model developed could mimic hyphal growth observed in the serpentine chip and provides a mechanism for the stable banding patterns seen in race tubes at the macroscopic scale and synchronization through molecules riding the growth wave in the device.

An integrate-and-fire approach to Ca2+ signaling. Part II: Cumulative refractoriness

Inositol 1,4,5-trisphosphate-induced Ca2+ signaling is a second messenger system used by almost all eukaryotic cells. The agonist concentration stimulating Ca2+ signals is encoded in the frequency of a Ca2+ concentration spike sequence. When a cell is stimulated, the interspike intervals (ISIs) often show a distinct transient during which they gradually increase, a system property we refer to as cumulative refractoriness. We extend a previously published stochastic model to include the Ca2+ concentration in the intracellular Ca2+ store as a slow adaptation variable. This model can reproduce both stationary and transient statistics of experimentally observed ISI sequences. We derive approximate expressions for the mean and coefficient of variation of the stationary ISIs. We also consider the response to the onset of a constant stimulus and estimate the length of the transient and the strength of the adaptation of the ISI. We show that the adaptation sets the coefficient of variation in agreement with current ideas derived from experiments. Moreover, we explain why, despite a pronounced transient behavior, ISI correlations can be weak, as often observed in experiments. Finally, we fit our model to reproduce the transient statistics of experimentally observed ISI sequences in stimulated HEK cells. The fitted model is able to qualitatively reproduce the relationship between the stationary interval correlations and the number of transient intervals, as well as the strength of the ISI adaptation. We also find positive correlations in the experimental sequence that cannot be explained by our model.

Making time and space for calcium control of neuron activity

Calcium directly controls or indirectly regulates numerous functions that are critical for neuronal network activity. Intracellular calcium concentration is tightly regulated by numerous molecular mechanisms because spatial domains and temporal dynamics (not just peak amplitude) are critical for calcium control of synaptic plasticity and ion channel activation, which in turn determine neuron spiking activity. The computational models investigating calcium control are valuable because experiments achieving high spatial and temporal resolution simultaneously are technically unfeasible. Simulations of calcium nanodomains reveal that specific calcium sources can couple to specific calcium targets, providing a mechanism to determine the direction of synaptic plasticity. Cooperativity of calcium domains opposes specificity, suggesting that the dendritic branch might be the preferred computational unit of the neuron.

T cell Ca2+ microdomains through the lens of computational modeling

Cellular Ca2+ signaling is highly organized in time and space. Locally restricted and short-lived regions of Ca2+ increase, called Ca2+ microdomains, constitute building blocks that are differentially arranged to create cellular Ca2+ signatures controlling physiological responses. Here, we focus on Ca2+ microdomains occurring in restricted cytosolic spaces between the plasma membrane and the endoplasmic reticulum, called endoplasmic reticulum-plasma membrane junctions. In T cells, these microdomains have been finely characterized. Enough quantitative data are thus available to develop detailed computational models of junctional Ca2+ dynamics. Simulations are able to predict the characteristics of Ca2+ increases at the level of single channels and in junctions of different spatial configurations, in response to various signaling molecules. Thanks to the synergy between experimental observations and computational modeling, a unified description of the molecular mechanisms that create Ca2+ microdomains in the first seconds of T cell stimulation is emerging.

An efficient reduced-lattice model of IP3R for probing Ca2+ dynamics

Numerous cellular processes are regulated by Ca2+ signals, and the endoplasmic reticulum (ER) membrane's inositol triphosphate receptor (IP3R) is critical for modulating intracellular Ca2+ dynamics. The IP3Rs are seen to be clustered in a variety of cell types. The combination of IP3Rs clustering and IP3Rs-mediated Ca2+-induced Ca2+ release results in the hierarchical organization of the Ca2+ signals, which challenges the numerical simulation given the multiple spatial and temporal scales that must be covered. The previous methods rather ignore the spatial feature of IP3Rs or fail to coordinate the conflicts between the real biological relevance and the computational cost. In this work, a general and efficient reduced-lattice model is presented for the simulation of IP3Rs-mediated multiscale Ca2+ dynamics. The model highlights biological details that make up the majority of the calcium events, including IP3Rs clustering and calcium domains, and it reduces the complexity by approximating the minor details. The model's extensibility provides fresh insights into the function of IP3Rs in producing global Ca2+ events and supports the research under more physiological circumstances. Our work contributes to a novel toolkit for modeling multiscale Ca2+ dynamics and advances knowledge of Ca2+ signals.

Modeling IP3-induced Ca2+ signaling based on its interspike interval statistics

Inositol 1,4,5-trisphosphate (IP3)-induced Ca2+ signaling is a second messenger system used by almost all eukaryotic cells. Recent research demonstrated randomness of Ca2+ signaling on all structural levels. We compile eight general properties of Ca2+ spiking common to all cell types investigated and suggest a theory of Ca2+ spiking starting from the random behavior of IP3 receptor channel clusters mediating the release of Ca2+ from the endoplasmic reticulum capturing all general properties and pathway-specific behavior. Spike generation begins after the absolute refractory period of the previous spike. According to its hierarchical spreading from initiating channel openings to cell level, we describe it as a first passage process from none to all clusters open while the cell recovers from the inhibition which terminated the previous spike. Our theory reproduces the exponential stimulation response relation of the average interspike interval Tav and its robustness properties, random spike timing with a linear moment relation between Tav and the interspike interval SD and its robustness properties, sensitive dependency of Tav on diffusion properties, and nonoscillatory local dynamics. We explain large cell variability of Tav observed in experiments by variability of channel cluster coupling by Ca2+-induced Ca2+ release, the number of clusters, and IP3 pathway component expression levels. We predict the relation between puff probability and agonist concentration and [IP3] and agonist concentration. Differences of spike behavior between cell types and stimulating agonists are explained by the different types of negative feedback terminating spikes. In summary, the hierarchical random character of spike generation explains all of the identified general properties.

Non-inositol 1,4,5-trisphosphate (IP3) receptor IP3-binding proteins

Conventionally, myo-D-inositol 1, 4,5-trisphosphate (IP₃) is thought to exert its second messenger effects through the gating of IP₃R Ca²⁺ release channels, located in Ca²⁺-storage organelles like the endoplasmic reticulum. However, there is considerable indirect evidence to support the concept that IP₃ might interact with other, non-IP₃R proteins within cells. To explore this possibility further, the Protein Data Bank was searched using the term "IP3". This resulted in the retrieval of 203 protein structures, the majority of which were members of the IP₃R/ryanodine receptor superfamily of channels. Only 49 of these structures were complexed with IP₃. These were inspected for their ability to interact with the carbon-1 phosphate of IP₃, since this is the least accessible phosphate group of its precursor, phosphatidylinositol 4,5-bisphosphate (PI(4,5)P₂). This reduced the number of structures retrieved to 35, of which 9 were IP₃Rs. The remaining 26 structures represent a diverse range of proteins, including inositol-lipid metabolizing enzymes, signal transducers, PH domain containing proteins, cytoskeletal anchor proteins, the TRPV4 ion channel, a retroviral Gag protein and fibroblast growth factor 2. Such proteins may impact on IP₃ signalling and its effects on cell-biology. This represents an area open for exploration in the field of IP₃ signalling.

Computational investigation of IP3 diffusion

Inositol 1,4,5-trisphosphate (IP₃) plays a key role in calcium signaling. After stimulation, it diffuses from the plasma membrane where it is produced to the endoplasmic reticulum where its receptors are localized. Based on in vitro measurements, IP₃ was long thought to be a global messenger characterized by a diffusion coefficient of ~ 280 μm²s^-1. However, in vivo observations revealed that this value does not match with the timing of localized Ca²⁺ increases induced by the confined release of a non-metabolizable IP₃ analog. A theoretical analysis of these data concluded that in intact cells diffusion of IP₃ is strongly hindered, leading to a 30-fold reduction of the diffusion coefficient. Here, we performed a new computational analysis of the same observations using a stochastic model of Ca²⁺ puffs. Our simulations concluded that the value of the effective IP₃ diffusion coefficient is close to 100 μm²s^-1. Such moderate reduction with respect to in vitro estimations quantitatively agrees with a buffering effect by non-fully bound inactive IP₃ receptors. The model also reveals that IP₃ spreading is not much affected by the endoplasmic reticulum, which represents an obstacle to the free displacement of molecules, but can be significantly increased in cells displaying elongated, 1-dimensional like geometries.

Characterizing spontaneous Ca2+ local transients in OPCs using computational modeling

Spontaneous Ca2+ local transients (SCaLTs) in isolated oligodendrocyte precursor cells are largely regulated by the following fluxes: store-operated Ca2+ entry (SOCE), Na+/Ca2+ exchange, Ca2+ pumping through Ca2+-ATPases, and Ca2+-induced Ca2+-release through ryanodine receptors and inositol-trisphosphate receptors. However, the relative contributions of these fluxes in mediating fast spiking and the slow baseline oscillations seen in SCaLTs remain incompletely understood. Here, we developed a stochastic spatiotemporal computational model to simulate SCaLTs in a homogeneous medium with ionic flow between the extracellular, cytoplasmic, and endoplasmic-reticulum compartments. By simulating the model and plotting both the histograms of SCaLTs obtained experimentally and from the model as well as the standard deviation of inter-SCaLT intervals against inter-SCaLT interval averages of multiple model and experimental realizations, we revealed the following: (1) SCaLTs exhibit very similar characteristics between the two data sets, (2) they are mostly random, (3) they encode information in their frequency, and (4) their slow baseline oscillations could be due to the stochastic slow clustering of inositol-trisphosphate receptors (modeled as an Ornstein-Uhlenbeck noise process). Bifurcation analysis of a deterministic temporal version of the model showed that the contribution of fluxes to SCaLTs depends on the parameter regime and that the combination of excitability, stochasticity, and mixed-mode oscillations are responsible for irregular spiking and doublets in SCaLTs. Additionally, our results demonstrated that blocking each flux reduces SCaLTs' frequency and that the reverse (forward) mode of Na+/Ca2+ exchange decreases (increases) SCaLTs. Taken together, these results provide a quantitative framework for SCaLT formation in oligodendrocyte precursor cells.

The macroscopic limit to synchronization of cellular clocks in single cells of Neurospora crassa

We determined the macroscopic limit for phase synchronization of cellular clocks in an artificial tissue created by a “big chamber” microfluidic device to be about 150,000 cells or less. The dimensions of the microfluidic chamber allowed us to calculate an upper limit on the radius of a hypothesized quorum sensing signal molecule of 13.05 nm using a diffusion approximation for signal travel within the device. The use of a second microwell microfluidic device allowed the refinement of the macroscopic limit to a cell density of 2166 cells per fixed area of the device for phase synchronization. The measurement of averages over single cell trajectories in the microwell device supported a deterministic quorum sensing model identified by ensemble methods for clock phase synchronization. A strong inference framework was used to test the communication mechanism in phase synchronization of quorum sensing versus cell-to-cell contact, suggesting support for quorum sensing. Further evidence came from showing phase synchronization was density-dependent.

Computational Approaches and Tools as Applied to the Study of Rhythms and Chaos in Biology

The temporal dynamics in biological systems displays a wide range of behaviors, from periodic oscillations, as in rhythms, bursts, long-range (fractal) correlations, chaotic dynamics up to brown and white noise. Herein, we propose a comprehensive analytical strategy for identifying, representing, and analyzing biological time series, focusing on two strongly linked dynamics: periodic (oscillatory) rhythms and chaos. Understanding the underlying temporal dynamics of a system is of fundamental importance; however, it presents methodological challenges due to intrinsic characteristics, among them the presence of noise or trends, and distinct dynamics at different time scales given by molecular, dcellular, organ, and organism levels of organization. For example, in locomotion circadian and ultradian rhythms coexist with fractal dynamics at faster time scales. We propose and describe the use of a combined approach employing different analytical methodologies to synergize their strengths and mitigate their weaknesses. Specifically, we describe advantages and caveats to consider for applying probability distribution, autocorrelation analysis, phase space reconstruction, Lyapunov exponent estimation as well as different analyses such as harmonic, namely, power spectrum; continuous wavelet transforms; synchrosqueezing transform; and wavelet coherence. Computational harmonic analysis is proposed as an analytical framework for using different types of wavelet analyses. We show that when the correct wavelet analysis is applied, the complexity in the statistical properties, including temporal scales, present in time series of signals, can be unveiled and modeled. Our chapter showcase two specific examples where an in-depth analysis of rhythms and chaos is performed: (1) locomotor and food intake rhythms over a 42-day period of mice subjected to different feeding regimes; and (2) chaotic calcium dynamics in a computational model of mitochondrial function.

High rates of calcium-free diffusion in the cytosol of living cells

Calcium (Ca²⁺) is a universal second messenger that participates in the regulation of innumerous physiological processes. The way in which local elevations of the cytosolic Ca²⁺ concentration spread in space and time is key for the versatility of the signals. Ca²⁺ diffusion in the cytosol is hindered by its interaction with proteins that act as buffers. Depending on the concentrations and the kinetics of the interactions, there is a large range of values at which Ca²⁺ diffusion can proceed. Having reliable estimates of this range, particularly of its highest end, which corresponds to the ions free diffusion, is key to understand how the signals propagate. In this work, we present the first experimental results with which the Ca²⁺-free diffusion coefficient is directly quantified in the cytosol of living cells. By means of fluorescence correlation spectroscopy experiments performed in Xenopus laevis oocytes and in cells of Saccharomyces cerevisiae, we show that the ions can freely diffuse in the cytosol at a higher rate than previously thought.

The Oscillation Amplitude, Not the Frequency of Cytosolic Calcium, Regulates Apoptosis Induction

Although a rising concentration of cytosolic Ca²⁺ has long been recognized as an essential signal for apoptosis, the dynamical mechanisms by which Ca²⁺ regulates apoptosis are not clear yet. To address this, we constructed a computational model that integrates known biochemical reactions and can reproduce the dynamical behaviors of Ca²⁺-induced apoptosis as observed in experiments. Model analysis shows that oscillating Ca²⁺ signals first convert into gradual signals and eventually transform into a switch-like apoptotic response. Via the two processes, the apoptotic signaling pathway filters the frequency of Ca²⁺ oscillations effectively but instead responds acutely to their amplitude. Collectively, our results suggest that Ca²⁺ regulates apoptosis mainly via oscillation amplitude, rather than frequency, modulation. This study not only provides a comprehensive understanding of how oscillatory Ca²⁺ dynamically regulates the complex apoptotic signaling network but also presents a typical example of how Ca²⁺ controls cellular responses through amplitude modulation.

Effects of store-operated and receptor-operated calcium channels on synchronization of calcium oscillations in astrocytes

Intercellular calcium signaling allows cells to communicate with each other and to interact with adjacent cells. Gap junction is the most common and important way for cellular communication. Recently, mathematical models have been widely used to gain a precise and quantitative understanding of the dynamics of intracellular calcium ions (Ca^2＋). In this paper, we establish a mathematical model considering the gap junction permeable to Ca^2＋ and to IP₃ for describing the calcium oscillations in coupled astrocytes. Store-operated calcium entry (SOCE) is viewed as the main process which controls the non-excitable cells, hence, we focus on the effect of store-operated calcium channel (SOCC) and receptor-operated calcium channel (ROCC) on the intercellular synchronization, respectively. By employing bifurcation analysis on this model, the dynamic behaviors of the coupled system with different physiological state cells is obtained with changes in the maximum capacity of the SOCC and the ROCC. The synchronization boundaries for different conditions are gained in the two parameters space of the channel parameters and the coupling strength. The results suggest that the variation of the maximum flow for different calcium channels determines the stable oscillations of the coupled system, as well as for the frequency and amplitude of oscillations. The SOCC has an expected effect on the change of the oscillatory interval while the ROCC demonstrated the influence on the amplitude modulation. Furthermore, the coupling strength and channel parameters could induce 1:1 locking of intercellular Ca^2＋ oscillations and the synchronization region like Arnol'd tongue is found.

Strong current response to slow modulation: A metabolic case-study

We study the current response to periodic driving of a crucial biochemical reaction network, namely, substrate inhibition. We focus on the conversion rate of substrate into product under time-varying metabolic conditions, modeled by a periodic modulation of the product concentration. We find that the system exhibits a strong nonlinear response to small driving frequencies both for the mean time-averaged current and for the fluctuations. For the first, we obtain an analytic formula by coarse-graining the original model to a solvable one. The result is nonperturbative in the modulation amplitude and frequency. We then refine the picture by studying the stochastic dynamics of the full system using a large deviation approach that allows us to show the resonant effect at the level of the time-averaged variance and signal-to-noise ratio. Finally, we discuss how this nonequilibrium effect may play a role in metabolic and synthetic networks.

A mathematical model of calcium dynamics: Obesity and mitochondria-associated ER membranes

Multiple cellular organelles tightly orchestrate intracellular calcium (Ca²⁺) dynamics to regulate cellular activities and maintain homeostasis. The interplay between the endoplasmic reticulum (ER), a major store of intracellular Ca²⁺, and mitochondria, an important source of adenosine triphosphate (ATP), has been the subject of much research, as their dysfunction has been linked with metabolic diseases. Interestingly, throughout the cell’s cytosolic domain, these two organelles share common microdomains called mitochondria-associated ER membranes (MAMs), where their membranes are in close apposition. The role of MAMs is critical for intracellular Ca²⁺ dynamics as they provide hubs for direct Ca²⁺ exchange between the organelles. A recent experimental study reported correlation between obesity and MAM formation in mouse liver cells, and obesity-related cellular changes that are closely associated with the regulation of Ca²⁺ dynamics. We constructed a mathematical model to study the effects of MAM Ca²⁺ dynamics on global Ca²⁺ activities. Through a series of model simulations, we investigated cellular mechanisms underlying the altered Ca²⁺ dynamics in the cells under obesity. We predict that, as the dosage of stimulus gradually increases, liver cells from obese mice will reach the state of saturated cytosolic Ca²⁺ concentration at a lower stimulus concentration, compared to cells from healthy mice.

It is well known that intracellular Ca²⁺ oscillations carry encoded signals in their amplitude and frequency to regulate various cellular processes, and accumulating evidence supports the importance of the interplay between the ER and mitochondria in cellular Ca²⁺ homeostasis. Miscommunications between the organelles may be involved in the development of metabolic diseases. Based on a recent experimental study that spotlighted a correlation between obesity and physical interactions of the ER and mitochondria in mouse hepatic cells, we constructed a mathematical model as a tool to probe the effects of the cellular changes linked with obesity on global cellular Ca²⁺ dynamics. Our model successfully reproduced the experimental study that observed a positive correlation between an increase in ER-mitochondrial junctions and the magnitude of mitochondrial Ca²⁺ responses. We postulate that hepatic cells from lean animals exhibit Ca²⁺ oscillations that are more robust under higher concentrations of stimulus, compared to cells from obese animals.

Dendritic spine geometry and spine apparatus organization govern the spatiotemporal dynamics of calcium

Dendritic spines are small subcompartments that protrude from the dendrites of neurons and are important for signaling activity and synaptic communication. These subcompartments have been characterized to have different shapes. While it is known that these shapes are associated with spine function, the specific nature of these shape-function relationships is not well understood. In this work, we systematically investigated the relationship between the shape and size of both the spine head and spine apparatus, a specialized endoplasmic reticulum compartment within the spine head, in modulating rapid calcium dynamics using mathematical modeling. We developed a spatial multicompartment reaction-diffusion model of calcium dynamics in three dimensions with various flux sources, including N-methyl-D-aspartate receptors (NMDARs), voltage-sensitive calcium channels (VSCCs), and different ion pumps on the plasma membrane. Using this model, we make several important predictions. First, the volume to surface area ratio of the spine regulates calcium dynamics. Second, membrane fluxes impact calcium dynamics temporally and spatially in a nonlinear fashion. Finally, the spine apparatus can act as a physical buffer for calcium by acting as a sink and rescaling the calcium concentration. These predictions set the stage for future experimental investigations of calcium dynamics in dendritic spines.