Identification of dynamic mass-action biochemical reaction networks using sparse Bayesian methods.

Identifying the reactions that govern a dynamical biological system is a crucial but challenging task in systems biology. In this work, we present a data-driven method to infer the underlying biochemical reaction system governing a set of observed species concentrations over time. We formulate the problem as a regression over a large, but limited, mass-action constrained reaction space and utilize sparse Bayesian inference via the regularized horseshoe prior to produce robust, interpretable biochemical reaction networks, along with uncertainty estimates of parameters. The resulting systems of chemical reactions and posteriors inform the biologist of potentially several reaction systems that can be further investigated. We demonstrate the method on two examples of recovering the dynamics of an unknown reaction system, to illustrate the benefits of improved accuracy and information obtained.

Interpretable polynomial neural ordinary differential equations.

Neural networks have the ability to serve as universal function approximators, but they are not interpretable and do not generalize well outside of their training region. Both of these issues are problematic when trying to apply standard neural ordinary differential equations (ODEs) to dynamical systems. We introduce the polynomial neural ODE, which is a deep polynomial neural network inside of the neural ODE framework. We demonstrate the capability of polynomial neural ODEs to predict outside of the training region, as well as to perform direct symbolic regression without using additional tools such as SINDy.

Epidemiological modeling in StochSS Live!

SUMMARY: We present StochSS Live!, a web-based service for modeling, simulation and analysis of a wide range of mathematical, biological and biochemical systems. Using an epidemiological model of COVID-19, we demonstrate the power of StochSS Live! to enable researchers to quickly develop a deterministic or a discrete stochastic model, infer its parameters and analyze the results. AVAILABILITY AND IMPLEMENTATION: StochSS Live! is freely available at https://live.stochss.org/. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Control over single-cell distribution of G1 lengths by WNT governs pluripotency.

The link between single-cell variation and population-level fate choices lacks a mechanistic explanation despite extensive observations of gene expression and epigenetic variation among individual cells. Here, we found that single human embryonic stem cells (hESCs) have different and biased differentiation potentials toward either neuroectoderm or mesendoderm depending on their G1 lengths before the onset of differentiation. Single-cell variation in G1 length operates in a dynamic equilibrium that establishes a G1 length probability distribution for a population of hESCs and predicts differentiation outcome toward neuroectoderm or mesendoderm lineages. Although sister stem cells generally share G1 lengths, a variable proportion of cells have asymmetric G1 lengths, which maintains the population dispersion. Environmental Wingless-INT (WNT) levels can control the G1 length distribution, apparently as a means of priming the fate of hESC populations once they undergo differentiation. As a downstream mechanism, global 5-hydroxymethylcytosine levels are regulated by G1 length and thereby link G1 length to differentiation outcomes of hESCs. Overall, our findings suggest that intrapopulation heterogeneity in G1 length underlies the pluripotent differentiation potential of stem cell populations.

Coordinating cell polarization and morphogenesis through mechanical feedback.

Many cellular processes require cell polarization to be maintained as the cell changes shape, grows or moves. Without feedback mechanisms relaying information about cell shape to the polarity molecular machinery, the coordination between cell polarization and morphogenesis, movement or growth would not be possible. Here we theoretically and computationally study the role of a genetically-encoded mechanical feedback (in the Cell Wall Integrity pathway) as a potential coordination mechanism between cell morphogenesis and polarity during budding yeast mating projection growth. We developed a coarse-grained continuum description of the coupled dynamics of cell polarization and morphogenesis as well as 3D stochastic simulations of the molecular polarization machinery in the evolving cell shape. Both theoretical approaches show that in the absence of mechanical feedback (or in the presence of weak feedback), cell polarity cannot be maintained at the projection tip during growth, with the polarization cap wandering off the projection tip, arresting morphogenesis. In contrast, for mechanical feedback strengths above a threshold, cells can robustly maintain cell polarization at the tip and simultaneously sustain mating projection growth. These results indicate that the mechanical feedback encoded in the Cell Wall Integrity pathway can provide important positional information to the molecular machinery in the cell, thereby enabling the coordination of cell polarization and morphogenesis.

Dynamic Phase Extraction: Applications in Pulse Rate Variability.

