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

## Abstract

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.

## INTRODUCTION

Twenty-first century cell biology has been transformed by rapid development of new technologies that have delivered our field into an era where scientists can now easily generate vast amounts of quantitative data, providing a broad and comprehensive view of problems that were previously accessible only by brief, laboriously achieved glances through tiny chinks. Northern blots have been replaced by RNA sequencing, Western blots are being replaced by proteomics, and the central cell biological tool of imaging has been revolutionized by a wealth of dramatic improvements in labeling methods, optical design, and digital imaging. While currently available data are just a tiny fraction of the amount that will eventually be needed to model entire cells, sufficient data are available for a number of cell biological processes to ask whether and how we can use it to understand the biological mechanisms underlying those processes.

Mathematical modeling can help (Cohen, 2004; Gunawardena, 2014). By generating a model of a process that we want to study, where we are able to define and control all the inputs and parameters, we can directly determine whether a mechanism we have hypothesized is actually sufficient to explain the phenomenon we have observed. While mathematical modeling of cell biological processes can take many forms, here we will focus exclusively on computational simulation. One of our goals is to discuss the current state of the art in computational cell biology and remaining open challenges, particularly regarding the spatial organization of signaling processes, the inclusion of stochastic effects, and the multiscale nature of many cell biological processes.

Computers are much better than humans at keeping track of large, complex systems with many interacting parts and also at mercilessly following predetermined rules. One subfield of cell biology that has been able to make great use of these abilities is the study of signal transduction, where typically many different molecular species interact with one another in branching and reticular networks. Biochemical and genetic experiments have been able to map out and characterize many individual pairwise interactions in multicomponent signaling pathways. However, the resulting “wiring diagrams” have provided little direct insight into how these networks of molecular interactions could give rise to the diverse and fascinating outputs of these systems, which might be able to generate oscillations, high-pass or low-pass filtering of receptor-mediated input signals, conversion of analog signals into switch-like binary outputs, and so forth. Computational simulations have been instrumental in bringing order and insight into this tangled web so that now it is possible to recognize recurrent motifs in the design of signal transduction systems and, sometimes, accurately predict cellular responses to external stimuli (Eungdamrong and Iyengar, 2004; Kestler *et al.*, 2008; Janes and Lauffenburger, 2013; Cao *et al.*, 2016).

Within this context, a very fertile ecosystem of computational tools for modeling cellular biochemistry has flourished (Bartocci and Lio, 2016).There are now also approaches that enable scientists who are not computational experts themselves to translate their hypotheses about cellular signaling mechanisms into formal models in compact and intuitive ways; these include using textual rules (Faeder, 2011; Maus *et al.*, 2011; Tiger *et al.*, 2012; Harris *et al.*, 2016; Boutillier *et al.*, 2018) or iconographic symbols (Zhang *et al.*, 2013; Schaff *et al.*, 2016; Sekar *et al.*, 2017) to specify molecular interactions. Based on those specifications, the tools generate the resulting computational representations of the signaling networks and allow modelers to easily modify their assumptions to explore the consequences of such simulated manipulations on cellular behavior (Lopez *et al.*, 2013). To allow for tool-independent formulation and sharing of computational models, the widely used “systems biology markup language” SBML was developed. (Hucka *et al.*, 2003) SBML is continuously evolving and is supported by a large number of software systems for simulation and data analysis (Keating *et al.*, 2020; SBML.org). While there are many variations, generally simulations of this kind are able to keep track of concentrations and interactions of many individual molecular species as they change over time, often the most interesting dimension for the study of signal transduction.

However, many cell biological processes, including some kinds of signaling, cannot be analyzed without the notion of space. Cell–cell communication is frequently based on the exchange of soluble messenger molecules, such as hormones, cytokines or chemokines, that diffuse through extracellular space before being captured by specific receptors at particular locations on cellular membranes. Direct cell–cell contacts also typically involve only a few of the receptors on a cell and generate localized signals that activate cascades of protein interactions and modifications to propagate from the membrane into the cytoplasm. Spatial simulation of cell biological phenomena is not new; in 1952, Hodgkin and Huxley simulated the propagation of an action potential down a neuronal axon (Hodgkin and Huxley, 1952; Hellander *et al.*, 2015), and in the same year, Alan Turing used computational simulation to demonstrate how chemical systems featuring both diffusion and reaction could generate regular spatial patterns from an initial uniform state (Turing, 1952). More recently, many researchers have developed computational simulations that employ state-of-the-art knowledge about the properties and interactions of individual molecular components to attempt spatially resolved reproduction of complex cell biological phenomena. Reaction–diffusion models of intracellular biochemistry have been used to explore a variety of cellular symmetry-breaking processes, including the establishment of cell polarity (Jilkine and Edelstein-Keshet, 2011). Although the exact mechanistic details vary, spatial simulation approaches have yielded insights for symmetry-breaking systems as diverse as yeast bud site selection (Wedlich-Soldner *et al.*, 2003), asymmetric cell division in early *Caenorhabditis elegans* embryos (Dawes and Munro, 2011) and neutrophil chemotaxis (Onsum and Rao, 2007). Establishment of spatial gradients that determine cell fate has been explored in cells ranging from giant syncytial *Drosophila* embryos (Gregor *et al.*, 2007) to tiny individual bacteria (Chen *et al.*, 2011). Simulations that explicitly consider spatial effects as one cell communicates with its neighbors have been used to understand the formation of regular stripes in *Drosophila* embryos (von Dassow *et al.*, 2000) and bizarre noncell-autonomous effects in the patterning of wing bristles in adult flies (Amonlirdviman *et al.*, 2005).

In many of the modeling efforts that we have mentioned so far, the computational model and/or simulation was formulated explicitly for the problem at hand. While this approach has enabled important scientific insights, we believe that spatial simulation for cell biological processes can become a much more widely used tool in the cell biologist’s toolbox if there was more general access to user-friendly implementations of general spatial modeling frameworks that do not require extensive computational expertise. Consider, for example, the wide variety of user-friendly open-source software packages now available for analysis of sequencing data (Rice *et al.*, 2000; Trapnell *et al.*, 2012). One particularly important benefit of more standardized approaches to spatial simulations in cell biology is that standardization may help to resolve whether conflicting conclusions arise because of fundamental scientific differences in model assumptions or because of details of numerical implementation.

Another aspect gaining importance as we zoom in closer on the building blocks of cellular structures is the fact that, at the molecular level, cellular biochemistry is governed by stochastic processes such as thermal Brownian motion and the collisions and interactions among individual particles (Schnoerr *et al.*, 2017). Only in the limit of high concentrations and homogenous spatial distributions can the behavior of the molecular components of cellular signaling pathways be described in terms of deterministic reaction rate equations. Many subcellular mechanisms operate far from this limit either due to highly nonhomogenous clustering of receptors and of the signaling components they recruit, such as studied in the MAPK signaling pathway (Takahashi *et al.*, 2010), or due to locally low copy numbers, for instance of multimolecular complexes regulating transcription in the nucleus (Cho *et al.*, 2018). To accurately capture the stochastic characteristics of such processes, computational models have to simulate the motion and interactions of individual molecules. A full simulation of Brownian dynamics (BD), following each single “Brownian hop” of all molecules of a cellular region would, however, in most cases be too computationally expensive and time consuming. Moreover, by choosing a particular time step for recreating Brownian hops on the computer, we would impose this timescale on our simulations, missing events, such as molecular encounters, that may occur ‘in between’ our time steps. Many approaches have been developed to deal with this problem and we will discuss several of them below. A common theme among all of them is, however, saving computational cost through temporal and spatial coarse graining without sacrificing too much accuracy in the simulation results. This challenge unites researchers looking at spatially resolved modeling (without a focus on stochastic effects) and those who try to capture the manifestations of stochastic fluctuations in cellular systems.

The question of how the field could best go about building and sharing broadly applicable computational tools for spatial and stochastic modeling of cell biological processes was the focus of our working group, organized by J.R.F. and R.F.M., which met with the support of the National Institute for Mathematical and Biological Synthesis (NIMBioS). We chose to focus specifically on the biochemical scale of molecular interactions. We did not attempt to include the enormous field of molecular dynamics (MD) simulations that explores forces and movements of individual atoms within proteins or other macromolecules. The reason is that, because of their intensive computational demands, MD simulations are currently limited to exploration of very small biological systems (a few macromolecules) over very short periods of time (typically in the microsecond range or below), too small and too fast to be incorporated into cell-scale computational simulations. Conversely, we also limited ourselves to considering simulations of interactions within systems of molecular complexes with specific stoichiometry rather than extending our analysis to mesoscopic-scale models that abstract the behaviors of complex molecular systems into continuum physical descriptions, such as those describing cytoskeletal filaments as elastic beams (Nedelec and Foethke, 2007; Odell and Foe, 2008) or the plasma membrane as a flexible thin film (Fowler *et al.*, 2016). While these kinds of models are enormously useful in cell biology, they rely on fundamental simplifying assumptions. In contrast, we are specifically interested in exploring whether detailed simulations of the behaviors and interactions among biomolecular complexes can succeed in predicting certain kinds of mesoscopic phenomena and hence may help determining under which conditions the simplifying assumptions are justified. As we will discuss in our conclusion, it will be an exciting future direction for the field of biological simulation when all these three levels of spatially resolved simulations can be seamlessly interconnected.

In this article, we will first briefly survey several existing approaches for spatial and stochastic cell simulations and then apply them to a series of “unit tests” and benchmark problems. The problems we chose are categorized to cover a variety of different aspects of spatially resolved and/or stochastic simulations of cellular behavior. Conceptionally, they are relatively simple and are meant to capture specific challenges related to essential aspects of biological processes. They must be accurately modeled by a simulation tool to ensure that that tool’s results are reproducible for the problem category covered by the unit tests or benchmark problems. With these examples in hand, we then summarize how features of a particular cell biological problem should guide selection of the appropriate modeling approach. One should note that, in practice, finding or developing the appropriate approach for a given cell biological phenomenon involves far more decisions than “just” whether the model should capture spatial or stochastic aspects (or both). A similarly important and related challenge is to find the right balance between biological realism and detail on the one hand and computational manageability and parameter estimation on the other hand. However, a systematic exploration of methods for identifying the right degree of model granularity would go beyond the scope of an article exploring the roles of space and stochasticity. Here, we will, therefore, limit our discussion of topics such as model reduction and parameter estimation to pointing out where we think these aspects become particularly important. Finally, we will present the results of our wide-ranging and, sometimes, highly opinionated discussions on future directions and challenges for the field. We all share an ambitious vision of the future power of spatial cell simulations both for exploring hypotheses about mechanisms and for coming to grips with the massive amounts of quantitative data now available to cell biologists. Our overall goal here is to map out where the field currently stands and propose a trajectory for the future.

## MODELING APPROACHES

### Overview of approaches for spatial and stochastic modeling and simulation of molecular reactions underlying cell biological phenomena

Computational models that simulate the biochemistry underlying cell biological processes need to be able to describe molecular players and their reactions. However, depending on the particular question at hand, taking into account spatial aspects and stochastic effects (Figure 1) may or may not be essential, as we show in *Results*. The addition of spatial resolution is computationally demanding, and stochastic simulations are usually more costly than their deterministic counterparts. To attempt a whole-cell simulation, for example, one must choose whether more components and a more complex reaction network are necessary, generally requiring the sacrifice of spatial resolution (Tomita *et al.*, 1999; Sanghvi *et al.*, 2013), or if spatial resolution is necessary, then the reaction network must be simplified (Ghaemi *et al.*, 2020). If spatial resolution is a priority, must species be resolved as individual particles, capturing fluctuations in copy numbers but at considerable extra expense (Supplemental Table S3), or is an efficient deterministic approach sufficient?

Whether to include spatial or stochastic resolution will also affect the kind of data and (a priori known) parameters that are required and the results that can be expected. Differential equation-based models can incorporate phenomenological elements such as Hill-type functions bridging different model elements whose interdependencies are either not well understood or whose details are considered less important for the overall quality of a modeling effort. Moreover, more abstract models sometimes permit identifying components (species or mechanisms) whose kinetics contribute little to the behavior of a model or that can be lumped with other species to simplify its computational representation (Rao *et al.*, 2014). This can be particularly useful for very large systems with well-defined constraints, such as metabolic network models (Masid *et al.*, 2020). With decreasing complexity of a model, it also becomes easier to perform robust parameter estimations and to determine how well the model is justified based on the available data (Raue *et al.*, 2009). However, model abstractions are also at risk of losing the ability to describe the behavior of biological entities (for instance molecular signaling components) considered relevant by experimental biologists and, ultimately, the choice between model simplicity and biological realism can be difficult (Meier-Schellersheim *et al.*, 2019).

### Nonspatial modeling approaches

In many situations, we can describe a biochemical system adequately in terms of the overall concentrations of interacting molecule types and complexes (collectively called “species”), while neglecting the spatial variations in these concentrations. Reaction rate equations (see Box 1) describe how the species concentrations evolve in time. The terms in these equations arise from the rates of the reactions that can occur in the system, which are often described by the Law of Mass Action. The rate of reaction between two interacting species can be given by the product of their concentrations and a rate constant, for example, *k*_{on} for the ligand-receptor binding described in Box 1. The bimolecular rate constants that appear in these equations are sometimes referred to as “macroscopic” rate constants because they describe the average rate of reaction assuming homogenous distribution of the reacting species. In contrast, “microscopic” rate constants, which govern reaction kinetics at the scale of interacting particles, may take into account more details about the way the molecules approach each other, as discussed below.

#### BOX 1: reaction rate equations

Technically speaking, reaction rate equations are ordinary differential equations (ODEs). Here, the “ordinary” refers to the fact that they involve only time (as opposed to, for instance, time *and* space). To describe the time evolution of multiple interacting molecule types, one uses coupled differential equations that express how the components’ concentration changes are linked (or coupled). For applications, see, for example, Aldridge *et al.* (2006) and Tyson *et al.* (2003). From a numerical/mathematical point of view, ODEs describing biochemical reactions are typically simple, and many tools exist that can solve them to obtain the temporal evolution of the concentrations in ODE models.

Consider a simple model of a receptor binding to a ligand. We call *R* the concentration of the receptor, *L* that of the ligand, and *RL* that of the complex formed by the binding of the two. The rate equations giving the time derivatives of *RL*, *R* and *L* for this reaction could be written as

Here, *k*_{on} and *k*_{off} are the association and dissociation constants, respectively. The time course of *RL* would look similar to the red curve in Figure 1a, whereas time courses of *R* and *L* would be similar to the blue curve. These equations can be solved analytically, but the additional complexity of most biologically relevant models generates equations that require numerical solution by computer.

