# Stochastic single-particle based simulations of cellular signaling embedded
into computational models of cellular morphology

Thorsten Prüstel  Computational Systems Biology Section
Laboratory of Immune System Biology
National Institute of Allergy and Infectious Diseases
National Institutes of Health, Bethesda, Maryland 20892, USA  Martin Meier-
Schellersheim  Computational Systems Biology Section
Laboratory of Immune System Biology
National Institute of Allergy and Infectious Diseases
National Institutes of Health, Bethesda, Maryland 20892, USA

###### Abstract

Cells exhibit a wide variety of different shapes. This diversity poses a
challenge for computational approaches that attempt to shed light on the role
cell geometry plays in regulating cell physiology and behavior. The simulation
platform Simmune is capable of embedding the computational representation of
signaling pathways into realistic models of cellular morphology. However,
Simmune’s current approach to account for the cell geometry is limited to
deterministic models of reaction-diffusion processes, thus providing a coarse-
grained description that ignores stochastic local fluctuations. Here we
present an extension of Simmune that removes these limitations by employing an
alternative computational representation of cellular geometry that is smooth
and grid-free. These features make it possible to incorporate a fully
stochastic, spatially resolved description of the cellular biochemistry. The
alternative computational representation is compatible with Simmune’s current
approach for specifying molecular interactions. This means that a modeler
using the approach needs to create a model of cellular biochemistry and
morphology only once to be able to use it for both, deterministic and
stochastic simulations.

11footnotetext: Email<EMAIL_ADDRESS><EMAIL_ADDRESS>
##  1 Introduction

Computational approaches that aim at gaining a mechanistic understanding of
cellular behavior by generating realistic predictive models of cellular
signaling networks face several challenges. First, signaling pathways
dynamically emerge from interactions among many signaling proteins and lipids.
These may form multi-molecular complexes, assembling and disassembling over
time in a highly organized manner [40, 13, 36]. Frequently, signal-transducing
complexes are characterized by several binding sites and internal states that,
in the course of molecular interactions, undergo post-translational
modifications, thereby modifying their biochemical behavior and generating new
interaction possibilities. This dynamically changing, combinatorial complexity
of signaling pathways may result in prohibitively large numbers of possible
functional states, rendering conventional approaches that describe a pathway
by explicitly describing its mass-action kinetics inapplicable [33].

Importantly, cellular signaling networks operate within given cellular
morphologies defining spatial contexts that can only be ignored if one assumes
that the biochemistry is well-stirred and evolves in a static, simple shaped
geometry, such as a cube. However, these assumptions are frequently not
justified, as is obvious when looking at the typically heterogeneous
spatiotemporal distributions of concentrations within cells [41, 45] and the
wide variety of intricate cellular morphologies. Because the specific cellular
shape can act on the biochemistry and vice versa, computational
representations of cellular shape are required that are flexible enough to
reflect this diverse biological reality.

Moreover, having models of cellular biochemistry and computational
representations of a cell’s morphology at our disposal, they have to be
integrated with each other. This step typically requires adapting algorithms
that were developed for simple flat geometries, such as planes and cubes.
Given the variety of cellular shapes and signaling pathways that can be
combined with each other this adaptation should occur automatically and be
left to the computer. An automatic incorporation requires, however, that the
computational models of cellular morphology and biochemistry are not
‘hardwired’ but can be freely interchanged.

The simulation software Simmune addresses these challenges and can be used to
generate predictive spatially-resolved models [5, 6]. Simmune is designed and
implemented in a highly modular fashion and provides software facilities that
are capable of creating computational models of signaling pathways and
cellular morphology independently from each other. The Simmune Modeler [66]
provides a visual language for molecular states and bimolecular interactions
that allows experimentalists without computational background to create
realistic models of cellular signaling pathways. The difficulty to specify
‘upfront’ all possible reactions within a reaction network is avoided by a
rule-based approach [48, 34, 44, 25, 43, 46, 24]: A user only has to specify
pairwise interactions among molecular binding sites. The Simmune Network
Viewer [15] can then visualize how these interactions are interconnected and
Simmune automatically generates the full resulting network of all possible
multi-molecular complexes. Computational models of cellular geometry can be
created with the Simmune CellDesigner. This software tool provides a graphical
user interface (GUI) enabling a user to populate the simulated system with
spatially discretized versions of cells, their intra-cellular organelles and
of the extracellular space. It also permits to assign certain regions to act
as the locations of the initial biochemistry.

Finally, the Simmune Simulator automatically integrates the representations
created by the Modeler and CellDesigner and computes the spatially-resolved
dynamics as time series and concentration profiles of the system [5, 6].

Currently, Simmune models the time evolution of the reaction network as a
deterministic process and molecular populations as concentrations. As a
consequence, it provides a coarse-grained description in terms of averaged
mean-field quantities and ignores stochastic fluctuations. For many cellular
signaling systems, e.g. those that involve large molecule copy numbers, a
deterministic description is justified and even warranted, because an
alternative stochastic description would be prohibitively inefficient without
contributing additional insights into the system’s behavior.

However, at a molecular scale, all processes underlying cell physiology are
fundamentally stochastic phenomena. Biochemical reactions appear as
probabilistic transformations of or interactions between stochastically moving
discrete molecules. The timing of individual reactions is nondeterministic and
the total number of molecules in a given state evolves in a discrete and
stochastic, or ’noisy’, fashion. As soon as the involved copy numbers are low
and one is interested in phenomena for which the discrete nature of molecules
and their location is essential, such as receptor clustering, this stochastic
aspect of biological reality can no longer be neglected and a deterministic
description fails.

In this manuscript, we present an extension of Simmune that removes its
limitation to deterministic simulations. In contrast to spatially discretized
grid-based approaches, we employ a computational representation of cellular
geometry that is smooth and grid-free. These features make it possible to
incorporate a fully stochastic, spatially resolved description of the cellular
biochemistry. The extension adheres to Simmune’s modular design: the new
computational representation is compatible with Simmune’s current approach,
thereby a user needs to create a model of cellular morphology only once to be
able to explore its behavior through deterministic or stochastic simulations.

