Reaction-diffusion Simulations Research Articles

Multi-condensate state as a functional strategy to optimize the cell signaling output

The existence of multiple biomolecular condensates inside living cells is a peculiar phenomenon not compatible with the predictions of equilibrium statistical mechanics. In this work, we address the problem of multiple condensates state (MCS) from a functional perspective. We combine Langevin dynamics, reaction-diffusion simulation, and dynamical systems theory to demonstrate that MCS can indeed be a function optimization strategy. Using Arp2/3 mediated actin nucleation pathway as an example, we show that actin polymerization is maximum at an optimal number of condensates. For a fixed amount of Arp2/3, MCS produces a greater response compared to its single condensate counterpart. Our analysis reveals the functional significance of the condensate size distribution which can be mapped to the recent experimental findings. Given the spatial heterogeneity within condensates and non-linear nature of intracellular networks, we envision MCS to be a generic functional solution, so that structures of network motifs may have evolved to accommodate such configurations.

Vesicle and reaction-diffusion hybrid modeling with STEPS

Vesicles carry out many essential functions within cells through the processes of endocytosis, exocytosis, and passive and active transport. This includes transporting and delivering molecules between different parts of the cell, and storing and releasing neurotransmitters in neurons. To date, computational simulation of these key biological players has been rather limited and has not advanced at the same pace as other aspects of cell modeling, restricting the realism of computational models. We describe a general vesicle modeling tool that has been designed for wide application to a variety of cell models, implemented within our software STochastic Engine for Pathway Simulation (STEPS), a stochastic reaction-diffusion simulator that supports realistic reconstructions of cell tissue in tetrahedral meshes. The implementation is validated in an extensive test suite, parallel performance is demonstrated in a realistic synaptic bouton model, and example models are visualized in a Blender extension module.

Observing the Dynamics of an Electrochemically Driven Active Material with Liquid Electron Microscopy.

Electrochemical liquid electron microscopy has revolutionized our understanding of nanomaterial dynamics by allowing for direct observation of their electrochemical production. This technique, primarily applied to inorganic materials, is now being used to explore the self-assembly dynamics of active molecular materials. Our study examines these dynamics across various scales, from the nanoscale behavior of individual fibers to the micrometer-scale hierarchical evolution of fiber clusters. To isolate the influences of the electron beam and electrical potential on material behavior, we conducted thorough beam-sample interaction analyses. Our findings reveal that the dynamics of these active materials at the nanoscale are shaped by their proximity to the electrode and the applied electrical current. By integrating electron microscopy observations with reaction-diffusion simulations, we uncover that local structures and their formation history play a crucial role in determining assembly rates. This suggests that the emergence of nonequilibrium structures can locally accelerate further structural development, offering insights into the behavior of active materials under electrochemical conditions.

Implications of different membrane compartmentalization models in particle-based in silico studies

Studying membrane dynamics is important to understand the cellular response to environmental stimuli. A decisive spatial characteristic of the plasma membrane is its compartmental structure created by the actin-based membrane-skeleton (fences) and anchored transmembrane proteins (pickets). Particle-based reaction–diffusion simulation of the membrane offers a suitable temporal and spatial resolution to analyse its spatially heterogeneous and stochastic dynamics. Fences have been modelled via hop probabilities, potentials or explicit picket fences. Our study analyses the different approaches’ constraints and their impact on simulation results and performance. Each of the methods comes with its own constraints; the picket fences require small timesteps, potential fences might induce a bias in diffusion in crowded systems, and probabilistic fences, in addition to carefully scaling the probability with the timesteps, induce higher computational costs for each propagation step.

Beyond analytic solution: Analysis of FRAP experiments by spatial simulation of the forward problem

Fluorescence redistribution after photobleaching is a commonly used method to understand the dynamic behavior of molecules within cells. Analytic solutions have been developed for specific, well-defined models of dynamic behavior in idealized geometries, but these solutions are inaccurate in complex geometries or when complex binding and diffusion behaviors exist. We demonstrate the use of numerical reaction-diffusion simulations using the Virtual Cell software platform to model fluorescence redistribution after photobleaching experiments. Multiple simulations employing parameter scans and varying bleaching locations and sizes can help to bracket diffusion coefficients and kinetic rate constants in complex image-based geometries. This approach is applied to problems in membrane surface diffusion as well as diffusion and binding in cytosolic volumes in complex cell geometries. In addition, we model diffusion and binding within phase-separated biomolecular condensates (liquid droplets). These are modeled as spherical low-affinity binding domains that also define a high viscosity medium for exchange of the free fluorescently labeled ligand with the external cytosol.

