Catalyst: Fast and flexible modeling of reaction networks.

We introduce Catalyst.jl, a flexible and feature-filled Julia library for modeling and high-performance simulation of chemical reaction networks (CRNs). Catalyst supports simulating stochastic chemical kinetics (jump process), chemical Langevin equation (stochastic differential equation), and reaction rate equation (ordinary differential equation) representations for CRNs. Through comprehensive benchmarks, we demonstrate that Catalyst simulation runtimes are often one to two orders of magnitude faster than other popular tools. More broadly, Catalyst acts as both a domain-specific language and an intermediate representation for symbolically encoding CRN models as Julia-native objects. This enables a pipeline of symbolically specifying, analyzing, and modifying CRNs; converting Catalyst models to symbolic representations of concrete mathematical models; and generating compiled code for numerical solvers. Leveraging ModelingToolkit.jl and Symbolics.jl, Catalyst models can be analyzed, simplified, and compiled into optimized representations for use in numerical solvers. Finally, we demonstrate Catalyst's broad extensibility and composability by highlighting how it can compose with a variety of Julia libraries, and how existing open-source biological modeling projects have extended its intermediate representation.

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.

How Retroactivity Affects the Behavior of Incoherent Feedforward Loops.

An incoherent feedforward loop (IFFL) is a network motif known for its ability to accelerate responses and generate pulses. It remains an open question to understand the behavior of IFFLs in contexts with high levels of retroactivity, where an upstream transcription factor binds to numerous downstream binding sites. Here we study the behavior of IFFLs by simulating and comparing ODE models with different levels of retroactivity. We find that increasing retroactivity in an IFFL can increase, decrease, or keep the network's response time and pulse amplitude constant. This suggests that increasing retroactivity, traditionally considered an impediment to designing robust synthetic systems, could be exploited to improve the performance of IFFLs. In contrast, we find that increasing retroactivity in a negative autoregulated circuit can only slow the response. The ability of an IFFL to flexibly handle retroactivity may have contributed to its significant abundance in both bacterial and eukaryotic regulatory networks.

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 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.

Modeling Genetic Circuit Behavior in Transiently Transfected Mammalian Cells.

Binning cells by plasmid copy number is a common practice for analyzing transient transfection data. In many kinetic models of transfected cells, protein production rates are assumed to be proportional to plasmid copy number. The validity of this assumption in transiently transfected mammalian cells is not clear; models based on this assumption appear unable to reproduce experimental flow cytometry data robustly. We hypothesize that protein saturation at high plasmid copy number is a reason previous models break down and validate our hypothesis by comparing experimental data and a stochastic chemical kinetics model. The model demonstrates that there are multiple distinct physical mechanisms that can cause saturation. On the basis of these observations, we develop a novel minimal bin-dependent ODE model that assumes different parameters for protein production in cells with low versus high numbers of plasmids. Compared to a traditional Hill-function-based model, the bin-dependent model requires only one additional parameter, but fits flow cytometry input-output data for individual modules up to twice as accurately. By composing together models of individually fit modules, we use the bin-dependent model to predict the behavior of six cascades and three feed-forward circuits. The bin-dependent models are shown to provide more accurate predictions on average than corresponding (composed) Hill-function-based models and predictions of comparable accuracy to EQuIP, while still providing a minimal ODE-based model that should be easy to integrate as a subcomponent within larger differential equation circuit models. Our analysis also demonstrates that accounting for batch effects is important in developing accurate composed models.

Biophysical assay for tethered signaling reactions reveals tether-controlled activity for the phosphatase SHP-1.