Moreover, the computational representations of the biochemistry created with
the Simmune Modeler can be used for the stochastic simulations as well. Thus,
the software now offers a flexible combination of computational tools to study
cellular behavior: The ’building blocks’ of a computational model, such as the
specifications of molecular interactions, representations of cellular geometry
and the type of dynamics (deterministic or stochastic) can be freely combined.
We expect that this flexibility will facilitate simulation studies comparing
differences between deterministic and stochastic dynamics that may emerge in
spite of shared biochemical rules and cellular morphology. Finally, we
anticipate that the close correspondence between the different computational
representations of cellular shapes will pave the way to develop hybrid
simulators capable of integrating deterministic with stochastic time
evolution.

We point out that there are other software tools available that offer
spatially-resolved and/or particle-based stochastic simulation possibilities,
such as, for example [16, 25, 7, 59, 4, 31, 55, 32, 47, 11].

Figure 1:  Schematic illustration of the technical difficulties to integrate a
discretized grid-like model of cellular morphology as used by Simmune with a
stochastic dynamics that generates molecule trajectories. These difficulties
become more severe for surface diffusion, because the stair-like structure
would require major modifications of the stochastic dynamics to adequately
define Brownian paths around corners and edges.

##  2 Computational models of cellular morphology

In Simmune, the geometry of a simulated system consists of the shapes of the
simulated cell(s) with their intracellular domains and organelles and the
structure of the extracellular space [5]. In a deterministic approach, a
system of coupled partial differential equations, known as reaction-diffusion
equations, governs the time evolution of the spatially resolved biochemistry.
Molecular populations are described in terms of local concentrations that are
coarse-grained, averaged quantities. Mathematically, concentrations are
continuous function of time and space–mathematical idealizations that can only
be represented at a finite number of spatial points by a computer. Hence, the
numerical solution of the reaction-diffusion equations requires a spatial
discretization of the system’s geometry. Simmune employs a finite volume
scheme [23] that amounts to a uniform division of the simulated space into
cubic volume elements, also referred to as voxels. As mentioned in the
Introduction, the Simmune CellDesigner and its graphical user interface (GUI)
allows users to generate this discretized space and populate it with
representations of cells, their membranes, cytosolic compartments and
intracellular organelles that are all composed of subsets of voxels.

While the reaction-diffusion equations can be adapted to a grid representation
of space in a natural way for cytosolic and extracellular diffusion, this is
not longer true for diffusion on the membrane. As discussed in Refs. [49, 5],
the surface elements of the voxels that enclose a simulated cell inadequately
approximate the cell’s surface. For instance, the circumference and area of a
discretized sphere does not converge to the correct values, even when the grid
constant tends to zero. Furthermore, a cubic discretization of the surface
Laplacian results in severe artifacts that have to be corrected for by
calculating the local membrane curvatures to adapt area, shape and distances
between surface elements of the cubic voxels. With these corrections being
taken into account, surface diffusion can be accurately described within the
finite volume scheme [49, 5].

Whereas cubic discretization artifacts can be controlled for the deterministic
case using straightforward numerical corrections, this is no longer the case
for stochastic particle-based dynamics. This is because the mathematical
formulation of stochastic dynamics [61] is quite different from the
corresponding deterministic formulation. The dynamic state of the biochemical
system at a point in time is completely specified by the finite set of
locations of the individual molecular complexes, rather than by local
concentrations. The corresponding equations of motion are stochastic
differential equations; their solutions are continuous trajectories (Brownian
paths) that describe the stochastic motion of individual molecular complexes
in time and space. It is difficult to see how a description based on
trajectories could be made compatible with a cubic discretization of space.
This becomes especially apparent for the case of surface diffusion, which
would require to profoundly adapt the stochastic differential equations in
such a way that Brownian paths around corners and edges of cubes can be
defined. These technical difficulties together with the above discussed
artifacts introduced by the surface discretization in terms of cubic voxels
motivate to look for alternative models of geometry.

Clearly, a smooth grid-free model of the system’s geometry is more suitable to
support a stochastic dynamics that generates the paths of individual molecular
complexes. Furthermore, models based on an explicit parametrization may be
excluded, because it would be difficult to construct appropriate parameter
maps efficiently for general surfaces. An alternative is offered by implicit
surfaces that refer to a widely used class of models, such as blobs, also
known as metaballs or soft objects [10, 28, 17], for representing smooth
objects. Instead of using an explicit parametrization for the object’s
surface, implicit surfaces specify certain constraints. For instance, a
metaball specifies a scalar function in three dimensional (3d) space

$\displaystyle\phi$ $\displaystyle:$
$\displaystyle\mathbb{R}^{3}\rightarrow\mathbb{R},$ (1)
$\displaystyle\mathbf{r}=(x,y,z)\rightarrow\phi(\mathbf{r}),$

also known as potential, influence field or density function. In addition, a
metaball defines a threshold value $s$ that serves to define the surface
$\mathcal{S}$ as isocontour

$\mathcal{S}=\\{\mathbf{r}\in\mathbb{R}^{3}:\phi(\mathbf{r})=s\\}.$ (2)

A major advantage of metaballs is that more complex geometries can readily be
generated by a process called blobby modeling [10] that refers to the smooth
blending of several metaballs. Here, the basic idea is to consider the
superposition of the individual influence fields and to define the surface of
the resulting object again as the corresponding isocontour

$\mathcal{S}=\\{\mathbf{r}\in\mathbb{R}^{3}:\sum^{\text{\\#
metaballs}}_{i}\phi_{i}(\mathbf{r})=s\\}.$ (3)