MitoTNT: Mitochondrial Temporal Network Tracking for 4D live-cell fluorescence microscopy data.

Mitochondria form a network in the cell that rapidly changes through fission, fusion, and motility. Dysregulation of this four-dimensional (4D: x,y,z,time) network is implicated in numerous diseases ranging from cancer to neurodegeneration. While lattice light-sheet microscopy has recently made it possible to image mitochondria in 4D, quantitative analysis methods for the resulting datasets have been lacking. Here we present MitoTNT, the first-in-class software for Mitochondrial Temporal Network Tracking in 4D live-cell fluorescence microscopy data. MitoTNT uses spatial proximity and network topology to compute an optimal tracking assignment. To validate the accuracy of tracking, we created a reaction-diffusion simulation to model mitochondrial network motion and remodeling events. We found that our tracking is >90% accurate for ground-truth simulations and agrees well with published motility results for experimental data. We used MitoTNT to quantify 4D mitochondrial networks from human induced pluripotent stem cells. First, we characterized sub-fragment motility and analyzed network branch motion patterns. We revealed that the skeleton node motion is correlated along branch nodes and is uncorrelated in time. Second, we identified fission and fusion events with high spatiotemporal resolution. We found that mitochondrial skeleton nodes near the fission/fusion sites move nearly twice as fast as random skeleton nodes and that microtubules play a role in mediating selective fission/fusion. Finally, we developed graph-based transport simulations that model how material would distribute on experimentally measured mitochondrial temporal networks. We showed that pharmacological perturbations increase network reachability but decrease network resilience through a combination of altered mitochondrial fission/fusion dynamics and motility. MitoTNT's easy-to-use tracking module, interactive 4D visualization capability, and powerful post-tracking analyses aim at making temporal network tracking accessible to the wider mitochondria research community.

Dissecting the lipidic determinants of inner mitochondrial membrane structure.

The inner mitochondrial membrane (IMM) is the primary site of ATP synthesis within eukaryotic cells. The IMM possesses unique, highly curved morphologies defined by cristae membranes (CM), whose large surface area functions to optimize respiration. While a plethora of studies have demonstrated the role of cristae-shaping proteins in generating characteristic CM morphologies, analysis of the roles mitochondrial phospholipids play in shaping CM structure is distinctly lacking. Here we combine lipidomic analysis with multi-scale modelling to demonstrate how multiple components of CMs shape IMM morphology and respiration in yeast. Tuning of phospholipid saturation in engineered strains revealed a surprising breakpoint in CM structure and respiration which phenocopied losses of CM-shaping proteins. Reaction-diffusion simulations on tomogram-derived mitochondrial morphologies showed that lower curvature CM structures can directly contribute to a loss of ATP synthesis. Molecular dynamics (MD) simulations of mitochondrial lipidomes revealed changes in membrane mechanical properties controlled by lipid saturation. We modelled CM tubule formation using continuum mechanics, which revealed that CM formation is defined by a snap-through instability, controlled by the mitochondrial lipidome. Our model also predicted a cooperative interaction between lipid saturation and cardiolipin (CL), which both contribute to curvature. While loss of CL has a minimal phenotype in aerobically-grown yeast, we found that it becomes essential in conditions of increased saturation, induced either genetically or by growth in microaerobic conditions. These results demonstrate how the formation of a highly curved CM is dependent on the mechanical contributions of multiple lipid species, which then dictate mitochondrial function.

4D dynamical whole-cell simulations of a growing minimal cell.

