Spatio-temporal correlations can drastically change the response of a MAPK pathway

Parallelization of Particle-Based Reaction-Diffusion Simulations Using MPI

Particle-based reaction-diffusion offers a high-resolution alternative to the continuum reaction-diffusion approach, capturing the volume-excluding nature of molecules undergoing stochastic dynamics. This is essential for simulating self-assembly into higher-order structures like filaments, lattices, or macromolecular complexes. Applications of self-assembly are ubiquitous in chemistry, biology, and materials science, but these higher-resolution methods increase computational cost. Here, we present a parallel implementation of the particle-based NERDSS software using the message passing interface (MPI), achieving close to linear scaling for up to 96 processors. By using a spatial decomposition of the system across processors, our approach extends to very large simulation volumes. The scalability of parallel NERDSS is evaluated for reversible reactions and several examples of higher-order self-assembly in 3D and 2D, with all test cases producing accurate solutions. Parallel efficiency depends on the system size, timescales, and reaction network, showing optimal scaling for smaller assemblies with slower timescales. We provide parallel NERDSS code open-source, supporting development and extension to other particle-based reaction-diffusion software.

Bistable Functions and Signaling Motifs in Systems Chemistry: Taking the Next Step Toward Synthetic Cells

A key challenge in modern chemistry research is to mimic life-like functions using simple molecular networks and the integration of such networks into the first functional artificial cell. Central to this endeavor is the development of signaling elements that can regulate the cell function in time and space by producing entities of code with specific information to induce downstream activity. Such artificial signaling motifs can emerge in nonequilibrium systems, exhibiting complex dynamic behavior like bistability, multistability, oscillations, and chaos. However, the de novo, bottom-up design of such systems remains challenging, primarily because the kinetic characteristics and energy aspects yielding bifurcation have not yet been globally defined. We herein review our recent work that focuses on the design and functional analysis of peptide-based networks, propelled by replication reactions and exhibiting bistable behavior. Furthermore, we rationalize and discuss their exploitation and implementation as variable signaling motifs in homogeneous and heterogeneous environments.The bistable reactions constitute reversible second-order autocatalysis as positive feedback to generate two distinct product distributions at steady state (SS), the low-SS and high-SS. Quantitative analyses reveal that a phase transition from simple reversible equilibration dynamics into bistability takes place when the system is continuously fueled, using a reducing agent, to keep it far from equilibrium. In addition, an extensive set of experimental, theoretical, and simulation studies highlight a defined parameter space where bistability operates.Analogous to the arrangement of protein-based bistable motifs in intracellular signaling pathways, sequential concatenation of the synthetic bistable networks is used for signal processing in homogeneous media. The cascaded network output signals are switched and erased or transduced by manipulating the order of addition of the components, allowing it to reach "on demand" either the low-SS or high-SS. The pre-encoded bistable networks are also useful as a programming tool for the downstream regulation of nanoscale materials properties, bridging together the Systems Chemistry and Nanotechnology fields. In such heterogeneous cascade pathways, the outputs of the bistable network serve as input signals for consecutive nanoparticle formation reaction and growth processes, which-depending on the applied conditions-regulate various features of (Au) nanoparticle shape and assembly.Our work enables the design and production of various signaling apparatus that feature higher complexity than previously observed in chemical networks. Future studies, briefly discussed at the end of the Account, will be directed at the design and analysis of more elaborate functionality, such as bistability under flow conditions, multistability, and oscillations. We propose that a profound understanding of the design principles facilitating the replication-based bistability and related functions bear implications for exploring the origin of protein functionality prior to the highly evolved replication-translation-transcription machinery. The integration of our peptide-based signaling motifs within future synthetic cells seems to be a straightforward development of the two alternating states as memory and switch elements for controlling cell growth and division and even communication among different cells. We furthermore suggest that such systems can be introduced into living cells for various biotechnology applications, such as switches for cell temporal and spatial manipulations.

Loss of bimolecular reactions in reaction-diffusion master equations is consistent with diffusion limited reaction kinetics in the mean field limit

We show that the resolution-dependent loss of bimolecular reactions in spatiotemporal Reaction-Diffusion Master Equations (RDMEs) is in agreement with the mean-field Collins-Kimball (C-K) theory of diffusion-limited reaction kinetics. The RDME is a spatial generalization of the chemical master equation, which enables studying stochastic reaction dynamics in spatially heterogeneous systems. It uses a regular Cartesian grid to partition space into locally well-mixed reaction compartments and treats diffusion as a jump reaction between neighboring grid cells. As the chance for reactants to be in the same grid cell decreases for smaller cell widths, the RDME loses bimolecular reactions in finer grids. We show that for a single homo-bimolecular reaction, the mesh spacing can be interpreted as the reaction radius of a well-mixed C-K rate. Then, the bimolecular reaction loss is consistent with diffusion-limited kinetics in the mean-field steady state. In this interpretation, the constant in a bimolecular reaction propensity is no longer the macroscopic reaction rate but the rate of the ballistic C-K step. For the same grid resolution, different diffusion models in RDME, such as those based on finite differences and Gaussian jumps, represent different reaction radii.

Parallelization of particle-based reaction-diffusion simulations using MPI

Particle-based reaction-diffusion models offer a high-resolution alternative to the continuum reaction-diffusion approach, capturing the discrete and volume-excluding nature of molecules undergoing stochastic dynamics. These methods are thus uniquely capable of simulating explicit self-assembly of particles into higher-order structures like filaments, spherical cages, or heterogeneous macromolecular complexes, which are ubiquitous across living systems and in materials design. The disadvantage of these high-resolution methods is their increased computational cost. Here we present a parallel implementation of the particle-based NERDSS software using the Message Passing Interface (MPI) and spatial domain decomposition, achieving close to linear scaling for up to 96 processors in the largest simulation systems. The scalability of parallel NERDSS is evaluated for bimolecular reactions in 3D and 2D, for self-assembly of trimeric and hexameric complexes, and for protein lattice assembly from 3D to 2D, with all parallel test cases producing accurate solutions. We demonstrate how parallel efficiency depends on the system size, the reaction network, and the limiting timescales of the system, showing optimal scaling only for smaller assemblies with slower timescales. The formation of very large assemblies represents a challenge in evaluating reaction updates across processors, and here we restrict assembly sizes to below the spatial decomposition size. We provide the parallel NERDSS code open source, with detailed documentation for developers and extension to other particle-based reaction-diffusion software.

Enzyme kinetics simulation at the scale of individual particles

Enzyme-catalyzed reactions involve two distinct timescales: a short timescale on which enzymes bind to substrate molecules to produce bound complexes and a comparatively long timescale on which the molecules of the complex are transformed into products. The uptake of the substrate in these reactions is the rate at which the product is made on the long timescale. Models often only consider the uptake to reduce the number of chemical species that need to be modeled and to avoid explicitly treating multiple timescales. Typically, the uptake rates cannot be described by mass action kinetics and are traditionally derived by applying singular perturbation theory to the system's governing differential equations. This analysis ignores short timescales by assuming that a pseudo-equilibrium between the enzyme and the enzyme-bound complex is maintained at all times. This assumption cannot be incorporated into current particle-based simulations of reaction-diffusion systems because they utilize proximity-based conditions to govern the instances of reactions that cannot maintain this pseudo-equilibrium for infinitely fast reactions. Instead, these methods must directly simulate the dynamics on the short timescale to accurately model the system. Due to the disparate timescales, such simulations require excessive amounts of computational time before the behavior on the long timescale can be observed. To resolve this problem, we use singular perturbation theory to develop a proximity-based reaction condition that enables us to ignore all fast reactions and directly reproduce non-mass action kinetics at long timescales. To demonstrate our approach, we implement simulations of a specific third order reaction with kinetics reminiscent of the prototypical Michaelis-Menten system.