Figure 2:  Schematic visualization of the local tangent plane spanned by the
tangent vectors $\mathbf{t}_{0}(\mathbf{r})$, $\mathbf{t}_{1}(\mathbf{r})$ and
the normal $\mathbf{n}(\mathbf{r})$ at a point $\mathbf{r}$ of the surface of
a smooth grid-free represented object. Exploitation of the locally Euclidean
structure reduces the degree of technical difficulties encountered in the
discrete case, cf. Fig. 1. A grid-free representation in terms of compositions
of metaballs offers the advantage that the unit normal and, in turn, the local
differential geometry such as the tangent plane can readily be determined.

While implicit surfaces can be used directly as computational models, they
have also been employed as auxiliary data structures to facilitate the
creation of discrete models of geometry, such as surface triangulations [17,
28]. In fact, Simmune’s generation of a discretized representation of the
systems’ geometry is implemented in terms of metaballs. In view of possible
applications to stochastic simulations, it is now natural to promote these
purely auxiliary data structures to the actual computational representation of
cellular geometry to enable stochastic simulations. With this extension, a
user needs to create a model of the system’s geometry only once and obtains
automatically both a discretized grid and smooth grid-free representation
amenable to deterministic and stochastic simulations. In addition, the close
relation between the two representations facilitates the comparison between
deterministic and stochastic simulation results.

Finally, let us point out the explicit mathematical form of the potential used
by Simmune’s metaball system

$\phi(\mathbf{r})=\left\\{\begin{array}[]{ll}&[1-(r/R)^{2})]^{2}\qquad\mathrm{if}\,r<R\\\
&0\qquad\qquad\qquad\quad\mathrm{otherwise},\end{array}\right.$ (4)

where $r=|\mathbf{r}|$ and $R$ refers to the metaball’s radius.

##  3 Integration of computational models of geometry and stochastic
integrators

Simmune uses a particle-based stochastic simulation algorithm [52, 53] that is
based on the Smoluchowski-Collins-Kimball (SCK) model [63, 29, 54] of
diffusion-influenced bimolecular reactions.

Different theoretical approaches exist to implement the SCK model’s account
for the diffusive component of bimolecular reactions [61, 27]. The Langevin
approach, also known as Brownian Dynamics (BD) in a numerical context,
accounts for the time evolution of the molecules’ positions directly as
specific solutions (trajectories) of stochastic differential equations. In
contrast, the Smoluchowski and Einstein diffusion equations govern the time
evolution of the conditional probability density function (PDF)
$p(\mathbf{r},t|\mathbf{r}_{0})$, also known as Green’s function (GF) that
provides the probability of finding a molecule at a position $\mathbf{r}$ at
time $t$, given it was sharply localized at $\mathbf{r}_{0}$ at time
$t_{0}=0$.

The SCK framework is widely used in stochastic simulation algorithms to
quantify the impact of diffusion on chemical reaction rates [20, 21, 42, 8,
62, 50, 39, 56]. Frequently, these type of algorithms employ GF to address the
problem of the high computational cost that single-particle stochastic models
based solely on BD typically entail.

The integration of computational representations of the cellular morphology
and the stochastic integrator that we will describe in the following does not
depend on the details of the stochastic approach. This ’loose coupling’ offers
the advantage that different stochastic integrators, even those that do not
use the SCK picture, such as Doi-based [18, 19, 64, 2, 53] approaches, could
be readily swapped with each other.

Typically, particle-based stochastic simulation algorithms assume that flat
Euclidean spaces, such as planes and cubes, describe the simulated system’s
geometry and tacitly employ special geometric features only provided by
Euclidean structures, such as the existence of a global coordinate system and
distance metric. Therefore, the transition to more general geometries
motivated by the wish for realistic models of cellular morphology necessitates
modifications to the stochastic description. The parts of the algorithms
describing reactions that involve a spatial aspect, such as bimolecular
associations and unimolecular decays (through the need to place the decay
products correctly) have to be adapted, while descriptions of ’non-spatial’
reactions of molecular complexes such as phosphorylations remain unaffected.
Among the different types of bimolecular reactions, such as associations
between two cytosolic complexes and between cytosolic and membrane complexes,
the most profound adaptions are necessary for associations between membrane
complexes.

###  3.1 Stochastic description of reaction-diffusion systems on curved
membranes

Already in the Euclidean case, reaction-diffusion systems in 2d are quite
different from corresponding 3d systems, due to distinct features that
diffusion exhibits in 2d. Two is the critical dimension with regard to
recurrence and transience of Brownian walks [60], the steady-state solution of
the diffusion equation is inconsistent with the boundary condition at infinity
[22], and the time-dependent reaction rate vanishes in the long-time limit.
These properties are related to each other and make it more complicated to
derive a mass-action rate from the underlying microscopic stochastic
description [65]. In addition, the mathematical treatment appears to be more
involved than in 1d and 3d [51, 30]. Finally, the 2d case is of particular
importance in cell biological applications, being a prerequisite for a
quantitative understanding of processes such as signal-induced inhomogeneities
and receptor clustering on cell membranes [9], where the receptors’ slow
lateral diffusion represents a less efficient mixing mechanism than in the
three dimensional cytoplasm.

The task at hand is now to incorporate the stochastic algorithms originally
developed for Euclidean planes with the computational representations of
general cellular shapes, thereby including curvature as a quantity that takes
part in the stochastic dynamics.

In mathematics, the corresponding problem is solved by exploiting the locally
flat structure of Riemannian manifolds: at any point $\mathbf{r}_{0}$ of its
surface a Riemann manifold gives rise to a local tangent plane, also known as
frame, that comes equipped with a locally Euclidean structure [38]. Equations
that govern the dynamics in a globally flat space keep their validity in the
local frame in a vicinity of $\mathbf{r}_{0}$. The size of surface patch
around $\mathbf{r}_{0}$ that is well approximated by the local frame is
determined by the local curvature.