Tethered enzymatic reactions are ubiquitous in signaling networks but are poorly understood. A previously unreported mathematical analysis is established for tethered signaling reactions in surface plasmon resonance (SPR). Applying the method to the phosphatase SHP-1 interacting with a phosphorylated tether corresponding to an immune receptor cytoplasmic tail provides five biophysical/biochemical constants from a single SPR experiment: two binding rates, two catalytic rates, and a reach parameter. Tether binding increases the activity of SHP-1 by 900-fold through a binding-induced allosteric activation (20-fold) and a more significant increase in local substrate concentration (45-fold). The reach parameter indicates that this local substrate concentration is exquisitely sensitive to receptor clustering. We further show that truncation of the tether leads not only to a lower reach but also to lower binding and catalysis. This work establishes a new framework for studying tethered signaling processes and highlights the tether as a control parameter in clustered receptor signaling.

Uniform asymptotic approximation of diffusion to a small target: Generalized reaction models.

The diffusion of a reactant to a binding target plays a key role in many biological processes. The reaction radius at which the reactant and target may interact is often a small parameter relative to the diameter of the domain in which the reactant diffuses. We develop uniform in time asymptotic expansions in the reaction radius of the full solution to the corresponding diffusion equations for two separate reactant-target interaction mechanisms: the Doi or volume reactivity model and the Smoluchowski-Collins-Kimball partial-absorption surface reactivity model. In the former, the reactant and target react with a fixed probability per unit time when within a specified separation. In the latter, upon reaching a fixed separation, they probabilistically react or the reactant reflects away from the target. Expansions of the solution to each model are constructed by projecting out the contribution of the first eigenvalue and eigenfunction to the solution of the diffusion equation and then developing matched asymptotic expansions in Laplace-transform space. Our approach offers an equivalent, but alternative, method to the pseudopotential approach we previously employed [Isaacson and Newby, Phys. Rev. E 88, 012820 (2013)PLEEE81539-375510.1103/PhysRevE.88.012820] for the simpler Smoluchowski pure-absorption reaction mechanism. We find that the resulting asymptotic expansions of the diffusion equation solutions are identical with the exception of one parameter: the diffusion-limited reaction rates of the Doi and partial-absorption models. This demonstrates that for biological systems in which the reaction radius is a small parameter, properly calibrated Doi and partial-absorption models may be functionally equivalent.

Imaging and characterizing cells using tomography.

We can learn much about cell function by imaging and quantifying sub-cellular structures, especially if this is done non-destructively without altering said structures. Soft X-ray tomography (SXT) is a high-resolution imaging technique for visualizing cells and their interior structure in 3D. A tomogram of the cell, reconstructed from a series of 2D projection images, can be easily segmented and analyzed. SXT has a very high specimen throughput compared to other high-resolution structure imaging modalities; for example, tomographic data for reconstructing an entire eukaryotic cell is acquired in a matter of minutes. SXT visualizes cells without the need for chemical fixation, dehydration, or staining of the specimen. As a result, the SXT reconstructions are close representations of cells in their native state. SXT is applicable to most cell types. The deep penetration of soft X-rays allows cells, even mammalian cells, to be imaged without being sectioned. Image contrast in SXT is generated by the differential attenuation soft X-ray illumination as it passes through the specimen. Accordingly, each voxel in the tomographic reconstruction has a measured linear absorption coefficient (LAC) value. LAC values are quantitative and give rise to each sub-cellular component having a characteristic LAC profile, allowing organelles to be identified and segmented from the milieu of other cell contents. In this chapter, we describe the fundamentals of SXT imaging and how this technique can answer real world questions in the study of the nucleus. We also describe the development of correlative methods for the localization of specific molecules in a SXT reconstruction. The combination of fluorescence and SXT data acquired from the same specimen produces composite 3D images, rich with detailed information on the inner workings of cells.

The influence of spatial variation in chromatin density determined by X-ray tomograms on the time to find DNA binding sites.