Efficient and scalable prediction of stochastic reaction-diffusion processes using graph neural networks

The dynamics of locally interacting particles that are distributed in space give rise to a multitude of complex behaviours. However the simulation of reaction-diffusion processes which model such systems is highly computationally expensive, the cost increasing rapidly with the size of space. Here, we devise a graph neural network based approach that uses cheap Monte Carlo simulations of reaction-diffusion processes in a small space to cast predictions of the dynamics of the same processes in a much larger and complex space, including spaces modelled by networks with heterogeneous topology. By applying the method to two biological examples, we show that it leads to accurate results in a small fraction of the computation time of standard stochastic simulation methods. The scalability and accuracy of the method suggest it is a promising approach for studying reaction-diffusion processes in complex spatial domains such as those modelling biochemical reactions, population evolution and epidemic spreading.

Compartmental exchange regulates steady states and stochastic switching of a phosphorylation network

The phosphoregulation of proteins with multiple phosphorylation sites is governed by biochemical reaction networks that can exhibit multistable behavior. However, the behavior of such networks is typically studied in a single reaction volume, while cells are spatially organized into compartments that can exchange proteins. In this work, we use stochastic simulations to study the impact of compartmentalization on a two-site phosphorylation network. We characterize steady states and fluctuation-driven transitions between them as a function of the rate of protein exchange between two compartments. Surprisingly, the average time spent in a state before stochastically switching to another depends nonmonotonically on the protein exchange rate, with the most frequent switching occurring at intermediate exchange rates. At sufficiently small exchange rates, the state of the system and mean switching time are controlled largely by fluctuations in the balance of enzymes in each compartment. This leads to negatively correlated states in the compartments. For large exchange rates, the two compartments behave as a single effective compartment. However, when the compartmental volumes are unequal, the behavior differs from a single compartment with the same total volume. These results demonstrate that exchange of proteins between distinct compartments can regulate the emergent behavior of a common signaling motif.

Kinetics of diffusion-influenced multisite phosphorylation with enzyme reactivation

The simplest way to account for the influence of diffusion on the kinetics of multisite phosphorylation is to modify the rate constants in the conventional rate equations of chemical kinetics. We have previously shown that this is not enough and new transitions between the reactants must also be introduced. Here we extend our results by considering enzymes that are inactive after modifying the substrate and need time to become active again. This generalization leads to a surprising result. The introduction of enzyme reactivation results in a diffusion-modified kinetic scheme with a new transition that has a negative rate constant. The reason for this is that mapping non-Markovian rate equations onto Markovian ones with time-independent rate constants is not a good approximation at short times. We then developed a non-Markovian theory that involves memory kernels instead of rate constants. This theory is now valid at short times, but is more challenging to use. As an example, the diffusion-modified kinetic scheme with new connections was used to calculate kinetics of double phosphorylation and steady-state response in a phosphorylation-dephosphorylation cycle. We have reproduced the loss of bistability in the phosphorylation-dephosphorylation cycle when the enzyme reactivation time decreases, which was obtained by particle-based computer simulations.

Phase plane dynamics of ERK phosphorylation

The extracellular signal-regulated kinase (ERK) controls multiple critical processes in the cell and is deregulated in human cancers, congenital abnormalities, immune diseases, and neurodevelopmental syndromes. Catalytic activity of ERK requires dual phosphorylation by an upstream kinase, in a mechanism that can be described by two sequential Michaelis-Menten steps. The estimation of individual reaction rate constants from kinetic data in the full mechanism has proved challenging. Here, we present an analytically tractable approach to parameter estimation that is based on the phase plane representation of ERK activation and yields two combinations of six reaction rate constants in the detailed mechanism. These combinations correspond to the ratio of the specificities of two consecutive phosphorylations and the probability that monophosphorylated substrate does not dissociate from the enzyme before the second phosphorylation. The presented approach offers a language for comparing the effects of mutations that disrupt ERK activation and function in vivo. As an illustration, we use phase plane representation to analyze dual phosphorylation under heterozygous conditions, when two enzyme variants compete for the same substrate.

Fluorescent biosensors illuminate the spatial regulation of cell signaling across scales

As cell signaling research has advanced, it has become clearer that signal transduction has complex spatiotemporal regulation that goes beyond foundational linear transduction models. Several technologies have enabled these discoveries, including fluorescent biosensors designed to report live biochemical signaling events. As genetically encoded and live-cell compatible tools, fluorescent biosensors are well suited to address diverse cell signaling questions across different spatial scales of regulation. In this review, methods of examining spatial signaling regulation and the design of fluorescent biosensors are introduced. Then, recent biosensor developments that illuminate the importance of spatial regulation in cell signaling are highlighted at several scales, including membranes and organelles, molecular assemblies, and cell/tissue heterogeneity. In closing, perspectives on how fluorescent biosensors will continue enhancing cell signaling research are discussed.

Endoplasmic reticulum network heterogeneity guides diffusive transport and kinetics

The endoplasmic reticulum (ER) is a dynamic network of interconnected sheets and tubules that orchestrates the distribution of lipids, ions, and proteins throughout the cell. The impact of its complex, dynamic morphology on its function as an intracellular transport hub remains poorly understood. To elucidate the functional consequences of ER network structure and dynamics, we quantify how the heterogeneity of the peripheral ER in COS7 cells affects diffusive protein transport. In vivo imaging of photoactivated ER membrane proteins demonstrates their nonuniform spreading to adjacent regions, in a manner consistent with simulations of diffusing particles on extracted network structures. Using a minimal network model to represent tubule rearrangements, we demonstrate that ER network dynamics are sufficiently slow to have little effect on diffusive protein transport. Furthermore, stochastic simulations reveal a novel consequence of ER network heterogeneity: the existence of "hot spots" where sparse diffusive reactants are more likely to find one another. ER exit sites, specialized domains regulating cargo export from the ER, are shown to be disproportionately located in highly accessible regions, further from the outer boundary of the cell. Combining in vivo experiments with analytic calculations, quantitative image analysis, and computational modeling, we demonstrate how structure guides diffusive protein transport and reactions in the ER.

The membrane surface as a platform that organizes cellular and biochemical processes

Membranes are essential for life. They act as semi-permeable boundaries that define cells and organelles. In addition, their surfaces actively participate in biochemical reaction networks, where they confine proteins, align reaction partners, and directly control enzymatic activities. Membrane-localized reactions shape cellular membranes, define the identity of organelles, compartmentalize biochemical processes, and can even be the source of signaling gradients that originate at the plasma membrane and reach into the cytoplasm and nucleus. The membrane surface is, therefore, an essential platform upon which myriad cellular processes are scaffolded. In this review, we summarize our current understanding of the biophysics and biochemistry of membrane-localized reactions with particular focus on insights derived from reconstituted and cellular systems. We discuss how the interplay of cellular factors results in their self-organization, condensation, assembly, and activity, and the emergent properties derived from them.