Following this mathematical lead, we now seek as a first step an algorithmic
description of free diffusion, i.e., Brownian motion on the surface of a
metaball or of a smooth blend of several of them.

As discussed in Sec. 2, this task is equivalent to find a numerical approach
to free Brownian motion on a surface defined by the general equation

$\phi(\mathbf{r})=0,$ (5)

see Eq. 3. This problem has been solved previously [35] in the spirit of the
mathematical constructions described above. Now, while Ref. [35] considers
free diffusion on a general surface, we use the same framework to describe
2d-2d interactions, such as (irr)reversible associations of membrane
molecules, and 2d-3d interactions, such as binding reactions between membrane
and extracellular or cytosolic complexes, for which the surface and bulk are
described by general shapes that can be represented by metaballs.

Figure 3:  Discrete and smooth representations of cellular morphology.
Top A discretized finite volume representation suitable for deterministic
simulations created via the Simmune CellDesigner. Colored patches specify
regions that allow to define the spatial distribution of the initial
biochemistry.
Bottom The automatically generated corresponding grid-free, smooth
representation of the same cellular morphology amenable to stochastic
simulations. For illustration, the grid is also displayed.

Following Ref. [35] and defining the surface as the points that satisfy Eq. 5,
one may characterize the local tangent plane at a point
$\mathbf{r}_{0}=(x_{0},y_{0},z_{0})$ as the set of points $\mathbf{r}=(x,y,z)$
that satisfy

$(\mathbf{r}-\mathbf{r}_{0})\cdot\nabla\phi(\mathbf{r}_{0})=0,$ (6)

that is, the vector $\mathbf{r}-\mathbf{r}_{0}$ is perpendicular to the
gradient $\nabla\phi(\mathbf{r}_{0})$, which is the surface’s normal at point
$\mathbf{r}_{0}$.

Now, we seek to construct an orthonormal basis of local tangent plane, that
is, two mutually perpendicular vectors $\mathbf{t}_{0}$, $\mathbf{t}_{1}$ that
span the tangent space at $\mathbf{r}_{0}$.

Knowledge of this basis enables one to define a free Brownian position update
analogously to the Euclidean plane case

$\mathbf{r}^{\ast}(t+\Delta t)=\mathbf{r}(t)+\sqrt{2D\Delta
t}\xi\mathbf{t}_{0}+\sqrt{2D\Delta t}\chi\mathbf{t}_{1}$ (7)

where $\xi$ and $\chi$ refer to samples from a Gaussian distribution with
vanishing mean and variance equal to unity. However, the updated position
given by Eq. 7 may not fulfill condition Eq. (5). Therefore,
$\mathbf{r}^{\ast}$ has to be projected back onto the surface as follows [35]

$\mathbf{r}(t+\Delta t)=\mathbf{r}^{\ast}(t+\Delta
t)-\frac{\phi(\mathbf{r}^{\ast}(t+\Delta
t))\nabla\phi(\mathbf{r}^{\ast}(t+\Delta
t))}{|\nabla\phi(\mathbf{r}^{\ast}(t+\Delta t))|^{2}}.$ (8)

Here, $\mathbf{r}(t+\Delta t)$ is the final location on the surface after one
completed simulation step.

It has been shown [35] that the algorithm represented by Eqs. (7), (8) is
stable and satisfies the detailed balance condition, provided that
$|\mathbf{r}(t+\Delta t)-\mathbf{r}^{\ast}(t+\Delta t)|$ is much smaller than
the local radius of curvature.

Therefore, the local curvature enforces an additional upper limit to the size
of the time step $\Delta t$; already in the Euclidean case there are upper
limits determined by the requirement that every molecule encounters maximal
one other molecule per time step on average [39] and by the unimolecular
reactions’ time scales.

A major advantage of the metaball representation is that it permits an
efficient and straightforward implementation of the algorithm just described.
The defining properties of metaballs are critical for a swift construction of
the local differential geometry, that is, the local tangent plane and the
curvature. Because of Eqs. (3), (5), one can immediately obtain the unit
normal vector at any point of the metaball’s surface

$\mathbf{n}(\mathbf{r}):=\nabla\phi(\mathbf{r})/|\nabla\phi(\mathbf{r})|$ (9)

Next, one can exploit that knowledge of $\mathbf{n}(\mathbf{r})$ is sufficient
to obtain the tangent plane, see [37, 26] for algorithms that construct a
complete orthonormal basis, given only a single vector. Also the Gaussian
curvature can be calculated from the unit normal [57, 35]

$K(\mathbf{r})=\frac{1}{R_{1}R_{2}}=\frac{1}{2}[-(\partial_{i}\phi_{j}(\mathbf{r}))^{2}+(\nabla\cdot\mathbf{n}(\mathbf{r}))^{2}].$
(10)

Here, $R_{1}$ and $R_{2}$ denote two principal radii of the curvature. In the
light of the central importance of the gradient, the simple mathematical form
of the potential Eq. (4) offers an additional advantage: even more complex
surfaces composed of several metaballs are represented by a linear
superposition of these potentials and hence the resulting gradient itself can
be calculated efficiently without the need to resort to numerical
differentiation.

In this way, the construction of a local coordinate system opens up the
possibility to reuse Euclidean geometry stochastic algorithms that not only
can handle free diffusion but also bimolecular reactions. All subsequent steps
of these algorithms, such as detection of encounters/reactions between two
molecules and correct two-body propagation of a non-reacting molecule pair,
can be extended to models involving general cellular geometries.

###  3.2 Associations between surface and bulk complexes

As mentioned, the ability to swiftly construct the local differential geometry
at any point of a metaball permits to export also techniques that deal with
2d-3d membrane-bulk binding reactions, such as associations between membrane
and cytosolic molecules, to general cellular morphologies specified by Eq.
(5).