The chemical master equation (CME; McQuarrie, 1967; Gillespie, 1992; Ge and Qian, 2013) considers the discrete and stochastic nature of the biochemical system, which can cause differences from deterministic rate equations (Samoilov and Arkin, 2006). The CME describes how the probability of the system being in a specific state evolves over time by using reaction probabilities (likelihood of occurrence per unit time) rather than the equivalent reaction rates. Just like reaction rate equations, the CME assumes well-stirred (homogeneous) systems. In most practical modeling applications, the CME cannot be solved analytically (i.e., with a closed-form expression), but simulations of the CME are conceptually straightforward and widely used (see Box 2).

Both the reaction rate approach and the CME approach simply cannot capture effects of inhomogeneous distributions of molecules in space, such as receptors clustered on membranes, or intrinsic time delays due to diffusion to localized targets such as membrane-bound receptors. For more realistic simulations of cellular process, we must turn to different computational approaches that explicitly include space.

### Spatial modeling approaches

The molecular components of living systems are not distributed homogeneously and the high spatial resolution of today’s fluorescence microscopy is continuously giving us more examples of biological phenomena where the spatial arrangement of the underlying biochemical processes is fundamentally important. To model such phenomena, we have to switch from nonspatial to spatial simulations. However, this switch is frequently not easy due to a growth in the number of model and system features that must be specified (Figure 2a). The most important difference between nonspatial and spatial simulations is that the latter take into account the translocation of the interacting molecules. In the simplest case, this means that, in addition to reactions, the diffusion of molecular species in space must be simulated. Along with diffusion, the system geometry must be specified and what happens at the “walls.” This is significantly more challenging to implement when the spatial system containing the interacting molecules is not simply a square box with rigid walls. Realistic spatially resolved models often aim to capture aspects such as particular cell morphologies; examples include synaptic structures with narrow regions connected to larger cell bodies (Rangamani *et al.*, 2016; Cugno *et al.*, 2019) or geometries that are flat and almost two-dimensional, like lamellipodial extensions in migrating eukaryotic cells (Nickaeen *et al.*, 2017) (Figure 2b). Geometries for cell simulations can be designed by hand, derived from microscope images (Schaff *et al.*, 2000), or generated from machine-learned cell models (Majarian *et al.*, 2019). Currently, there are only a few simulation tools that can model cellular biochemistry within dynamic morphologies (Angermann *et al.*, 2012; Tanaka *et al.*, 2015) and the computational treatment of reaction–diffusion processes within domains that exhibit moving boundaries is still a very active field of research (Wolgemuth and Zajac, 2010; Novak and Slepchenko, 2014). We note that models allowing moving boundaries (which usually represent the membrane) do not necessarily capture the biophysics of membrane dynamics (Figure 2d). They can be decoupled methodologically, although fundamentally, they are not.

Spatial models have the capacity to build in many additional features such as mechanics, electrostatics, or coarse molecular structure that may be particularly important for useful simulations of cell biological processes (Figure 2d). However, modeling efforts that take these additional features into account to render their simulations of cellular behavior more realistic are exceedingly rare. One reason for this is that they introduce further challenges to both the numerical implementation and the mathematical descriptions of the physical model. Another reason for shying away from this level of detail is that, frequently, it would be difficult or (currently) impossible to measure sufficient parameters to be able to estimate the remaining, unknown, parameters through computational “fitting.” Limiting themselves to smaller regions of cellular membranes, several studies have integrated (simple) biochemical processes with membrane biophysics to explore how actin dynamics drive membrane protrusions (Mogilner and Rubinstein, 2005; Atilgan *et al.*, 2006). There is reason for hope that such efforts may, at some point, be scaled up to comprise larger membrane domains or even entire cells, given the ever increasing resolution and quantitative information contents in imaging-based measurements of the underlying molecular diffusion processes (Saha *et al.*, 2016; Swaminathan *et al.*, 2017).

*Reaction–diffusion equations* represent the most straightforward extension of reaction rate equations for the inclusion of spatial aspects. Instead of depending just on time as a variable, the behavior of molecular species now additionally depends on spatial coordinates. Equations describing these reactions must therefore be formulated as partial differential equations (PDEs) rather than as ODEs (Figure 1c). As is frequently the case for mathematical models of biological phenomena, only very simple situations can be described through equations that can be explicitly solved in such a way that the solution describes the behavior of the modeled system as a continuous function of space and time (Lipkow and Odde, 2008). In most cases, one has to explore reaction–diffusion equations through numerical simulations that divide space into subvolumes (frequently called ‘voxels’) and calculate how diffusion leads to exchange among the subvolumes. Reaction–diffusion equations have been widely used to model spatially resolved biomolecular dynamics and interactions of cell biological systems (Loew and Schaff, 2001). The spatial dynamics can be extended beyond pure diffusion (e.g., to include advection), and reactions can be defined phenomenologically (Hill-type or Michaelis–Menten). Like reaction rate equations, deterministic PDE reaction–diffusion equations do not capture stochastic fluctuations in species numbers and thus cannot, for example, capture pattern formation driven by a system’s sensitivity to low copy numbers (Howard and Rutenberg, 2003).

One way that stochastic reaction–diffusion equations can be formulated is the spatial extension of the CME, known as the reaction–diffusion master equation (RDME) (Figure 1d). In the RDME, instead of only defining how a system switches from one set of numbers of molecules in particular states to another (for instance when a complex in the system decays into two molecules), the RDME includes “hops” from one location to another. Importantly, just like the CME, the RDME describes discrete changes. That means, it requires a spatial discretization into subvolumes within each of which well-mixed conditions are assumed to prevail. Diffusion events of molecules are tracked only when they occur between adjacent subvolumes, not within an individual subvolume (Fange *et al.*, 2010). Similar to the nonspatial case, it is usually not possible to solve the RDME analytically, and instead it is standard practice to compute solutions by simulating a particular stochastic time evolution of the system with, for instance, the SSA (see Box 2) that adds diffusional hops to the list of events that can occur, as is done with Lattice Microbes (Roberts *et al.*, 2013), StochSS (Drawert *et al.*, 2016), and STEPS (Wils and De Schutter, 2009; Chen and De Schutter, 2017). However, care must be taken to choose the right degree of spatial resolution that strikes the appropriate balance between capturing spatial details and avoiding subvolume sizes that are so small that discretization dilutes the molecules to a point where they essentially do not “see” their potential reaction partners anymore because the molecules are spread out over distinct subvolumes (see Box 3). A recent review (Smith and Grima, 2019) discusses the relationship between the RDME model and the particle-based models of reaction–diffusion.

#### BOX 2: the Chemical Master Equation (CME)

The CME is a set of coupled linear ODEs that describe the time dependence of the probability to occupy each of a set of discrete states, with each state defined by the copy numbers of molecular species. The CME may be written as

where the composition vector *N* is composed of copy numbers of each molecular species, and thus one has a set of equations, one for each possible instantiation of *N*. Here, the sum runs over all possible reactions R. The *υ _{r}* is the stoichiometric vector of reaction

*r*that describes how this reaction changes the number of molecules in the composition vector

*N*and

*α*(

_{r}*N*) is the probability per unit time that reaction

*r*occurs, given that the system is in the state described by

*N*. The CME is important from a conceptual point of view as it represents a framework to describe probabilistic transitions and thus captures the stochasticity underlying all molecular interactions (Grima and Schnell, 2008; Schnoerr

*et al.*, 2017). The computational cost of solving the CME equations scales exponentially with the number of chemical species, and, although clever approaches have extended the size of systems for which the CME can be solved (Munsky and Khammash, 2006), high computational cost still limits biological applications. An intuitively simple way to calculate a solution of the CME would be to set up a simulation where time ticks forward in small, discrete intervals (time steps). However, the fixed time step in this integration scheme has finite error that is only eliminated in the limit

*Δt*→ 0, since reactions may occur even during shorter time steps than the one chosen to propagate the system in time.

One popular and precise method used to generate trajectories through the state space sampled by the CME without the need to choose a discrete time step is the stochastic simulation algorithm (SSA) (Gillespie, 1976; see Gillespie *et al.*, 2013, for a detailed review). In an SSA simulation, the time interval until the next reaction occurs is itself sampled, as is the type of reaction that will occur (Gillespie, 1976). Molecular species that can react quickly and have many possible interaction partners will be selected frequently, while rarer molecules associated with slower reactions will be selected rarely. Since the simulations proceed with one reaction at a time, the computational cost depends strongly on the number of particles and reaction rates in the system. In contrast, the effort required to integrate (or solve) reaction rate equations depends mostly on how many molecule types are involved and whether their interactions occur on different or similar timescales. Various approaches have been developed to increase the efficiency of both exact (Gupta and Mendes, 2018) and approximate (Schnoerr *et al.*, 2017) stochastic simulations of the CME. In addition, efficient methods have been developed to compute distributions and moments directly from the CME itself (Hasenauer *et al.*, 2014; Hellander *et al.*, 2017). A stochastic simulation of the example system from Box 1 would look similar to Figure 1b. Note, however, that such a trajectory represents only one possible time course compatible with the underlying CME. This means that many stochastic simulation trajectories must be collected to determine probability distributions and moments of the CME.

#### BOX 3: spatial discretization of reaction–diffusion equations

To numerically solve reaction–diffusion processes modeled as PDEs, it can be challenging to choose the appropriate spatial discretization of the modeled biological geometry. A discretization that is too coarse will suppress many spatial details and will represent a poor approximation of the underlying biology. However, keeping track of the contents of many very small voxels will not only be very expensive computationally but may also lead to situations where the assumptions of mass-action kinetics no longer strictly hold since the fraction of the molecules in the system that populate a single voxel becomes so small that the very concept of an average concentration becomes problematic. Furthermore, the accuracy and efficiency of PDE solvers is not just sensitive to the resolution of the spatial discretization (sometimes called lattice or mesh) but also to the discretization scheme as manifested, for instance, in the shape of the voxels. The practice of designing adaptive meshes, that is, combining voxels of different shape and size in one simulation to capture small-scale spatial details where needed while keeping the total number of voxels as low as possible, is a field of active research. The structure of the mesh also has to be adjusted to the numerical method chosen to solve the reaction–diffusion PDEs. Finite volume methods directly simulate diffusional exchange between voxels. In contrast, finite element algorithms optimize the coefficients of interpolation functions at the nodes of the mesh to achieve good approximations of concentration profile resulting from the combination of reactions and diffusion. See, for example, Richmond *et al.* (2005). We note that mesh-free approaches to solving PDEs provide an alternate to spatial discretization methods.

Similar to PDEs, the accuracy and cost of the RDME is sensitive to the spatial mesh; this problem is inherent to all spatially discretized simulations. Computational costs grow rapidly as the mesh resolution increases. Importantly, the accuracy of an RDME model does not always increase with a finer mesh. A very small mesh size violates the assumption that species are dilute and their own molecular volume is small relative to the voxel (Erban and Chapman, 2009; Isaacson, 2009; Wolf *et al.*, 2010; Isaacson and Zhang, 2018). For specific nonfundamental reaction types, RDME has an additional limitation in that it does not always converge to the CME solutions in the limit of fast diffusion, as expected (Smith and Grima, 2016). Hence, the RDME may be viewed as a *nonconvergent* approximation of more microscopic spatially continuous models, such as the Smoluchowski model discussed below. We note that a variety of lattice methods have been designed to overcome the small voxel size issue (Chew *et al.*, 2018) and to address convergence to a more microscopic model (Isaacson, 2013; Isaacson and Zhang, 2018).

*Particle-based spatial simulation methods* take into account the stochastic motion and interactions of individual molecules in continuous time and space and are thus capable of modeling biochemical processes that involve low copy numbers and strongly heterogeneous molecular spatial distributions (Figure 1e). These methods have the highest resolution (Figure 2c), but they come with a high computational cost. Importantly, simulations that treat each particle as an individual also offer the possibility of building in more detailed molecular features (Figure 2d). Typically, particle-based approaches resolve a bimolecular reaction A+B→C of a pair of molecules as two physically distinct stochastic processes. First, the molecules’ diffusive (Brownian) motion leads to their encounter. Then, the molecules either form a bond with a reaction probability determined by the reaction rate constant and their current separation or else they diffuse away from each other. Because the macroscopic kinetics of association must depend on both diffusion and reaction rate constants, this generally results in a distinction between a microscopic and macroscopic rate (see Box 4). Numerical approaches to single-particle reaction–diffusion calculations for biological systems can perhaps be best categorized into classes based on whether they model reactions to occur on collisions (von Smoluchowski, 1917; Collins and Kimball, 1949), or whether they model reactions to occur within a reactive volume (Doi, 1976; Erban and Chapman, 2009). Within both of these classes, different algorithms introduce approximations that affect accuracy in recovering the underlying physical model, flexibility, and connection to experimental (macroscopic) rates (Figure 2c).

#### BOX 4: microscopic versus macroscopic rates, and the sensitivity of strong binding to diffusion

For all nonspatial models, as well as PDE and RDME models, bimolecular association reactions are parameterized by the macroscopic rate constants, k_{on}, corresponding to the rates one would measure from a binding experiment in bulk solution. This is because in all these spatial models, species that are localized in a small volume are assumed well mixed, thus obeying the same mass-action kinetics used in nonspatial models. In the single-particle methods, however, molecular interaction kinetics is split into two steps, as described in the text. This results in a purely diffusive contribution to the bimolecular encounter and then what is effectively an energetic contribution defined by a microscopic on-rate. In the Smoluchowski model, the encounter occurs on collision at a specific binding radius *σ* with microscopic rate k_{a}. In 3D, this results in the long-known relationship:

where *D* is the sum of both species’ diffusion coefficients. With this relationship we can directly assess the impact of diffusion on controlling macroscopic kinetics. As shown in the image on the left below, for large macroscopic rates, the macroscopic kinetics of the A+B→∅ reaction is noticeably dependent on diffusion. For smaller rates, as shown in the Figure on the right, the effect of diffusion is negligible, despite D_{A} = D_{B} dropping from 100 to 1 µm^{2}/s. Here, A_{0} = B_{0} = 1000 particles (here corresponding to 62 μM). Hence, large k_{a}, or strong binding, is diffusion-limited, and small k_{a} is rate-limited.

In two and one dimensions (e.g., on surfaces and filaments), the relation between microscopic and macroscopic parameters is more complicated, due to the properties of diffusion in lower dimensions smaller than three. No single relationship exists between microscopic and macroscopic rates (see, e.g., Yogurtcu and Johnson, 2015a), but meaningful theoretical relationships can be defined if the system size is considered (Szabo *et al.*, 1980), where in 2D we further correct for system density using (Yogurtcu and Johnson, 2015a):

*N _{A}* and

*N*are the copy numbers of reactants A,B and the system size is S.