Colloidal Physics Modeling Reveals How Per-Ribosome Productivity Increases with Growth Rate in Escherichia coli

Faster-growing cells must synthesize proteins more quickly. Increased ribosome abundance only partly accounts for increases in total protein synthesis rates. The productivity of individual ribosomes must increase too, almost doubling by an unknown mechanism. Prior models point to diffusive transport as a limiting factor but raise a paradox: faster-growing cells are more crowded, yet crowding slows diffusion. We suspected that physical crowding, transport, and stoichiometry, considered together, might reveal a more nuanced explanation. To investigate, we built a first-principles physics-based model of Escherichia coli cytoplasm in which Brownian motion and diffusion arise directly from physical interactions between individual molecules of finite size, density, and physiological abundance. Using our microscopically detailed model, we predicted that physical transport of individual ternary complexes accounts for ~80% of translation elongation latency. We also found that volumetric crowding increases during faster growth even as cytoplasmic mass density remains relatively constant. Despite slowed diffusion, we predicted that improved proximity between ternary complexes and ribosomes wins out, illustrating a simple physics-based mechanism for how individual elongating ribosomes become more productive. We speculate that crowding imposes a physical limit on growth rate and undergirds cellular behavior more broadly. Unfitted colloidal-scale modeling offers systems biology a complementary "physics engine" for exploring how cellular-scale behaviors arise from physical transport and reactions among individual molecules. IMPORTANCE Ribosomes are the factories in cells that synthesize proteins. When cells grow faster, there are not enough ribosomes to keep up with the demand for faster protein synthesis without individual ribosomes becoming more productive. Yet, faster-growing cells are more crowded, seemingly making it harder for each ribosome to do its work. Our computational model of the physics of translation elongation reveals the underlying mechanism for how individual ribosomes become more productive: proximity and stoichiometry of translation molecules overcome crowding. Our model also suggests a universal physical limitation of cell growth rates.

Extracting, quantifying, and comparing dynamical and biomechanical properties of living matter through single particle tracking

A panoply of new tools for tracking single particles and molecules has led to an explosion of experimental data, leading to novel insights into physical properties of living matter governing cellular development and function, health and disease. In this Perspective, we present tools to investigate the dynamics and mechanics of living systems from the molecular to cellular scale via single-particle techniques. In particular, we focus on methods to measure, interpret, and analyse complex data sets that are associated with forces, materials properties, transport, and emergent organisation phenomena within biological and soft-matter systems. Current approaches, challenges, and existing solutions in the associated fields are outlined in order to support the growing community of researchers at the interface of physics and the life sciences. Each section focuses not only on the general physical principles and the potential for understanding living matter, but also on details of practical data extraction and analysis, discussing limitations, interpretation, and comparison across different experimental realisations and theoretical frameworks. Particularly relevant results are introduced as examples. While this Perspective describes living matter from a physical perspective, highlighting experimental and theoretical physics techniques relevant for such systems, it is also meant to serve as a solid starting point for researchers in the life sciences interested in the implementation of biophysical methods.

A Gaussian jump process formulation of the reaction–diffusion master equation enables faster exact stochastic simulations

We propose a Gaussian jump process model on a regular Cartesian lattice for the diffusion part of the Reaction–Diffusion Master Equation (RDME). We derive the resulting Gaussian RDME (GRDME) formulation from analogy with a kernel-based discretization scheme for continuous diffusion processes and quantify the limits of its validity relative to the classic RDME. We then present an exact stochastic simulation algorithm for the GRDME, showing that the accuracies of GRDME and RDME are comparable, but exact simulations of the GRDME require only a fraction of the computational cost of exact RDME simulations. We analyze the origin of this speedup and its scaling with problem dimension. The benchmarks suggest that the GRDME is a particularly beneficial model for diffusion-dominated systems in three dimensional spaces, often occurring in systems biology and cell biology.

Algebra, Geometry and Topology of ERK Kinetics