Approaches to stochastic algorithms based on Green’s function assisted
Brownian Dynamics involve the key step to rewrite the two-body diffusion
equation describing an isolated pair of reacting molecules into two diffusion
equations [62, 39, 56]. The resulting equations describe free diffusion of the
center-of-diffusion vector and diffusion of the inter-molecule vector subject
to a boundary condition that implements the reaction.

This alternative form of the two-body diffusion equation is enabled by a
transformation to a center-of-diffusion coordinate system. While this
transformation is straightforward for a pair of two bulk molecules, it becomes
more complicated for a pair of one molecule diffusing on a plane and the other
moving in the bulk.

In this scenario, the transformation that separates the two-body diffusion
equation reads [56]

$\mathbf{r}=\begin{pmatrix}x_{3d}-x_{2d}\\\ y_{3d}-y_{2d}\\\
\epsilon(z_{3d}-z_{2d})\end{pmatrix},\quad\epsilon=\sqrt{1+\frac{D_{2d}}{D_{3d}}}$
(11)

and

$\mathbf{R}=\frac{1}{D_{2d}+D_{3d}}\begin{pmatrix}D_{2d}x_{3d}+D_{3d}x_{2d}\\\
D_{2d}y_{3d}+D_{3d}y_{2d}\\\ z_{2d}\end{pmatrix},$ (12)

provided the 2d plane is described by $z=\text{const}$. Here,
$\mathbf{r}_{3d}=(x_{3d},y_{3d},z_{3d})$, $D_{3d}$ and
$\mathbf{r}_{2d}=(x_{2d},y_{2d},z_{2d})$, $D_{2d}$ refer to the position and
diffusion constant describing the bulk and surface molecule, respectively.
Now, using the construction of the local differential geometry as described in
Sec. 3.1, the center-of-diffusion transformation and all subsequent steps of
the stochastic algorithm can be applied with respect to the local frame. In
this way, 2d-3d binding reactions can be extended to general models of
cellular morphology.

##  4 Numerical Tests

We performed several tests of the presented simulation algorithm.

First, we simulated free diffusion of a single molecule and the diffusion-
controlled irreversible reaction of an isolated pair on the surface of a
sphere. The aim was to numerically construct the corresponding probability
density function, i.e., the Green’s function $p(\theta,t|\theta_{0})$ [30],
where $\theta$ denotes the polar angle of standard spherical coordinates. The
simulation set up was as follows. For free diffusion, the molecule was
initially placed at the North Pole of the sphere, that is, $\theta_{0}=0$. The
simulation was run for a time $t_{\text{sim}}=1\text{s}$ during which molecule
underwent a diffusive motion and the simulation time step is $\Delta
t=0.001\text{s}$. The final position at time $t_{\text{sim}}$ was recorded and
a new simulation was initiated. Fig. 4 shows the simulation results for two
different sets of the diffusion constant $D$ and of the sphere radius $R$. The
used sets are $D=1\mu\text{m}^{2}\text{s}^{-1}$, $R=1\mu\text{m}$ and
$D=1.25\mu\text{m}^{2}\text{s}^{-1}$, $R=1.5\mu\text{m}$. We found excellent
agreement between the simulation results and the exact analytical expression
of the GF describing free diffusion on a sphere (Eq. 13).

In addition, we considered an isolated pair of molecules $A$ and $B$. Molecule
$A$ was held fixed at the sphere’s South Pole $\theta=\pi$, while molecule $B$
was placed at several initial positions. The simulation was run for a time
$t_{\text{sim}}=2\text{s}$, during which molecule $B$ diffused and potentially
associated with molecule $A$ that was modeled as a perfect sink. The radius of
the sphere, the effective encounter radius [3, 30] and the diffusion constant
are $R=1\mu\text{m}$, $a_{\text{eff}}=0.1\mu\text{m}$ and
$D=1\mu\text{m}^{2}\text{s}^{-1}$, respectively. After each run, we recorded
$B$’s final position, unless it was bound to $A$. The resulting histogram was
normalized to account for the number of bound states at $t_{\text{sim}}$. Fig.
5 shows the results for various initial positions together with the
corresponding graphs of the exact analytical representation (Eq. (14), see
Refs. [14, 30]). Again, we found excellent agreement.

Figure 4:  Numerical construction of the GF describing the free diffusive
motion of a molecule on the surface of a sphere. The figure shows the $\theta$
dependence of the GF for two sets of the diffusion constant and sphere radius
$R$, as indicated, at time $t_{\text{sim}}=1\text{s}$. The solid lines
correspond to the exact analytical expression (Eq. (13)). The various markers
refer to the stochastic simulation results.  Figure 5:  Numerical construction
of the GF describing the diffusion-controlled irreversible binding of an
isolated pair on the surface of a sphere. The figure shows the $\theta$
dependence of the GF satisfying a completely absorbing boundary condition for
different values of the initial position $\theta_{0}$, as indicated, at time
$t_{\text{sim}}=2\text{s}$. The solid lines correspond to the exact analytical
expressions (Eq. (14)). The various markers refer to simulation results.

Next, we turned to many-molecule systems and simulated the reversible binding
reaction $A+B\leftrightarrow C$ on a surface of a sphere with surface area
$A=0.994\mu\text{m}^{2}$. Initially, $N_{A}(t=0)=N_{B}(t=0)=994$ copies of $A$
and $B$ were uniformly distributed over the surface. The simulation was
performed with $D_{\text{eff}}=D_{A}+D_{B}=2\mu\text{m}^{2}\text{s}^{-1}$ and
$a_{\text{eff}}=a_{A}+a_{B}=10^{-3}\mu\text{m}$. Using the dimensionless time
step $\Delta\tau=\Delta tD_{\text{eff}}/a^{2}_{\text{eff}}=0.2$., the
simulation was executed until $t_{\text{sim}}=1\text{s}$ and the copy numbers
$N_{A}(t)$, $N_{B}(t)$, $N_{C}(t)$ over time $t$ were recorded. The further
parameters are $\kappa_{a}=1\mu\text{m}^{2}\text{s}^{-1}$ and
$\kappa_{d}=1\text{s}^{-1}$, where $\kappa_{a}$ and $\kappa_{d}$ refer to the
intrinsic association and dissociation constant, respectively, which
characterize the so-called backreaction boundary condition (BC) [3]. The
backreaction BC generalizes the classic Smoluchowski BC to account for
diffusion-influenced reversible reactions [54, 3]. Fig. 6 shows the simulation
result and the numerical solution of the modified rate equation (MRE) approach
[58, 65].