_{B}For our test cases below, we thus always derive the microscopic rates k_{a} to reproduce the macroscopic rates using Eq. A or Eq. B. This is because it is already clear from nonspatial models that changes to the macroscopic rates will necessarily alter the reaction kinetics. Since we are not focused on probing the influence of kinetic parameters on molecular behavior, but rather the role of explicit spatial representations in controlling species distributions and encounter times, we preserve all k_{on} values. However, it is worth noting that for a reaction pair with a large k_{a} value, if diffusion slows throughout the simulation due to, for instance, formation of large complexes, then the macroscopic kinetics will also slow down, an effect that is naturally captured in GF-based methods.

The definition of the reaction probability is the primary challenge and distinguishing feature of different single-particle algorithms. For the first class of collision-based numerical methods, the reaction probabilities are derived or matched to the Smoluchowski model of diffusion-influenced reactions, which naturally captures excluded volume (Figure 2c). We note that excluded volume is a critical feature for dense systems and for single-particle simulations of clustering or assembly where molecular structure/volume impacts interactions between species. For either class of models, the Green’s function (GF) approaches provide the most accurate solution by predicting the encounter probability for pairs of particles based on the particles’ initial positions. GF methods allow for much larger time steps, although they can only calculate the encounter probability for two particles at a time, meaning that the simulated system has to be segmented into two-particle subsystems (see Box 5). In practice, this turns out to be feasible for many interesting biological problems. Frameworks that have been developed to take advantage of this GF approach include FPR (Johnson and Hummer, 2014), NERDSS (Varga *et al.*, 2020), GFRD (van Zon and ten Wolde, 2005), and eGFRD (Sokolowski *et al.*, 2019). The SpatioCyte methods are performed on a lattice and recover the correct kinetics (beyond short times) and equilibrium of the Smoluchowski model (Chew *et al.*, 2018, 2019). Smoldyn is also derived to use large time steps (albeit without excluded volume), and it is simpler to implement than GF approaches (Andrews and Bray, 2004a; Andrews *et al.*, 2010; Andrews, 2017). However, the reaction parameters are coupled to the time step size, rather than representing independent model features (e.g., binding radii and microscopic rates), resulting in a time dependence (and 2D equilibrium) that is not as rigorously correct.

#### BOX 5. single-particle reaction–diffusion simulations

Single-particle methods simulate the stochastic behavior of individual and (in spite of the name) pairs of interacting molecules. Any biochemical network whose description does not include ad-hoc phenomenological processes (such as, for example, Hill coefficients describing nonlinear dose–response characteristics) can be described as composed of uni- and bimolecular reactions. As unimolecular reactions are only time-dependent, they are typically modeled as Poisson processes. For bimolecular reactions, the distance between a pair of particles influences the probability that they will diffuse to either collision or their reactive volume and react with one another in a time step. The time evolution of the molecules’ positions is described by a stochastic differential equation, the overdamped Langevin equation (Van Kampen, 2007). Its numerical implementation, known as BD (Ermak and Mccammon, 1978; Northrup *et al.*, 1984), requires tiny time steps to accurately resolve molecular encounters, which renders the BD scheme highly inefficient. It is, however, possible to derive reaction probabilities for pairs of molecules that are nonetheless propagated using BD updates but use larger steps. GFRD is the only method that does not use BD updates for reactive pairs. To calculate the distance-dependent reaction probabilities for pairs of molecules, the most accurate approach is to use the GF defined below. This can enhance the efficiency of BD simulations by resolving bimolecular reactions within one large time step, Δ*t*, without approximation.

The GF *p*(*r,* Δ*t*|*r*_{0}) can be obtained as the solution of the diffusion equation that describes a pair of molecules *A, B* that diffuse with diffusion constants *D _{A},D_{B}*, respectively, and may undergo a reaction

*A*+

*B*→ as follows

where *D* = *D _{A}* +

*D*and

_{B}*r*,

*r*

_{0}refer to the distance between

*A*and

*B*after and before the time step, respectively. In accordance with the two-step picture described in the main text, reactions are incorporated by imposing boundary conditions that specify the physics at or within an encounter distance

*r*= σ. In the collision-based Smoluchowski model,

*r*is always ≥ σ, and the Collins-Kimball Boundary Conditions (Collins and Kimball, 1949) in 3D is written as:

where *k _{a}* refers to the intrinsic reaction constant. In the volume reactive, or Doi model, reactions occur whenever

*r*≤ σ, with intrinsic rate

*λ*, creating a reactive sink between the two particles (Doi, 1976). It is worth noting that finding the appropriate GF to capture the desired properties of the molecular interactions can be a challenge. For the case of

*reversible*diffusion-influenced bimolecular reactions of an isolated molecule pair in 2D, the GF was derived only in 2012 (Prüstel and Meier-Schellersheim, 2012).

For the second class of volume-reactive numerical methods, the reaction probabilities are derived based on a distance cutoff between particles (often called Doi model), which thus naturally lacks excluded volume (Figure 2c). There are no GF-based algorithms for this model. Erban and Chapman derived reaction probabilities and corresponding microscopic rates for this model in the limit of small time steps (Erban and Chapman, 2009), which is the basis for the implementations ReaDDy (Schoneberg and Noe, 2013), and SpringSaLaD (Michalski and Loew, 2016). Both these implementations introduce methods to capture excluded volume via, for example, repulsive short-range forces, which ultimately require adjustments to properly recover reversible reactions, as done in ReaDDY 2 (Hoffmann *et al.*, 2019). With these methods, the ability to reproduce the underlying model will depend on using small time steps and assessment of the extent to which the modeled forces introduce error into the kinetics of many-body systems. Last, the widely used implementation MCell (Kerr *et al.*, 2008) is not based on either the Smoluchowski or the Doi model, and instead derives reaction probabilities that quantify collisions within a volume, where this instantaneous volume depends on time step and diffusion constants. MCell lacks excluded volume, but it can take large time steps and recover the proper equilibrium in reversible reactions.

The future extensibility of all methods hinges on the feasibility of finding mathematical expressions for the crucial reaction probabilities that incorporate additional features and details, such as curved surfaces, intramolecular constraints, and external and internal deterministic forces (Figure 2, b and d). Reaction probabilities for multisite molecules have been derived, for example, by assuming rigid molecules and simplifying the dependence on orientation (Johnson, 2018). The addition of any interaction potentials (and therefore forces) between particles can significantly alter reaction kinetics (Zhou, 1990), and quantifying reaction probabilities has required either substantial computational overhead (Johnson and Hummer, 2014) or steady-state assumptions in dilute systems (Dibak *et al.*, 2019). Careful validation of these additional features is critical for producing models that can be quantitatively reproduced across multiple simulation platforms.

## RESULTS

We present a series of test problems relevant to spatial modeling of cellular and subcellular processes. This list is not meant to be exhaustive, but rather to permit a manageably sized survey of the kinds of problems that different simulation programs may be challenged to solve within a larger biological study. Our selection thus includes very simple problems that can be solved by many simulation tools as well as complex problems that can be solved (at present) by only a few. By framing the problems explicitly, we facilitate a direct comparison among different simulation packages both with respect to accuracy of execution and how they encode these particular scenarios. While all the models presented have been simulated previously, they have not been subject to the quantitative comparative analysis performed here across multiple model and method types. This comparison provides us with distinctive insight into the sensitivity of quantitative and qualitative behavior that emerges with specific biologically relevant features, which we summarize at the end of *Results*.