In this work, we examine how volume exclusion caused by regions of high chromatin density might influence the time required for proteins to find specific DNA binding sites. The spatial variation of chromatin density within mouse olfactory sensory neurons is determined from soft X-ray tomography reconstructions of five nuclei. We show that there is a division of the nuclear space into regions of low-density euchromatin and high-density heterochromatin. Volume exclusion experienced by a diffusing protein caused by this varying density of chromatin is modeled by a repulsive potential. The value of the potential at a given point in space is chosen to be proportional to the density of chromatin at that location. The constant of proportionality, called the volume exclusivity, provides a model parameter that determines the strength of volume exclusion. Numerical simulations demonstrate that the mean time for a protein to locate a binding site localized in euchromatin is minimized for a finite, nonzero volume exclusivity. For binding sites in heterochromatin, the mean time is minimized when the volume exclusivity is zero (the protein experiences no volume exclusion). An analytical theory is developed to explain these results. The theory suggests that for binding sites in euchromatin there is an optimal level of volume exclusivity that balances a reduction in the volume searched in finding the binding site, with the height of effective potential barriers the protein must cross during the search process.

Uniform asymptotic approximation of diffusion to a small target.

The problem of the time required for a diffusing molecule, within a large bounded domain, to first locate a small target is prevalent in biological modeling. Here we study this problem for a small spherical target. We develop uniform in time asymptotic expansions in the target radius of the solution to the corresponding diffusion equation. Our approach is based on combining expansions of a long-time approximation of the solution, involving the first eigenvalue and eigenfunction of the Laplacian, with expansions of a short-time correction calculated by a pseudopotential approximation. These expansions allow the calculation of corresponding expansions of the first passage time density for the diffusing molecule to find the target. We demonstrate the accuracy of our method in approximating the first passage time density and related statistics for the spherically symmetric problem where the domain is a large concentric sphere about a small target centered at the origin.

A convergent reaction-diffusion master equation.

The reaction-diffusion master equation (RDME) is a lattice stochastic reaction-diffusion model that has been used to study spatially distributed cellular processes. The RDME is often interpreted as an approximation to spatially continuous models in which molecules move by Brownian motion and react by one of several mechanisms when sufficiently close. In the limit that the lattice spacing approaches zero, in two or more dimensions, the RDME has been shown to lose bimolecular reactions. The RDME is therefore not a convergent approximation to any spatially continuous model that incorporates bimolecular reactions. In this work we derive a new convergent RDME (CRDME) by finite volume discretization of a spatially continuous stochastic reaction-diffusion model popularized by Doi. We demonstrate the numerical convergence of reaction time statistics associated with the CRDME. For sufficiently large lattice spacings or slow bimolecular reaction rates, we also show that the reaction time statistics of the CRDME may be approximated by those from the RDME. The original RDME may therefore be interpreted as an approximation to the CRDME in several asymptotic limits.

Reaction-diffusion master equation, diffusion-limited reactions, and singular potentials.

To model biochemical systems in which both noise in the chemical reaction process and spatial movement of molecules is important, both the reaction-diffusion master equation (RDME) and Smoluchowski diffusion-limited reaction (SDLR) partial differential equation (PDE) models have been used. In previous work we showed that the solution to the RDME may be interpreted as an asymptotic approximation in the reaction radius to the solution of the SDLR PDE [S. A. Isaacson, SIAM J. Appl. Math. 70, 77 (2009)]. The approximation was shown to be divergent in the limit that the lattice spacing in the RDME approached zero. In this work we expand upon these results for the special case of the two-molecule annihilation reaction, A+B-->Ø. We first introduce a third stochastic reaction-diffusion PDE model that incorporates a pseudopotential based bimolecular reaction mechanism. The solution to the pseudopotential model is then shown to be an asymptotic approximation to the solution of the SDLR PDE for small reaction radii. We next illustrate how the RDME may be obtained by a formal discretization of the pseudopotential model, motivating why the RDME is itself an asymptotic approximation of the SDLR PDE. Finally, we give a more detailed numerical analysis of the difference between solutions to the RDME and SDLR PDE models as a function of both the reaction-radius and the lattice spacing (in the RDME).