The MEK/ERK signalling pathway is involved in cell division, cell specialisation, survival and cell death (Shaul and Seger in Biochim Biophys Acta (BBA)-Mol Cell Res 1773(8):1213-1226, 2007). Here we study a polynomial dynamical system describing the dynamics of MEK/ERK proposed by Yeung et al. (Curr Biol 2019, https://doi.org/10.1016/j.cub.2019.12.052 ) with their experimental setup, data and known biological information. The experimental dataset is a time-course of ERK measurements in different phosphorylation states following activation of either wild-type MEK or MEK mutations associated with cancer or developmental defects. We demonstrate how methods from computational algebraic geometry, differential algebra, Bayesian statistics and computational algebraic topology can inform the model reduction, identification and parameter inference of MEK variants, respectively. Throughout, we show how this algebraic viewpoint offers a rigorous and systematic analysis of such models.

Detailed balance for particle models of reversible reactions in bounded domains

In particle-based stochastic reaction–diffusion models, reaction rates and placement kernels are used to decide the probability per time a reaction can occur between reactant particles and to decide where product particles should be placed. When choosing kernels to use in reversible reactions, a key constraint is to ensure that detailed balance of spatial reaction fluxes holds at all points at equilibrium. In this work, we formulate a general partial-integral differential equation model that encompasses several of the commonly used contact reactivity (e.g., Smoluchowski-Collins-Kimball) and volume reactivity (e.g., Doi) particle models. From these equations, we derive a detailed balance condition for the reversible A + B ⇆ C reaction. In bounded domains with no-flux boundary conditions, when choosing unbinding kernels consistent with several commonly used binding kernels, we show that preserving detailed balance of spatial reaction fluxes at all points requires spatially varying unbinding rate functions near the domain boundary. Brownian dynamics simulation algorithms can realize such varying rates through ignoring domain boundaries during unbinding and rejecting unbinding events that result in product particles being placed outside the domain.

Dephosphorylation accelerates the dissociation of ZAP70 from the T cell receptor

Protein-protein binding domains are critical in signaling networks. Src homology 2 (SH2) domains are binding domains that interact with sequences containing phosphorylated tyrosines. A subset of SH2 domain-containing proteins has tandem domains, which are thought to enhance binding affinity and specificity. However, a trade-off exists between long-lived binding and the ability to rapidly reverse signaling, which is a critical requirement of noise-filtering mechanisms such as kinetic proofreading. Here, we use modeling to show that the unbinding rate of tandem, but not single, SH2 domains can be accelerated by phosphatases. Using surface plasmon resonance, we show that the phosphatase CD45 can accelerate the unbinding rate of zeta chain-associated protein kinase 70 (ZAP70), a tandem SH2 domain-containing kinase, from biphosphorylated peptides from the T cell receptor (TCR). An important functional prediction of accelerated unbinding is that the intracellular ZAP70-TCR half-life in T cells will not be fixed but rather, dependent on the extracellular TCR-antigen half-life, and we show that this is the case in both cell lines and primary T cells. The work highlights that tandem SH2 domains can break the trade-off between signal fidelity (requiring long half-life) and signal reversibility (requiring short half-life), which is a key requirement for T cell antigen discrimination mediated by kinetic proofreading.

Morphology and Transport in Eukaryotic Cells

Transport of intracellular components relies on a variety of active and passive mechanisms, ranging from the diffusive spreading of small molecules over short distances to motor-driven motion across long distances. The cell-scale behavior of these mechanisms is fundamentally dependent on the morphology of the underlying cellular structures. Diffusion-limited reaction times can be qualitatively altered by the presence of occluding barriers or by confinement in complex architectures, such as those of reticulated organelles. Motor-driven transport is modulated by the architecture of cytoskeletal filaments that serve as transport highways. In this review, we discuss the impact of geometry on intracellular transport processes that fulfill a broad range of functional objectives, including delivery, distribution, and sorting of cellular components. By unraveling the interplay between morphology and transport efficiency, we aim to elucidate key structure-function relationships that govern the architecture of transport systems at the cellular scale. Expected final online publication date for the Annual Review of Biophysics, Volume 51 is May 2022. Please see http://www.annualreviews.org/page/journal/pubdates for revised estimates.

Multiscale molecular kinetics by coupling Markov state models and reaction-diffusion dynamics

A novel approach to simulate simple protein-ligand systems at large time and length scales is to couple Markov state models (MSMs) of molecular kinetics with particle-based reaction-diffusion (RD) simulations, MSM/RD. Currently, MSM/RD lacks a mathematical framework to derive coupling schemes, is limited to isotropic ligands in a single conformational state, and lacks multiparticle extensions. In this work, we address these needs by developing a general MSM/RD framework by coarse-graining molecular dynamics into hybrid switching diffusion processes. Given enough data to parameterize the model, it is capable of modeling protein-protein interactions over large time and length scales, and it can be extended to handle multiple molecules. We derive the MSM/RD framework, and we implement and verify it for two protein-protein benchmark systems and one multiparticle implementation to model the formation of pentameric ring molecules. To enable reproducibility, we have published our code in the MSM/RD software package.

Influence of molecular rebinding on the reaction rate of complex formation

We simulated Brownian diffusion and reaction-diffusion processes to study the influence of molecular rebinding on the reaction rates of bimolecular reactions. We found that the number of rebinding events, N_reb, is proportional to the target's size and inversely proportional to the diffusion coefficient D and simulation time-step Δt. We found the proportionality constant close to π^-1/2. We confirmed that N_reb is defined as a ratio of the activation-limited rate constant k_a to the diffusion-limited rate constant, k_D. We provide the formula describing the reactivity coefficient κ, modelling the transient-native complex transition for the activation-controlled reaction rates. We show that κ is proportional to (D/Δt)^1/2. Finally, we apply our rebinding-including reaction rate model to the real reactions of photoacid dissociation and protein association. Based on literature data for both types of reactions, we found the Δt time-scale. We show that for the photodissociation of a proton, the Δt is equal to 171 ± 18 fs and the average number of rebinding events is approximately equal to 40. For proteins, Δt is of the order of 100 ps with around 20 rebinding events. In both cases the timescale is similar to the timescale of fluctuation of the solvent molecules surrounding the reactants; vibrations and bending in the case of photoacid dissociation and diffusional motion for proteins.

A multiscale compartment-based model of stochastic gene regulatory networks using hitting-time analysis

Spatial stochastic models of single cell kinetics are capable of capturing both fluctuations in molecular numbers and the spatial dependencies of the key steps of intracellular regulatory networks. The spatial stochastic model can be simulated both on a detailed microscopic level using particle tracking and on a mesoscopic level using the reaction-diffusion master equation. However, despite substantial progress on simulation efficiency for spatial models in the last years, the computational cost quickly becomes prohibitively expensive for tasks that require repeated simulation of thousands or millions of realizations of the model. This limits the use of spatial models in applications such as multicellular simulations, likelihood-free parameter inference, and robustness analysis. Further approximation of the spatial dynamics is needed to accelerate such computational engineering tasks. We here propose a multiscale model where a compartment-based model approximates a detailed spatial stochastic model. The compartment model is constructed via a first-exit time analysis on the spatial model, thus capturing critical spatial aspects of the fine-grained simulations, at a cost close to the simple well-mixed model. We apply the multiscale model to a canonical model of negative-feedback gene regulation, assess its accuracy over a range of parameters, and demonstrate that the approximation can yield substantial speedups for likelihood-free parameter inference.

Space-time histories approach to fast stochastic simulation of bimolecular reactions

Computational models of reaction-diffusion systems involving low copy numbers or strongly heterogeneous molecular spatial distributions, such as those frequently found in cellular signaling pathways, require approaches that account for the stochastic dynamics of individual particles, as opposed to approaches representing them through their average concentrations. Efforts to remedy the high computational cost associated with particle-based stochastic approaches by taking advantage of Green's functions are hampered by the need to draw random numbers from complicated, and therefore costly, non-standard probability distributions to update particle positions. Here, we introduce an approach that permits the reconstruction of entire molecular trajectories, including bimolecular encounters, in retrospect, after a simulated time step, while avoiding inefficient draws from non-standard distributions. This means that highly accurate stochastic simulations can be performed for system sizes that would be prohibitively costly to simulate with conventional Green's function based methods. The algorithm applies equally well to one, two, and three dimensional systems and can be readily extended to include deterministic forces specified by an interaction potential, such as the Coulomb potential.

Diffusive search and trajectories on tubular networks: a propagator approach

Several organelles in eukaryotic cells, including mitochondria and the endoplasmic reticulum, form interconnected tubule networks extending throughout the cell. These tubular networks host many biochemical pathways that rely on proteins diffusively searching through the network to encounter binding partners or localized target regions. Predicting the behavior of such pathways requires a quantitative understanding of how confinement to a reticulated structure modulates reaction kinetics. In this work, we develop both exact analytical methods to compute mean first passage times and efficient kinetic Monte Carlo algorithms to simulate trajectories of particles diffusing in a tubular network. Our approach leverages exact propagator functions for the distribution of transition times between network nodes and allows large simulation time steps determined by the network structure. The methodology is applied to both synthetic planar networks and organelle network structures, demonstrating key general features such as the heterogeneity of search times in different network regions and the functional advantage of broadly distributing target sites throughout the network. The proposed algorithms pave the way for future exploration of the interrelationship between tubular network structure and biomolecular reaction kinetics.

Speed limits of protein assembly with reversible membrane localization

Protein assembly is often studied in a three-dimensional solution, but a significant fraction of binding events involve proteins that can reversibly bind and diffuse along a two-dimensional surface. In a recent study, we quantified how proteins can exploit the reduced dimensionality of the membrane to trigger complex formation. Here, we derive a single expression for the characteristic timescale of this multi-step assembly process, where the change in dimensionality renders rates and concentrations effectively time-dependent. We find that proteins can accelerate dimer formation due to an increase in relative concentration, driving more frequent collisions, which often win out over slow-downs due to diffusion. Our model contains two protein populations that dimerize with one another and use a distinct site to bind membrane lipids, creating a complex reaction network. However, by identifying two major rate-limiting pathways to reach an equilibrium steady-state, we derive an excellent approximation for the mean first passage time when lipids are in abundant supply. Our theory highlights how the "sticking rate" or effective adsorption coefficient of the membrane is central in controlling timescales. We also derive a corrected localization rate to quantify how the geometry of the system and diffusion can reduce rates of membrane localization. We validate and test our results using kinetic and particle-based reaction-diffusion simulations. Our results establish how the speed of key assembly steps can shift by orders-of-magnitude when membrane localization is possible, which is critical to understanding mechanisms used in cells.

Concentration sensing in crowded environments

Signal transduction within crowded cellular compartments is essential for the physiological function of cells. Although the accuracy with which receptors can probe the concentration of ligands has been thoroughly investigated in dilute systems, the effect of macromolecular crowding on the inference of concentration remains unclear. In this work, we develop an algorithm to simulate reversible reactions between reacting Brownian particles. Our algorithm facilitates the calculation of reaction rates and correlation times for ligand-receptor systems in the presence of macromolecular crowding. Using this method, we show that it is possible for crowding to increase the accuracy of estimated ligand concentration based on receptor occupancy. In particular, we find that crowding can enhance the effective association rates between small ligands and receptors to a degree sufficient to overcome the increased chance of rebinding due to caging by crowding molecules. For larger ligands, crowding decreases the accuracy of the receptor's estimate primarily by decreasing the microscopic association and dissociation rates.

Statistics of Reaction Flux and Dynamical Activity Associated with a Diffusion-Influenced Ligand-Binding Reaction

In this paper, we consider a macromolecule with two competitive binding sites where a ligand can bind to and gives rise to a unicyclic reaction network consisting of four states-(i) a single state with both binding sites vacant, (ii) two states with one bound site and one free binding site, and (iii) an another single state with both sites occupied. We obtain probability densities of the time-integrated current along the clockwise direction and the dynamical activity or mean number of jumps between different states for finite times at a fast diffusion limit. On the other hand, in the diffusion-limited case, ligand diffusion between the two binding sites directly connects the mono-ligated states-changing the reaction scheme. Addition of the new reaction channel alters the precision of ligand occupancy to a single site, the mean dynamical activity, and the mean entropy production rate. All of these quantities are calculated with varying degrees of competition between the two sites for ligands, and we find that increase in the competition between the two sites decreases all above-mentioned additive functionals. The upper bound of precision associated with a single-site ligand occupancy for a diffusion-influenced reaction network is set by the mean dynamical activity (mean entropy production rate) at small (large) ligand concentrations at the steady-state limit.

Quantifying the roles of space and stochasticity in computer simulations for cell biology and cellular biochemistry

Most of the fascinating phenomena studied in cell biology emerge from interactions among highly organized multimolecular structures embedded into complex and frequently dynamic cellular morphologies. For the exploration of such systems, computer simulation has proved to be an invaluable tool, and many researchers in this field have developed sophisticated computational models for application to specific cell biological questions. However, it is often difficult to reconcile conflicting computational results that use different approaches to describe the same phenomenon. To address this issue systematically, we have defined a series of computational test cases ranging from very simple to moderately complex, varying key features of dimensionality, reaction type, reaction speed, crowding, and cell size. We then quantified how explicit spatial and/or stochastic implementations alter outcomes, even when all methods use the same reaction network, rates, and concentrations. For simple cases, we generally find minor differences in solutions of the same problem. However, we observe increasing discordance as the effects of localization, dimensionality reduction, and irreversible enzymatic reactions are combined. We discuss the strengths and limitations of commonly used computational approaches for exploring cell biological questions and provide a framework for decision making by researchers developing new models. As computational power and speed continue to increase at a remarkable rate, the dream of a fully comprehensive computational model of a living cell may be drawing closer to reality, but our analysis demonstrates that it will be crucial to evaluate the accuracy of such models critically and systematically.

Crosstalk and ultrasensitivity in protein degradation pathways

Protein turnover is vital to cellular homeostasis. Many proteins are degraded efficiently only after they have been post-translationally “tagged” with a polyubiquitin chain. Ubiquitylation is a form of Post-Translational Modification (PTM): addition of a ubiquitin to the chain is catalyzed by E3 ligases, and removal of ubiquitin is catalyzed by a De-UBiquitylating enzyme (DUB). Nearly four decades ago, Goldbeter and Koshland discovered that reversible PTM cycles function like on-off switches when the substrates are at saturating concentrations. Although this finding has had profound implications for the understanding of switch-like behavior in biochemical networks, the general behavior of PTM cycles subject to synthesis and degradation has not been studied. Using a mathematical modeling approach, we found that simply introducing protein turnover to a standard modification cycle has profound effects, including significantly reducing the switch-like nature of the response. Our findings suggest that many classic results on PTM cycles may not hold in vivo where protein turnover is ubiquitous. We also found that proteins sharing an E3 ligase can have closely related changes in their expression levels. These results imply that it may be difficult to interpret experimental results obtained from either overexpressing or knocking down protein levels, since changes in protein expression can be coupled via E3 ligase crosstalk. Understanding crosstalk and competition for E3 ligases will be key in ultimately developing a global picture of protein homeostasis.

Previous work has shown that substrates of Post-Translational Modification (PTM) cycles can have coupled responses if those substrates share enzymes. This implies that modifications leading to substrate degradation (e.g. ubiquitylation by an E3 ligase) could introduce coupling in concentrations of substrates sharing a ligase. Using mathematical models, we found adding protein turnover to a PTM cycle diminishes both sensitivity and ultrasensitivity, particularly in models admitting long ubiquitin chains. We also found that proteins sharing an E3 ligase can indeed have coupled changes in both expression and sensitivity to signals. These results imply that accounting for crosstalk in protein degradation networks is crucial for the interpretation of results from a wide variety of common experimental perturbations to living systems.

Strong intracellular signal inactivation produces sharper and more robust signaling from cell membrane to nucleus

For a chemical signal to propagate across a cell, it must navigate a tortuous environment involving a variety of organelle barriers. In this work we study mathematical models for a basic chemical signal, the arrival times at the nuclear membrane of proteins that are activated at the cell membrane and diffuse throughout the cytosol. Organelle surfaces within human B cells are reconstructed from soft X-ray tomographic images, and modeled as reflecting barriers to the molecules’ diffusion. We show that signal inactivation sharpens signals, reducing variability in the arrival time at the nuclear membrane. Inactivation can also compensate for an observed slowdown in signal propagation induced by the presence of organelle barriers, leading to arrival times at the nuclear membrane that are comparable to models in which the cytosol is treated as an open, empty region. In the limit of strong signal inactivation this is achieved by filtering out molecules that traverse non-geodesic paths.

The inside of cells is a complex spatial environment, filled with organelles, filaments and proteins. It is an open question how cell signaling pathways function robustly in the presence of such spatial heterogeneity. In this work we study how organelle barriers influence the most basic of chemical signals; the diffusive propagation of an activated protein from the cell membrane to nucleus. Three-dimensional B cell organelle and membrane geometries reconstructed from soft X-ray tomographic images are used in building mathematical models of the signal propagation process. Our models demonstrate that organelle barriers significantly increase the time required for a diffusing protein to traverse from the cell membrane to nucleus when compared to a cell with an empty cytosolic space. We also show that signal inactivation, a fundamental component of all signaling pathways, can provide robustness in the signal arrival time in two ways. Increasing rates of signal inactivation reduce variability in the arrival time, while also dramatically reducing the degree to which organelle barriers increase the arrival time (in comparison to a cell with an empty cytosol).

Getting around the cell: physical transport in the intracellular world

Eukaryotic cells face the challenging task of transporting a variety of particles through the complex intracellular milieu in order to deliver, distribute, and mix the many components that support cell function. In this review, we explore the biological objectives and physical mechanisms of intracellular transport. Our focus is on cytoplasmic and intra-organelle transport at the whole-cell scale. We outline several key biological functions that depend on physically transporting components across the cell, including the delivery of secreted proteins, support of cell growth and repair, propagation of intracellular signals, establishment of organelle contacts, and spatial organization of metabolic gradients. We then review the three primary physical modes of transport in eukaryotic cells: diffusive motion, motor-driven transport, and advection by cytoplasmic flow. For each mechanism, we identify the main factors that determine speed and directionality. We also highlight the efficiency of each transport mode in fulfilling various key objectives of transport, such as particle mixing, directed delivery, and rapid target search. Taken together, the interplay of diffusion, molecular motors, and flows supports the intracellular transport needs that underlie a broad variety of biological phenomena.

NERDSS: A Nonequilibrium Simulator for Multibody Self-Assembly at the Cellular Scale

Currently, a significant barrier to building predictive models of cellular self-assembly processes is that molecular models cannot capture minutes-long dynamics that couple distinct components with active processes, whereas reaction-diffusion models cannot capture structures of molecular assembly. Here, we introduce the nonequilibrium reaction-diffusion self-assembly simulator (NERDSS), which addresses this spatiotemporal resolution gap. NERDSS integrates efficient reaction-diffusion algorithms into generalized software that operates on user-defined molecules through diffusion, binding and orientation, unbinding, chemical transformations, and spatial localization. By connecting the fast processes of binding with the slow timescales of large-scale assembly, NERDSS integrates molecular resolution with reversible formation of ordered, multisubunit complexes. NERDSS encodes models using rule-based formatting languages to facilitate model portability, usability, and reproducibility. Applying NERDSS to steps in clathrin-mediated endocytosis, we design multicomponent systems that can form lattices in solution or on the membrane, and we predict how stochastic but localized dephosphorylation of membrane lipids can drive lattice disassembly. The NERDSS simulations reveal the spatial constraints on lattice growth and the role of membrane localization and cooperativity in nucleating assembly. By modeling viral lattice assembly and recapitulating oscillations in protein expression levels for a circadian clock model, we illustrate the adaptability of NERDSS. NERDSS simulates user-defined assembly models that were previously inaccessible to existing software tools, with broad applications to predicting self-assembly in vivo and designing high-yield assemblies in vitro.

The interplay of spatial organization and biochemistry in building blocks of cellular signalling pathways

Biochemical pathways and networks are central to cellular information processing. While a broad range of studies have dissected multiple aspects of information processing in biochemical pathways, the effect of spatial organization remains much less understood. It is clear that space is central to intracellular organization, plays important roles in cellular information processing and has been exploited in evolution; additionally, it is being increasingly exploited in synthetic biology through the development of artificial compartments, in a variety of ways. In this paper, we dissect different aspects of the interplay between spatial organization and biochemical pathways, by focusing on basic building blocks of these pathways: covalent modification cycles and two-component systems, with enzymes which may be monofunctional or bifunctional. Our analysis of spatial organization is performed by examining a range of 'spatial designs': patterns of localization or non-localization of enzymes/substrates, theoretically and computationally. Using these well-characterized in silico systems, we analyse the following. (i) The effect of different types of spatial organization on the overall kinetics of modification, and the role of distinct modification mechanisms therein. (ii) How different information processing characteristics seen experimentally and studied from the viewpoint of kinetics are perturbed, or generated. (iii) How the activity of enzymes (bifunctional enzymes in particular) may be spatially manipulated, and the relationship between localization and activity. (iv) How transitions in spatial organization (encountered either through evolution or through the lifetime of cells, as seen in multiple model organisms) impacts the kinetic module (and pathway) behaviour, and how transitions in chemistry may be impacted by prior spatial organization. The basic insights which emerge are central to understanding the role of spatial organization in biochemical pathways in both bacteria and eukaryotes, and are of direct relevance to engineering spatial organization of pathways in bottom-up synthetic biology in cellular and cell-free systems.

Multisite reversible association in membranes and solutions: From non-Markovian to Markovian kinetics

The role of diffusion on the kinetics of reversible association to a macromolecule with two inequivalent sites is studied. Previously, we found that, in the simplest possible description, it is not sufficient to just renormalize the rate constants of chemical kinetics, but one must introduce direct transitions between the bound states in the kinetic scheme. The physical reason for this is that a molecule that just dissociated from one site can directly rebind to the other rather than diffuse away into the bulk. Such a simple description is not valid in two dimensions because reactants can never diffuse away into the bulk. In this work, we consider a variety of more sophisticated implementations of our recent general theory that are valid in both two and three dimensions. We compare the predicted time dependence of the concentrations for a wide range of parameters and establish the range of validity of various levels of the general theory.

Presence or Absence of Ras Dimerization Shows Distinct Kinetic Signature in Ras-Raf Interaction

Initiations of cell signaling pathways often occur through the formation of multiprotein complexes that form through protein-protein interactions. Therefore, detecting their presence is central to understanding the function of a cell signaling pathway, aberration of which often leads to fatal diseases, including cancers. However, the multiprotein complexes are often difficult to detect using microscopes due to their small sizes. Therefore, currently, their presence can be only detected through indirect means. In this article, we propose to investigate the presence or absence of protein complexes through some easily measurable kinetic parameters, such as activation rates. As a proof of concept, we investigate the Ras-Raf system, a well-characterized cell signaling system. It has been hypothesized that Ras dimerization is necessary to create activated Raf dimers. Although there are circumstantial evidences supporting the Ras dimerization hypothesis, direct proof of Ras dimerization is still inconclusive. In the absence of conclusive direct experimental proof, this hypothesis can only be examined through indirect evidences of Ras dimerization. In this article, using a multiscale simulation technique, we provide multiple criteria that distinguishes an activation mechanism involving Ras dimerization from another mechanism that does not involve Ras dimerization. The provided criteria will be useful in the investigation of not only Ras-Raf interaction but also other two-protein interactions.

Inference of Multisite Phosphorylation Rate Constants and Their Modulation by Pathogenic Mutations

Multisite protein phosphorylation plays a critical role in cell regulation [1-3]. It is widely appreciated that the functional capabilities of multisite phosphorylation depend on the order and kinetics of phosphorylation steps, but kinetic aspects of multisite phosphorylation remain poorly understood [4-6]. Here, we focus on what appears to be the simplest scenario, when a protein is phosphorylated on only two sites in a strict, well-defined order. This scenario describes the activation of ERK, a highly conserved cell-signaling enzyme. We use Bayesian parameter inference in a structurally identifiable kinetic model to dissect dual phosphorylation of ERK by MEK, a kinase that is mutated in a large number of human diseases [7-12]. Our results reveal how enzyme processivity and efficiencies of individual phosphorylation steps are altered by pathogenic mutations. The presented approach, which connects specific mutations to kinetic parameters of multisite phosphorylation mechanisms, provides a systematic framework for closing the gap between studies with purified enzymes and their effects in the living organism.

pSpatiocyte: a high-performance simulator for intracellular reaction-diffusion systems

Background

Studies using quantitative experimental methods have shown that intracellular spatial distribution of molecules plays a central role in many cellular systems. Spatially resolved computer simulations can integrate quantitative data from these experiments to construct physically accurate models of the systems. Although computationally expensive, microscopic resolution reaction-diffusion simulators, such as Spatiocyte can directly capture intracellular effects comprising diffusion-limited reactions and volume exclusion from crowded molecules by explicitly representing individual diffusing molecules in space. To alleviate the steep computational cost typically associated with the simulation of large or crowded intracellular compartments, we present a parallelized Spatiocyte method called pSpatiocyte.

Results

The new high-performance method employs unique parallelization schemes on hexagonal close-packed (HCP) lattice to efficiently exploit the resources of common workstations and large distributed memory parallel computers. We introduce a coordinate system for fast accesses to HCP lattice voxels, a parallelized event scheduler, a parallelized Gillespie’s direct-method for unimolecular reactions, and a parallelized event for diffusion and bimolecular reaction processes. We verified the correctness of pSpatiocyte reaction and diffusion processes by comparison to theory. To evaluate the performance of pSpatiocyte, we performed a series of parallelized diffusion runs on the RIKEN K computer. In the case of fine lattice discretization with low voxel occupancy, pSpatiocyte exhibited 74% parallel efficiency and achieved a speedup of 7686 times with 663552 cores compared to the runtime with 64 cores. In the weak scaling performance, pSpatiocyte obtained efficiencies of at least 60% with up to 663552 cores. When executing the Michaelis-Menten benchmark model on an eight-core workstation, pSpatiocyte required 45- and 55-fold shorter runtimes than Smoldyn and the parallel version of ReaDDy, respectively. As a high-performance application example, we study the dual phosphorylation-dephosphorylation cycle of the MAPK system, a typical reaction network motif in cell signaling pathways.

Conclusions

pSpatiocyte demonstrates good accuracies, fast runtimes and a significant performance advantage over well-known microscopic particle methods in large-scale simulations of intracellular reaction-diffusion systems. The source code of pSpatiocyte is available at https://spatiocyte.org.

Hierarchical algorithm for the reaction-diffusion master equation

We have developed an algorithm coupling mesoscopic simulations on different levels in a hierarchy of Cartesian meshes. Based on the multiscale nature of the chemical reactions, some molecules in the system will live on a fine-grained mesh, while others live on a coarse-grained mesh. By allowing molecules to transfer from the fine levels to the coarse levels when appropriate, we show that we can save up to three orders of magnitude of computational time compared to microscopic simulations or highly resolved mesoscopic simulations, without losing significant accuracy. We demonstrate this in several numerical examples with systems that cannot be accurately simulated with a coarse-grained mesoscopic model.

Potential based, spatial simulation of dynamically nested particles

BACKGROUND

To study cell biological phenomena which depend on diffusion, active transport processes, or the locations of species, modeling and simulation studies need to take space into account. To describe the system as a collection of discrete objects moving and interacting in continuous space, various particle-based reaction diffusion simulators for cell-biological system have been developed. So far the focus has been on particles as solid spheres or points. However, spatial dynamics might happen at different organizational levels, such as proteins, vesicles or cells with interrelated dynamics which requires spatial approaches that take this multi-levelness of cell biological systems into account.

RESULTS

Based on the perception of particles forming hollow spheres, ML-Force contributes to the family of particle-based simulation approaches: in addition to excluded volumes and forces, it also supports compartmental dynamics and relating dynamics between different organizational levels explicitly. Thereby, compartmental dynamics, e.g., particles entering and leaving other particles, and bimolecular reactions are modeled using pair-wise potentials (forces) and the Langevin equation. In addition, forces that act independently of other particles can be applied to direct the movement of particles. Attributes and the possibility to define arbitrary functions on particles, their attributes and content, to determine the results and kinetics of reactions add to the expressiveness of ML-Force. Its implementation comprises a rudimentary rule-based embedded domain-specific modeling language for specifying models and a simulator for executing models continuously. Applications inspired by cell biological models from literature, such as vesicle transport or yeast growth, show the value of the realized features. They facilitate capturing more complex spatial dynamics, such as the fission of compartments or the directed movement of particles, and enable the integration of non-spatial intra-compartmental dynamics as stochastic events.

CONCLUSIONS

By handling all dynamics based on potentials (forces) and the Langevin equation, compartmental dynamics, such as dynamic nesting, fusion and fission of compartmental structures are handled continuously and are seamlessly integrated with traditional particle-based reaction-diffusion dynamics within the cell. Thereby, attributes and arbitrary functions allow to flexibly describe diverse spatial phenomena, and relate dynamics across organizational levels. Also they prove crucial in modeling intra-cellular or intra-compartmental dynamics in a non-spatial manner, and, thus, to abstract from spatial dynamics, on demand which increases the range of multi-compartmental processes that can be captured.

Revisiting the Reduction of Stochastic Models of Genetic Feedback Loops with Fast Promoter Switching

Propensity functions of the Hill type are commonly used to model transcriptional regulation in stochastic models of gene expression. This leads to an effective reduced master equation for the mRNA and protein dynamics only. Based on deterministic considerations, it is often stated or tacitly assumed that such models are valid in the limit of rapid promoter switching. Here, starting from the chemical master equation describing promoter-protein interactions, mRNA transcription, protein translation, and decay, we prove that in the limit of fast promoter switching, the distribution of protein numbers is different than that given by standard stochastic models with Hill-type propensities. We show the differences are pronounced whenever the protein-DNA binding rate is much larger than the unbinding rate, a special case of fast promoter switching. Furthermore, we show using both theory and simulations that use of the standard stochastic models leads to drastically incorrect predictions for the switching properties of positive feedback loops and that these differences decrease with increasing mean protein burst size. Our results confirm that commonly used stochastic models of gene regulatory networks are only accurate in a subset of the parameter space consistent with rapid promoter switching.

The Influence of Molecular Reach and Diffusivity on the Efficacy of Membrane-Confined Reactions

Signaling by surface receptors often relies on tethered reactions whereby an enzyme bound to the cytoplasmic tail of a receptor catalyzes reactions on substrates within reach. The overall length and stiffness of the receptor tail, the enzyme, and the substrate determine a biophysical parameter termed the molecular reach of the reaction. This parameter determines the probability that the receptor-tethered enzyme will contact the substrate in the volume proximal to the membrane when separated by different distances within the membrane plane. In this work, we develop particle-based stochastic reaction-diffusion models to study the interplay between molecular reach and diffusion. We find that increasing the molecular reach can increase reaction efficacy for slowly diffusing receptors, whereas for rapidly diffusing receptors, increasing molecular reach reduces reaction efficacy. In contrast, if reactions are forced to take place within the two-dimensional plasma membrane instead of the three-dimensional volume proximal to it or if molecules diffuse in three dimensions, increasing molecular reach increases reaction efficacy for all diffusivities. We show results in the context of immune checkpoint receptors (PD-1 dephosphorylating CD28), a standard opposing kinase-phosphatase reaction, and a minimal two-particle model. The work highlights the importance of the three-dimensional nature of many two-dimensional membrane-confined interactions, illustrating a role for molecular reach in controlling biochemical reactions.

A Model for the Spatiotemporal Design of Gene Regulatory Circuits †

Mathematical modeling assists the design of synthetic regulatory networks by providing a detailed mechanistic understanding of biological systems. Models that can predict the performance of a design are fundamental for synthetic biology since they minimize iterations along the design-build-test lifecycle. Such predictability depends crucially on what assumptions (i.e., biological simplifications) the model considers. Here, we challenge a common assumption when it comes to the modeling of bacterial-based gene regulation: considering negligible the effects of intracellular physical space. It is commonly assumed that molecules, such as transcription factors (TF), are homogeneously distributed inside a cell, so there is no need to model their diffusion. We describe a mathematical model that accounts for molecular diffusion and show how simulations of network performance are decisively affected by the distance between its components. Specifically, the model focuses on the search by a TF for its target promoter. The combination of local searches, via one-dimensional sliding along the chromosome, and global searches, via three-dimensional diffusion through the cytoplasm, determine TF-promoter interplay. Previous experimental results with engineered bacteria in which the distance between TF source and target was minimized or enlarged were successfully reproduced by the spatially resolved model we introduce here. This suggests that the spatial specification of the circuit alone can be exploited as a design parameter in synthetic biology to select programmable output levels.

The Linear Noise Approximation for Spatially Dependent Biochemical Networks

An algorithm for computing the linear noise approximation (LNA) of the reaction-diffusion master equation (RDME) is developed and tested. The RDME is often used as a model for biochemical reaction networks. The LNA is derived for a general discretization of the spatial domain of the problem. If M is the number of chemical species in the network and N is the number of nodes in the discretization in space, then the computational work to determine approximations of the mean and the covariances of the probability distributions is proportional to [Formula: see text] in a straightforward implementation. In our LNA algorithm, the work is proportional to [Formula: see text]. Since N usually is larger than M, this is a significant reduction. The accuracy of the approximation in the algorithm is estimated analytically and evaluated in numerical experiments.

Spatial Stochastic Intracellular Kinetics: A Review of Modelling Approaches

Models of chemical kinetics that incorporate both stochasticity and diffusion are an increasingly common tool for studying biology. The variety of competing models is vast, but two stand out by virtue of their popularity: the reaction-diffusion master equation and Brownian dynamics. In this review, we critically address a number of open questions surrounding these models: How can they be justified physically? How do they relate to each other? How do they fit into the wider landscape of chemical models, ranging from the rate equations to molecular dynamics? This review assumes no prior knowledge of modelling chemical kinetics and should be accessible to a wide range of readers.

Sharing of Phosphatases Promotes Response Plasticity in Phosphorylation Cascades

Sharing of positive or negative regulators between multiple targets is frequently observed in cellular signaling cascades. For instance, phosphatase sharing between multiple kinases is ubiquitous within the MAPK pathway. Here we investigate how such phosphatase sharing could shape robustness and evolvability of the phosphorylation cascade. Through modeling and evolutionary simulations, we demonstrate that 1) phosphatase sharing dramatically increases robustness of a bistable MAPK response, and 2) phosphatase-sharing cascades evolve faster than nonsharing cascades. This faster evolution is particularly pronounced when evolving from a monostable toward a bistable phenotype, whereas the transition speed of a population from a bistable to monostable response is not affected by phosphatase sharing. This property may enable the phosphatase-sharing design to adapt better in a changing environment. Analysis of the respective mutational landscapes reveal that phosphatase sharing reduces the number of limiting mutations required for transition from monostable to bistable responses, hence facilitating a faster transition to such response types. Taken together, using MAPK cascade as an example, our study offers a general theoretical framework to explore robustness and evolutionary plasticity of signal transduction cascades.