JCVI-syn3A is a genetically minimal bacterial cell with a single chromosome consisting of 493 genes that has retained few regulatory proteins or small RNAs. We recently developed a 3D fully dynamical whole-cell model that simulates the first 20 minutes of the cell cycle for this minimal cell before replication and growth. This model, based on experimental characterization of Syn3A and measurements from numerous biochemical experiments in related organisms, consisted of complete reaction networks for genetic information processing and metabolism totaling over 7200 chemical reactions. We now extend this 4D model to simulate the trajectory of a complete cell cycle while we dynamically update the shape of the membrane, chromosome configuration (including replication state), and polysome structures. To capture the broad range of length- and time-scales of the various dynamic processes, we periodically communicate cell state information among several coupled simulation techniques, including: reaction-diffusion and chemical master equation simulations for stochastic reactions, numerical integration of ODE's for metabolic reactions, and coarse-grained Brownian dynamics for the chromosome and membrane. Membrane shapes during division are inferred from experimental fluorescent images of dividing cells. The chromosome configurations are placed using physically realistic energy functions. From the model, we determine time-dependent concentrations, spatial distributions, and reaction fluxes that offer insight into the principles of life for this minimal cell, for example: prediction of doubling time matching the experimental value, realistic distributions of mRNA half-lives, and time-dependent balance of ATP costs and production for all reactions in the cell.

Quantitative characterization of membrane-protein reversible association using FCS

Functionally meaningful reversible protein-membrane interactions mediate many biological events. Fluorescence correlation spectroscopy (FCS) is increasingly used to quantitatively study the non-reversible binding of proteins to membranes using lipid vesicles in solution. However, the lack of a complete description of the phase and statistical equilibria in the case of reversible protein-membrane partitioning has hampered the application of FCS to quantify the partition coefficient (Kx). In this work, we further extend the theory that describes membrane-protein partitioning to account for spontaneous protein-membrane dissociation and reassociation to the same or a different lipid vesicle. We derive the probability distribution of proteins on lipid vesicles for reversible binding and demonstrate that FCS is a suitable technique for accurate Kx quantification of membrane-protein reversible association. We also establish the limits to Kx determination by FCS studying the Cramer-Rao bound on the variance of the retrieved parameters. We validate the mathematical formulation against reaction-diffusion simulations to study phase and statistical equilibria and compare the Kx obtained from a computational FCS titration experiment with the experimental ground truth. Finally, we demonstrate the application of our methodology studying the association of anti-HIV broadly neutralizing antibody (10E8-3R) to the membrane.

Investigating thrombin-loaded fibrin in whole blood clot microfluidic assay via fluorogenic peptide

During clotting under flow, thrombin rapidly generates fibrin, whereas fibrin potently sequesters thrombin. This co-regulation was studied using microfluidic whole blood clotting on collagen/tissue factor, followed by buffer wash, and a start/stop cycling flow assay using the thrombin fluorogenic substrate, Boc-Val-Pro-Arg-AMC. After 3 min of clotting (100 s−1) and 5 min of buffer wash, non-elutable thrombin activity was easily detected during cycles of flow cessation. Non-elutable thrombin was similarly detected in plasma clots or arterial whole blood clots (1000 s−1). This thrombin activity was ablated by Phe-Pro-Arg-chloromethylketone (PPACK), apixaban, or Gly-Pro-Arg-Pro to inhibit fibrin. Reaction-diffusion simulations predicted 108 nM thrombin within the clot. Heparin addition to the start/stop assay had little effect on fibrin-bound thrombin, whereas addition of heparin-antithrombin (AT) required over 6 min to inhibit the thrombin, indicating a substantial diffusion limitation. In contrast, heparin-AT rapidly inhibited thrombin within microfluidic plasma clots, indicating marked differences in fibrin structure and functionality between plasma clots and whole blood clots. Addition of GPVI-Fab to blood before venous or arterial clotting (200 or 1000 s−1) markedly reduced fibrin-bound thrombin, whereas GPVI-Fab addition after 90 s of clotting had no effect. Perfusion of AF647-fibrinogen over washed fluorescein isothiocyanate (FITC)-fibrin clots resulted in an intense red layer around, but not within, the original FITC-fibrin. Similarly, introduction of plasma/AF647-fibrinogen generated substantial red fibrin masses that did not penetrate the original green clots, demonstrating that fibrin cannot be re-clotted with fibrinogen. Overall, thrombin within fibrin is non-elutable, easily accessed by peptides, slowly accessed by average-sized proteins (heparin/AT), and not accessible to fresh fibrinogen.