For these test cases, we contrast results from stochastic, deterministic, spatial, and nonspatial modeling approaches. We use Virtual Cell software for all nonspatial simulations and spatial deterministic simulations. For spatial stochastic simulations, we use the single-particle softwares NERDSS, Smoldyn, MCell, and eGFRD. We provide executable model inputs and numerical outputs for each model in our publicly accessible repository (https://github.com/spatialnimbios/testcases/__).__ We summarize a broader range of actively developed tools and their features in Supplemental Table S2, as distinct software tools have introduced selected complex features of RD systems (also discussed previously; Takahashi *et al.*, 2005; Schoneberg *et al.*, 2014), and this is an additional consideration for users when selecting a tool for their biological problem.

### Category 1: “unit test” cases

This category of problems represents fundamental building blocks for which there is a known correct answer (at least at steady state). We emphasize that because these reactions form the basis of much more complex models and geometries, it is essential that they be carefully tested for accuracy with regard to reaction kinetics and, in the case of bimolecular reactions, reversibility. For all of these, we initialize simulations with well-mixed components, and thus one may expect that any modeling approach would give the same outcome. However, due to differences in both approaches and algorithmic choices, we find that differences do in fact emerge in specific parameter regimes, particularly at short times before the system approaches steady state (Figure 3).

#### 1A: bimolecular association in 3D, 2D, and from 3D to 2D.

Although seemingly simple, bimolecular association events require both a diffusional encounter and a reactive event; thus, the rate-constants and the kinetics are dependent on the dimensionality of the systems, and even for well-mixed systems, spatial details can cause deviations from nonspatial models. For reversible bimolecular association of well-mixed reactants in a closed system, the equilibrium state is analytically solveable and the kinetics of nonspatial rate equations also have analytical solutions. The effects of diffusion (Agmon and Szabo, 1990; Zhou and Szabo, 1996; Gopich *et al.*, 2001; Gopich and Szabo, 2002a, b), electrostatics (Zhou, 1993; Schreiber *et al.*, 2009), orientation (Shoup *et al.*, 1981; Zhou, 1990), and dimensionality (Szabo *et al.*, 1980; Torney and McConnell, 1983; Prüstel and Meier-Schellersheim, 2012; Yogurtcu and Johnson, 2015) on reaction kinetics have also received considerable theoretical study, providing a rich basis for understanding spatial effects. One may note that when bimolecular association is reversible, recovering the proper equilibrium is a simple test that can nonetheless be challenging for single-particle methods. Reaction probabilities and the placement of reactants on un/binding events must be derived to ensure equilibrium is reached (Box 5).

In 3D, all models and tools produce nearly identical results, even for this strongly diffusion-influenced reaction (k_{on} = 1.48 10^{7} M^{–1}s^{–1}). This is as expected (Figure 3, top panel). Although at short times the kinetics is slightly faster for Smoluchowski-type simulations (NERDSS), the kinetics rather rapidly converges to the macroscopic rate equations (Johnson and Hummer, 2014). Differences between single-particle and nonspatial methods can also emerge for reversible reactions as they approach equilibrium (Mattis and Glasser, 1998; Tauber *et al.*, 2005), but these can only be effectively observed with high numerical precision and statistics—usually they are dwarfed by the copy number fluctuations.

Unlike in 3D, macroscopic rate equations in 2D only approximate the dynamics captured in Smoluchowski-type approaches at all times (Fange *et al.*, 2010; Hellander *et al.*, 2012; Yogurtcu and Johnson, 2015) (see Box 4). All macroscopic rate-based methods produce the same kinetics as each other (Figure 3, middle panel). Here we see distinctions between the spatial PDE and the spatial single-particle methods. Although species diffuse in the PDE, because they are present at all positions in space (due to uniform initial conditions), association is not dependent on their spatial distribution. For single-particle methods, there is always a distribution of starting separations between species that leads to some very fast reactions initially and at long times produces slower reactive collisions as particles that started off close to each other have already been consumed in the reaction. The Smoldyn method uses the steady-state solution to the Smoluchowski model to derive reaction parameters (Andrews and Bray, 2004b), but in 2D there is no steady state, and thus the reaction parameters are approximate. Because of this, Smoldyn can generate inaccurate kinetics in certain parameter regimes, with deviations being typically small in 3D but significant in 2D.

For binding between 3D particles and 2D particles (relevant for biological cases where soluble cytoplasmic proteins bind to membrane proteins or lipids), all models produce the same equilibrium, but the spatial models have slower kinetics delayed by diffusion to the surface (Figure 3, bottom panel). The extent of divergence between the nonspatial and the spatial models is driven by three factors, the “height” the solution volume stretches from the membrane surface, the diffusion coefficient of the 3D particles, and the speed of the binding. Here we simulated a fast, strongly diffusion-influenced reaction (8.4 × 10^{7} M^{–1}s^{–1}), meaning nearly every collision results in a reaction. For a simulation box with *h* = 0.2 μm, one can estimate an average time to diffuse to a surface particle would be ∼60 μs (D^{3D} = 30 μm^{2}/s), which is relatively fast. However, with the numerous 2D particles mixed in the solution volume for a nonspatial simulation, the time to diffuse to a “surface” particle drops to ∼6 μs. Thus we find a mean relaxation time of ∼200 μs without space versus ∼700 μs with space (Figure 3, bottom panel). By dropping the reaction rate to more moderate protein–protein interaction levels, the spatial and nonspatial results begin to converge. Smoldyn shows excellent agreement for larger steps, here 10^{–6} s, although the kinetics shift slower for shorter steps. We note that when particles can only collide with one another from one side (because one is embedded in a surface, for example), this reduces the binding by a half, and solvers should explicitly account for this so that user-defined rates produce the equilibrium expected from a nonspatial model.

Last, transitions to the surface can be modeled using adsorption, which uses an effectively 1D rate. This is more efficient but importantly, it need not reach the same equilibrium as explicit particle simulations, because the occupancy of surface binding sites is not accounted for. Modeling explicit particles thus gives more control over the surface properties, and algorithms for binding to surfaces while accounting for site occupancy using implicit sites rival adsorption models in speed (Fu *et al.*, 2019). Not all tools allow for all types of surface binding; hence, it is important to recognize these distinctions between adsorption versus single-site binding.

#### 1B: crowding.

Inside of living cells, the extremely high density of macromolecules (with typical spacing on the order of nanometers) can alter the speed of molecular diffusion (Ando and Skolnick, 2010) and kinetics of intermolecular reactions, either increasing or decreasing biochemical reaction rates as compared with rates in dilute solution, depending on the size and mobility of the crowders (Zimmerman and Minton, 1993; Schreiber *et al.*, 2009). As crowders become larger and less mobile (e.g., vesicles), they act more like barriers to encounters, slowing down rates of association (Minton, 2006; Zhou *et al.*, 2008). However, as we see here in our test-case, when crowders are of comparable sizes and similar mobility to the reactant species, they drive up rates of association. This rate increase is due to a reduction in the total volume available to the reactants, effectively concentrating them without providing substantial barriers to encounters (Minton, 2006; Zhou *et al.*, 2008). To quantify how increasing concentrations of crowding agents alter bimolecular association rates, here we simulated the bimolecular reaction A+B→B+C in the presence of additional, inert, C crowders, where all species are mobile, the same size, and exclude volume (Figure 4a). The analytical solution for no crowding/no excluded volume, *A* (*t*) = *A*_{0} exp (–*k _{macro} B_{tot}t*), provides a convenient baseline and fit function for interpreting deviations due to crowding/excluded volume (Figure 4b).

For the single-particle algorithms that capture excluded volume (GFRD, van Zon and ten Wolde, 2005; and FPR, Johnson and Hummer, 2014), two main results emerge. First, the overall kinetics of association increase with increasing crowding fraction, up from *k*_{on} = 63.5 nm^{3}/µs at zero crowding/no excluded volume (3.8 × 10^{7} M^{–1}s^{–1}) to ∼85–100 nm^{3}/µs with 25% crowding fraction. This result is consistent across both algorithms, indicating that for GF-based solutions to the Smoluchowski model applied to a strongly diffusion-influenced reaction, small mobile crowders will enhance reaction rates for concentrated reactants (here A_{0} = 13.3 mM). This same trend was observed for simulations at lower reactant concentrations but comparable rate constants (Kim and Yethiraj, 2009). In contrast, simulations that immobilized the crowders caused them to act as a rigid barrier, leading to a reduction in reaction rates despite using the same reactant concentrations studied here (Andrews, 2020). Second, for high crowding regimes, we find that the kinetics is not described by a single rate constant (Figure 4d), whereas at low crowding, the results fit extremely well to the nonspatial analytical solution, with a new rate-constant (Figure 4c). This is perhaps not surprising; for diffusion-influenced reactions, short-timescale kinetics is dominated by reactants that are already close together where we expect crowding agents to promote their repeated collisions. Then, as the reactant populations decrease either to zero or toward an equilibrium state, the kinetics slows relative to nonspatial predictions, although this shift may be hard to detect (Mattis and Glasser, 1998; Tauber *et al.*, 2005; Johnson and Hummer, 2014; Yogurtcu and Johnson, 2015a). Our results suggest that crowding agents exacerbate this slow search for the final reactants, similar to what happens in 2D, making the deviations from a single-rate constant kinetics easier to detect.

Although we show here that mobile crowders impact observed reaction rates, the changes are often quite modest. For the strongly diffusion-influenced reaction simulated here, we measure clear increases in rates, but for more rate-limited reactions (*k*_{on} = 6 × 10^{5} M^{–1}s^{–1}), the effect of crowders on the reaction rate is minimal (data not shown). Finally, FPR and GFRD are not in perfect agreement in terms of the quantitative size of the change in kinetics, although qualitatively they both predict a higher rate. At these high densities, GFRD converts to a brute-force BD algorithm, rather than its exact event-driven method, because the overhead cost of propagating GFRD is more expensive in dense systems. Because FPR combines the GF approach with simple Brownian updates and reweighting corrections, it performs well in dense systems, but still must also take extremely small time steps (10^{–10} to 10^{–12} s) to prevent particle overlap. It is not clear which method is more correct, as neither will preserve exactly two-body problems at each step. Further, we note that over the span of a picosecond time step (10^{–12} s), particle dynamics is not truly diffusive but is inertial. Capturing these dynamics would require a different model (e.g., generalized Langevin Dynamics; Van Kampen, 2007) that tracked both positions and velocities of particles and would be a valuable comparison in future work.

### Category 2: “intermediate” cases

This category includes slightly more realistic biological cases with interesting emergent properties. These are particularly useful for illustrating the fundamental conceptual differences among the different modeling approaches.

#### 2A: exploiting membrane localization to stabilize protein–protein interactions.

In a variety of biological processes, including clathrin-mediated endocytosis and initiation of signal transduction, multivalent proteins localize to membranes and assemble into larger, multiprotein complexes. Reducing dimensionality from a 3D search to a 2D search can accelerate kinetics of receptor binding (Adam and Delbruck, 1968; Berg and Purcell, 1977; Axelrod and Wang, 1994) and stabilize interactions for macromolecules restricted to the surface (Minton, 1995; Kholodenko *et al.*, 2000; Abel *et al.*, 2012). For soluble cytoplasmic protein binding partners, however, localization to the 2D membrane can further drive dramatic increases in protein complex stability, as was recently quantified using a simple model of two proteins that bind to one another and to a specific membrane lipid (Yogurtcu and Johnson, 2018) (Figure 5, a and b). The origin of the increased stability is largely a concentration effect, where the proteins collide with one another more frequently on the surface than in solution. While a change also must occur in 3D versus 2D equilibrium constants (Wu *et al.*, 2011), the magnitude is usually much smaller than the change in 3D versus 2D search space (V/A > K^{3D}/K^{2D}) (Yogurtcu and Johnson, 2018).

With the ability to localize to the membrane and reduce their search space, over 90% of proteins end up bound to one another (Figure 5c), which contrasts markedly with the purely 3D solution binding result calculated at equilibrium where only 4.6% would end up bound (for K_{D} = 20 μM and [A]_{0} = [B]_{0} = 1 μM). Although all simulation methods are expected to produce the same equilibrium, and in this example stochastic effects are minimal, we find here that spatial effects arise because diffusion slows the localization of A and B molecules to the membrane (Figure 5d). The height of the box in the spatial simulations is 5 μm, and so despite efficient solution diffusion (50 μm^{2}/s), the proteins are recruited more slowly to the membrane in spatial simulations than in the well-mixed simulations (compare red and green curves in Figure 5d with purple and blue curves). The protein–lipid complexes plotted in Figure 5d peak before dropping, as they bind to one another to equilibrate, and the peak is lower when they do not localize rapidly to the membrane. Errors in the results calculated using Smoldyn arise due to approximate treatment of purely 2D interactions, which lead to quantitative deviations from expected equilibria and kinetics for these reactions. This current limitation is being actively addressed in Smoldyn software (S. Andrews, personal communication).

Overall, these types of reactions form a critical component of more complex models of membrane-mediated assembly. Due to the quantitative differences we observe here in time dependence of species numbers, if the model is coupled to reactions that drive it out of an equilibrium steady state, this could then drive qualitative changes in the biological outcomes (as in case 2B below).

#### 2B: increasing stochastic fluctuations in a system with multiple steady states.

Positive feedback in combination with other interactions can give rise to systems with multiple, distinct steady states. Such systems can exhibit large differences in their dynamics depending on whether they are simulated deterministically or stochastically. A simple model that illustrates such effects is the autophosphorylating kinase model first introduced by Lisman (Lisman, 1985) and studied more recently by Agarwal *et al.* (Agarwal *et al.*, 2012), who analyzed a stochastic version of the model. Figure 6a shows a diagram of the model, which consists of a kinase that can activate itself through phosphorylation (reactions 2 and 3) and a phosphatase that can bind and dephosphorylate the active kinase (reactions 4 and 5). Because the production rate of active kinase (Ap) increases with the amount of Ap starting at low concentrations (Figure 6b, red curve), the system exhibits positive feedback. The rate balance plot shown in Figure 6b illustrates a case where the production and degradation rates of Ap as a function of Ap intersect multiple times to give rise to multiple steady states. This model has three steady states, also called “fixed points,” two of which are stable (Figure 6b, filled circles). For a stable fixed point, the balance of production and decay returns the system back to the fixed point following any slight changes in Ap away from the steady state value, whereas at an unstable fixed point, the balance of rates carries the system away from the fixed point as a result of any tiny fluctuation.

Deterministic simulations of the system, whether spatial or not, converge rapidly to one of the stable steady states depending on the initial level of Ap (Figure 6c, black dashed lines). Simulations starting from a state with low initial Ap will reach the steady state with lower Ap and vice versa for high Ap. Stochastic trajectories, on the other hand, may initially stay in the vicinity of the closer steady state, but fluctuations due to noise occasionally induce switches between states. In both trajectories shown in Figure 6c, the system starts in the lower steady state (as indicated by gray shading) but after a minute or so switches to the higher steady state where it continues to display fluctuations, some of which lead to short-lived excursions back to the lower state. For the set of rate parameters we used here, the system spends more total time in the higher state than the lower state (88% vs. 12%). In addition, the length of each segment in the higher state, called the residence time, is longer on average (2.6 s vs. 0.34 s) (see *Materials and Methods* and Supplemental Figure S1).

Despite the dramatic differences between the deterministic and the stochastic simulations for this model, we find the addition of space has relatively modest effects, with spatial models still producing bistable switching as diffusion slows. Figure 6c shows a nonspatial SSA simulation in the top panel (SSA) and an explicitly spatial simulation in the lower panel (NERDSS). For the geometries and diffusion constant values chosen here, these two simulations yield small differences between the probabilities and the residence times for each state (Supplemental Table S4). When the diffusion constant is reduced by a factor of 10, the probabilities of the lower state drop significantly although not dramatically, and the residence times of both states drop (Supplemental Figure S2 and Supplemental Table S4). The shorter residence times in spatial simulations relative to the nonspatial likely result from larger fluctuations in copy numbers due to transient spatial inhomogeneities, making both steady-states slightly less stable. Reduction of the diffusion coefficient by another factor of 10 makes simulation incompatible with the macroscopic reaction rates specified in the model. In other words, the reaction rate is limited by diffusion (see Box 4). We found that the effect of varying diffusion constant was quantitatively similar in the same model, but with all rates reduced by a factor of 10, reducing the fastest rate from 8 × 10^{8} to 8 × 10^{7} M^{–1}s^{–1} (Supplemental Table S5). A study on a similar model of bistability found that slow diffusion limited the parameter regimes where one observed bistable switching due to fluctuations in local concentration, but also by changing the effective rates (Endres, 2015). Here we kept the macroscopic rates the same as the nonspatial model even as diffusion constants slow, and our results therefore indicate that timescales are sensitive to spatial inhomogeneities that persist longer with slower diffusion.

### Category 3: “application” cases

Proper implementation of the test cases in this category requires the building blocks described above. These well-defined problems yield rich, biologically interpretable outputs that are sensitive to parameter choices.

#### 3A: stochastic effects in gene expression.

Proteins that are important for controlling circadian clocks exhibit regular oscillations in expression levels. These oscillations are quite robust even in the presence of stochastic copy-number fluctuations. In fact, stochastic fluctuations can support oscillations in regimes that a purely deterministic model cannot (Supplemental Figure S3). Using a simple model of circadian oscillations (Vilar *et al.*, 2002) (Figure 7a), we simulated the behavior of an activator protein A and repressor protein R that are produced from mRNA transcribed from a single copy of a gene (one for each protein). Coupling of A and R expression is driven by positive feedback of the activator A, which binds to each gene’s promoters to enhance transcription. Protein R also binds to A to effectively degrade it, and all proteins and mRNA are also degraded spontaneously at a constant rate. If the spontaneous degradation rate of protein R is slow, the oscillations will quench in the deterministic model, but persist in the stochastic solutions, as demonstrated in the original work and reproduced in Supplemental Figure S3.

We find that with the addition of space to the model, with all species diffusing at *D* = 10 μm^{2}/s, the oscillation times show no real significant differences in single-particle or deterministic solutions relative to the nonspatial model (Supplemental Table S6). To quantitatively compare the kinetics across all models, we cannot use simple steady-state values due to the oscillations. These time-dependent oscillations are nearly perfectly regular in the deterministic models, but are quite imperfect (although recurring) in all stochastic and single-particle methods (Figure 7b). We therefore measure the average time interval between peaks in the expression of A and the lag time between the appearance of a peak in A expression followed by a peak in R (see *Materials and Methods*). The similar results across all methods show how with relatively small spatial dimensions (sphere of R = 1 μm), purely 3D reactions, and all species diffusing at the same *D* = 10 μm^{2}/s, no spatial dependence was distinguishable.

The lack of any significant spatial dependence is somewhat surprising because we pushed the reaction rates into the strongly diffusion-influenced regime. The model was actually formulated to describe the slow oscillations of gene expression regulation, with all rates reported in units of hr^{–1}, rendering diffusion times irrelevant. Here, we chose to accelerate the rates by a factor of 3600 (from hr^{–1} to s^{–1}) due to the computational expense of simulating spatial models with explicit diffusion. This result overall indicates that the reactants still mixed faster than any spatial inhomogeneities could emerge, preventing deviations between spatial and nonspatial simulations.

We were finally able to measure a significant dependence on diffusion when we simulated the same set of reactions but we localized and immobilized the promoters for A and R to a small nucleus in the center of the volume, with unrestricted diffusion in and out of this nucleus (Figure 7c). For the same diffusion constants of *D* = 10 μm^{2}/s for all other species, this immobilization of promoters was not enough to have an effect, and oscillation times remained not significantly different (Supplemental Figure S4 and Supplemental Table S7). However, when the diffusion constant of the activator A was slowed to 2 μm^{2}/s, we observed a persistently higher concentration of A and its mRNA near the cell center relative to the cell periphery (Figure 7 and Supplemental Figure S4). This localization of A near the gene promoters had the effect of shortening the oscillation period in the PDE from 25 to 22.5 s (Supplemental Table S7). This change was independent of the diffusion constant of R, which does not bind the genes. Further, this same trend was observed in the single-particle simulations, which also produced faster oscillations as D_{A} slowed to 2 μm^{2}/s (Supplemental Table S7). Overall, we thus found that while slowing the search time for the Activator for both genes does modulate the oscillation times, the oscillations are quite robust and strongly dictated by total copy numbers and the unimolecular decay reactions, which are inherently independent of diffusion and spatial dimensions. This fairly complicated scenario is one where intuitively we had expected to see strong spatial effects, so it is rather interesting that the actual spatial effects are minor in this case. Certainly this outcome might vary depending on the parameter regimes selected.

#### 3B: spatial and temporal oscillations in MinCDE.

When coupled to reactions, diffusion that is sufficiently slow or occurs over large enough length scales can establish spatial gradients of concentrations and subsequent oscillations in space and time. By construction, spatial oscillations are lost in nonspatial models, and no temporal oscillations occur either (Supplemental Figure S5). In a simplified model of bacterial cell division, the MinD and MinE proteins (Huang *et al.*, 2003) spatially control the location of cell division by creating an oscillating spatial gradient of both proteins in the cytosol and on the membrane (Figure 8). The site of cell division is determined by assembly of a ring constructed from the bacterial tubulin homolog FtsZ. Because the Min proteins inhibit the assembly of the FtsZ ring, the spatial oscillation of the Min proteins from one end to the other in the cylindrical bacterial cell results in the localization of the FtsZ ring to a site very close to the geometrical center (Raskin and de Boer, 1999). This model involves 3D and 3D→2D bimolecular reactions and unimolecular reactions with no 2D reactions. Instead of forming explicit polymers, the membrane-bound proteins act to locally increase protein density through recruitment.

We are able to produce very similar spatial and temporal oscillations of the MinD-ATP protein on the membrane in both a stochastic single-particle model (Smoldyn; Andrews, 2017) and a deterministic PDE solution (using Vcell; Moraru *et al.*, 2008; Schaff *et al.*, 2016). The major distinction between the two models is that the deterministic PDE is able to support symmetric or striped oscillations in the protein, as would be generated by two traveling waves in opposite directions. However, we find that these dynamics are not stable even in the deterministic solution—the accumulation of numerical precision error eventually breaks the symmetry and transitions to pole-to-pole oscillations. We show that by tightening the error tolerance on the numerical integration (by a factor of 10) or by increasing the PDE mesh resolution (by 2), one can delay the onset of this transition, clearly demonstrating the dependence on numerical precision (Figure 8e and Supplemental Figure S6). The stochastic single-particle method is therefore not able to support the symmetric oscillations due to copy number fluctuations and always produces pole-to-pole oscillations (Figure 8b). This shows quite clearly that the pole-to-pole oscillations are more robust to stochastic fluctuations in copy numbers for this geometry, and in the wild-type biological systems, these are the types of oscillations always observed for cells of this length (Raskin and de Boer, 1999).

A challenging aspect of model comparison for any simulations tracking the spatial and time-dependent concentration of species is deciding how best to quantify the results. Analysis is often defined in a problem-specific way. For the MinCDE model, depending on the cell geometry and model parameters, it may not even reach a steady state over long timescales, but can continue to oscillate and change oscillation amplitude or frequency. For symmetric geometries (like cylinders), we reduced the analysis to a function of time by tracking the MinD-ATP^{2D} molecules at a specific point on the membrane (Supplemental Figure S7). Then, similar to the analysis for Figure 6, we calculated the average period between oscillations (*Materials and Methods*), finding values of 43.0 ± 0.7 s for Smoldyn and 41.3 ± 0.3 s for the PDE.

The oscillations in this model are fairly sensitive to initial concentrations of species (Supplemental Figure S6) and their relative stoichiometry. If the initial concentrations were cut in half, we found the oscillations disappeared in both models. The MinCDE model has been extensively studied, with the major determining features of oscillations being initial concentrations and cell geometry (see Halatek and Frey, 2012 and reference therein). The oscillations are relatively robust even to small changes in the reaction network, as long as membrane recruitment (3D→2D) and ATP hydrolysis are included (Halatek and Frey, 2012). We note that for Smoldyn, the oscillations were sensitive to the diffusion coefficients of the membrane-bound particles, set here to D = 0.05 μm^{2}/s. If proteins did not diffuse on the membrane, the oscillations disappeared in Smoldyn, whereas they were retained in the PDE.

### Summary of test case outcomes

By simulating the same biological cases using a range of spatial, nonspatial, stochastic, and deterministic methods, we have shown here how specific biological features and in some cases algorithm selection (e.g., classes of single-particle approaches) can alter quantitative and and even qualitative outcomes. As summarized in Figure 9, quantitative differences can be small, medium, or large based on system geometry and parameter regimes. Our results with fast versus slow diffusion are consistent with a recent review on spatial simulations that contain a more detailed analysis of the RDME method (Smith and Grima, 2019). We find that qualitative differences can be minor, but can also be major, producing entirely new behavior patterns that emerge. For example, major impacts could be observed in application scenario A2, which combines an elongated geometry with a comparatively slow diffusion (in comparison to the reaction rate), and includes both 3D and 2D dynamics. In general, we find that effects of introducing stochasticity to a system may be orthogonal to its spatial properties, as wrongly applying a deterministic method is likely to dominate all other choices. This can be most clearly seen for the results as shown in I2. Overall, we found that, while stochasticity played a minor role in models with reversible interactions (U1, I1), it was capable of driving major changes in models with irreversible interactions and feedback (I2, A1). Quantitatively, crowding might have a large impact, in particular if combined with slow diffusion and high reaction rates. Also, the impact on dynamics that involve reactions in 2D, be this 3D→2D or 2D dynamics, is significant and cannot be ignored in many cases. The suitability of specific simulation methods and their limitations must also be taken into account when any investigator is making choices about which methods to employ for a particular biological problem (see Supplemental Table S2). For example, if crowding and sterics need to be considered, particle-based approaches, in particular those with excluded volume, should be applied. However, in the case of very dense crowding, these approaches do not scale well because of their computational cost.

The spatial dimension in our test cases primarily drives quantitative changes in outcomes of varying magnitude. However, once multiple features are combined, as in the A2 MinCDE model, major qualitative differences emerge. We expect that our results will extrapolate to informing model selection in a broader range of biological processes. For example, the effect of clustered molecules acts to spatially localize interactions, which similar to the membrane localization (I1) can drive dramatic increases in complex stability and changes to kinetics. Compartmentalized or highly elongated and narrow geometries also can act to locally alter species concentrations, not unlike the role of dimensionality in driving quantitative changes (U1, I1) or qualitative changes (A2). Last, in all the models studied here, all species were initially well mixed in their system volumes, albeit in U1 (3D→2D) and I1 there was both a membrane and a solution volume. By introducing an explicit spatial concentration gradient via localized sources or sinks of molecules, spatial effects are inevitable and can only be captured with explicit spatial representations.

## DISCUSSION

Since the use of GFP as a genetically encoded fluorescent tag came into widespread use over a quarter-century ago (Chalfie *et al.*, 1994), methods for live-cell light microscopy have advanced spectacularly, and the spatial resolution and temporal resolution now possible for imaging experiments in living cells is truly awe-inspiring (Chen *et al.*, 2014; Guo *et al.*, 2018; Liu *et al.*, 2018). At the same time, a great deal of technical development has gone into the generation of live-cell fluorescent reporters for biological readouts including signaling events, enzymatic activity, and force generation (Greenwald *et al.*, 2018; Yasunaga *et al.*, 2019). Many of the cell biological phenomena illuminated by direct observation of dynamic molecular events in living cells reveal a spatial component that could not be fully appreciated during the era when the most sensitive assays for molecular events were biochemical in nature (Machacek *et al.*, 2009; Tay *et al.*, 2010).

At this point, the quality of time-resolved live cell imaging data for cell biological events, and its quantitative accuracy, is far ahead of our ability to model these events computationally. As a community we could avoid doing the work to develop detailed, spatially resolved models when the data quality was relatively poor, relying on images of fixed cells at poor spatial resolution, but now that the data are excellent there is no excuse. Quantitative data demand quantitative models.

Accurate, spatially resolved 3D modeling of subcellular processes is inherently challenging, particularly as compared with the larger field of time-resolved cellular modeling where space is not taken into account. In exploring the currently available tools for 3D modeling, and dreaming about what the future might hold, we have identified several sticking points that will require focused community effort to overcome.

### Spatially accurate modeling is intrinsically computationally expensive

Most interesting cell biological events unfold over minutes, hours, or days, requiring that relevant computational models be able to simulate comparable lengths of real time. For spatially resolved models, keeping track of the positions as well as the numbers and identities of relevant molecular species requires substantially greater computational overhead than simpler models that include only numbers (or concentrations) changing over time. In practice, this difficulty often limits spatially resolved models to keeping track of only a few molecular species.

An interesting analogy can be drawn to the field of MD simulations, which aim to model relative movements of atoms within proteins (Dror *et al.*, 2012), while models of the scales that we are discussing typically aim to simulate the relative movements and transformations of molecules within cells. The equations of motion for atoms within protein molecules are relatively well understood, but accurately calculating their trajectories over even the millisecond timescales relevant for protein folding or conformational changes requires massive computational power (Shaw *et al.*, 2010). Generally, increases in complexity and timescale for simulations of atoms in proteins have benefited from two types of efforts, improvements in hardware including purpose-built computers and distributed computing architectures for parallelization (Snow *et al.*, 2002), accompanied by improvements in statistical sampling approaches to capture rare but significant events (Zuckerman, 2011). Even so, most MD simulations still cannot simulate processes central to the lives of proteins such as enzymatic catalysis, because they essentially treat atoms as hard balls subject to defined forces, but ignore the quantum mechanics that would be required to describe chemical transformations. Some hybrid methods have shown promise by including quantum-mechanical detail only in the enzyme’s active site (Karelina and Kulik, 2017), but naturally this improvement in chemical detail comes at substantial computational cost.

Similarly, for spatially resolved models at the cell biological scale, there is room for substantial improvement in efficient design of hardware and computational parallelization, development of appropriate statistical sampling methods, and hybrid methods that could smoothly integrate full detail in some particular locations or at some particular times with computationally efficient approximations elsewhere (Hellander *et al.*, 2017). GPU architectures offer many advantages for efficient parallelization and have been effectively used by several groups to push cell biological models to much larger scales (Hallock *et al.*, 2014; Ghaemi *et al.*, 2020), though these approaches have not yet been widely adopted. There are also some kinds of modeling problems that intrinsically cannot be parallelized, where each step in time depends on the state of the system as a whole (Enculescu *et al.*, 2010).

Several features of cell-scale modeling render this a massively more difficult computational challenge as compared with MD. One, of course, is the scale; a protein may include thousands or tens of thousands of atoms, while a typical cell will include 2–4 million proteins per cubic micrometer of volume (Milo, 2013), or over 10^{10} in an average mammalian cell, not to mention lipids, nucleic acids, carbohydrates, small metabolites and ions that may all be relevant to propagation of cellular signals. At present, all spatially resolved models for cell biology therefore necessarily include only a very small sampling of the molecular species that may influence the process. Furthermore, most parameters as fundamental as molecular concentrations and reaction rate constants are not known, or known only approximately (Gutenkunst *et al.*, 2007; Schillings *et al.*, 2015). For that reason, cell biological models must typically also perform a broad sampling of “parameter space.” There are increasing efforts to accelerate and optimize this parameter search (Mitra *et al.*, 2018; Shockley *et al.*, 2018; Mitra *et al.*, 2019), but the size of parameter space nonetheless limits the depth to which any particular instantiation of a spatially resolved model can be explored.

### There is no common framework for sharing and reproducing spatially resolved models

Over the past decade, there has been a humbling realization that a large fraction of published results in experimental biology, including many that have led directly to significant efforts toward drug development, cannot be reproduced (Prinz *et al.*, 2011; Begley and Ellis, 2012). This has led to widespread reform in both the performance and the publication of experimental research, notably with the US National Institutes of Health (NIH) now requiring explicit training in experimental rigor and reproducibility and altering the structure of research grant applications to address these issues directly (Health, 2018). In order for the modeling enterprise to advance, it will be necessary for our community to establish standards for sharing and reproducing spatially resolved computational models as well.

Sharing and reproducing models are only possible if every quantitative parameter is strictly defined in a meaningful way. This will require a common “language” that is flexible enough to describe a wide variety of complex cell biological problems, and it is particularly challenging to define parameters that can be implemented in an equivalent way for both PDE-based and particle-based simulations of the same process. Even within the confines of our small working group, where we had all agreed to work together closely and were in constant communication throughout the process, we frequently found it challenging to run the “same” model on different platforms because of inconsistent parameter definitions. We also found that published models often lack critical parameters even if they appear at first glance to be completely described.

For time-resolved cell biological models that do not include an explicit spatial component, the SBML represents a large-scale community effort to improve standardization of model representation. An extension to SBML explicitly designed for spatially resolved modeling, SBML-spatial, is nearing completion (SBML.org) and is supported to varying degrees by some modeling systems, although the specification has not yet been formally released and is not widely adopted. Some leaders in the field have begun to organize resources to aid in efforts toward verifying simulation results (https://reproduciblebiomodels.org). Several related fields that require highly complex computational models, such as finite element modeling in biomechanics, have also begun to define community standards for reporting and verifying model parameters and results (Erdemir *et al.*, 2012).

Minimal steps toward ensuring reproducibility in cell biological spatial modeling must include a general expectation that all code will be readily accessible (via GitHub or equivalent) and directly runnable by outside users, not dependent on local files or libraries. All code for spatial modeling should include explicitly defined unit tests, such as those we have included here. Furthermore, all parameters used in any model must include meaningful units. This is not a trivial consideration; for example, if the problem being studied involves a transition between a soluble form and a membrane-associated form of a particular molecular species, we have found that rate constants and their units are usually not carefully documented in most published models, and in some cases it is not even clear that different researchers agree on what the “right” units are for some typical problems.

### It is difficult to compare outcomes of different types of models and to compare model results to experimental data

A typical goal for cell biological modeling is to determine whether a particular underlying mechanistic hypothesis is consistent with experimental observation. For processes with a spatially dependent stochastic component, which includes many cell biological problems of interest, each computational run (and indeed each independent experiment) will give a slightly different result. To compare experimental results with simulations, and to compare outcomes of different simulations with one another, the simplest approach is to use some kind of summary statistics that described probabilities of various outcomes. However, it is not often trivial to determine which summary statistics are most appropriate. For example, one influential early model illustrating the importance of stochasticity in genetic regulation leading to cell fate determination was an exploration of the events underlying the switch for bacteriophage lambda from a lysogenic to a lytic state, published in 1998 (Arkin *et al.*, 1998). However, more than a decade passed before a useful method for calculating summary statistics on these extremely rare switching events was developed (Morelli *et al.*, 2009). Ironically, it may in many cases be easier to repeat an experiment for thousands or tens of thousands of cells, for example using automated videomicroscopy (Cai *et al.*, 2018) than to run a computationally intensive spatial simulation the same number of times. In general, the field would benefit from greater attention to scalability in model construction.

Some statistical strategies that were originally developed in other fields have been creatively adapted for applications in spatially defined stochastic simulations. These include extensions to the singular spectrum analysis methods that are widely used for analysis of time-series data (Shlemov *et al.*, 2015) and the weighted ensemble sampling method that is useful when rare events are of particular interest (Donovan *et al.*, 2016).

Visualization of model outcomes should also be standardized. This is much more conceptually challenging for spatially resolved models than for models where the only output is a change in concentration of molecular species over time. Here recent developments in the field of visual analytics, which combines statistical analysis with visualization and user interactions to cope with large complex data sets, may offer some promise (Matković *et al.*, 2018). Visual analytics take the human into the loop for exploring data and the space of simulation experiments. Other approaches aim at testing requirements and confirming behavioral expectations automatically. One example of a promising avenue for automatic analysis would be spatio-temporal model checking approaches, which rely on explicitly stating spatio-temporal properties to be checked on the produced trajectories of the spatial simulation (Bartocci *et al.*, 2015).

### It may not be possible to address distinct model features required for different kinds of cell biological problems within a single modeling framework

In our comparative work, we emphasized simple biochemical interactions because we were able to implement these examples in comparable ways across multiple different computational platforms. We chose to use only very simple geometries for our example cases, but some of the stochastic simulation platforms we considered are capable of incorporating complex cell geometries, and complex geometries are also standard in most PDE solvers. However, none of the packages we used here are designed to incorporate correct subcellular physics for cytoskeletal mechanics, force generation, fluid flow, or electrostatics. Clearly these physical elements are important in accurate simulations of cellular behaviors. It has been widely recognized for over a century that cellular organization is tightly coupled to cellular structure (Abbot, 1916).

In our opinion, the highest priorities for integration of correct subcellular physics with stochastic molecule-based methods of the kind described here are accurate mechanical models for cytoskeletal filaments, molecular motors, and membranes. Interactions of the cytoskeleton and its associated motors with membranes determine overall cell shape and organelle localization, and these also underlie interactions of cells with one another and with their extracellular matrix. Actions of cytoskeletal filaments and motors directly deform membranes to produce invaginations or protrusions, and membrane surfaces influence cytoskeletal filament growth both because of their mechanical resistance (Keren *et al.*, 2008) and because of their ability to accumulate and spatially organize key regulators (Mullins *et al.*, 2018). Simulation frameworks that could properly deal with these kinds of mechanical and biochemical interactions should be adaptable for a wide variety of related problems, for example, exploring the assembly of viral “replication factories” on intracellular membranes of host cells, a process which requires large-scale oligomerization of viral polymerase (Spagnolo *et al.*, 2010).

While we recognize that all models are necessarily approximations, we believe that some approximations are better than others. For membrane mechanics, the Helfrich model (Helfrich, 1973) is the most widely used framework for calculation of energies associated with membrane bending and deformation (Guckenberger and Gekle, 2017). While the theory itself is simple and elegant, calculation of energies for highly complex geometries in a continuum mechanics model can become complicated. A few interesting approaches have been used to enable cellular geometries to adapt in response to biochemical reactions (Angermann *et al.*, 2012; Tanaka *et al.*, 2015), but so far none of them have attempted to incorporate correct Helfrich-derived bending energies. PDE-based models for membrane–protein interactions can better approximate the correct underlying physics, but have so far been limited to fairly simple geometries (Rangamani *et al.*, 2014; Wu *et al.*, 2018). Ultimately, the membrane mechanics must be coupled to the proteins and protein assemblies that provide most of the force necessary to induce membrane bending and remodeling.

Capturing the dynamics of protein assemblies, both as highly ordered filaments or lattices (e.g., microtubules or clathrin-coated vesicles) and as disordered condensates, will be essential for bridging biochemical interactions with mechanics. Recent advances in extending single-particle methods to include rigid (Varga *et al.*, 2020) or polymer-like structures (Michalski and Loew, 2016; Hoffmann *et al.*, 2019) have enabled simulations of highly ordered clathrin lattices and viral shells (Varga *et al.*, 2020), or disordered condensate-like assemblies (Chattaraj *et al.*, 2019), respectively. For these simulations, excluded volume is a critical feature to reproduce realistic molecular features and prevent unphysical overlap between multicomponent assemblies. For cytoskeletal mechanics, in many cases a simple approximation of cytoskeletal filaments as elastic rods and motor proteins as discrete localized generators of quantized directional force can give impressively realistic simulation outcomes (Nedelec and Foethke, 2007). This approach has been widely used for studies of dynamic microtubule and actin structures in eukaryotic cells (Karsenti *et al.*, 2006; Odell and Foe, 2008; Rubinstein *et al.*, 2009). It has recently been demonstrated that mechanically accurate modeling of cytoskeletal filament dynamics can be fruitfully combined with a continuum model for membrane bending in order to simulate endocytosis (Akamatsu *et al.*, 2020).

More generally, it is clear that much more work needs to be done in the development of multiscale simulations and hybrid computational approaches (Groen *et al.*, 2019). For example, in some intermediate regimes, approximate subvolume-based simulations may be able to bridge the gap between well-mixed and particle-based simulations. It will require careful work to understand what approximations are generally useful when jumping across scales in hybrid models and what approximations give rise to propagating errors.

Overall, the members of our working group have thoroughly enjoyed this opportunity to critically evaluate the state of tools that are currently available for simulation of cell biological processes. Although we have of course identified important gaps, we are generally optimistic about recent developments and look forward to incorporation of more sophisticated computational methods to improve the efficiency and accuracy of cell biological modeling as these methods are developed. Spatially accurate cellular simulation is a field still in its infancy, but its future is bright.

## MATERIALS AND METHODS

Request a protocol through Bio-protocol.

### Simulation solvers

ODE, PDE, and SSA simulations were all run using Virtual Cell, version 7 (Moraru *et al.*, 2008). Unless otherwise noted, the PDE solver used a fully implicit finite volume, regular grid (variable time step), max step of 0.1 s. Absolute error tol: 10^{–9}, relative error: 10^{–7}. The mesh used the default values, which typically resulted in mesh side lengths of ∼0.02 μm. We note that for the PDE, the mesh, volume, or concentrations sometimes had to be adjusted slightly to ensure the copy numbers were the same for accurate comparisons. The ODE used the Combined Stiff Solver (IDA/CVODE), abs error tol: 10^{–9}, relative error tol: 10^{–9}. SSA simulations used the Gibson–Bruck method, an exact SSA, implemented in Virtual cell, or for the autophosphorylation model, RuleBender (Smith *et al.*, 2012).

All single-particle solvers propagated BD, with no additional forces added. Smoldyn simulations were run either through Virtual Cell, version 7, or using the stand-alone Smoldyn software, version 2.60/2.61 (Andrews, 2017). FPR or NERDSS simulations were run using the NERDSS software, version 1 (Varga *et al.*, 2020), which uses the FPR algorithms (Johnson and Hummer, 2014; Yogurtcu and Johnson, 2015; Johnson, 2018) to solve the single-particle reaction–diffusion model. GFRD simulations were run using eGFRD software (Sokolowski *et al.*, 2019) from August 3, 2019. MCell simulations were run using MCell version 3.3 (Kerr *et al.*, 2008).

### Models and simulation parameters

Unless otherwise noted, species in spatial simulations were initialized as uniformly distributed.

#### 1. Units tests.

*1a. 3D Reversible binding A+B C.* Model parameters: V = 3.1934 μm^{3} cube, [A]_{0} = [B]_{0} = 52 μM, or 100,000 copies, [C]_{0} = 0. k_{on} = 1.476 × 10^{7} M^{–1}s^{–1}, k_{off} = 0.02451 s^{–1}. Microscopic rates, k_{a} = 1000 nm^{3}μs^{–1}, k_{b} = 1s^{–1}. K_{D} = 0.00166 μM. For FPR/NERDSS, *σ* = 1 nm. D_{A} = D_{B} = D_{C} = 1 μm^{2}/s. [A]_{eq} = 0.293 μM, or 563.5 copies, reached after ∼1 s. *Sim parameters:* single-particle methods ∆t = 10^{–7}-10^{–6} s. FPR/NERDSS N_{traj} = 18, SSA N_{traj} = 10, and Smoldyn N_{traj} = 1.

*1b. 2D Reversible binding A+B ⇌ C.* Model parameters: A = 1 μm^{2}, flat surface. [A]_{0} = [B]_{0} = 1000 μm^{–2}, or 1000 copies, [C]_{0} = 0. k_{on} = 3.07 μm^{2}s^{–1}, k_{off} = 3.07 s^{–1}. Microscopic rates, k_{a} = 10 μm^{2}s^{–1}, k_{b} = 10 s^{–1}. K_{D} = 0.00166 μM. For FPR/NERDSS, *σ* = 1 nm. D_{A} = D_{B} = 1 μm^{2}/s, D_{C} = 0.5 μm^{2}/s. [A]_{eq} = 31.1267 μm^{–2}, reached after ∼0.1 s. *Sim parameters:* single-particle methods ∆t = 10^{–7}s. FPR/NERDSS N_{traj} = 5, SSA N_{traj} = 5, Smoldyn N_{traj} = 2.

*1c. 3D→2D Reversible binding A+B ⇌ C.* Model parameters: V = 1 μm^{3} cube, flat membrane surface of 2.2 μm × 2.2 μm. [A]_{0} = 1 μM, or 602 copies, [B]_{0} = 6045 μm^{–2}, or 29,258 copies, [C]_{0} = 0. k_{on} = 8.44 × 10^{7} M^{–1}s^{–1}, k_{off} = 70.085 s^{–1}. Microscopic rates, k_{a}^{3D} = 500 nm^{3}μs^{–1}, k_{b} = 250 s^{–1}. K_{D} = 0.83 μM. *σ* = 1 nm. D_{A} = 30 μm^{2}/s, D_{B} = 1 μm^{2}/s, D_{C} = 0.97 μm^{2}/s. [A]_{eq} = 10.3 copies, or 0.017 μM, reached after ∼0.02 s. *Sim parameters:* single-particle methods ∆t = 10^{–7}s. FPR/NERDSS N_{traj} = 10, SSA N_{traj} = 5, and Smoldyn N_{traj} = 3. Smoldyn ∆t = 10^{–6} s, with the longer steps converging to the PDE solution.

*1d. Crowding A+B→B+C.* Model parameters: V = 12,500 nm^{3} cube (side length of 23.21 nm), periodic boundary conditions enforced. [B] = 100 copies (13.3 mM). [A]_{0} = 100 copies. [C]_{0} varied across 6 systems: [0, 994, 2187, 3381, 4575, 5768], corresponding to total volume fractions occupied of [0.00838, 0.05, 0.1, 0.15, 0.2, 0.25]. Reactions: 1. A+B→B+C k_{a} = 85 nm^{3}/μs, 2. A+A k_{a} = 0, 3. A+C k_{a} = 0, 4. B+B k_{a} = 0, 5. B+C k_{a} = 0, 6. C+C k_{a} = 0. With *no volume exclusion,* Rxn 1 has k_{on} = 63.5 nm^{3}/μs = 3.82 × 10^{7} M^{–1}s^{–1}. *σ* = 1 nm, and thus each particle is effectively modeled as a volume-excluding sphere of diameter = 1 nm. D_{A} = D_{A} = D_{c} = 10 μm^{2}/s.

*Sim parameters:* for FPR/NERDSS(Varga *et al.*, 2020), a maximal ∆t for each system was limited by the density, values for each crowding fraction: [10^{–11} s, 10^{–11} s, 10^{–11} s, 5 × 10^{–12} s, 5 × 10^{–12} s, 10^{–12} s]. N_{traj} = 40–80 for each crowding fraction. Another N_{traj} = 40–80 were run at another time step (5 times faster) to verify quantitatively similar results. For eGFRD (Sokolowski *et al.*, 2019), N_{traj} = 10 for each crowding fraction. Algorithmic adjustments to standard FPR: first, when an A+B→B+C reaction occurred, the products were kept at their same coordinates as the reactants. Second, at each time step, particle positions are updated, and they must avoid overlapping (r ≥ *σ*) all other particles, considering all particles within collision distance (r < *R _{max}*) (Johnson and Hummer, 2014). Due to the high crowding, we used a more stringent criteria to update particle positions. We defined clusters that contained all the particles with a capacity to possibly overlap all their partners. The positions of all particles in this cluster were simultaneously updated to avoid any overlap. For large clusters, to reduce the number of rejected updates, the time step was reduced by a factor of 10 and then updated 10 times to the full ∆t. Third, initial configurations for high crowded simulations were generated on a grid. This crystalline configuration was then “melted” by simply setting the rate of reaction 1 to zero to prevent binding events. After 0.1 s, the now disordered but well-mixed configuration was used to start simulations.

#### 2. Intermediate tests.

*2a. Membrane localization model.* Model parameters: V = 1.1045 μm^{3}, A = 0.2209 μm^{2}, rectangular cube, 0.47 μm × 0.47 μm surface, height = 5 μm. [A]_{0} = [B]_{0} = 1 μM or 665 copies, [M]_{0} = 5.6474 μM = 17,000/μm^{2}, or 3755 copies. [AM]_{0} = [BM]_{0} = [AB]_{0} = [ABM]_{0} = . [MAB]_{0} = [MABM]_{0} = 0. Reactions: 1. A+B *⇌* AB k_{on} = 0.05 μM^{–1}s^{–1}. 2. A+M *⇌* AM k_{on} = 2 μM^{–1}s^{–1}. 3. B+M *⇌* BM k_{on} = 2 μM^{–1}s^{–1}*.* 4. A+BM *⇌* ABM k_{on} = 0.05 μM^{–1}s^{–1}. 5. B+MA *⇌* MAB k_{on} = 0.05 μM^{–1}s^{–1}. 6. AB+M *⇌* ABM k_{on} = 2 μM^{–1}s^{–1}. 7. AB+M *⇌* MAB k_{on} = 2 μM^{–1}s^{–1}. 8. MA+BM *⇌* MABM 0.0415 μm^{2}s^{–1}. 9. MAB+M *⇌* MABM 1.6693 μm^{2}s^{–1}. 10. M+ABM *⇌* MABM 1.6693 μm^{2}s^{–1}. For all reactions, k_{off} = 1 s^{–1}. At equilibrium, 90.3% of proteins were in complex with one another on the membrane (MABM), or 600 copies, and 1.16% were remaining in solution, or 7.7 copies. *σ* = 1 nm for all rxns, D_{A} = D_{B} = 50 μm^{2}/s. D_{AB} = 25 μm^{2}/s. D_{AM} = D_{BM} = 0.495 μm^{2}/s. D_{ABM} = D_{MAB} = 0.49 μm^{2}/s. D_{MABM} = 0.248 μm^{2}/s D_{M} = 0.5 μm^{2}/s. All membrane bound had D_{z} = 0.

*Sim parameters:* single-particle methods ∆t = 10^{–7}s. FPR/NERDSS N_{traj} = 10. SSA N_{traj} = 10 Smoldyn N_{traj} = 10. We note that the model includes 10 coupled pairwise reactions (three are 3D, four are 3D→2D, and three are 2D), although only 6 equilibrium constants are unique. This is because there are only 3 reactants and multiple thermodynamic cycles that can result in the final assemblies (see, e.g., MABM in Figure 5a). To preserve a thermodynamic equilibrium at steady state, not all rates can be defined independently.

*2b. Autophosphorylation:* model parameters A: V = 0.003 μm^{3}. [A]_{0} = 103 copies, [Ap]_{0} = 5 copies, [P] _{0} = 9 copies. [A_Ap]_{0} = 0, [Ap_P]_{0} = 0. 1. A→Ap at 21.2 s^{–1}. 2. A+Ap→A_Ap at 10^{7}M^{–1}s^{–1}, k_{a} = 16.7 nm^{3}μs^{–1}, *σ* = 1 nm. 3.A_Ap→Ap+Ap at 200 s^{–1}. 4. Ap+P→Ap_P at 8*10^{8} M^{–1}s^{–1}, k_{a} = 2820 nm^{3}μs^{–1}, *σ* = 1 nm and 5. Ap_P→A+P at 539 s^{–1}. D_{A} = D_{P} = 100 μm^{2}/s. A second spatial stochastic simulation had D_{A} = D_{P} = 10 μm^{2}/s, which required new bimolecular parameters for rxn 2: k_{a} = 17.8 nm^{3}μs^{–1}, *σ* = 1 nm, rxn 4: k_{a} = 11191 nm^{3}μs^{–1}, *σ* = 6 nm. Note: the binding radius had to be increased to recover the large macroscopic binding rates for rxn 4. A third spatial stochastic simulation had D_{A} = D_{P} = 20 μm^{2}/s, which required new bimolecular parameters for rxn 2: k_{a} = 17.2 nm^{3}μs^{–1}, *σ* = 1 nm, rxn 4: k_{a} = 3919 nm^{3}μs^{–1}, *σ* = 4 nm.

*Sim parameters:* single particle methods: ∆t = 2*10^{–8}s, N_{traj} = 2. For D = 10, ∆t = 1*10^{–7}s, N_{traj} = 2. For D = 20, ∆t = 1*10^{–7}s, N_{traj} = 2. All simulations were ≥200 s to obtain statistics on switching times.

B version of model with all rates reduced by 10: 1. 2.12 s^{–1}. 2. 10^{6} M^{–1}s^{–1}, k_{a} = 1.67 nm^{3}μs^{–1}, *σ* = 1 nm. 3. 20 s^{–1}. 4. 8*10^{7} M^{–1}s^{–1}, k_{a} = 282.0 nm^{3}μs^{–1}, *σ* = 1 nm and 5. 53.9 s^{–1}. D_{A} = D_{P} = 10 μm^{2}/s. A second spatial stochastic simulation had D_{A} = D_{P} = 100 μm^{2}/s, which required for rxn 4: k_{a} = 140.31 nm^{3}μs^{–1}, *σ* = 1 nm. A third spatial stochastic simulation had D_{A} = D_{P} = 1 μm^{2}/s, which required for rxn 2: k_{a} = 1.78 nm^{3}μs^{–1}, *σ* = 1 nm. rxn 4: k_{a} = 1119.14 nm^{3}μs^{–1}, *σ* = 6 nm.

*Sim parameters:* ∆t = 2*10^{–7}s, N_{traj} = 2. For D = 100, ∆t = 2*10^{–8}s, N_{traj} = 2. For D = 1, ∆t = 1*10^{–6} s, N_{traj} = 2. All statistics were collected for at least 2000 s, except the fast diffusion constant. With only 100 s for D = 100, we observed transitions but insufficient statistics to determine trends relative to the other models.

#### 3. Application tests.

*3a. Circadian clock oscillator model.* Model parameters A: 4.189 μm^{3}, cube, or a sphere with R = 1 μm. [prmA]_{0} = [prmR]_{0} = 1 copy, or 0.000397 μM. [A]_{0} = [R]_{0} = [C]_{0} = [mRNAA]_{0} = [mRNAR]_{0} = [prmA.A]_{0} = [prmR.A]_{0} = 0. Reactions: 1. A + R → C k_{on} = 1204 µM^{–1}s^{–1}, k_{a} = 356,263 nm^{3}μs^{–1}, s = 8 nm. 2. prmA → prmA + mRNAA k_{2} = 50 s^{–1}. 3. prmA.A → prmA.A + mRNAA k_{3} = 500 s^{–1}. 4. prmR → prmR + mRNAR k_{4} = 0.01 s^{–1}. 5. prmR.A → prmR.A + mRNAR k_{5} = 50 s^{–1}. 6. prmA + A ↔ prmA.A k_{6f} = 602 µM^{–1}s^{–1}, k_{6b} = 50 s^{–1}, k_{a} = 4889 nm^{3}μs^{–1}, k_{b} = 244.5 s^{–1} s = 5 nm. 7. prmR + A↔ prmR.A k_{7f} = 602 µM^{–1}s^{–1}, k_{7b} = 100 s^{–1}, k_{a} = 4889 nm^{3}μs^{–1}, k_{b} = 489 s^{–1}*σ* = 5 nm. 8. mRNAA → mRNAA + A k_{8} = 50 s^{–1}. 9. mRNA_R → mRNA_R + R k_{9} = 5 s^{–1}. 10. mRNA_A → NULL k_{10} = 10 s^{–1}. 11. mRNA_R → NULL k_{11} = 0.5 s^{–1}. 12. A → NULL k_{12} = 1 s^{–1}. 13. R→ NULL k_{13} = 0.2 s^{–1}. 14. C→ R k_{14} = 1 s^{–1}. D = 10 μm^{2}/s (for all 9 species).

B version of the model was run with one modification from the original A model: 13. R→ NULL k_{13} = 0.05 s^{–1}, eliminating oscillations in the deterministic model Supplemental Figure S3.

C version of this model was created with modifications from original A model: a) promoters were spatially localized to a sphere at the center of the volume of R = 0.1 μm. b) D_{prmA} = D_{prmR} = D_{prmA.A} = D_{prmA.R} = 0. c) [prmA]_{0} = [prmR]_{0} = 1 copy, or here 0.3964 μM. No barrier existed for proteins to access this ‘nucleus’ in the cell center. Full details are in the Supplemental Supporting Methods.