Kinetics and equilibrium constants of oligonucleotides at low concentrations. Hybridization and melting study

Although DNA hybridization/melting is one of the most important biochemical reactions, the non-trivial kinetics of the process is not yet fully understood. In this work, we use Förster resonance energy transfer (FRET) to investigate the influence of temperature, ionic strength, and oligonucleotide length on the kinetic and equilibrium constants of DNA oligonucleotide binding and dissociation. We show that at low reagent concentrations and ionic strength, the time needed to establish equilibrium between single and double strand forms may be of the order of days, even for simple oligonucleotides of a length of 20 base pairs or less. We also identify and discuss the possible artifacts related to fluorescence-based experiments conducted in extremely dilute solutions. The results should prove useful for the judicious design of technologies based on DNA-matching, including sensors, DNA multiplication, sequencing, and gene manipulation.

eGFRD in all dimensions

Biochemical reactions often occur at low copy numbers but at once in crowded and diverse environments. Space and stochasticity therefore play an essential role in biochemical networks. Spatial-stochastic simulations have become a prominent tool for understanding how stochasticity at the microscopic level influences the macroscopic behavior of such systems. While particle-based models guarantee the level of detail necessary to accurately describe the microscopic dynamics at very low copy numbers, the algorithms used to simulate them typically imply trade-offs between computational efficiency and biochemical accuracy. eGFRD (enhanced Green's Function Reaction Dynamics) is an exact algorithm that evades such trade-offs by partitioning the N-particle system into M ≤ N analytically tractable one- and two-particle systems; the analytical solutions (Green's functions) then are used to implement an event-driven particle-based scheme that allows particles to make large jumps in time and space while retaining access to their state variables at arbitrary simulation times. Here we present "eGFRD2," a new eGFRD version that implements the principle of eGFRD in all dimensions, thus enabling efficient particle-based simulation of biochemical reaction-diffusion processes in the 3D cytoplasm, on 2D planes representing membranes, and on 1D elongated cylinders representative of, e.g., cytoskeletal tracks or DNA; in 1D, it also incorporates convective motion used to model active transport. We find that, for low particle densities, eGFRD2 is up to 6 orders of magnitude faster than conventional Brownian dynamics. We exemplify the capabilities of eGFRD2 by simulating an idealized model of Pom1 gradient formation, which involves 3D diffusion, active transport on microtubules, and autophosphorylation on the membrane, confirming recent experimental and theoretical results on this system to hold under genuinely stochastic conditions.