STEPS 4.0: Fast and memory-efficient molecular simulations of neurons at the nanoscale

Recent advances in computational neuroscience have demonstrated the usefulness and importance of stochastic, spatial reaction-diffusion simulations. However, ever increasing model complexity renders traditional serial solvers, as well as naive parallel implementations, inadequate. This paper introduces a new generation of the STochastic Engine for Pathway Simulation (STEPS) project (http://steps.sourceforge.net/), denominated STEPS 4.0, and its core components which have been designed for improved scalability, performance, and memory efficiency. STEPS 4.0 aims to enable novel scientific studies of macroscopic systems such as whole cells while capturing their nanoscale details. This class of models is out of reach for serial solvers due to the vast quantity of computation in such detailed models, and also out of reach for naive parallel solvers due to the large memory footprint. Based on a distributed mesh solution, we introduce a new parallel stochastic reaction-diffusion solver and a deterministic membrane potential solver in STEPS 4.0. The distributed mesh, together with improved data layout and algorithm designs, significantly reduces the memory footprint of parallel simulations in STEPS 4.0. This enables massively parallel simulations on modern HPC clusters and overcomes the limitations of the previous parallel STEPS implementation. Current and future improvements to the solver are not sustainable without following proper software engineering principles. For this reason, we also give an overview of how the STEPS codebase and the development environment have been updated to follow modern software development practices. We benchmark performance improvement and memory footprint on three published models with different complexities, from a simple spatial stochastic reaction-diffusion model, to a more complex one that is coupled to a deterministic membrane potential solver to simulate the calcium burst activity of a Purkinje neuron. Simulation results of these models suggest that the new solution dramatically reduces the per-core memory consumption by more than a factor of 30, while maintaining similar or better performance and scalability.

A Kinetic Model for Predicting Trace Gas Uptake and Reaction.

A model is developed to describe trace gas uptake and reaction with applications to aerosols and microdroplets. Gas uptake by the liquid is formulated as a coupled equilibria that links gas, surface, and bulk regions of the droplet or solution. Previously, this framework was used in explicit stochastic reaction-diffusion simulations to predict the reactive uptake kinetics of ozone with droplets containing aqueous aconitic acid, maleic acid, and sodium nitrite. With the use of prior data and simulation results, a new equation for the uptake coefficient is derived, which accounts for both surface and bulk reactions. Lambert W functions are used to obtain closed form solutions to the integrated rate laws for the multiphase kinetics; similar to previous expressions that describe Michaelis-Menten enzyme kinetics. Together these equations couple interface and bulk processes over a wide range of conditions and do not require many of the limiting assumptions needed to apply resistor model formulations to explain trace gas uptake and reaction.

Reaction–diffusion simulations of hydrogen isotope trapping and release from cavities in tungsten, II: Array of cavities

We present a reaction–diffusion model of hydrogen (H) trapping and release from a metal containing an array of cavities. The model is based on a kinetic description of H trapping in chemisorption sites at the cavity surface and H2 gas precipitation in the cavity volume. Significant differences compared to the case of a single cavity (Zibrov and Schmid, 2022) are observed at high cavity number densities due to retrapping effects. They affect the shape of thermal desorption spectroscopy (TDS) curves, the variation of the TDS peak position with the heating rate, and the desorption at a fixed temperature. The TDS spectra from cavities exhibit features that are very difficult to describe with conventional models but are readily observed in experiments: Sharp rising edges of desorption peaks and constant desorption at a fixed temperature. Also the dynamics of H transport and trapping in the material containing cavities differs from that for point defects, especially at low temperatures. The non-equilibrium effects in the interaction of H with cavities are particularly important for metals with low H solubility, such as tungsten and beryllium.

Analytical model of hydrogen inventory saturation in the subsurface of the wall material and comparison to Reaction-Diffusion simulations

We herein introduce an analytical model of hydrogen inventory saturation in the subsurface of materials (several micrometers depth) under plasma implantation. This model is valid for materials for which the desorption process is not limited by hydrogen recombination at the surface. It is based on a simplified assumption for the implantation of both hydrogen ions and atoms (point sources) and on a stationary approach. The model provides an approximation of the density profile of mobile/interstitial hydrogen and of the density profile of hydrogen trapped at materials defects in the subsurface layer. The analytical model shows good agreement with Reaction-Diffusion simulations of deuterium implantation in tungsten at different material temperatures. For the fusion relevant materials tungsten and beryllium, it is shown that most of the inventory is found in traps. The model gives the filling ratio of traps in the subsurface at steady-state fstat,iBULK. This simple parameter indicates how the total subsurface inventory builds up during plasma exposure and provides a simple way to understand the retention dynamics observed during non-linear Reaction-Diffusion simulations.

An automated near-real time computational method for induction and treatment of scar-related ventricular tachycardias.

Catheter ablation is currently the only curative treatment for scar-related ventricular tachycardias (VTs). However, not only are ablation procedures long, with relatively high risk, but success rates are punitively low, with frequent VT recurrence. Personalized in-silico approaches have the opportunity to address these limitations. However, state-of-the-art reaction diffusion (R-D) simulations of VT induction and subsequent circuits used for in-silico ablation target identification require long execution times, along with vast computational resources, which are incompatible with the clinical workflow. Here, we present the Virtual Induction and Treatment of Arrhythmias (VITA), a novel, rapid and fully automated computational approach that uses reaction-Eikonal methodology to induce VT and identify subsequent ablation targets. The rationale for VITA is based on finding isosurfaces associated with an activation wavefront that splits in the ventricles due to the presence of an isolated isthmus of conduction within the scar; once identified, each isthmus may be assessed for their vulnerability to sustain a reentrant circuit, and the corresponding exit site automatically identified for potential ablation targeting. VITA was tested on a virtual cohort of 7 post-infarcted porcine hearts and the results compared to R-D simulations. Using only a standard desktop machine, VITA could detect all scar-related VTs, simulating activation time maps and ECGs (for clinical comparison) as well as computing ablation targets in 48 minutes. The comparable VTs probed by the R-D simulations took 68.5 hours on 256 cores of high-performance computing infrastructure. The set of lesions computed by VITA was shown to render the ventricular model VT-free. VITA could be used in near real-time as a complementary modality aiding in clinical decision-making in the treatment of post-infarction VTs.

Modernizing the NEURON Simulator for Sustainability, Portability, and Performance.

The need for reproducible, credible, multiscale biological modeling has led to the development of standardized simulation platforms, such as the widely-used NEURON environment for computational neuroscience. Developing and maintaining NEURON over several decades has required attention to the competing needs of backwards compatibility, evolving computer architectures, the addition of new scales and physical processes, accessibility to new users, and efficiency and flexibility for specialists. In order to meet these challenges, we have now substantially modernized NEURON, providing continuous integration, an improved build system and release workflow, and better documentation. With the help of a new source-to-source compiler of the NMODL domain-specific language we have enhanced NEURON's ability to run efficiently, via the CoreNEURON simulation engine, on a variety of hardware platforms, including GPUs. Through the implementation of an optimized in-memory transfer mechanism this performance optimized backend is made easily accessible to users, providing training and model-development paths from laptop to workstation to supercomputer and cloud platform. Similarly, we have been able to accelerate NEURON's reaction-diffusion simulation performance through the use of just-in-time compilation. We show that these efforts have led to a growing developer base, a simpler and more robust software distribution, a wider range of supported computer architectures, a better integration of NEURON with other scientific workflows, and substantially improved performance for the simulation of biophysical and biochemical models.

Virtual induction and treatment of arrhythmias (VITA): a fast automated computational tool to induce scar-related tachycardia and identify ablation targets

Abstract Funding Acknowledgements Type of funding sources: Public grant(s) – National budget only. Main funding source(s): British Heart Foundation, Wellcome Trust Background Catheter ablation therapy of post-infarction ventricular tachycardia (VT) is often a lengthy procedure, with a high risk profile, and VT recurrence remains common. While computational modeling may aid pre-procedure planning and target identification, state-of-the-art reaction diffusion (R-D) simulations require access to high performance computing facilities, making them incompatible with clinical workflows. Purpose To present the Virtual Induction and Treatment of Arrhythmias (VITA), a novel, fast and fully automated computational tool to induce VT and identify subsequent ablation targets. Methods VITA employs multisite pacing to automatically find split activation wavefronts associated with the presence of channel isthmuses of viable myocytes within the inexcitable scar. The algorithm probes for all possible associated reentrant pathways and refines these to obtain a set of unique circuits sustaining potential clinical VTs. VT sustenance is assessed by comparing the electrical wavelength with the anatomical circuit path length. Corresponding exit sites of susceptible circuits are automatically identified, and a minimal lesion set constructed for ablation targeting. VITA also allows for rapid simulation of ECGs for clinical comparison. A summary of the VITA is shown in the Figure. VITA was tested on a virtual cohort of 7 post-infarcted porcine hearts and results compared to R-D simulations. Results Using only a desktop computer, VITA successfully identified all possible scar-related VT circuits and optimal ablation targets in just 48mins (per model), compared to 68.5hrs using a high-performance computer required by the R-D approach. Furthermore, VITA consistently identified more VTs associated with intra-scar circuits than R-D (Figure C). Conclusion VITA provides the first near real-time computational tool to aid treatment of post-infarction VTs.

Sliding across a surface: Particles with fixed and mobile ligands.

A quantitative model of the mobility of ligand-presenting particles at the interface is pivotal to understanding important systems in biology and nanotechnology. In this work, we investigate the emerging dynamics of particles featuring ligands that selectively bind receptors decorating an interface. The formation of a ligand-receptor complex leads to a molecular bridge anchoring the particle to the surface. We consider systems with reversible bridges in which ligand-receptor pairs bind/unbind with finite reaction rates. For a given set of bridges, the particle can explore a tiny fraction of the surface as the extensivity of the bridges is finite. We show how, at timescales longer than the bridges' lifetime, the average position of the particle diffuses away from its initial value. We distill our findings into two analytic equations for the sliding diffusion constant of particles carrying mobile and fixed ligands. We quantitatively validate our theoretical predictions using reaction-diffusion simulations. We compare our findings with results from recent literature studies and discuss the molecular parameters that likely affect the particle's mobility most. Our results, along with recent literature studies, will allow inferring the microscopic parameters at play in complex biological systems from experimental trajectories.

Large self-assembled clathrin lattices spontaneously disassemble without sufficient adaptor proteins.

Clathrin-coated structures must assemble on cell membranes to internalize receptors, with the clathrin protein only linked to the membrane via adaptor proteins. These structures can grow surprisingly large, containing over 20 clathrin, yet they often fail to form productive vesicles, instead aborting and disassembling. We show that clathrin structures of this size can both form and disassemble spontaneously when adaptor protein availability is low, despite high abundance of clathrin. Here, we combine recent in vitro kinetic measurements with microscopic reaction-diffusion simulations and theory to differentiate mechanisms of stable vs unstable clathrin assembly on membranes. While in vitro conditions drive assembly of robust, stable lattices, we show that concentrations, geometry, and dimensional reduction in physiologic-like conditions do not support nucleation if only the key adaptor AP-2 is included, due to its insufficient abundance. Nucleation requires a stoichiometry of adaptor to clathrin that exceeds 1:1, meaning additional adaptor types are necessary to form lattices successfully and efficiently. We show that the critical nucleus contains ~25 clathrin, remarkably similar to sizes of the transient and abortive structures observed in vivo. Lastly, we quantify the cost of bending the membrane under our curved clathrin lattices using a continuum membrane model. We find that the cost of bending the membrane could be largely offset by the energetic benefit of forming curved rather than flat structures, with numbers comparable to experiments. Our model predicts how adaptor density can tune clathrin-coated structures from the transient to the stable, showing that active energy consumption is therefore not required for lattice disassembly or remodeling during growth, which is a critical advance towards predicting productive vesicle formation.

Large self-assembled clathrin lattices spontaneously disassemble without sufficient adaptor proteins

Clathrin-coated structures must assemble on cell membranes to internalize receptors. These structures can grow surprisingly large, containing over 20 clathrin seemingly locked in a hexagonal lattice, yet often fail to form productive vesicles. We show that these abortive events happen spontaneously when adaptor availability is low, despite abundant clathrin. Here, we combine recent in vitro kinetic measurements with stochastic reaction-diffusion simulations and theory to differentiate mechanisms of stable vs unstable clathrin assembly on membranes.