*Sim parameters:* NERDSS ∆t = 0.5, 2, and 10 μs, and Smoldyn ∆t = 100 μs. SSA N_{traj} = 10. Smoldyn N_{traj} = 2, NERDSS N_{traj} = 2 (600 s each); see Supplemental Tables S6 and S7. Time steps were limited by the fast unbinding reaction of the single promoter, Rxn 7, rather than the particle density. With only a single event possible due to *N* = 1, a larger time step results in slightly fewer events occuring relative to the theoretical value.

*3b. Model: MinCDE cell division model.* Model parameters: cylindrical volume, h = 4 μm, R = 0.5 μm. Five species: [MinD-ATP]_{0} = 2.1143 μM, or 4000 copies. [MinE]_{0} = 0.74 μM, or 1400 copies. [MinD-ADP]_{0} = [MinD-ATP2D]_{0} = [MinD-ATP.MinE2D]_{0} = 0, where 2D indicates membrane bound. Six irreversible reactions. 1. MinD-ATP→MinD-ATP2D, *k*_{1} = 0.025 µm.s^{–1}. 2. MinD-ATP + MinD-ATP2D→2MinD-ATP2D, k_{2} = 0.903 s^{–1}.µM^{–1}. 3. MinD-ATP + MinD-ATP.MinE2D→MinD-ATP2D + MinD-ATP.MinE2D, k_{3} = 0.903 s^{–1}.µM^{–1}. 4. MinE+MinD-ATP2D→MinD-ATP.MinE2D, k_{4} = 56.0 s^{–1}.µM^{–1}. 5. MinD-ATP.MinE2D→MinD-ADP + MinE, k_{5} = 0.7 s^{–1}. 6. MinD-ADP→MinD-ATP, k_{6} = 1 s^{–1}. *D* = 2.5 μm^{2}/s for MinD-ATP, MinD-ADP, and MinE. D = 0 μm^{2}/s for MinD-ATP2D and MinD-ATP.MinE2D in the PDE. Initial species all well-mixed in PDE. In Smoldyn, MinD-ATP copies split between solution and membrane. Necessary modification for Smoldyn to produce oscillations: *D* = 0.05 μm^{2}/s for MinD-ATP2D and MinD-ATP.MinE2D.