Finally, we considered reversible associations between membrane and cytosolic
molecules. To this end, we simulated the reversible binding reaction
$A+B\leftrightarrow C$, initially uniformly placing $N_{A}(t=0)=1007$,
$N_{B}(t=0)=991$ copies inside a sphere with radius $R=0.249\mu\text{m}$ and
on its surface, respectively. The further parameters are $\kappa_{a}=5\cdot
10^{-4}\mu\text{m}^{3}\text{s}^{-1}$, $\kappa_{d}=1\text{s}^{-1}$ and
$a_{\text{eff}}=10^{-3}\mu\text{m}$. The simulation was run for a time
$t_{\text{sim}}=10\text{s}$ and the copy numbers $N_{A}(t)$, $N_{B}(t)$,
$N_{C}(t)$ over time $t$ were recorded. The used dimensionless time step is
$\Delta\tau=\Delta tD_{\text{eff}}/a^{2}_{\text{eff}}=10.1$. Fig. 7 compares
the simulation result with the numerical solution of the standard mass-action
equations.

Figure 6:  Simulation results vs. MRE solution for the reversible
$A+B\leftrightarrow C$ reaction taking place on the surface of a sphere.
Simulation results and numerical MRE solution are shown in dashed red and
solid black curves, respectively, for
$D_{A}=D_{B}=1\mu\text{m}^{2}\text{s}^{-1}$ and
$\kappa_{a}=1\mu\text{m}^{2}\text{s}^{-1}$ and $\kappa_{d}=1\text{s}^{-1}$.
Figure 7:  Simulation results vs. numerical solution of mass-action kinetics
for the reversible $A+B\leftrightarrow C$ reaction taking place between $A$
and $B$ molecules, diffusing inside the sphere and on its surface,
respectively. Simulation results and numerical solutions are shown in dashed
red and solid black for $D_{A}/10=D_{B}=1\mu\text{m}^{2}s^{-1}$ and
$\kappa_{a}=5\cdot 10^{-4}\mu\text{m}^{3}s^{-1}$ and $\kappa_{d}=1s^{-1}$.

##  5 Appendix

The Green’s function for free diffusion on a sphere with radius $R$ is given
by [12]

$p_{\text{free}}(\theta,t|\theta_{0}=0)=\frac{1}{4\pi}\sum_{l=0}^{\infty}(2l+1)\exp\bigg{(}-\frac{Dtl(l+1)}{R^{2}}\bigg{)}P_{l}(\cos\theta),$
(13)

where $P_{l}$ refers to the Legendre polynomials [1].

Quite recently, an exact expression has been obtained [30] that gives the
survival probability of an isolated pair undergoing a diffusion-influenced
reversible reaction on the surface of a sphere. Here, we use the GF rather
than the survival probability; as pointed out in Ref. [30], the corresponding
expression can be readily inferred from the results presented in [30]. For the
limiting case of a perfect sink one has [14, 30]

$p_{\text{abs}}(\theta,t|\theta_{0})=\frac{1}{2\pi}\sum_{n=0}^{\infty}b^{2}_{n}\exp\bigg{(}-\frac{Dt\nu_{n}(\nu_{n}+1)}{R^{2}}\bigg{)}P_{\nu_{n}}(\cos\theta)P_{\nu_{n}}(\cos\theta_{0})$
(14)

and $p_{\text{abs}}(\theta,t|\theta_{0})$ satisfies the completely absorbing
BC

$p_{\text{abs}}(\theta_{a},t|\theta_{0})=0,$ (15)

where $\theta_{a}$ is given by

$\cos\theta_{a}=-\sqrt{1-(a/R)^{2}}.$ (16)

Here, $a,R$ denote the encounter radius of the isolated pair and the radius of
the sphere, respectively. Note that the Legendre functions of fractional order
are given by Gauss’ hypergeometric function ${}_{2}F(a,b;c;x)$, more precisely
one has [1]

$P_{\nu}(x):=\,_{2}F(-\nu,\nu+1;1;(1-x)/2)).$ (17)

The $\nu_{n}$ are the solutions of the equation [30]

$P_{\nu}(\cos(\theta_{a}))=0.$ (18)

Then, the $b_{n}$ terms can be calculated according to [30]

$b_{n}=\bigg{(}\int^{1}_{a}[P_{\nu_{n}}(x)]^{2}\,dx\bigg{)}^{-1/2}.$ (19)

To apply the expression given in Eq. (14), we numerically calculated the first
eight roots $\nu_{n}$ and the corresponding coefficients $b_{n}$.

## Acknowledgments

This research was supported by the Intramural Research Program of the NIH,
NIAID.