Unbiased Atomistic Insight into the Mechanisms and Solvent Role for Globular Protein Dimer Dissociation

Association and dissociation of proteins are fundamental processes in nature. Although simple to understand conceptually, the details of the underlying mechanisms and role of the solvent are poorly understood. Here, we investigate the dissociation of the hydrophilic β-lactoglobulin dimer by employing transition path sampling. Analysis of the sampled path ensembles reveals a variety of mechanisms: (1) a direct aligned dissociation (2) a hopping and rebinding transition followed by unbinding, and (3) a sliding transition before unbinding. Reaction coordinate and transition-state analysis predicts that, besides native contact and neighboring salt-bridge interactions, solvent degrees of freedom play an important role in the dissociation process. Bridging waters, hydrogen-bonded to both proteins, support contacts in the native state and nearby lying transition-state regions, whereas they exhibit faster dynamics in further lying transition-state regions, rendering the proteins more mobile and assisting in rebinding. Analysis of the structure and dynamics of the solvent molecules reveals that the dry native interface induces enhanced populations of both disordered hydration water near hydrophilic residues and tetrahedrally ordered hydration water nearby hydrophobic residues. Although not exhaustive, our sampling of rare unbiased reactive molecular dynamics trajectories enhances the understanding of protein dissociation via complex pathways including (multiple) rebinding events.