*Sim parameters:* Smoldyn: N_{traj} = 4, ∆t = 2 ms.

### Analysis methods

*Equilibria and error:* for all stochastic methods, we calculated averages and variances of copy numbers in time over *N* trajectories. We used the standard error of the mean at time points: SEM(*t*) = . For equilibrium values, we used differences from the mean to determine when equilibrium was reached. We note that for all models, we selected parameters such that the equilibrium was not too close to zero, which is important for assessing clear deviations from the proper equilibrium values.

*Clock oscillator timescales*: the average time separation between A peaks, between R peaks, and the lag time between A and R peaks was calculated in two ways. First, by finding the peak maxima and simply calculating the separation between neighboring peaks, this was then averaged over all peaks. Second, to calculate time period of peaks, we used a discrete FFT that was zero-padded out to 5000 s to increase the resolution of the frequencies sampled. The frequency with the largest coefficient, f_{max}, was used to define the period of the oscillation (1/f_{max}). The cross-correlation between the A and R time-series was used to identify the lag-time, again keeping the time with the largest correlation coefficient. Both methods produced similar results, and were performed using MATLAB functions. We report in the Supplemental tables the periods determined from the find peaks method, with SEMs reporting the calculated .

*Fit to kinetics of crowding simulations*: to define the k_{fit} values for the crowding simulations as a function of packing fraction, we fit the copies of A versus time *A*(*t*) = *A*_{0} exp (–*k _{fit} B_{tot}t*), where B