## References

  * [1] M. Abramowitz and I. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, Dover, New York, 1965.
  * [2] I. Agbanusi and S. Isaacson, Bull. Math. Biol., 76 (2014), p. 922.
  * [3] N. Agmon and A. Szabo, Theory of reversible diffusion-influenced reactions, The Journal of Chemical Physics, 92 (1990), pp. 5270–5284.
  * [4] S. S. Andrews, N. J. Addy, R. Brent, and A. P. Arkin, Detailed simulations of cell biology with smoldyn 2.1, PLoS computational biology, 6 (2010).
  * [5] B. R. Angermann, F. Klauschen, A. D. Garcia, T. Prustel, F. Zhang, R. N. Germain, and M. Meier-Schellersheim, Computational modeling of cellular signaling processes embedded into dynamic spatial contexts, Nat Methods, 9 (2012), pp. 283–9.
  * [6] B. R. Angermann and M. Meier-Schellersheim, Using python for spatially resolved modeling with simmune, Methods Mol Biol, 1945 (2019), pp. 161–177.
  * [7] S. Arjunan and M. Tomita, A new multicompartmental reaction-diffusion modeling method links transient membrane attachment of e. coli mine to e-ring formation, Nature Precedings, (2009), pp. 1–1.
  * [8] T. M. Barenbrug, E. Peters, and J. D. Schieber, Accurate method for the brownian dynamics simulation of spherical particles with hard-body interactions, The Journal of chemical physics, 117 (2002), pp. 9202–9214.
  * [9] I. Bethani, S. Skanland, I. Dikic, and A. Acker-Palmer, EMBO J., 29 (2010), p. 2677.
  * [10] J. F. Blinn, A generalization of algebraic surface drawing, ACM transactions on graphics (TOG), 1 (1982), pp. 235–256.
  * [11] M. L. Blinov, J. C. Schaff, D. Vasilescu, I. I. Moraru, J. E. Bloom, and L. M. Loew, Compartmental and spatial rule-based modeling with virtual cell, Biophysical journal, 113 (2017), pp. 1365–1372.
  * [12] D. R. Brillinger, A particle migrating randomly on a sphere, Journal of Theoretical Probability, 10 (1997), pp. 73–87.
  * [13] M. D. Brown and D. B. Sacks, Protein scaffolds in map kinase signalling, Cell Signal, 21 (2009), pp. 462–9.
  * [14] N. Chao, S. H. Young, and M. Poo, Localization of cell membrane components by surface diffusion into a” trap”, Biophysical journal, 36 (1981), pp. 139–153.
  * [15] H.-C. Cheng, B. R. Angermann, F. Zhang, and M. Meier-Schellersheim, Networkviewer: visualizing biochemical reaction networks with embedded rendering of molecular interaction rules, BMC Syst Biol, 8 (2014), p. 70.
  * [16] J. Czech, M. Dittrich, and J. R. Stiles, Rapid creation, monte carlo simulation, and visualization of realistic 3d cell models, in Systems Biology, Springer, 2009, pp. 237–287.
  * [17] B. R. De Araújo, D. S. Lopes, P. Jepp, J. A. Jorge, and B. Wyvill, A survey on implicit surface polygonization, ACM Computing Surveys (CSUR), 47 (2015), pp. 1–39.
  * [18] M. Doi, J. Phys. A: Math. Gen., 9 (1976), p. 1465.
  * [19] , J. Phys. A: Math. Gen., 9 (1976), p. 1479.
  * [20] A. L. Edelstein and N. Agmon, Brownian dynamics simulations of reversible reactions in one dimension, The Journal of chemical physics, 99 (1993), pp. 5396–5404.
  * [21] , Brownian simulation of many-particle binding to a reversible receptor array, Journal of Computational Physics, 132 (1997), pp. 260–275.
  * [22] C. Emeis and P. Fehder, J. Am. Chem. Soc., 92 (1970), p. 2246.
  * [23] R. Eymard, T. Gallouët, and R. Herbin, Finite volume methods, Handbook of numerical analysis, 7 (2000), pp. 713–1018.
  * [24] J. R. Faeder, M. L. Blinov, and W. S. Hlavacek, Rule-based modeling of biochemical systems with bionetgen, Methods Mol Biol, 500 (2009), pp. 113–67.
  * [25] J. Feret, V. Danos, J. Krivine, R. Harmer, and W. Fontana, Internal coarse-graining of molecular systems, Proc Natl Acad Sci U S A, 106 (2009), pp. 6453–8.
  * [26] J. R. Frisvad, Building an orthonormal basis from a 3d unit vector without normalization, Journal of Graphics Tools, 16 (2012), pp. 151–159.
  * [27] C. Gardiner, Stochastic methods, vol. 4, Springer Berlin, 2009.
  * [28] A. Gomes, I. Voiculescu, J. Jorge, B. Wyvill, and C. Galbraith, Implicit curves and surfaces: mathematics, data structures and algorithms, Springer Science & Business Media, 2009.
  * [29] U. Gösele, Prog. React. Kinet., 13 (1984), p. 63.
  * [30] D. S. Grebenkov, Reversible reactions controlled by surface diffusion on a sphere, J Chem Phys, 151 (2019), p. 154103.
  * [31] G. Gruenert, B. Ibrahim, T. Lenser, M. Lohel, T. Hinze, and P. Dittrich, Rule-based spatial modeling with diffusing, geometrically constrained molecules, BMC bioinformatics, 11 (2010), p. 307.
  * [32] I. Hepburn, W. Chen, S. Wils, and E. De Schutter, Steps: efficient simulation of stochastic reaction–diffusion models in realistic morphologies, BMC systems biology, 6 (2012), p. 36.
  * [33] W. S. Hlavacek, J. R. Faeder, M. L. Blinov, A. S. Perelson, and B. Goldstein, The complexity of complexes in signal transduction, Biotechnol Bioeng, 84 (2003), pp. 783–94.
  * [34] W. S. Hlavacek, J. R. Faeder, M. L. Blinov, R. G. Posner, M. Hucka, and W. Fontana, Rules for modeling signal-transduction systems, Sci STKE, 2006 (2006), p. re6.
  * [35] R. Hołyst, D. Plewczyński, A. Aksimentiev, and K. Burdzy, Diffusion on curved, periodic surfaces, Phys Rev E Stat Phys Plasmas Fluids Relat Interdiscip Topics, 60 (1999), pp. 302–7.
  * [36] S. Hong, R. B. Troyanovsky, and S. M. Troyanovsky, Spontaneous assembly and active disassembly balance adherens junction homeostasis, Proc Natl Acad Sci U S A, 107 (2010), pp. 3528–33.
  * [37] J. F. Hughes and T. Möller, Building an orthonormal basis from a unit vector, Journal of Graphics Tools, 4 (1999), pp. 33–35.
  * [38] N. Ikeda and S. Watanabe, Stochastic differential equations and diffusion processes, Kodansha scientific books, North-Holland, 1989.
  * [39] M. E. Johnson and G. Hummer, Free-propagator reweighting integrator for single-particle dynamics in reaction-diffusion models of heterogeneous protein-protein interaction systems, Phys Rev X, 4 (2014).
  * [40] R. B. Jones, A. Gordus, J. A. Krall, and G. MacBeath, A quantitative protein interaction network for the erbb receptors using protein microarrays, Nature, 439 (2006), pp. 168–74.
  * [41] B. N. Kholodenko, Four-dimensional organization of protein kinase signaling cascades: the roles of diffusion, endocytosis and molecular motors, J Exp Biol, 206 (2003), pp. 2073–82.
  * [42] H. Kim, M. Yang, and K. J. Shin, Dynamic correlation effect in reversible diffusion-influenced reactions: Brownian dynamics simulation in three dimensions, The Journal of chemical physics, 111 (1999), pp. 1068–1075.
  * [43] M. Koschorreck and E. D. Gilles, Alc: automated reduction of rule-based models, BMC Syst Biol, 2 (2008), p. 91.
  * [44] L. Lok and R. Brent, Automatic generation of cellular reaction networks with moleculizer 1.0, Nat Biotechnol, 23 (2005), pp. 131–6.
  * [45] C. I. Maeder, M. A. Hink, A. Kinkhabwala, R. Mayr, P. I. H. Bastiaens, and M. Knop, Spatial regulation of fus3 map kinase activity through a reaction-diffusion mechanism in yeast pheromone signalling, Nat Cell Biol, 9 (2007), pp. 1319–26.
  * [46] A. Mallavarapu, M. Thomson, B. Ullian, and J. Gunawardena, Programming with models: modularity and abstraction provide powerful capabilities for systems biology, J R Soc Interface, 6 (2009), pp. 257–70.
  * [47] T. Masaru and P. K. Dhar, E-cell system: basic concepts and applications., Springer, 2013.
  * [48] M. Meier-Schellersheim, X. Xu, B. Angermann, E. J. Kunkel, T. Jin, and R. N. Germain, Key role of local regulation in chemosensing revealed by a new molecular interaction-based modeling method, PLoS Comput Biol, 2 (2006), p. e82.
  * [49] I. L. Novak, F. Gao, Y.-S. Choi, D. Resasco, J. C. Schaff, and B. M. Slepchenko, Diffusion on a curved surface coupled to diffusion in the volume: Application to cell biology, J Comput Phys, 226 (2007), pp. 1271–1290.
  * [50] T. Oppelstrup, V. V. Bulatov, A. Donev, M. H. Kalos, G. H. Gilmer, and B. Sadigh, First-passage kinetic monte carlo method, Physical Review E, 80 (2009), p. 066701.
  * [51] T. Prüstel and M. Meier-Schellersheim, J. Chem. Phys., 137 (2012), p. 054104.
  * [52] T. Prüstel and M. Meier-Schellersheim.  arXiv:1509.06810v1 [q-bio.QM], 2015.
  * [53] T. Prüstel and M. Meier-Schellersheim, Unified path integral approach to theories of diffusion-influenced reactions, Physical Review E, 96 (2017), p. 022151.
  * [54] S. Rice, Diffusion Limited Reactions, Elsevier, New York, 1985.
  * [55] M. W. Sneddon, J. R. Faeder, and T. Emonet, Efficient modeling, simulation and coarse-graining of biological complexity with nfsim, Nat Methods, 8 (2011), pp. 177–83.
  * [56] T. R. Sokolowski, J. Paijmans, L. Bossen, T. Miedema, M. Wehrens, N. B. Becker, K. Kaizu, K. Takahashi, M. Dogterom, and P. R. ten Wolde, egfrd in all dimensions, The Journal of chemical physics, 150 (2019), p. 054108.
  * [57] M. Spivak, A Comprehensive Introduction to Differential Geometry, no. v. 3 in A Comprehensive Introduction to Differential Geometry, Publish or Perish, Incorporated, 1975.
  * [58] A. Szabo, Theoretical approaches to reversible diffusion-influenced reactions: Monomer–excimer kinetics, The Journal of chemical physics, 95 (1991), pp. 2481–2490.
  * [59] D. P. Tolle and N. Le Novère, Meredys, a multi-compartment reaction-diffusion simulator using multistate realistic molecular complexes, BMC systems biology, 4 (2010), p. 24.
  * [60] D. Toussaint and F. Wilczek, J. Chem. Phys., 78 (1983), p. 2642.
  * [61] N. G. Van Kampen, Stochastic processes in physics and chemistry, vol. 1, Elsevier, 1992.
  * [62] J. S. van Zon and P. R. ten Wolde, Simulating biochemical networks at the particle level and in time and space: Green’s function reaction dynamics, Phys Rev Lett, 94 (2005), p. 128103.
  * [63] M. von Smoluchowski, Z. Phys. Chem., 92 (1917), p. 129.
  * [64] G. Wilemski and M. Fixman, J. Chem. Phys., 58 (1973), p. 4009.
  * [65] O. N. Yogurtcu and M. E. Johnson, Theory of bi-molecular association dynamics in 2d for accurate model and experimental parameterization of binding rates, J Chem Phys, 143 (2015), p. 084117.
  * [66] F. Zhang, B. R. Angermann, and M. Meier-Schellersheim, The simmune modeler visual interface for creating signaling networks based on bi-molecular interactions, Bioinformatics, 29 (2013), pp. 1229–30.