Pulse rate variability is a physiological parameter that has been extensively studied and correlated with many physical ailments. However, the phase relationship between inter-beat interval, IBI, and breathing has very rarely been studied. Develop a technique by which the phase relationship between IBI and breathing can be accurately and efficiently extracted from photoplethysmography (PPG) data. A program based on Lock-in Amplifier technology was written in Python to implement a novel technique, Dynamic Phase Extraction. It was tested using a breath pacer and a PPG sensor on 6 subjects who followed a breath pacer at varied breathing rates. The data were then analyzed using both traditional methods and the novel technique (Dynamic Phase Extraction) utilizing a breath pacer. Pulse data was extracted using a PPG sensor. Dynamic Phase Extraction (DPE) gave the magnitudes of the variation in IBI associated with breathing [Formula: see text] measured with photoplethysmography during paced breathing (with premature ventricular contractions, abnormal arrhythmias, and other artifacts edited out). [Formula: see text] correlated well with two standard measures of pulse rate variability: the Standard Deviation of the inter-beat interval (SDNN) (ρ = 0.911) and with the integrated value of the Power Spectral Density between 0.04 and 0.15 Hz (Low Frequency Power or LF Power) (ρ = 0.885). These correlations were comparable to the correlation between the SDNN and the LF Power (ρ = 0.877). In addition to the magnitude [Formula: see text], Dynamic Phase Extraction also gave the phase between the breath pacer and the changes in the inter-beat interval (IBI) due to respiratory sinus arrythmia (RSA), and correlated well with the phase extracted using a Fourier transform (ρ = 0.857). Dynamic Phase Extraction can extract both the phase between the breath pacer and the changes in IBI due to the respiratory sinus arrhythmia component of pulse rate variability ([Formula: see text], but is limited by needing a breath pacer.

Accelerated regression-based summary statistics for discrete stochastic systems via approximate simulators.

BACKGROUND: Approximate Bayesian Computation (ABC) has become a key tool for calibrating the parameters of discrete stochastic biochemical models. For higher dimensional models and data, its performance is strongly dependent on having a representative set of summary statistics. While regression-based methods have been demonstrated to allow for the automatic construction of effective summary statistics, their reliance on first simulating a large training set creates a significant overhead when applying these methods to discrete stochastic models for which simulation is relatively expensive. In this τ work, we present a method to reduce this computational burden by leveraging approximate simulators of these systems, such as ordinary differential equations and τ-Leaping approximations. RESULTS: We have developed an algorithm to accelerate the construction of regression-based summary statistics for Approximate Bayesian Computation by selectively using the faster approximate algorithms for simulations. By posing the problem as one of ratio estimation, we use state-of-the-art methods in machine learning to show that, in many cases, our algorithm can significantly reduce the number of simulations from the full resolution model at a minimal cost to accuracy and little additional tuning from the user. We demonstrate the usefulness and robustness of our method with four different experiments. CONCLUSIONS: We provide a novel algorithm for accelerating the construction of summary statistics for stochastic biochemical systems. Compared to the standard practice of exclusively training from exact simulator samples, our method is able to dramatically reduce the number of required calls to the stochastic simulator at a minimal loss in accuracy. This can immediately be implemented to increase the overall speed of the ABC workflow for estimating parameters in complex systems.

Associations of longitudinal D-Dimer and Factor II on early trauma survival risk.

BACKGROUND: Trauma-induced coagulopathy (TIC) is a disorder that occurs in one-third of severely injured trauma patients, manifesting as increased bleeding and a 4X risk of mortality. Understanding the mechanisms driving TIC, clinical risk factors are essential to mitigating this coagulopathic bleeding and is therefore essential for saving lives. In this retrospective, single hospital study of 891 trauma patients, we investigate and quantify how two prominently described phenotypes of TIC, consumptive coagulopathy and hyperfibrinolysis, affect survival odds in the first 25 h, when deaths from TIC are most prevalent. METHODS: We employ a joint survival model to estimate the longitudinal trajectories of the protein Factor II (% activity) and the log of the protein fragment D-Dimer ([Formula: see text]g/ml), representative biomarkers of consumptive coagulopathy and hyperfibrinolysis respectively, and tie them together with patient outcomes. Joint models have recently gained popularity in medical studies due to the necessity to simultaneously track continuously measured biomarkers as a disease evolves, as well as to associate them with patient outcomes. In this work, we estimate and analyze our joint model using Bayesian methods to obtain uncertainties and distributions over associations and trajectories. RESULTS: We find that a unit increase in log D-Dimer increases the risk of mortality by 2.22 [1.57, 3.28] fold while a unit increase in Factor II only marginally decreases the risk of mortality by 0.94 [0.91,0.96] fold. This suggests that, while managing consumptive coagulopathy and hyperfibrinolysis both seem to affect survival odds, the effect of hyperfibrinolysis is much greater and more sensitive. Furthermore, we find that the longitudinal trajectories, controlling for many fixed covariates, trend differently for different patients. Thus, a more personalized approach is necessary when considering treatment and risk prediction under these phenotypes. CONCLUSION: This study reinforces the finding that hyperfibrinolysis is linked with poor patient outcomes regardless of factor consumption levels. Furthermore, it quantifies the degree to which measured D-Dimer levels correlate with increased risk. The single hospital, retrospective nature can be understood to specify the results to this particular hospital's patients and protocol in treating trauma patients. Expanding to a multi-hospital setting would result in better estimates about the underlying nature of consumptive coagulopathy and hyperfibrinolysis with survival, regardless of protocol. Individual trajectories obtained with these estimates can be used to provide personalized dynamic risk prediction when making decisions regarding management of blood factors.

An arbitrary Lagrangian Eulerian smoothed particle hydrodynamics (ALE-SPH) method with a boundary volume fraction formulation for fluid-structure interaction.

We present a new weakly-compressible smoothed particle hydrodynamics (SPH) method capable of modeling non-slip fixed and moving wall boundary conditions. The formulation combines a boundary volume fraction (BVF) wall approach with the transport-velocity SPH method. The resulting method, named SPH-BVF, offers detection of arbitrarily shaped solid walls on-the-fly, with small computational overhead due to its local formulation. This simple framework is capable of solving problems that are difficult or infeasible for standard SPH, namely flows subject to large shear stresses or at moderate Reynolds numbers, and mass transfer in deformable boundaries. In addition, the method extends the transport-velocity formulation to reaction-diffusion transport of mass in Newtonian fluids and linear elastic solids, which is common in biological structures. Taken together, the SPH-BVF method provides a good balance of simplicity and versatility, while avoiding some of the standard obstacles associated with SPH: particle penetration at the boundaries, tension instabilities and anisotropic particle alignments, that hamper SPH from being applied to complex problems such as fluid-structure interaction in a biological system.

Ontogeny of Circadian Rhythms and Synchrony in the Suprachiasmatic Nucleus.

In mammals, the suprachiasmatic nucleus (SCN) of the hypothalamus coordinates daily rhythms including sleep-wake, hormone release, and gene expression. The cells of the SCN must synchronize to each other to drive these circadian rhythms in the rest of the body. The ontogeny of circadian cycling and intercellular coupling in the SCN remains poorly understood. Recent in vitro studies have recorded circadian rhythms from the whole embryonic SCN. Here, we tracked the onset and precision of rhythms in PERIOD2 (PER2), a clock protein, within the SCN isolated from embryonic and postnatal mice of undetermined sex. We found that a few SCN cells developed circadian periodicity in PER2 by 14.5 d after mating (E14.5) with no evidence for daily cycling on E13.5. On E15.5, the fraction of competent oscillators increased dramatically corresponding with stabilization of their circadian periods. The cells of the SCN harvested at E15.5 expressed sustained, synchronous daily rhythms. By postnatal day 2 (P2), SCN oscillators displayed the daily, dorsal-ventral phase wave in clock gene expression typical of the adult SCN. Strikingly, vasoactive intestinal polypeptide (VIP), a neuropeptide critical for synchrony in the adult SCN, and its receptor, VPAC2R, reached detectable levels after birth and after the onset of circadian synchrony. Antagonists of GABA or VIP signaling or action potentials did not disrupt circadian synchrony in the E15.5 SCN. We conclude that endogenous daily rhythms in the fetal SCN begin with few noisy oscillators on E14.5, followed by widespread oscillations that rapidly synchronize on E15.5 by an unknown mechanism.SIGNIFICANCE STATEMENT We recorded the onset of PER2 circadian oscillations during embryonic development in the mouse SCN. When isolated at E13.5, the anlagen of the SCN expresses high, arrhythmic PER2. In contrast, a few cells show noisy circadian rhythms in the isolated E14.5 SCN and most show reliable, self-sustained, synchronized rhythms in the E15.5 SCN. Strikingly, this synchrony at E15.5 appears before expression of VIP or its receptor and persists in the presence of blockers of VIP, GABA or neuronal firing. Finally, the dorsal-ventral phase wave of PER2 typical of the adult SCN appears ∼P2, indicating that multiple signals may mediate circadian synchrony during the ontogeny of the SCN.

Mechanical feedback coordinates cell wall expansion and assembly in yeast mating morphogenesis.

The shaping of individual cells requires a tight coordination of cell mechanics and growth. However, it is unclear how information about the mechanical state of the wall is relayed to the molecular processes building it, thereby enabling the coordination of cell wall expansion and assembly during morphogenesis. Combining theoretical and experimental approaches, we show that a mechanical feedback coordinating cell wall assembly and expansion is essential to sustain mating projection growth in budding yeast (Saccharomyces cerevisiae). Our theoretical results indicate that the mechanical feedback provided by the Cell Wall Integrity pathway, with cell wall stress sensors Wsc1 and Mid2 increasingly activating membrane-localized cell wall synthases Fks1/2 upon faster cell wall expansion, stabilizes mating projection growth without affecting cell shape. Experimental perturbation of the osmotic pressure and cell wall mechanics, as well as compromising the mechanical feedback through genetic deletion of the stress sensors, leads to cellular phenotypes that support the theoretical predictions. Our results indicate that while the existence of mechanical feedback is essential to stabilize mating projection growth, the shape and size of the cell are insensitive to the feedback.

The effect of cell geometry on polarization in budding yeast.

The localization (or polarization) of proteins on the membrane during the mating of budding yeast (Saccharomyces cerevisiae) is an important model system for understanding simple pattern formation within cells. While there are many existing mathematical models of polarization, for both budding and mating, there are still many aspects of this process that are not well understood. In this paper we set out to elucidate the effect that the geometry of the cell can have on the dynamics of certain models of polarization. Specifically, we look at several spatial stochastic models of Cdc42 polarization that have been adapted from published models, on a variety of tip-shaped geometries, to replicate the shape change that occurs during the growth of the mating projection. We show here that there is a complex interplay between the dynamics of polarization and the shape of the cell. Our results show that while models of polarization can generate a stable polarization cap, its localization at the tip of mating projections is unstable, with the polarization cap drifting away from the tip of the projection in a geometry dependent manner. We also compare predictions from our computational results to experiments that observe cells with projections of varying lengths, and track the stability of the polarization cap. Lastly, we examine one model of actin polarization and show that it is unlikely, at least for the models studied here, that actin dynamics and vesicle traffic are able to overcome this effect of geometry.

Analysis of the role of thrombomodulin in all-trans retinoic acid treatment of coagulation disorders in cancer patients.

BACKGROUND: Clinical studies have shown that all-trans retinoic acid (RA), which is often used in treatment of cancer patients, improves hemostatic parameters and bleeding complications such as disseminated intravascular coagulation (DIC). However, the mechanisms underlying this improvement have yet to be elucidated. In vitro studies have reported that RA upregulates thrombomodulin (TM) expression on the endothelial cell surface. The objective of this study was to investigate how and to what extent the TM concentration changes after RA treatment in cancer patients, and how this variation influences the blood coagulation cascade. RESULTS: In this study, we introduced an ordinary differential equation (ODE) model of gene expression for the RA-induced upregulation of TM concentration. Coupling the gene expression model with a two-compartment pharmacokinetic model of RA, we obtained the time-dependent changes in TM and thrombomodulin-mRNA (TMR) concentrations following oral administration of RA. Our results indicated that the TM concentration reached its peak level almost 14 h after taking a single oral dose (110 [Formula: see text]) of RA. Continuous treatment with RA resulted in oscillatory expression of TM on the endothelial cell surface. We then coupled the gene expression model with a mechanistic model of the coagulation cascade, and showed that the elevated levels of TM over the course of RA therapy with a single daily oral dose (110 [Formula: see text]) of RA, reduced the peak thrombin levels and endogenous thrombin potential (ETP) up to 50 and 49%, respectively. We showed that progressive reductions in plasma levels of RA, observed in continuous RA therapy with a once-daily oral dose (110 [Formula: see text]) of RA, did not affect TM-mediated reduction of thrombin generation significantly. This finding prompts the hypothesis that continuous RA treatment has more consistent therapeutic effects on coagulation disorders than on cancer. CONCLUSIONS: Our results indicate that the oscillatory upregulation of TM expression on the endothelial cells over the course of RA therapy could potentially contribute to the treatment of coagulation abnormalities in cancer patients. Further studies on the impacts of RA therapy on the procoagulant activity of cancer cells are needed to better elucidate the mechanisms by which RA therapy improves hemostatic abnormalities in cancer.

Functional network inference of the suprachiasmatic nucleus.

In the mammalian suprachiasmatic nucleus (SCN), noisy cellular oscillators communicate within a neuronal network to generate precise system-wide circadian rhythms. Although the intracellular genetic oscillator and intercellular biochemical coupling mechanisms have been examined previously, the network topology driving synchronization of the SCN has not been elucidated. This network has been particularly challenging to probe, due to its oscillatory components and slow coupling timescale. In this work, we investigated the SCN network at a single-cell resolution through a chemically induced desynchronization. We then inferred functional connections in the SCN by applying the maximal information coefficient statistic to bioluminescence reporter data from individual neurons while they resynchronized their circadian cycling. Our results demonstrate that the functional network of circadian cells associated with resynchronization has small-world characteristics, with a node degree distribution that is exponential. We show that hubs of this small-world network are preferentially located in the central SCN, with sparsely connected shells surrounding these cores. Finally, we used two computational models of circadian neurons to validate our predictions of network structure.

Identification of influential proteins in the classical retinoic acid signaling pathway.

BACKGROUND: In the classical pathway of retinoic acid (RA) mediated gene transcription, RA binds to a nuclear hormone receptor dimer composed of retinoic acid receptor (RAR) and retinoid X receptor (RXR), to induce the expression of its downstream target genes. In addition to nuclear receptors, there are other intracellular RA binding proteins such as cellular retinoic acid binding proteins (CRABP1 and CRABP2) and cytochrome P450 (CYP) enzymes, whose contributions to the RA signaling pathway have not been fully understood. The objective of this study was to compare the significance of various RA binding receptors, i.e. CRABP1, CRABP2, CYP and RAR in the RA signaling pathway. In this regard, we developed a mathematical model of the RA pathway, which is one of the few models, if not the only one, that includes all main intracellular RA binding receptors. We then performed a global sensitivity analysis (GSA) to investigate the contribution of the RA receptors to RA-induced mRNA production, when the cells were treated with a wide range of RA levels, from physiological to pharmacological concentrations. RESULTS: Our results show that CRABP2 and RAR are the most and the least important proteins, respectively, in controlling the model performance at physiological concentrations of RA (1-10 nM). However, at higher concentrations of RA, CYP and RAR are the most sensitive parameters of the system. Furthermore, we found that depending on the concentrations of all RA binding proteins, the rate of metabolism of RA can either change or remain constant following RA therapy. The cellular levels of CRABP1 are more important than that of CRABP2 in controlling RA metabolite formation at pharmacological conditions (RA = 0.1-1 µM). Finally, our results indicate a significant negative correlation between total mRNA production and total RA metabolite formation at pharmacological levels of RA. CONCLUSIONS: Our simulations indicate that the significance of the RA binding proteins in the RA pathway of gene expression strongly depends on intracellular concentration of RA. This study not only can explain why various cell types respond to RA therapy differently, but also can potentially help develop pharmacological methods to increase the efficacy of the drug.

Bayesian inference of elastic properties with resonant ultrasound spectroscopy.

Bayesian modeling and Hamiltonian Monte Carlo (HMC) are utilized to formulate a robust algorithm capable of simultaneously estimating anisotropic elastic properties and crystallographic orientation of a specimen from a list of measured resonance frequencies collected via Resonance Ultrasound Spectroscopy (RUS). Unlike typical optimization procedures which yield point estimates of the unknown parameters, computing a Bayesian posterior yields probability distributions for the unknown parameters, and HMC is an efficient way to compute this posterior. The algorithms described are demonstrated on RUS data collected from two parallelepiped specimens of structural metal alloys. First, the elastic constants for a specimen of fine-grain polycrystalline Ti-6Al-4 V with random crystallographic texture and isotropic elastic symmetry are estimated. Second, the elastic constants and crystallographic orientation for a single crystal Ni-based superalloy CMSX-4 specimen are accurately determined, using only measurements of the specimen geometry, mass, and resonance frequencies. The unique contributions of this paper are as follows: the application of HMC for sampling the Bayesian posterior of a probabilistic RUS model, and the procedure for simultaneous estimation of elastic constants and lattice-specimen misorientation. Compared to previous approaches these algorithms demonstrate superior convergence behavior, particularly when the initial parameterization is unknown, and enable substantially simplified experimental procedures.

Macromolecular Crowding Regulates the Gene Expression Profile by Limiting Diffusion.

We seek to elucidate the role of macromolecular crowding in transcription and translation. It is well known that stochasticity in gene expression can lead to differential gene expression and heterogeneity in a cell population. Recent experimental observations by Tan et al. have improved our understanding of the functional role of macromolecular crowding. It can be inferred from their observations that macromolecular crowding can lead to robustness in gene expression, resulting in a more homogeneous cell population. We introduce a spatial stochastic model to provide insight into this process. Our results show that macromolecular crowding reduces noise (as measured by the kurtosis of the mRNA distribution) in a cell population by limiting the diffusion of transcription factors (i.e. removing the unstable intermediate states), and that crowding by large molecules reduces noise more efficiently than crowding by small molecules. Finally, our simulation results provide evidence that the local variation in chromatin density as well as the total volume exclusion of the chromatin in the nucleus can induce a homogenous cell population.

Stochastic Simulation Service: Bridging the Gap between the Computational Expert and the Biologist.

We present StochSS: Stochastic Simulation as a Service, an integrated development environment for modeling and simulation of both deterministic and discrete stochastic biochemical systems in up to three dimensions. An easy to use graphical user interface enables researchers to quickly develop and simulate a biological model on a desktop or laptop, which can then be expanded to incorporate increasing levels of complexity. StochSS features state-of-the-art simulation engines. As the demand for computational power increases, StochSS can seamlessly scale computing resources in the cloud. In addition, StochSS can be deployed as a multi-user software environment where collaborators share computational resources and exchange models via a public model repository. We demonstrate the capabilities and ease of use of StochSS with an example of model development and simulation at increasing levels of complexity.

Reaction rates for reaction-diffusion kinetics on unstructured meshes.

The reaction-diffusion master equation is a stochastic model often utilized in the study of biochemical reaction networks in living cells. It is applied when the spatial distribution of molecules is important to the dynamics of the system. A viable approach to resolve the complex geometry of cells accurately is to discretize space with an unstructured mesh. Diffusion is modeled as discrete jumps between nodes on the mesh, and the diffusion jump rates can be obtained through a discretization of the diffusion equation on the mesh. Reactions can occur when molecules occupy the same voxel. In this paper, we develop a method for computing accurate reaction rates between molecules occupying the same voxel in an unstructured mesh. For large voxels, these rates are known to be well approximated by the reaction rates derived by Collins and Kimball, but as the mesh is refined, no analytical expression for the rates exists. We reduce the problem of computing accurate reaction rates to a pure preprocessing step, depending only on the mesh and not on the model parameters, and we devise an efficient numerical scheme to estimate them to high accuracy. We show in several numerical examples that as we refine the mesh, the results obtained with the reaction-diffusion master equation approach those of a more fine-grained Smoluchowski particle-tracking model.

Mesoscopic-microscopic spatial stochastic simulation with automatic system partitioning.

The reaction-diffusion master equation (RDME) is a model that allows for efficient on-lattice simulation of spatially resolved stochastic chemical kinetics. Compared to off-lattice hard-sphere simulations with Brownian dynamics or Green's function reaction dynamics, the RDME can be orders of magnitude faster if the lattice spacing can be chosen coarse enough. However, strongly diffusion-controlled reactions mandate a very fine mesh resolution for acceptable accuracy. It is common that reactions in the same model differ in their degree of diffusion control and therefore require different degrees of mesh resolution. This renders mesoscopic simulation inefficient for systems with multiscale properties. Mesoscopic-microscopic hybrid methods address this problem by resolving the most challenging reactions with a microscale, off-lattice simulation. However, all methods to date require manual partitioning of a system, effectively limiting their usefulness as "black-box" simulation codes. In this paper, we propose a hybrid simulation algorithm with automatic system partitioning based on indirect a priori error estimates. We demonstrate the accuracy and efficiency of the method on models of diffusion-controlled networks in 3D.