_{tot}is 0.008/nm

^{3}, and we performed nonlinear fitting to the exponential rather than the logarithm of

*A*(

*t*). We fit each individual trajectory and report the mean of these values , where N is total trajectories. The SEM is calculated based on the variance around this mean value. Fits to the averaged trajectory gave similar values, and fits using the ln(

*A*(

*t*)) produced the same trends, but with slightly lower

*k*

_{fit}values. The same approach was used for both FPR and GFRD trajectories.

*MinCDE timescales:* for the MinCDE models, we calculated the time-scales of oscillations by choosing fixed points in space at either end of the cylinder (Supplemental Figure S7). Once the oscillations were pole-to-pole, we quantified the period of oscillations in time using the same approach as applied to the clock oscillator model above by finding the peak maxima and taking the average of time separations between adjacent peaks. This produced periods of 43.0 ± 0.7 s for Smoldyn and 41.3 ± 0.3 s for the PDE.

*Autophosphorylation state assignments:* to assign stochastic trajectories to states, 1 = low phosphorylation, 2 = high phosphorylation, we report results based on thresholding the values of *A*(*t*) and Ap(t). We created a transition region, where points could be in either 1 or 2 depending on their previous state, to minimize rapid recrossings between states (see Supplemental Figure S1). If [*A*(*t*) > 62 + Ap(t)], assign to state 1, else if [*A*(*t*) < 35 + Ap(t)] assign state 2, else [if(State(t-1) = *i*, State(t) = *i*], where *i* = [1,2]. We then counted the number of times in each state, length of intervals spent in each state, and transitions between points separated by a time *δ*t to construct . We calculated , and , where *N* is the number of intervals spent in state i, each of duration τ_{i}. Distributions of these residence times are shown in Supplemental Figure S1. To calculate errors on the probability of each state, each trajectory was split into 10 chunks (length at least the larger residence time), and these 10 state probabilities were used to calculate a mean and SD/SEM.

## FOOTNOTES

This article was published online ahead of print in MBoC in Press (http://www.molbiolcell.org/cgi/doi/10.1091/mbc.E20-08-0530) on November 25, 2020.

**Abbreviations used:**

BD | Brownian dynamics |

CME | chemical master equation |

GF | Green’s function |

MD | molecular dynamics |

ODE | ordinary differential equation |

PDE | partial differential equation |

RDME | reaction diffusion master equation |

SBML | systems biology markup language |

SSA | stochastic simulation algorithm. |

## ACKNOWLEDGMENTS

This work was funded by NIMBioS. We thank Carlos Lopez for stimulating discussions. M.E.J. also gratefully acknowledges computing resources from Johns Hopkins University’s MARCC supercomputer and XSEDE through XRAC MCB150059 and support from an NIH MIRA R35GM133644. M.M.S. and T.P. acknowledge support by the intramural research program of the National Institute of Allergy and Infectious Diseases, NIH. P.H. and A.M.U. acknowledge support from the DFG (German Research Foundation) Collaborative Research Center 1270/1 - 299150580 ELAINE.

## REFERENCES

**Am J Psychol**

*27*, 245–250. Crossref, Google Scholar (

**J Phys Chem B**

*116*, 3630–3640. Crossref, Medline, Google Scholar (

**Structural Chemistry and Molecular Biology**, San Francisco: Freeman, 198–215. Google Scholar (

**J Chem Phys**

*137*, 044105. Crossref, Medline, Google Scholar (

**J Chem Phys**

*92*, 5270–5284. Crossref, Google Scholar (

**Elife**

*9*. Crossref, Medline, Google Scholar (

**Nat Cell Biol**

*8*, 1195–1203. Crossref, Medline, Google Scholar (

**Science**

*307*, 423–426. Crossref, Medline, Google Scholar (

**Proc Natl Acad Sci USA**

*107*, 18457–18462. Crossref, Medline, Google Scholar (

**Bioinformatics**

*33*, 710–717. Crossref, Medline, Google Scholar (

**Phys Biol**

*17*, 045001. Crossref, Medline, Google Scholar (

**PLoS Comput Biol**

*6*, e1000705. Crossref, Medline, Google Scholar (

**Phys Biol**

*1*, 137–151. Crossref, Medline, Google Scholar (

**Phys Biol**

*1*, 137–151. Crossref, Medline, Google Scholar (

**Nature Methods**

*9*, 283–289. Crossref, Medline, Google Scholar (

**Genetics**

*149*, 1633–1648. Medline, Google Scholar (

**Biophys J**

*90*, 65–76. Crossref, Medline, Google Scholar (

**Biophys J**

*66*, 588–600. Crossref, Medline, Google Scholar (

**Hybrid Systems Biology. HSB 2015. Lecture Notes in Computer Science, tructural Chemistry and Molecular Biology**, ed. A AbateD Šafránek, Springer. Crossref, Google Scholar (

**PLoS Comput Biol**

*12*, e1004591. Crossref, Medline, Google Scholar (

**Nature**

*483*, 531–533. Crossref, Medline, Google Scholar (

**Biophys J**

*20*, 193–219. Crossref, Medline, Google Scholar (

**Bioinformatics**

*34*, i583–i592. Crossref, Medline, Google Scholar (

*et al.*(2018). Experimental and computational framework for a dynamic protein atlas of human cell division.

**Nature**

*561*, 411–415. Crossref, Medline, Google Scholar ,

**Bull Math Biol**

*78*, 617–661. Crossref, Medline, Google Scholar (

**Science**

*263*, 802–805. Crossref, Medline, Google Scholar (

**Biophys J**

*116*, 560–572. Crossref, Medline, Google Scholar (

**Front Neuroinform**

*11*, 13. Crossref, Medline, Google Scholar (

*et al.*(2014). Lattice light-sheet microscopy: Imaging molecules to embryos at high spatiotemporal resolution.

**Science**

*346*, 1257998. Crossref, Medline, Google Scholar ,

**Proc Natl Acad Sci USA**

*108*, 1052–1057. Crossref, Medline, Google Scholar (

**Phys Rev E**

*98*, 032418. Crossref, Google Scholar (

**Phys Rev E**

*99*, 042411. Crossref, Medline, Google Scholar (

**Science**

*361*, 412–415. Crossref, Medline, Google Scholar (

**PLoS Biol**

*2*, e439. Crossref, Medline, Google Scholar (

**J Colloid Sci**

*4*, 425–437. Crossref, Google Scholar (

**Sci Rep**

*9*, 11676. Crossref, Medline, Google Scholar (

**Biophys J**

*101*, 1412–1422. Crossref, Medline, Google Scholar (

**J Chem Phys**

*151*. Crossref, Medline, Google Scholar (

**J Phys a-Math Gen**

*9*, 1479–1495. Crossref, Google Scholar (

**PLoS Comput Biol**

*12*, e1004611. Crossref, Medline, Google Scholar (

*et al.*(2016). Stochastic simulation service: bridging the gap between the computational expert and the biologist.

**PLoS Comput Biol**

*12*, 1005220. Crossref, Google Scholar ,

**Annu Rev Biophys**

*41*, 429–452. Crossref, Medline, Google Scholar (

**Biophys J**

*98*, 1571–1581. Crossref, Medline, Google Scholar (

**PLoS One**

*10*, e0121681. Crossref, Medline, Google Scholar (

**Phys Biol**

*6*, 046001. Crossref, Medline, Google Scholar (

**J Biomech**

*45*, 625–633. Crossref, Medline, Google Scholar (

**J Chem Phys**

*69*, 1352–1360. Crossref, Google Scholar (

**Trends Cell Biol**

*14*, 661–669. Crossref, Medline, Google Scholar (

**BMC Biol**

*9*, 68. Crossref, Medline, Google Scholar (

**Proc Natl Acad Sci USA**

*107*, 19820–19825. Crossref, Medline, Google Scholar (

**Soft Matter**

*12*, 7792–7803. Crossref, Medline, Google Scholar (

**J Chem Phys**

*151*, 124115. Crossref, Medline, Google Scholar (

**Chemical Master Equation,**New York: Springer. Crossref, Google Scholar (

**PLoS Comput Biol**

*16*, e1007717. Crossref, Medline, Google Scholar (

**J Comput Phys**

*22*, 403–434. Crossref, Google Scholar (

**Physica A**

*188*, 404–425. Crossref, Google Scholar (

**J Chem Phys**

*138*, 170901. Crossref, Medline, Google Scholar (

**Chemical Physics**

*284*, 91–102. Crossref, Google Scholar (

**J Chem Phys**

*117*, 507–517. Crossref, Google Scholar (

**Phys Rev Lett**

*86*, 922–925. Crossref, Medline, Google Scholar (

**Chem Rev**

*118*, 11707–11794. Crossref, Medline, Google Scholar (

**Cell**

*130*, 141–152. Crossref, Medline, Google Scholar (

**Essays Biochem**

*45*, 41–56. Crossref, Medline, Google Scholar (

**Philos Trans A Math Phys Eng Sci**

*377*, 20180147. Medline, Google Scholar (

**J Phys Condens Matter**

*29*, 203001. Crossref, Medline, Google Scholar (

**BMC Biol**

*12*, 29. Crossref, Medline, Google Scholar (

*et al.*(2018). Visualizing intracellular organelle and cytoskeletal interactions at nanoscale resolution on millisecond timescales.

**Cell**

*175*, 1430–1442.e1417. Crossref, Medline, Google Scholar ,

**Computation (Basel)**

*6*, 9. Crossref, Medline, Google Scholar (

**PLoS Comput Biol**

*3*, e189. Crossref, Google Scholar (

**Cell Rep**

*1*, 741–752. Crossref, Medline, Google Scholar (

**Parallel Comput**

*40*, 86–99. Crossref, Medline, Google Scholar (

**Bioinformatics**

*32*, 3366–3368. Crossref, Medline, Google Scholar (

**J Math Biol**

*69*, 687–735. Crossref, Medline, Google Scholar (

**Z Naturforsch C J Biosci**

*28*, 693–703. Crossref, Medline, Google Scholar (

**Phys Rev E**

*85*, 042901. Crossref, Google Scholar (

**Phys Rev E Stat Nonlin Soft Matter Phys**

*91*, 023312. Crossref, Medline, Google Scholar (

**J Chem Phys**

*147*, 234101. Crossref, Medline, Google Scholar (

**J Physiol**

*117*, 500–544. Crossref, Medline, Google Scholar (

**PLoS Comput Biol**

*15*, e1006830. Crossref, Medline, Google Scholar (

**Phys Rev Lett**

*90*, 128102. Crossref, Medline, Google Scholar (

**Proc Natl Acad Sci USA**

*100*, 12724–12728. Crossref, Medline, Google Scholar (

*et al.*(2003). The systems biology markup language (SBML): a medium for representation and exchange of biochemical network models.

**Bioinformatics**

*19*, 524–531. Crossref, Medline, Google Scholar ,

**SIAM J Appl Math**

*70*, 77–111. Crossref, Google Scholar (

**J Chem Phys**

*139*. Crossref, Medline, Google Scholar (

**J Comput Phys**

*374*, 954–983. Crossref, Google Scholar (

**J Cell Sci**

*126*, 1913–1921. Crossref, Medline, Google Scholar (

**PLoS Comput Biol**

*7*, e1001121. Crossref, Medline, Google Scholar (

**J Phys Chem B**

*122*, 11771–11783. Crossref, Medline, Google Scholar (

**Physical Review X**

*4*, 031037. Crossref, Medline, Google Scholar (

**J Chem Theory Comput**

*13*, 563–576. Crossref, Medline, Google Scholar (

**Nat Cell Biol**

*8*, 1204–1211. Crossref, Medline, Google Scholar (

*et al.*(2020). SBML Level 3: an extensible format for the exchange and reuse of biological models.

**Mol Syst Biol**

*16*. Crossref, Medline, Google Scholar ,

**Nature**

*453*, 475–480. Crossref, Medline, Google Scholar (

**SIAM J Sci Comput**

*30*, 3126–3149. Crossref, Medline, Google Scholar (

**Bioessays**

*30*, 1110–1125. Crossref, Medline, Google Scholar (

**Trends Cell Biol**

*10*, 173–178. Crossref, Medline, Google Scholar (

**Biophys J**

*96*, 1333–1340. Crossref, Medline, Google Scholar (

**Cell Mol Bioeng**

*1*, 84–92. Crossref, Medline, Google Scholar (

**Proc Natl Acad Sci USA**

*82*, 3055–3057. Crossref, Medline, Google Scholar (

*et al.*(2018). Observing the cell in its native state: Imaging subcellular dynamics in multicellular organisms.

**Science**

*360*, eaaq1392. Crossref, Medline, Google Scholar ,

**Trends Biotechnol**

*19*, 401–406. Crossref, Medline, Google Scholar (

**Mol Syst Biol**

*9*, 646. Crossref, Medline, Google Scholar (

**Nature**

*461*, 99–103. Crossref, Medline, Google Scholar (

**PLoS Comput Biol**

*15*, e1006199. Crossref, Medline, Google Scholar (

**Nat Commun**

*11*, 2821. Crossref, Medline, Google Scholar (

**Rev Mod Phys**

*70*, 979–1001. Crossref, Google Scholar (

**BMC Syst Biol**

*5*, 166. Crossref, Medline, Google Scholar (

**Stochastic Approach to Chemical Kinetics**, London: Methuen. Google Scholar (

**Front Immunol**

*10*, 2268. Crossref, Medline, Google Scholar (

**Biophys J**

*110*, 523–529. Crossref, Medline, Google Scholar (

**BioEssays**

*35*, 1050–1055. Crossref, Medline, Google Scholar (

**Biophys J**

*68*, 1311–1322. Crossref, Medline, Google Scholar (

**J Cell Sci**

*119*, 2863–2869. Crossref, Medline, Google Scholar (

**Nat Commun**

*9*, 3901. Crossref, Medline, Google Scholar (

**iScience**

*19*, 1012–1036. Crossref, Medline, Google Scholar (

**Biophys J**

*89*, 782–795. Crossref, Medline, Google Scholar (

**Iet Systems Biology**

*2*, 352–362. Crossref, Medline, Google Scholar (

**Proc Natl Acad Sci**

*106*, 8101–8106. Crossref, Medline, Google Scholar (

**Biophys Rev**

*10*, 1537–1551. Crossref, Medline, Google Scholar (

**J Chem Phys**

*124*, 044104. Crossref, Medline, Google Scholar (

**New J Physics**

*9*, 427–427. Crossref, Google Scholar (

**PLoS Comput Biol**

*13*, e1005862. Crossref, Medline, Google Scholar (

**J Chem Phys**

*80*, 1517–1526. Crossref, Google Scholar (

**J Comput Phys**

*270*, 203–213. Crossref, Medline, Google Scholar (

**J Cell Biol**

*183*, 471–483. Crossref, Medline, Google Scholar (

**PLoS Comput Biol**

*3*, e36. Crossref, Medline, Google Scholar (

**Nature Reviews Drug Discovery**

*10*, 712–712. Crossref, Medline, Google Scholar (

**J Chem Phys**

*137*, 054104. Crossref, Medline, Google Scholar (

**Proc Natl Acad Sci**

*113*, E5298–E5307. Crossref, Medline, Google Scholar (

**Biophys J**

*107*, 751–762. Crossref, Medline, Google Scholar (

**BMC Syst Biol**

*8*, 52. Crossref, Medline, Google Scholar (

**Proc Natl Acad Sci USA**

*96*, 4971–4976. Crossref, Medline, Google Scholar (

**Bioinformatics**

*25*, 1923–1929. Crossref, Medline, Google Scholar (

**Trends Genet**

*16*, 276–277. Crossref, Medline, Google Scholar (

*283a*, 259–274. Google Scholar (

**J Comput Chem**

*34*, 245–255. Crossref, Medline, Google Scholar (

**Phys Biol**

*6*, 016005. Crossref, Medline, Google Scholar (

**Mol Biol Cell**

*27*, 3616–3626. Link, Google Scholar (

**Nat Biotechnol**

*24*, 1235–1240. Crossref, Medline, Google Scholar (

**Nat Methods**

*10*, 1192–1195. Crossref, Medline, Google Scholar (

**Numerical Computer Methods, Part C, Vol. 321**, New York: Elsevier, 1–23. Crossref, Google Scholar (

**Bioinformatics**

*32*, 2880–2882. Crossref, Medline, Google Scholar (

**PLoS Comput Biol**

*11*, e1004457. Crossref, Medline, Google Scholar (

*50*, 1–60. Google Scholar (

**Plos One**

*8*, e74261. Crossref, Medline, Google Scholar (

**Bmc Biophysics**

*7*, https://doi.org/10.1186/s13628-014-0011-5. Crossref, Medline, Google Scholar (

**Chem Rev**

*109*, 839–860. Crossref, Medline, Google Scholar (

**PLoS Comput Biol**

*13*, e1005857. Crossref, Medline, Google Scholar (

*et al.*(2010). Atomic-level characterization of the structural dynamics of proteins.

**Science**

*330*, 341–346. Crossref, Medline, Google Scholar ,

**BioMed Res Int**

*2015*, Article ID 689745. Google Scholar (

**Bioinformatics**

*34*, 695–697. Crossref, Medline, Google Scholar (

**Biophys J**

*36*, 697–714. Crossref, Medline, Google Scholar (

**Phys Rev E**

*93*, 052135. Crossref, Medline, Google Scholar (

**Bull Math Biol**

*81*, 2960–3009. Crossref, Medline, Google Scholar (

**BMC Bioinformatics**

*13*(Suppl 8), S3. Crossref, Medline, Google Scholar (

**Nature**

*420*, 102–106. Crossref, Medline, Google Scholar (

**J Chem Phys**

*150*, 054108. Crossref, Medline, Google Scholar (

**RNA**

*16*, 382–393. Crossref, Medline, Google Scholar (

*et al.*(2017). Actin retrograde flow actively aligns and orients ligand-engaged integrins in focal adhesions.

**Proc Natl Acad Sci USA**

*114*, 10648–10653. Crossref, Medline, Google Scholar ,

**J Chem Phys**

*72*, 4350–4357. Crossref, Google Scholar (

**FEBS Lett**

*579*, 1783–1788. Crossref, Medline, Google Scholar (

**Proc Natl Acad Sci USA**

*107*, 2473–2478. Crossref, Medline, Google Scholar (

**Bioinformatics**

*31*, 2340–2347. Crossref, Medline, Google Scholar (

**J Phys A-Math Gen**

*38*, R79–R131. Crossref, Google Scholar (

**Nature**

*466*, 267–271. Crossref, Medline, Google Scholar (

**Mol Syst Biol**

*8*, 578. Crossref, Medline, Google Scholar (

*et al.*(1999). E-CELL: software environment for whole-cell simulation.

**Bioinformatics**

*15*, 72–84. Crossref, Medline, Google Scholar ,

*387*, 147–170. Google Scholar (

**Nat Protoc**

*7*, 562–578. Crossref, Medline, Google Scholar (

*237*, 37–72. Google Scholar (

**Curr Opin Cell Biol**

*15*, 221–231. Crossref, Medline, Google Scholar (

**Stochastic Processes in Physics and Chemistry**, 3rd ed., Amsterdam: Elsevier. Google Scholar (

**Phys Rev Lett**

*94*, 128103. Crossref, Medline, Google Scholar (

**Biophys J**

*118*, P3026–P3040. Crossref, Medline, Google Scholar (

**Proc Natl Acad Sci USA**

*99*, 5988–5992. Crossref, Medline, Google Scholar (

**Nature**

*406*, 188–192. Crossref, Medline, Google Scholar (

**Z Phys Chem**

*92*, 129. Google Scholar (

**Science**

*299*, 1231–1235. Crossref, Medline, Google Scholar (

**Front Neuroinform**

*3*, 15. Crossref, Medline, Google Scholar (

**BMC Syst Biol**

*4*, 42. Crossref, Medline, Google Scholar (

**J Comput Phys**

*229*, 7287–7308. Crossref, Medline, Google Scholar (

**Nat Commun**

*9*, 136. Crossref, Medline, Google Scholar (

**Nature**

*475*, 510–513. Crossref, Medline, Google Scholar (

**Phys Biol**

*17*, 011001. Crossref, Medline, Google Scholar (

**J Chem Phys**

*143*, 084117. Crossref, Medline, Google Scholar (

**PLoS Comput Biol**

*14*, e1006031. Crossref, Medline, Google Scholar (

**Bioinformatics**

*29*, 1229–1230. Crossref, Medline, Google Scholar (

**J Phys Chem**

*94*, 8794–8800. Crossref, Google Scholar (

**Biophys J**

*64*, 1711–1726. Crossref, Medline, Google Scholar (

**Annu Rev Biophys**

*37*, 375–397. Crossref, Medline, Google Scholar (

**Biophys J**

*71*, 2440–2457. Crossref, Medline, Google Scholar (

**Annu Rev Biophys Biomol Struct**

*22*, 27–65. Crossref, Medline, Google Scholar (

**Annu Rev Biophys**

*40*, 41–62. Crossref, Medline, Google Scholar (