# Enzyme kinetics simulation at the scale of individual particles

Taylor Kearney111Monash University, Clayton, Victoria, Australia,
<EMAIL_ADDRESS>Mark B. Flegg222Monash University, Clayton,
Victoria, Australia<EMAIL_ADDRESS>
###### Abstract

Enzyme-catalysed reactions involve two distinct timescales. There is a short
timescale on which enzymes bind to substrate molecules to produce bound
complexes, and a comparatively long timescale on which the complex is
transformed into a product. The rate at which the substrate is converted into
product is characteristically non-linear and is traditionally derived by
applying of singular perturbation theory to the system’s governing equations.
Central to this analysis is the assumption that complex formation is
effectively instantaneous on the timescale over which significant substrate
degradation occurs. This prevents many particle-based simulations of reaction-
diffusion systems from accurately modelling enzyme kinetics since they rely on
proximity based reaction conditions that do not correctly model the fast
reactions associated with the complex on the long timescale. In this paper we
derive a new proximity based reaction condition that correctly incorporates
the reactions that occur on the short timescale for a specific enzymatic
system. We present proof of concept particle-based simulations that
demonstrate that non-linear reaction rates typical of enzyme kinetics can be
reproduced without needing to explicitly simulate reactions on the short
timescale.

Keywords— Enzyme kinetics, diffusion controlled reactions, Smoluchowski
kinetics, particle-based simulation

##  1 Introduction

Whole cell models promise to revolutionize systems biology. The capability to
simulate the integrated function of every gene and molecule in a cell would
assist clinicians in individualizing therapy [18, 32, 22] and enable computer-
aided designs in synthetic biology [33, 29]. Their development would encourage
the unification of our currently disconnected and heterogeneous biological
datasets [6, 24] and facilitate the discovery of emergent phenomena [46]. This
worthy goal has been touted as a ‘grand challenge’ for $21^{\text{st}}$
century systems biology [47] and will require extensive interdisciplinary
collaboration if we are to be successful [25].

Early models attempted to describe the entirety of a cell’s function using a
single mathematical technique; namely ordinary differential equations (ODEs)
[38]. ODE models alone are not sufficient to describe all robust cellular
behaviours. Over time, the disparate spatial and temporal scales of the
involved phenomena inspired the development of hybrid models that are a
conglomerate of many submodels that each target a specific biological module
within the cell [19, 43, 23]. Despite this progress, ODE models still abound
and are typically used to describe the chemical reactions that dictate a
cell’s function [7]. Such models are underpinned by the presumption that the
chemical species involved can be accurately represented as well-mixed,
deterministic, time-varying continuous concentrations. However, in biological
systems molecules can occur in very low numbers and in highly localised
distributions. For example, an entire cell may only contain a single molecule
of mRNA for a particular gene [45]. When molecules become so sparsely
distributed, it is impossible to select a neighborhood about a point that
contains sufficiently many molecules to define a meaningful concentration
[14]. Moreover, cellular functions exhibit inherent stochasticity [28, 35]
owing to their origin in molecular interactions. In this regime,
concentrations must be replaced by a collection of individual molecules
undergoing reaction-diffusion processes.

Ideally, we explicitly model the molecular dynamics (MD) [20] that govern
biochemical reactions, but such an approach results in simulations too
computationally intensive for current computing hardware [16]. Failing this,
we are forced to make simplifying assumptions about the involved physical
processes in the hopes of obtaining a model that represents individual
molecules, but abstains from explicit calculation of the intricate molecular
interactions. Smoluchowski proposed such a model in $1917$ that describes the
interaction of chemicals as diffusive point particles on a continuous domain,
and has become one of the most widely accepted idealized models for reaction-
diffusion systems of molecules [39]. Characteristic of this approach is the
presumption that the system is sparse and that the relevant molecules can be
treated as individuals that undergo isotropic diffusion as a result of their
collisions with implicit solvent molecules. Bimolecular reactions are modelled
by imposing that two molecules undergo a reaction if they become separated by
less than a predefined distance $\sigma$. When a reaction occurs, the reactant
molecules are removed from the system and replaced with a single molecule of
the product of the reaction. We note for the sake of completeness that in
Smoluchowski based frameworks unimolecular reactions are assumed to occur
instantaneously and can be modelled as Poissonian processes that are
independent of molecular diffusion [4, 40]. Higher order reactions can also be
incorporated by way of an extension developed by Flegg [17] that we will
review briefly in subsection 2.1.

Smoluchowski’s original model has been criticized for neglecting several
important physical mechanisms that can influence reaction rates. These
critiques have inspired the development of many derivative models that attempt
to account for additional mechanisms including: activation energies [8],
intermolecular forces [10, 26] and hydrodynamic effects [21, 51]. Despite this
apparent diversity and in some cases the introduction of additional reaction
parameters - see for instance models by Collins and Kimball [8], or Doi [11,
14] \- all current derivatives of Smoluchowski’s theory still describe
bimolecular reactions as a proximity-dependent interaction between diffusing
point particles. In each case, the definition of $\sigma$ can be altered to
include extra information relating to the additional physical mechanisms
considered, but its role as a parameter that summarizes the molecular
interaction remains unchanged [17]. This underlying commonality is a testament
to the robustness of Smoluchowski’s original theory, and it serves as the
foundation for many prominent software packages for particle based simulation
of reaction-diffusion systems, including: MCell [15, 42], Smoldyn [4], Green’s
function reaction dynamics (GFRD) [49, 48], ReaDDy [36] and enhanced Green’s
function reaction dynamics (eGFRD) [44, 41]. All of these packages are capable
of accurately simulating elementary unimolecular and bimolecular reactions in
isolation. Thus, it seems reasonable to conclude that the same software can
accurately simulate reaction networks composed of multiple elementary
reactions, but for many biochemical networks this approach overlooks a
fundamental assumption of Smoluchowski’s theory.

Let us apply Smolcuhowski’s model to a reaction between two chemical species
$A$ and $B$ which produces a product $C$. We assume that initially the
molecules of $A$ and $B$ are distributed uniformly at random within a volume
$V$. The central result of Smoluchowski’s theory states that two molecules
(originally modelled as hard spheres) diffusing in a sufficiently large volume
$V$, will - after an initial transient ($t_{s}$) - come into contact (at
distance $\sigma$) at a constant rate per unit time $K$ given by,

$K=\frac{k}{V}=\frac{4\pi\sigma\hat{D}_{2}}{V}.$ (1)

Here $\hat{D}_{2}=D_{A}+D_{B}$ is the relative diffusion coefficient, and
$D_{A}$ and $D_{B}$ are the diffusion coefficients associated with a molecule
of $A$ and $B$ respectively. Equation Eq. 1 allows us to select $\sigma$ so
that the reaction rate of our model matches the reaction rate of our
bimolecular reaction. Crucially, this relation is only valid once the
distribution of $B$ about $A$ molecules (and vice versa) reaches a steady
state, which usually happens very quickly; within $10$ns for a typical system
[34]. During the transient $t_{s}$ before the steady state is established, the
reaction rate is artificially inflated as any molecules of $A$ and $B$ that
are initialized within a distance $\sigma$ of each other undergo an
instantaneous reaction and others in the neighbourhood of $\sigma$ undergo an
inflated reaction rate temporarily. Due to this, the validity of
Smoluchowski’s model hinges on the assumption that a negligible number of
reactions occur during $t_{s}$. For this to be true, we require that the
expected separation of reactant pairs is large in comparison to $\sigma$ upon
initialization. In other words, the concentrations of $A$ and $B$ need to be
sufficiently small, or the reactants must have a sufficiently low affinity for
one another so that the associated reaction radius $\sigma$ is small.

These requirements quickly become problematic when examining biochemical
networks, since their action is often facilitated by biological catalysts
called enzymes. Enzymes bind selectively to compounds known as substrates to
enable essential biochemical reactions. They are fundamental to life and play
a critical role in metabolic processes, cell regulation and signal
transduction [31, 2]. Enzymatic systems are typically characterized by two
distinct timescales. There is a short timescale $t_{c}$ on which enzymes bind
to substrate molecules to form bound complexes, and a comparatively long
timescale $t_{p}$ on which molecules of the complex are converted to products,
degrading the substrate molecule and freeing the bound enzyme in the process.
The classical analysis of these systems is based on the pioneering work by
Michaelis and Menten who derived the degradation rate of the substrate for a
prototypical enzymatic system [30]. Although Michaelis and Menten did not
employ the theory themselves, it has long been known that the
characteristically non-linear degradation rate of the substrate on the long
timescale can be obtained by applying singular perturbation theory to the
system of ODEs that governs the reaction [37]. Central to the application of
the theory is the assumption that the timescale $t_{c}$ is effectively
instantaneous in comparison to $t_{p}$. In essence, it is assumed that free
enzymes bind so quickly to any available substrate molecules that the amount
of bound complex is always in an instantaneous steady state with the amount of
substrate; a so called quasi or pseudo-steady-state.

The pseudo-steady-state assumption poses a serious problem for Smoluchowski’s
framework. In order for the model to be valid, we require that relatively few
reactions occur during the initial transient. Conversely, if we are to mimic
the results of the singular perturbation analysis we require the reactions
between the enzyme and substrate molecules to be effectively instantaneous. If
we attempt to rectify this by explicitly simulating the dynamics on the fast
timescale then we are faced with the prospect of waiting for what in theory
would be an infinitely long time (and in practice a prohibitively long time)
before the substrate concentration has degraded appreciably. The issue is
further exacerbated because reactions involving enzymes are reversible,
meaning that the enzyme and substrate are free to unbind or disassociate from
each other before the substrate is converted to a product. Through this
mechanism, substrate and enzyme molecules are reintroduced into the system at
positions uncontrolled by the modeler. Several initialisation methods exist
that allow us to avoid artificial geminate recombination events, which arise
when the reactants are initialised too close together and so immediately react
[4, 27]. However, these methods consider the newly created pair of molecules
in isolation and do nothing to account for the artificial reactions that can
occur if the molecules are placed too close to other reactants in the system.
The related Collins Kimball model avoids the geminate recombination problem by
altering the Smoluchoswki reaction condition so that reactions only occur with
a fixed probability once reactants are deemed close enough [8, 1]. This
reduces the time it takes the system to reach a steady state following
initialisation, but the problematic transient is not removed entirely, and the
fundamental issue remains the same.

To address these issues, we propose a modification to Smoluchowski’s original
framework that allows us to reproduce the non-linear reaction rates
characteristic of enzymatic systems. We explore a new reaction condition
informed by the singular perturbation analysis of Michaelis-Menten kinetics
that incorporates the influence of the complex formation without requiring
events on $t_{c}$ to be explicitly modelled. This completely circumvents all
of the issues posed to by the disparate timescales, and enables us to directly
embed the results of singular perturbation theory within Smoluchowski’s
framework. The generalised reaction conditions can be easily implemented in
current software packages, and we demonstrate this by presenting proof of
concept simulations.

##  2 Mathematical Background

###  2.1 Generalized Smoluchowski Theory

In his work described in the introduction, Smoluchowski only considered
bimolecular reactions, but his theory can be viewed as a special case of a
more general framework that incorporates reactions of any order and was
developed by Flegg [17]. This extension is fundamental to our results in this
paper so we provide a brief summary, although for more details readers are
directed to the paper by Flegg.

Consider a system of $N$ distinct molecules that are initially well-mixed
(distributed uniformly at random) within a domain $\Omega$ of volume $V$,
where $V$ is finite but very large. In addition, let $D_{i}$ and
$\mathbf{x}_{i}$ to denote the respective diffusion constant and the
$3$-dimensional position of the i-th molecule. Since our reaction conditions
will be functions of the molecules’ relative proximity, it is convenient to
transform the coordinate system into diffusive Jacobi coordinates or
separation coordinates defined by,

$\displaystyle\bm{\eta}_{1}$ $\displaystyle=\bar{\bm{x}}_{N}$ (2)
$\displaystyle\bm{\eta}_{i}$
$\displaystyle=\bm{x}_{i}-\bar{\bm{x}}_{i-1},\quad
i=2,3,...,N,\quad\text{where,}$ (3) $\displaystyle\bar{\bm{x}}_{i}$
$\displaystyle=\frac{\sum^{i}_{j=1}\bm{x}_{j}D_{j}^{-1}}{\sum^{i}_{k=1}D_{k}^{-1}},$
(4)

is the ‘centre of diffusion’ of the first $i$ molecules (analogous to the
center of mass except that the positions are weighted by their inverse
diffusion coefficients rather than their masses). In our new coordinate system
$\bm{\eta}_{1}$ is the center of diffusion of the $N$ molecules and
$\bm{\eta}_{i}$ for $i>1$ describe the separation between $\bm{x}_{i}$ and the
center of diffusion, $\bar{\bm{x}}_{i-1}$, of the previous $i-1$ molecules.
The state vector $\bm{\eta}_{i}$ can be shown to undergo independent linear
diffusion with diffusion constant,

$\hat{D}_{i}=\begin{cases}\bar{D}_{N}\quad\text{ when }i=1,\\\
D_{i}+\bar{D}_{i-1}\quad\text{ when }i>1,\end{cases}$ (5)

where $\bar{D}_{j}$ is the diffusion constant associated with
$\bar{\bm{x}}_{j}$,

$\bar{D}_{j}=\frac{1}{\sum_{i=1}^{j}D_{i}^{-1}}.$ (6)

Under this specific transformation, the molecules, and by extension the state
vectors, diffuse independently. This framework is useful to work in since the
physical constraint that ‘translations of the whole system should not cause a
reaction’ can be simply translated as ‘reaction conditions must be independent
of $\bm{\eta}_{1}$’. As a result, for a well-mixed system, the joint
probability density, $P\left(,\bm{\eta},t\right)$, to find the reactants in an
unreacted state is also independent of $\bm{\eta}_{1}$ and is given by the
diffusion equation,

$\frac{\partial P\left(,\bm{\eta},t\right)}{\partial
t}=\left[\sum_{i=2}^{N}\hat{D}_{i}\hat{\nabla}^{2}_{i}\right]P\left(,\bm{\eta},t\right),$
(7)

where
$\bm{\eta}=\left\\{\bm{\eta}_{2},\bm{\eta}_{3},...,\bm{\eta}_{N}\right\\}$
describes the state of the system and $\hat{\nabla}^{2}_{i}$ denotes the
Laplacian with respect to the coordinates of $\bm{\eta}_{i}$. Initially, no
reactions have occurred and the probability density can be found by
normalization,

$P\left(\bm{\eta},0\right)=P_{\infty}=\frac{1}{V^{N-1}},$ (8)

noting that the power here is $N-1$ since the state does not depend on the
coordinates of $\bm{\eta}_{1}$ and so $P$ is the probability density to find
the other coordinates only. The inner boundary corresponds to the reaction
condition and defines a region, $\Omega_{R}$, upon whose boundary,
$\partial\Omega_{R}$, the system is absorbed and reacts,

$P\left(\bm{\eta}\in\partial\Omega_{R},t\right)=0.$ (9)

Sufficiently far from the origin we expect $P$ to be unperturbed by the
absorption of states on the inner boundary, and require,

$\text{lim}_{\bm{\eta}\rightarrow\infty}P\left(\bm{\eta},t\right)=P_{\infty}=\frac{1}{V^{N-1}}.$
(10)

To recover Smoluchowski’s original result we restrict ourselves to $N=2$ and
adopt the inner boundary,

$\Omega_{R}=\left\\{\bm{\eta}=\bm{\eta}_{2}:r_{2}\leq\sigma\right\\},\quad\text{where
}r_{2}=||\bm{\eta}_{2}||.$ (11)

Assuming the system is initially well-mixed and that $\sigma$ is comparatively
small, we can solve Equation Eq. 7 via Laplace transform which yields the
radially symmetric solution [34],

$P\left(r_{2},t\right)=\frac{1}{V}\left[1-\frac{\sigma}{r_{2}}\text{erfc}\left(\frac{r_{2}-\sigma}{\sqrt{4\hat{D}_{2}t}}\right)\right].$
(12)

The corresponding reaction rate is given by the total flux of the probability
density over the inner boundary,

$K=\frac{4\pi\sigma^{2}\hat{D}_{2}}{V}\left.\frac{\partial
P\left(r_{2},t\right)}{\partial
r_{2}}\right|_{\sigma}=\frac{4\pi\sigma\hat{D}_{2}}{V}\left(1+\frac{\sigma}{\sqrt{4\hat{D}_{2}t}}\right).$
(13)

We see that in general the probability density and hence Smoluchowski’s
reaction rate, is time dependent. The time dependent rate in Equation Eq. 13
is only well approximated by the more commonly stated steady-state rate in
Equation Eq. 1 once the probability density has converged sufficiently to the
pseudo-stationary distribution,

$\text{lim}_{t\rightarrow\infty}P\left(r_{2},t\right)=\frac{1}{V}\left(1-\frac{\sigma}{r_{2}}\right).$
(14)

Before the system reaches this steady state the reaction rate can be
arbitrarily high, diverging to infinity as we approach $t=0$. This transient
behavior is an artifact of the initial condition, where two reactants have a
non-zero probability of being initialized such that $r_{2}\leq\sigma$ (or even
an elevated chance of being initialized $r_{2}\sim\sigma$). Any such pair will
undergo an instantaneous (or increased rate of) reaction that is not diffusion
limited, i.e. the reaction rate is not determined by reactant diffusion. If
the expected separation between reactants is sufficiently large in comparison
to $\sigma$ then reactions that occur during the transient can be neglected
without incurring significant errors. This is the antithesis of the conditions
required for the pseudo-steady-state assumption to be valid. Instead, we
expect that a significant number of reactions will occur during this
transient, which would cause the reaction rate to exceed what is predicted by
Equation Eq. 1.

##  3 Generalized reaction conditions

To illustrate the aforementioned issues with Smoluchowski’s theory and our
proposed solution, we consider a network that consists of three reactions and
involves three chemical species $A$, $B$ and $C$. In the first reaction,
molecules of $A$ and $B$ are able to bind to form a chemical complex that we
will denote $X$. Once formed, a molecule of $X$ is able to disassociate into
its constituent $A$ and $B$ molecules, or it can react with a molecule of $C$,
converting the complex to some product $P$. The production of $P$ does not
consume $C$ molecules and instead frees them to participate in other
reactions. Our reaction network can be summarized succinctly in chemical
shorthand as,

$\ce{A+B<=>[k_{1}][k_{-1}]X}\quad\text{and}\quad\ce{X+C->[k_{2}]P+C},$ (15)

where $k_{1}$ and $k_{2}$ are the second-order rate constants that govern the
formation rate of $X$ and $P$ respectively, while $k_{-1}$ is the first-order
rate constant that controls the rate of dissociation of $X$ [9].

The Law of Mass Action states that the rate of an elementary reaction is
proportional to the product of the reactant concentrations, with the constant
of proportionality being the rate constant associated with the reaction.
Applying this to Reaction Eq. 15 yields a system of five ODEs that can be
reduced to just two equations,

$\displaystyle\frac{da}{dt}$ $\displaystyle=-k_{1}ab+k_{-1}x,\quad\text{and}$
(16) $\displaystyle\frac{dx}{dt}$ $\displaystyle=k_{1}ab-k_{-1}x-k_{2}xc,$
(17)

where the lower case letters denote the concentrations of the corresponding
chemical species. The traditional approach - more rigorously justified by
Briggs and Haladane [5] than Michaelis and Menten [30] \- now proceeds by
applying the pseudo-steady-state approximation to the system Eq. 16-Eq. 17.
The approximation is often stated verbatim, but it is more instructive to view
it as the result of singular perturbation theory. To obtain non-linear
kinetics we require the complex concentration to be in a pseudo-steady-state
on the long timescale. For this to occur the complex needs to be short-lived
so will we consider a situation where it dissociates much more quickly than it
is formed. In addition, we assume the concentration of $C$ is much larger than
that of $A$, or $B$, so that reactions between $X$ and $C$ occur much faster
than those between $A$ and $B$. That is, we assume that both $k_{-1}$ and the
initial concentration of $C$, $c_{0}$, are
$O\left(\frac{1}{\varepsilon}\right)$ for a sufficiently small positive
dimensionless parameter $\varepsilon$ such that,

$k_{-1}=\frac{\bar{k}_{-1}}{\varepsilon}\quad\text{and}\quad
c_{0}=\frac{\bar{c}_{0}}{\varepsilon},$ (18)

where $\bar{k}_{-1}$ and $\bar{c}_{0}$ are both $O(1)$.

To proceed we define dimensionless variables,

$\bar{a}=\frac{a}{a_{0}},\quad\bar{b}=\frac{b}{b_{0}},\quad\bar{c}=\frac{c}{c_{0}},\quad\bar{x}=\frac{x}{a_{0}},\quad\text{and}\quad
T=k_{1}b_{0}t,$ (19)

where $a_{0}$ and $b_{0}$ are the initial concentrations of $A$ and $B$
respectively. Applying our change of variables to the system Eq. 16-Eq. 17
gives,

$\displaystyle\varepsilon\frac{d\bar{a}}{dT}$
$\displaystyle=-\varepsilon\bar{a}\bar{b}+\mu\bar{x}$ (20)
$\displaystyle\varepsilon\frac{d\bar{x}}{dT}$
$\displaystyle=\varepsilon\bar{a}\bar{b}-\mu\bar{x}-\nu\bar{x}\bar{c}$ (21)

where we have defined dimensionless parameters,

$\mu=\frac{\bar{k}_{-1}}{k_{1}b_{0}},\quad\text{and}\quad\nu=\frac{k_{2}\bar{c}_{0}}{k_{1}b_{0}}.$
(22)

From Equation Eq. 20 we find that $\bar{x}$ is $O(\varepsilon)$ in which case
Equation Eq. 21 yields,

$\bar{x}\left(T\right)=\frac{\varepsilon\bar{a}\bar{b}}{\mu+\nu\bar{c}}+O\left(\varepsilon^{2}\right).$
(23)

Substituting Equation Eq. 23 into Equation Eq. 20 and redimensionalizing gives
to leading order in $\varepsilon$,

$\frac{da}{dt}=-\frac{k_{1}abc}{\Gamma+c},$ (24)

where,

$\Gamma=\frac{k_{-1}}{k_{2}},$ (25)

is the Michaelis constant.

The preceding analysis shows that on the long timescale Reaction Eq. 15 can be
reduced to a single trimolecular reaction,

$\ce{A+B+C->[k_{3}\left(c\right)]P+C},$ (26)

where $k_{3}$ is a third-order rate constant,

$k_{3}\left(c\right)=\frac{k_{1}}{\Gamma+c}.$ (27)

This reduction is valid so long as $\varepsilon$ is sufficiently small. This
requires that the rate at which a single molecule of $X$ forms is slow when
compared to the rate at which it dissociates and/or the rate at which it is
converted into the product. In our analysis we assumed that the reaction
between $X$ and $C$ was fast due to $c_{0}$ being very large in comparison to
$a_{0}$ and $b_{0}$, however we could have assumed instead that $k_{2}$ was
much larger than $k_{1}$ and arrived at the same result. Since Smoluchowski’s
theory is unable to incorporate the disparate timescales required to arrive at
Reaction Eq. 26, we seek to design a new trimolecular reaction condition that
reproduces the reaction rate in Equation Eq. 27 directly.

Suppose we alter our system slightly so that it now contains $N_{c}=cV$
molecules of $C$ \- but still just a single molecule each of $A$ and $B$ \-
where, as before, $c$ is a well-mixed concentration of $C$ molecules. Since
there are now multiple molecules of $C$ we must consider $N_{c}$ distinct
states; one for each combination of $A$, $B$ and $C$ molecules. We recall that
$\bm{\eta}_{2}$ and $\bm{\eta}_{3}$ diffuse independently so that,
$P\left(\bm{\eta},t\right)=P_{2}\left(\bm{\eta}_{2},t\right)P_{3}\left(\bm{\eta}_{3},t\right)$,
and,

$\displaystyle\left[\frac{\partial}{\partial
t}-\hat{D}_{2}\hat{\nabla}^{2}_{2}\right]P_{2}\left(\bm{\eta}_{2},t\right)$
$\displaystyle\equiv\mathcal{L}_{2}P_{2}\left(\bm{\eta}_{2},t\right)=0,\quad\text{and}$
(28) $\displaystyle\left[\frac{\partial}{\partial
t}-\hat{D}_{3}\hat{\nabla}^{2}_{3}\right]P_{3}\left(\bm{\eta}_{3},t\right)$
$\displaystyle\equiv\mathcal{L}_{3}P_{3}\left(\bm{\eta}_{3},t\right)=0,$ (29)

where $\mathcal{L}_{2}$ and $\mathcal{L}_{3}$ are diffusion operators on the
$3$-dimensional spaces defined by $\bm{\eta}_{2}$ and $\bm{\eta}_{3}$
respectively. Since there is just one molecule of $B$, all states (which
consist of one $A$, one $B$ and one $C$ molecule) lie on manifolds of constant
$\bm{\eta}_{2}$, whilst the specific instance of $\bm{\eta}_{2}$ also diffuses
according to Equation Eq. 28. On the manifold, states diffuse independently in
the $\bm{\eta}_{3}$ space in accordance with Equation Eq. 29. The first state
incident on the inner boundary $\partial\Omega$ will cause a reaction, and
therefore we wish to know the dynamics of the state with the minimum magnitude
$|\bm{\eta}_{3}|$.

We need to understand the well-mixed steady state of this system. In the well-
mixed state, we can assume that $\bm{\eta}_{3}$ of the $N_{c}$ states are
uniformly and independently distributed in $\Omega$. Consider now a single
particular molecule of $C$ and let $\bm{\mathcal{H}}_{3}$ denote the event
that this particular molecule is associated with the minimum $|\bm{\eta}_{3}|$
when compared with any other $C$ molecule in the system. Moreover, let the
probability distribution function,
$P\left(,\bm{\eta}_{3},t|\bm{\mathcal{H}}_{3}\right)$, denote the probability
density for finding this molecule of $C$ at $\bm{\eta}_{3}$ at time $t$, given
the fact that has the smallest $|\bm{\eta}_{3}|$ of any molecule of $C$. Since
the molecules are well-mixed, the probability that an arbitrary $C$ molecule
has the minimum $|\bm{\eta}_{3}|$ is,

$P\left(,\bm{\mathcal{H}}_{3},t\right)=\frac{1}{N_{c}}.$ (30)

The probability of $\bm{\mathcal{H}}_{3}$, that a particular molecule of $C$
is closest to the origin, given it has a known $\bm{\eta}_{3}$, is equal to
the probability that all the other $N_{c}-1$ molecules of $C$ lie outside a
sphere, $V_{3}$, of radius $r_{3}=||\bm{\eta}_{3}||$ centered on the origin,
i.e.,

$P\left(,\bm{\mathcal{H}}_{3},t|\bm{\eta}_{3}\right)=\left[1-\int_{V_{3}}P\left(,\bm{\eta^{\prime}}_{3},t\right)dV_{3}^{\prime}\right]^{N_{c}-1},$
(31)

where $dV_{3}^{\prime}$ is an elemental volume for coordinates of
$\bm{\eta^{\prime}}_{3}$. Bayes Theorem then yields,

$P\left(,\bm{\eta}_{3},t|\bm{\mathcal{H}}_{3}\right)=N_{c}\left[1-\int_{V_{3}}P\left(,\bm{\eta^{\prime}}_{3},t\right)dV_{3}^{\prime}\right]^{N_{c}-1}P\left(,\bm{\eta}_{3},t\right).$
(32)

The system is very large and in the limit that $V$ \- and hence $N_{c}$ \-
tends to infinity, $P\left(,\bm{\eta}_{3},t|\bm{\mathcal{H}}_{3}\right)$ goes
to zero in accordance with Equation Eq. 30. To ensure we take the appropriate
limit in Equation Eq. 32 we define the scaled probability distributions,

which when substituted into Equation Eq. 32 give,

$\Phi\left(\bm{\eta}_{3},t\right)=\phi\left(\bm{\eta}_{3},t\right)\left[1-\frac{c}{N_{c}}\int_{V_{3}}\phi\left(\bm{\eta^{\prime}}_{3},t\right)dV_{3}^{\prime}\right]^{N_{c}-1}.$
(33)

Taking the limit $N_{c}\rightarrow\infty$ we obtain,

$\Phi\left(\bm{\eta}_{3},t\right)=\phi\left(\bm{\eta}_{3},t\right)\text{exp}\left(-c\int_{V_{3}}\phi\left(\bm{\eta^{\prime}}_{3},t\right)dV_{3}^{\prime}\right).$
(34)

That is, using (28) and (29), $\Phi\left(\bm{\eta}_{3},t\right)$ evolves
according the diffusion-advection equation,

$\frac{\partial\Phi\left(\bm{\eta}_{3},t\right)}{\partial
t}=\hat{D}_{3}\hat{\nabla}^{2}_{3}\Phi\left(\bm{\eta}_{3},t\right)+\hat{D}_{3}\hat{\nabla}_{3}\cdot\left(4\pi
r_{3}^{2}c\phi\Phi\left(\bm{\eta}_{3},t\right)\hat{\bm{r}_{3}}\right),$ (35)

as shown in Appendix A. We note here that the isotropic linear diffusion term
describes the independent Brownian motion of $\bm{\eta}_{2}$ and
$\bm{\eta}_{3}$ whilst advection towards $\bm{\eta}_{3}=\mathbf{0}$ represents
the flux of the likelihood that $C$ molecule with the second smallest
$|\bm{\eta}_{3}|$ value diffuses over the sphere of radius $r_{3}$ set by the
current $C$ molecule.

Typically unless a boundary-free steady state in $\mathbb{R}^{6}$ is sought
(see Equation (34)), this PDE is very difficult to solve. This is because of
the intrinsic relationship between $\phi$ and $\Phi$. In our particular case
however, we will be assuming that there is a very thin absorbing boundary
which is long in $r_{3}=|\bm{\eta}_{3}|$ but thin in $r_{2}=|\bm{\eta}_{2}|$.
We expect therefore that $\phi$ is equal to its well-mixed value of $\phi=1$
with a small perturbation caused by undulations of the absorbing surface. As
we will only concern ourselves with the leading order solution of this PDE
with a thin absorbing boundary, using $\phi=1$ and Equation
LABEL:eq:Phi_definition we arrive at the governing equation for
$P\left(,\bm{\eta}_{3},t|\bm{\mathcal{H}}_{3}\right)$,

$\frac{\partial P\left(,\bm{\eta}_{3},t|\bm{\mathcal{H}}_{3}\right)}{\partial
t}=\hat{D}_{3}\hat{\nabla}^{2}_{3}P\left(,\bm{\eta}_{3},t|\bm{\mathcal{H}}_{3}\right)+\hat{D}_{3}\hat{\nabla}_{3}\cdot\left(4\pi
cr_{3}^{2}P\left(,\bm{\eta}_{3},t|\bm{\mathcal{H}}_{3}\right)\hat{\bm{r}_{3}}\right).$
(36)

As $\bm{\eta}_{2}$ diffuses independently of $\bm{\eta}_{3}$ the evolution of
the joint probability density,
$P\left(,\bm{\eta},t|\bm{\mathcal{H}}_{3}\right)=P\left(\bm{\eta}_{2},\bm{\eta}_{3},t|\bm{\mathcal{H}}_{3}\right)$,
for finding the separation of the state with the minimum value of
$\bm{\eta}_{3}$ is governed by,

$\frac{\partial P\left(,\bm{\eta},t|\bm{\mathcal{H}}_{3}\right)}{\partial
t}=\left[\sum_{i=2}^{3}\hat{D}_{i}\hat{\nabla}^{2}_{i}\right]P\left(,\bm{\eta},t|\bm{\mathcal{H}}_{3}\right)+\hat{D}_{3}\hat{\nabla}_{3}\cdot\left(4\pi
cr_{3}^{2}P\left(,\bm{\eta},t|\bm{\mathcal{H}}_{3}\right)\hat{\bm{r}_{3}}\right).$
(37)

To find the correct boundary conditions for the probability
$P\left(,\bm{\eta},t|\bm{\mathcal{H}}_{3}\right)$ we need to remind ourselves
that this probability is zero (absorbing boundary condition) if the state
$\bm{\eta}$ ever reaches a manifold on which a reaction condition is met. We
note that if the condition for a reaction occurs on an absorbing boundary
extending small distances $r_{2}$ and $r_{3}$ then the advection term in (36)
becomes negligible as was the case for the trimolecular generalization of
Smolochowski reaction condition presented by Flegg [17] (although never
directly addressed in that paper). That is, any small absorbing boundary will
reach a pseudo-equilibrium that returns trimolecular mass action. Instead, in
order to obtain nonlinear reaction kinetics in the concentration of $C$ we
propose a long thin boundary $\partial\Omega_{R}$ where

$\Omega_{R}=\left\\{(\bm{\eta}_{2},\bm{\eta}_{3}):r_{2}<f(r_{3})\right\\},\quad
r_{2}=||\bm{\eta}_{2}||,$ (38)

where $0<f(r_{3})<\sigma$ is a monotonically decreasing function of $r_{3}$
and $\sigma$ is a small positive constant relative to the diffusion
coefficients $\hat{D}_{i}$ and the characteristic scale of the domain of $f$.
Since this boundary is thin, to leading order we have that
$\mathbf{\hat{n}}=\bm{\hat{\eta}}_{2}$ and that the steady state solution to
(37) with $P\left(,\bm{\eta},t|\bm{\mathcal{H}}_{3}\right)=0$ on
$\partial\Omega_{R}$ where $\Omega_{R}$ is given by (38) is weakly dependent
on $\bm{\eta}_{3}$ compared to $\bm{\eta}_{2}$. That is, the problem reduces
approximately to solving for the steady state of the problem

$\frac{\partial P\left(,\bm{\eta}_{2},t\right)}{\partial
t}=\hat{D}_{2}\hat{\nabla}^{2}_{2}P\left(,\bm{\eta}_{2},t\right),\quad\text{with
}P\left(,\bm{\eta}_{2},t\right)=0\text{ when }r_{2}=f(r_{3}),$ (39)

and the value of $P\left(,\bm{\eta}_{2},t\right)=Q(r_{3})$ at infinity is
given by the steady state solution of (36) on $\mathbb{R}^{3}$ – that is,
$Q(r_{3})=c\exp(-4\pi r_{3}^{3}c/3)$.

Finding the flux in the $\mathbf{\hat{n}}$ direction over $\partial\Omega_{R}$
and matching it to reaction rate $K$ is the same as solving the bimolecular
Smoluchowski reaction boundary problem for the reaction rate at each $r_{3}$
where the reaction radius is $r_{2}=f(r_{3})$, call this $K_{2}(f(r_{3})))$
where $K(\rho)=4\pi\hat{D}_{2}\rho$ is the well known Smoluchowski result,
multiplying this rate by the probability of finding a reaction with a given
$r_{3}$, $Q(r_{3})$, and integrating over all spheres of radius $r_{3}$ to
find the total flux. That is,

$K=\int^{\infty}_{0}K_{2}\left(f(r_{3})\right)Q\left(r_{3}\right)4\pi
r_{3}^{2}dr_{3}.$ (40)

Making the substitution $\lambda=4\pi r_{3}^{3}/3$ and recalling that $K\equiv
K(c)$ here is a rate that depends on the concentration $c$ for each molecule
of $A$ and $B$,

$\frac{K(c)}{4\pi c\hat{D}_{2}}=\int^{\infty}_{0}F(\lambda)\exp(-\lambda
c)d\lambda=\mathcal{L}(F(\lambda)).$ (41)

where $\mathcal{L}$ is the Laplace transform and $F(4\pi
r_{3}^{3}/3)=f(r_{3})$. Recalling that for the Michaelis-Menten system (15),
$K(c)=ck_{3}(c)$ where $k_{3}(c)$ is given by Eq. 27. In the case of
Michaelis-Menten kinetics it is easy therefore to take the inverse Laplace
transform to find $F$ and therefore $f$ (the unknown function that we require
to construct a reaction condition),

$f(r_{3})=\frac{k_{1}}{4\pi\hat{D}_{2}}\text{exp}\left(\frac{-4\pi\Gamma
r_{3}^{3}}{3}\right).$ (42)

Finally, substituting this result into Equation Eq. 38 we obtain a reaction
boundary $\partial\Omega_{MM}$,

$\Omega_{MM}=\left\\{\bm{\eta}:r_{2}<\frac{k_{1}}{4\pi\hat{D}_{2}}\text{exp}\left(\frac{-4\pi\Gamma
r_{3}^{3}}{3}\right)\right\\},$ (43)

that reproduces Michaelis-Menten kinetics.

##  4 Simulation of non-Linear kinetics

Particle based simulations of Smoluchowski’s framework typically use either a
time-driven (TD), or an event-driven (ED) methodology. TD approaches progress
through time using finite preset timesteps and the state of the system - the
position of each molecule - is updated during each timestep. Following the
position updates, the distance between relevant reactants can be calculated
and reactions are performed according to the imposed reaction condition. In
contrast, ED algorithms calculate the first passage times associated with the
movement of molecules, enabling accurate sampling of the exact event times.
For instance, the distribution for the first passage time for a molecule to
leave a section of the domain and the time at which two molecules first
approach within a set distance are usually of interest. Event times are
sampled from the appropriate distribution and the system is updated according
to a time ordered queue of events that is dynamically updated as each event is
processed. Inspired by the approach taken by Vijaykumar, Bolhuis and Ten Wolde
[50] we adopt a hybrid methodology that contains both TD and ED components.
The simulation switches between these two modes - referred to as TD mode and
ED mode henceforth - in order to efficiently track molecule positions and
apply reaction conditions.

The TD components of our simulation are based on a popular constant timestep
algorithm developed by Andrews and Bray in $2004$ [4]. Their algorithm is
simple yet accurate, making it easily adaptable to modifications of
Smoluchowski’s original framework. The original algorithm has been implemented
in the Smolydn software package, which has been used widely in literature [3].
Smolydn simulates diffusion and bimolecular reactions with single molecule
detail and can be extended in a straightforward manner to higher order
reactions [17].

In TD mode, the simulation progresses through time via a discrete timestep
$\Delta t$. The position, $\bm{x}_{i}\left(t\right)$, of each molecule $i$ is
updated randomly each step according to [12],

$\bm{x}_{i}(t+\Delta t)=\bm{x}_{i}(t)+\sqrt{2D_{i}\Delta t}\bm{\xi}_{i},$ (44)

where $\bm{\xi}_{i}$ is a three-dimensional vector of independent, normally
distributed, random numbers with unit variance and zero mean. Once the
positions have been updated, the separation vectors $\bm{\eta}_{2}$ and
$\bm{\eta}_{3}$ are calculated according to Equation Eq. 3 and Condition Eq.
43 is tested to determine if any reaction events took place during the last
timestep. Equation Eq. 44 provides an exact simulation of molecular diffusion,
but does not account for the fact that reactants should instantaneously react
once their separations satisfy the reaction condition. Molecules do not
undergo continuous motion and can ‘jump’ through the reaction boundary,
artificially skipping a reaction during the timestep. The issue can be
circumvented by taking the timestep to be sufficiently small, but this usually
requires excessive computational time. Instead, it is common to make use of a
numerically derived reaction boundary that is slightly larger than the
continuous-time theoretical boundary. The required numerical radius is
determined by the desired reaction rate and $\Delta t$ allowing it to be
precomputed. Corrections of this form are straightforward for fixed reaction
boundaries, but for condition Eq. 43 the size of the boundary is dictated by
the current separation vectors. In this case, we need to compute a series of
corrections so that the appropriate value may be retrieved for any combination
of separation vectors.

Reaction Eq. 26 is trimolecular, but it is convenient to imagine that it is
the result of two bimolecular reactions so that the original Smoldyn protocol
can be used. We treat molecules of $A$ and $B$ as undergoing a bimolecular
reaction to form a complex $AB$ according to the reaction condition,

$\Omega_{AB}=\left\\{\bm{\eta}_{2}:r_{2}\leq f(0)\right\\},$ (45)

where $f$ is given by Equation Eq. 42. Molecules of $AB$ can then react with
molecules of $C$ to produce molecules of $P$ so that the system evolves
according to,

$\ce{A+B<=>AB}\quad\text{and}\quad\ce{AB+C->P+C}.$ (46)

A molecule of $AB$ is initialized at the center of diffusion of the reacting
$A$ and $B$ molecules, which is defined in Equation Eq. 4 as,

$\bar{\bm{x}}_{2}=\frac{D_{1}^{-1}}{D_{1}^{-1}+D_{2}^{-1}}\bm{x}_{1}+\frac{D_{2}^{-1}}{D_{1}^{-1}+D_{2}^{-1}}\bm{x}_{2},$
(47)

where the labels $1$ and $2$ have been used to refer to the molecule $A$ and
$B$ respectively. Similarly, the corresponding separation vector can be
calculated according to Equation Eq. 3,

$\bm{\eta}_{2}=\bm{x}_{2}-\bm{x}_{1}.$ (48)

Both $\bar{\bm{x}}_{2}$ and $\bm{\eta}_{2}$ undergo independent isotropic
linear diffusion with the respective diffusion constants $\bar{D}_{2}$
(Equation Eq. 6) and $\hat{D}_{2}$ (Equation Eq. 5) [17], and can be updated
in TD mode using Equation Eq. 44.

Reaction Eq. 46 is deceptively similar to Reaction Eq. 15 with the complexes
$AB$ and $X$ appearing interchangeable. However, $AB$ is not actually a true
chemical species. Instead a molecule of $AB$ is constructed purely as a
convenient way to track a pair of $A$ and $B$ molecules that are close enough
that a reaction with a molecule of $C$ might be possible. Pairs of $A$ and $B$
molecules that are not close enough to form one of these virtual complexes are
treated as being nonreactive with $C$, which avoids unnecessary testing of the
reaction condition. If the constituents of a molecule of $AB$ move far enough
apart that reaction with $C$ becomes impossible the complex is dissolved and
the $A$ and $B$ molecules are reintroduced at the positions,

$\bm{x}_{1}=\bar{\bm{x}}_{2}-\frac{D_{2}^{-1}}{D_{1}^{-1}+D_{2}^{-1}}\bm{\eta}_{2}\quad\text{and}\quad\bm{x}_{2}=\bar{\bm{x}}_{2}+\frac{D_{1}^{-1}}{D_{1}^{-1}+D_{2}^{-1}}\bm{\eta}_{2},$
(49)

respectively. When viewed in this way Reaction Eq. 26 can be thought of as a
bimolecular reaction between $A$ and $B$ where the reaction boundary for each
pair of molecules is determined by the proximity of the associated complex to
the closest molecule of $C$. The position of the complex is the center of
diffusion of the $A$ and $B$ molecules allowing for simple calculation of the
separation to each molecule of $C$ using Equation Eq. 3,

$\bm{\eta}_{3}=\bm{x}_{3}-\bar{\bm{x}}_{2},$ (50)

where the subscript $3$ is associated with a particular molecule of $C$. Once
the minimum value of $\bm{\eta}_{3}$ has been found, Condition Eq. 43 can be
checked to determine if $A$ and $B$ are close enough for a reaction to occur,
as shown in Figure 1. The relevant reaction radius is selected from a
continuum of radii that need to be corrected to account for the use of finite
timesteps. To enable this, we precompute a table of corrections following the
protocol described in [4] so that the numerical reaction radius that
corresponds to any value of $f(\bm{\eta}_{3})$ can be retrieved efficiently.

(a)

(b)

Figure 1: (a) A pair of molecules $A$ and $B$ \- depicted as a blue and green
point respectively - are simulated as undergoing a bimolecular reaction to
form a molecule of $AB$ whenever their separation $\bm{\eta}_{2}$ is such that
$||\bm{\eta}_{2}||\leq f(0)$. This reaction condition (see Equation Eq. 45)
defines a spherical boundary of radius $f(0)$ centered on $A$, which is shown
as a dotted circle. The radius of the boundary represents the greatest
separation between the two molecules at which it is possible for Reaction 26
to occur according to Condition Eq. 43 and is illustrated in the graph on the
right as a dotted line. (b) When $A$ and $B$ react, a molecule of $AB$ (the
gray point) is initialised at their center of diffusion, $\bar{\bm{x}}_{2}$.
The separation, $\bm{\eta}_{3}$, between $AB$ and the closest molecule of $C$
(the red point) defines the radius of the reaction boundary that surrounds $A$
as shown by the dotted circle on the left and the dotted line on the graph of
Condition Eq. 43 on the right. $\bar{\bm{x}}_{2}$ and $\bm{\eta}_{2}$ diffuse
independently and can be updated without explicitly tracking $A$ and $B$
although their positions can be recovered at any time as indicated by the
dotted outlines of the blue and green points. If $||\bm{\eta}_{2}||\leq
f(||\bm{\eta}_{3}||)$ during any timestep a reaction is performed. Similarly,
if $||\bm{\eta}_{2}||>f(0)$ the molecule of $AB$ is dissolved and $A$ and $B$
are reintroduced at their respective positions.

Through the use of $AB$ molecules we are able to avoid unnecessary testing of
the reaction condition, but a significant amount of time can still be wasted
propagating $A$ and $B$ molecules that are not close enough to be considered
reactive. This is a common criticism of TD methods and the issue is mitigated
by switching the simulation of such molecules to ED mode. The exact position
of a molecule only needs to be known if it could be involved in a reaction in
the next timestep. Molecules that are isolated from other relevant reactants
do not need to be tracked explicitly, and instead can be placed in single
particle domains equivalent to those defined in eGFRD [44, 41]. The
introduction of protective domains means that at any time the simulation may
contain a mixture of molecules in TD and ED mode, as shown in Figure 2. The
escape time of each ED molecule is placed in a queue, which is used to
determine if an escape event will occur within the next TD timestep. If an
event is scheduled to occur between $\left(t,t+\Delta t\right]$ the
corresponding event time, $t_{\text{event}}$, is removed from the queue. The
event is then processed, and the position of all molecules in TD mode are
updated using the timestep, $\Delta t^{\prime}=t_{\text{event}}-t$. Since the
positions of molecules in TD mode are updated at every step a molecule may
approach a protective domain before the scheduled escape time. If a molecule
approaches close enough - maybe even crossing the boundary of the domain - it
is possible for a reaction to occur between it and the molecule within the
domain. However, we cannot be sure of the exact position of the molecule
within the protective domain, so to avoid missing a potential reaction we must
prematurely dissolve, or burst the domain. Upon bursting a domain the position
of the enclosed molecule is sampled as described in [41]. If the newly sampled
position is too close to a neighbouring domain this domain is also burst and
this process continues until all domains have been burst or are sufficiently
isolated from any of the molecules in TD mode. Any molecule that has escaped
its domain or had it burst, is placed in TD mode until it becomes sufficiently
isolated from the other reactants to be placed in a new protective domain.

Protective domains allow the simulation to make large adaptive jumps forward
in time during uninteresting periods where all the reactants are too far from
each other for reaction events to be possible. If two or more reactants are
close enough that reaction conditions need to be checked, then only the
relevant molecules need be switched to TD mode, greatly reducing the number of
reactant combinations that have to be considered. Protective domains can be
applied to $AB$ molecules, but they are more like the pair domains used in
eGFRD than the single particle domains described thus far. Similar to a
molecule of $A$, $B$, or $C$, an $AB$ molecule might escape its protective
domain by simply reaching its boundary, or the domain might be burst because a
$C$ molecule moved close enough that a reaction is possible. However, an
additional domain needs to be constructed for $\bm{\eta}_{2}$ since the
complex should be dissolved if $r_{2}>f(0)$ and a reaction should occur if the
closest $C$ molecule is such that $r_{2}\leq f(r_{3})$. This means that each
molecule of $AB$ has two protective domains from which it can escape; one for
each of $\bar{\bm{x}}_{2}$ and $\bm{\eta}_{2}$. The protective domain for
$\bar{\bm{x}}_{2}$ is equivalent to the single particle domain, and the
problem is identical to that of the ‘center of motion’ considered in eGFRD.
Similarly, the problem for $\bm{\eta}_{2}$ is analogous to that of the ‘inter-
particle vector’ in eGFRD, but differs crucially in that the inner boundary
condition that corresponds to the reaction of the $A$ and $B$ molecules is no
longer static. Instead, this boundary is determined by the proximity of the
nearest molecule of $C$ preventing reuse of the standard pair domain methods.
If we ignore this inner boundary it is possible to find the greens function
for an $\bm{\eta}_{2}$ domain that only has an outer boundary at $r_{2}=f(0)$.
This additional domain means that both $\bar{\bm{x}}_{2}$ and $\bm{\eta}_{2}$
have to be sampled when an escape event occurs or if the domain is burst, but
otherwise the domain functions similarly to a single particle domain.
Unfortunately, the utility of such domains is diminished by the fact that
reactions can occur even when $r_{3}$ is large so long as $r_{2}$ is
sufficiently small.

Figure 2: The molecule labeled $B_{1}$ is in close proximity to the reactive
molecule $A$ and its position must be tracked explicitly and updated each
timestep $\Delta t$ to enable the reaction condition - shown as a dotted
circle - to be checked. However, the molecules $B_{2}$ and $B_{3}$ are too far
from $A$ to react within the next timestep and their positions do not need to
be known exactly until $A$ approaches close enough that the reaction condition
needs to be tested. To avoid computing position updates for $B_{2}$ and
$B_{3}$ they are placed in single particle protective domains, labeled
$\mathcal{D}^{R_{2}}_{2}$ and $\mathcal{D}^{R_{3}}_{3}$ respectively, until
they either escape, or $A$ approaches close enough to the domain to burst it.

##  5 Numerical Results

Our modified Smoluchowski framework is not restricted to Michaelis-Menten
kinetics, but in order to demonstrate its validity and explore its
implications we will restrict our attention to three numerical tests of
Condition Eq. 43. The first test considers a well-mixed trimolecular system
and is designed to validate the simulation framework presented in section 4.
The second test examines the robustness of our approximation to the reaction
boundary given in Equation Eq. 42 by varying its shape. While the third and
final test explores whether Condition Eq. 43 reproduces the non-linear
behaviour described in Equation Eq. 27 as the concentration $c$ is varied.

###  5.1 Well-mixed trimolecular test system

To verify our simulation method we consider the trimolecular system,

$\ce{A+B+C->[k_{3}\left(c\right)]C},\quad\ce{\emptyset->[k_{0}]A},$ (51)

where $k_{0}$ is a zeroth order rate constant that controls the production of
$A$ and $k_{3}\left(c\right)$ is given in Equation Eq. 27. The population of
$A$ molecules in the system is governed by the Law of Mass Action,

$\frac{dN_{A}}{dt}=k_{0}^{\prime}-\frac{k_{1}^{\prime}}{\Gamma^{\prime}+N_{C}}N_{A}N_{B}N_{C},$
(52)

where $N_{A}$, $N_{B}$ and $N_{C}$ are the respective number of $A$, $B$ and
$C$ molecules at time $t$, and we have defined the parameters,

$k_{0}^{\prime}=Vk_{0},\quad
k_{1}^{\prime}=\frac{k_{1}}{V},\quad\text{and}\quad\Gamma^{\prime}=V\Gamma.$
(53)

The simulation is conducted in a well-mixed domain that is a dimensionless
unit cube with periodic boundary conditions. Within the domain we randomly
place a single immortal molecule of $B$, so that $N_{B}=1$ for the duration of
the simulation. Accounting for this single molecule of $B$ and defining
dimensionless variables,

$\bar{N}_{A}=\frac{k_{1}^{\prime}N_{B}}{k_{0}^{\prime}}N_{A},\quad\bar{N}_{c}=\frac{N_{C}}{\Gamma^{\prime}},\quad\text{and}\quad\tau=k_{1}^{\prime}N_{B}t,$
(54)

Equation Eq. 52 becomes,

$\frac{d\bar{N}_{A}}{d\tau}=1-\frac{\bar{N}_{A}\bar{N}_{C}}{1+\bar{N}_{C}}.$
(55)

We set $\bar{N}_{C}=5$ by choosing $\Gamma^{\prime}=1$ and placing $N_{C}=5$
immortal $C$ molecules uniformly at random within the domain. Using
dimensionless diffusion coefficients $D_{1}=D_{2}=D_{3}=1$ and time steps of
$\Delta\tau=2\times 10^{-6}$ we simulate the system for a non-dimensional
duration of $\tau=10$. The steady state distribution of $\bar{N}_{A}$ is
governed by the master equation associated with this birth-death process, from
which we expect a Poisson distribution [13],

$\bar{N}_{A}=\text{Pois}\left(\frac{1+\bar{N}_{C}}{\bar{N}_{C}}\right).$ (56)

The scaled rate constants $k_{0}^{\prime}$ and $k_{1}^{\prime}$ are chosen
such that $\bar{N}_{A}=\frac{1}{5}N_{A}$ and at the beginning of each
simulation, we sample $N_{A}$ from the expected steady state population,
$N_{A}=\text{Pois}\left(6\right)$, and uniformly distribute the molecules
throughout the volume. The simulation was repeated $2\times 10^{4}$ times and
the population of $A$ molecules was sampled at fixed intervals. In Figure 3,
the sampled distribution of $N_{A}$ (blue bars) is compared with the expected
Poisson distribution (red dots). While the mean of the initial Poisson
distribution has remained unchanged $6.00\pm 0.01$, the variance has been
reduced slightly from $6$ to $5.96\pm 0.02$. This indicates a slight increase
in the reaction rate, which may be the result of the corrections made to the
reaction radii required by the TD portion of the simulation. Despite these
discrepancies, the simulated distribution remains approximately Poissonian and
the simulation demonstrates good agreement with the theoretical kinetics for
these choices of reaction parameters considered.

Figure 3: The simulated distribution of the number of molecules of $A$
($N_{A}$) present in System Eq. 51. The value of $N_{A}$ in the simulation is
sampled with a dimensionless time interval of $1$ over a dimensionless
duration of $\tau=10$. The blue bars depict the distribution sampled from
$2\times 10^{4}$ simulations of the generalised reaction condition, while the
red dots represent the Poisson distribution predicted by the law of mass
action. The simulated reaction rate slightly exceeds the expected rate causing
it to favour lower tail of the distribution.

###  5.2 Non-linear kinetics

System Eq. 51 degrades $A$ at a rate that grows non-linearly in $N_{C}$ in
accordance with Equation Eq. 52. The dimensionless reaction rate is given by
the inverse of the mean of the steady state distribution for $\bar{N}_{A}$ in
Equation Eq. 56,

$K_{E}=\frac{\bar{N}_{C}}{1+\bar{N}_{C}}.$ (57)

Similarly, a simulated value for this rate can be calculated from the inverse
of the mean of the steady-state distribution of $\bar{N}_{A}$ obtained from a
series of simulations.

To demonstrate that our framework reproduces this expected non-linear
behaviour we use the same methodology described in subsection 5.1 to simulate
System Eq. 51 with $\Gamma^{\prime}=2$ and $k_{0}^{\prime}$ and
$k_{1}^{\prime}$ chosen such that $\bar{N_{A}}=\frac{3}{10}N_{A}$. These
reaction parameters were selected since they yield convenient dimensional
reaction rates for several of the $\bar{N}_{C}$ values considered but are of
no particular importance otherwise. By altering the value of $N_{C}$
accordingly we perform $1\times 10^{4}$ simulations for each value in
$\bar{N}_{C}\in\left\\{1,2,3,4,5,7.5,10,15,20\right\\}$ and calculate the
corresponding dimensionless reaction rates. The simulated rate values are
shown by the blue points in Figure 4, where they are compared against a plot
of Equation Eq. 57 shown by the red dashed line. We can see that although the
reaction parameters have been altered slightly from subsection 5.1, the
simulation again reproduces the theoretical reaction rate almost exactly for
$\bar{N}_{C}=5$. For values of $\bar{N}_{C}$ above this point we see that the
simulated rate diverges from the theoretical rate, with the ratio of the
simulated to theoretical values decreasing as $\bar{N}_{C}$ increases. This is
likely a result of the fact that the chance a reaction event is missed during
the simulation increases as the number of $C$ molecules is increased due to
crowding.

Initially for values below $\bar{N}_{C}=5$ we observed the opposite trend and
the ratio of the simulated rates to the theoretical rates increased
significantly as $\bar{N}_{C}$ decreased. We attributed the majority of this
discrepancy to the fact that the simulated system becomes an increasingly poor
reproduction of the situation considered in section 3. The derivation of
condition Eq. 43 is only valid in the limit that $V\rightarrow\infty$ so that
an effectively infinite number of independent $C$ molecules are present
regardless of the concentration $c$. In the simulation we approximate this
infinite population by placing periodic boundary conditions on our volume,
however the ‘image’ position of any molecule within the volume is perfectly
correlated with the molecule’s actual position. If we consider $\bar{N}_{C}=1$
as an example, the each simulation is initialised with $N_{C}=2$, which means
that the simulation only contains two independent $C$ molecules within the
volume. To obtain the points in Figure 4 we have attempted to correct for this
effect by expanding the volume and increasing $N_{C}$ in such a way that
$\bar{N}_{C}$ is kept constant. Figure 5 shows that the same correction can be
made by scaling $\Gamma^{\prime}$ and keeping the simulation volume constant.
Here $k^{\prime}_{0}$ and $k^{\prime}_{1}$ are chosen as before, but
$\Gamma^{\prime}$ and $N_{C}$ are varied such that $\bar{N}_{C}=4$ for each
simulation. The dimensionless simulated rate for each value of
$\Gamma^{\prime}$ is shown in blue and theoretical rate given by Equation Eq.
57 is shown the red dotted line. We can see that as $\Gamma^{\prime}$
decreases the reaction rate increases and significantly exceeds the
theoretical reaction rate below $\Gamma^{\prime}=1$ ($N_{C}=4$). We also note
that as $\Gamma^{\prime}$, and hence $N_{C}$, become large the simulated rate
falls below the theoretical rate and this is likely a consequence of the
crowding observed earlier.

Figure 4: The dimensionless reaction rate for System Eq. 51 plotted as a
function of $\bar{N}_{C}$. The theoretical reaction rate given in Equation Eq.
57 is plotted as a red dashed line and compared to the simulated rate shown by
the blue points for $\bar{N_{C}}\in\left\\{1,2,3,4,5,7.5,10,15,20\right\\}$.
Figure 5: The dimensionless reaction rate for System Eq. 51 with
$\bar{N}_{C}=4$ plotted for $\Gamma^{\prime}\in\\{0.25,0.5,1,2,4,8\\}$. The
theoretical reaction rate given in Equation Eq. 57 is plotted as a red dashed
line and compared to the simulated rate shown by the blue points. We can see
that as $\Gamma^{\prime}$ and hence the number of independent molecules of $C$
in the simulation is decreased the simulated reaction rate increases. There is
good agreement between the simulated and theoretical reaction rates for
$\Gamma^{\prime}=2$ (the volume contains $8$ molecules of $C$ for this test),
while for values higher than the simulated rate falls below the theoretical
rate due to crowding effects.

##  6 Conclusions

We proposed a modification to Smoluchowski’s model of reaction diffusion
systems that enables us to reproduce non-linear reaction rates characteristic
of enzyme kinetics. While Smoluchowski’s model is unable to correctly
incorporate the fast reactions associated with the formation of enzyme-
substrate complexes, the reaction condition in Equation Eq. 43 allows us to
reproduce the non-linear kinetics of Reaction Eq. 15 without needing to
explicitly simulate these reactions. Our reaction condition differs from the
static bimolecular conditions traditionally used in derivatives of
Smoluchowski’s framework as is size is determined by the relative proximity of
the reactants. Although Condition Eq. 43 is trimolecular, it can be thought of
as governing a bimolecular reaction where the size of the reaction boundary is
determined by the proximity of the third molecule. This view allows
trimolecular conditions of this form to be easily and efficiently incorporated
into time driven simulations of Smoluchowski’s framework. In addition, we have
identified several components of eGFRD that can be easily extended to our
modified framework. Leveraging these ideas we have conducted proof of concept
simulations that demonstrate our framework reproduces the expected non-linear
kinetics. While the theory presented here is limited to systems that can be
reduced to a single trimolecular reaction, in a future publication we hope to
extend our framework to a wider variety of enzymatic systems.

##  Appendix A Evolution of the closest molecule

The scaled probability distribution function,
$\Phi\left(\bm{\eta}_{3},t\right)$, defined originally in Equation
LABEL:eq:Phi_definition,

$\Phi\left(\bm{\eta}_{3},t\right)=\frac{P\left(,\bm{\eta}_{3},t|\bm{\mathcal{H}}_{3}\right)}{c},$
(58)

is proportional to the probability that a particular molecule of $C$ is
associated with $\bm{\eta}_{3}$ at time $t$, given that we know it is the
closest $C$ to the origin. We assume our system is large and consider the
limit, $N_{c}\rightarrow\infty$, so that $\Phi\left(\bm{\eta}_{3},t\right)$ is
given by Equation Eq. 34,

$\Phi\left(\bm{\eta}_{3},t\right)=\phi\left(\bm{\eta}_{3},t\right)\text{exp}\left(-c\int_{V_{3}}\phi\left(\bm{\eta^{\prime}}_{3},t\right)dV_{3}^{\prime}\right),$
(59)

where $\phi\left(\bm{\eta}_{3},t\right)$ is the scaled probability that the
molecule of $C$ is associated with $\bm{\eta}_{3}$ originally defined in
Equation section 3,

$\phi\left(\bm{\eta}_{3},t\right)=\frac{N_{c}P\left(,\bm{\eta}_{3},t\right)}{c}.$
(60)

To derive the governing equation for $\Phi\left(\bm{\eta}_{3},t\right)$ we
consider $\mathcal{L}_{3}\Phi\left(\bm{\eta}_{3},t\right)$ where,

$\mathcal{L}_{3}\equiv\frac{\partial}{\partial
t}-\hat{D}_{3}\hat{\nabla}^{2}_{3},$ (61)

is the diffusion operator on the $3$-dimensional space spanned by
$\bm{\eta}_{3}$, originally defined in Equation Eq. 29. The time derivative of
$\Phi\left(\bm{\eta}_{3},t\right)$ is given by,

$\Phi_{t}=g\phi_{t}-c\phi
g\left(\int_{V_{3}}\phi_{t}\left(\bm{\eta^{\prime}}_{3},t\right)dV_{3}^{\prime}\right),$
(62)

where we have used the $t$ subscript to denote differentiation with respect to
time and defined $g\equiv
g\left(\bm{\eta}_{3},t\right)=\text{exp}\left(-c\int_{V_{3}}\phi\left(\bm{\eta^{\prime}}_{3},t\right)dV_{3}^{\prime}\right)$
for notational convenience. We can also take Laplacian of $\Phi$ with respect
to the coordinates of $\bm{\eta}_{3}$,

$\hat{\nabla}^{2}_{3}\Phi=g\hat{\nabla}^{2}_{3}\phi+2\left(\hat{\nabla}_{3}\phi\right)\cdot\left(\hat{\nabla}_{3}g\right)+\phi\hat{\nabla}^{2}_{3}g,$
(63)

which when combined with Equation Eq. 62 yields,

$\mathcal{L}_{3}\Phi=g\phi_{t}-c\phi
g\left(\int_{V_{3}}\phi_{t}\left(\bm{\eta^{\prime}}_{3},t\right)dV_{3}^{\prime}\right)-\hat{D}_{3}\left[g\hat{\nabla}^{2}_{3}\phi-2\left(\hat{\nabla}_{3}\phi\right)\cdot\left(\hat{\nabla}_{3}g\right)-\phi\hat{\nabla}^{2}_{3}g\right].$
(64)

Recalling that $\mathcal{L}_{3}\phi=0$ from Equation Eq. 29 and applying the
Divergence theorem we find,

$\begin{split}\mathcal{L}_{3}\Phi&=\hat{D}_{3}\left[-c\phi
g\left(\int_{V_{3}}\hat{\nabla}^{2}_{3}\phi\left(\bm{\eta^{\prime}}_{3},t\right)dV_{3}^{\prime}\right)-2\left(\hat{\nabla}_{3}\phi\right)\cdot\left(\hat{\nabla}_{3}g\right)-\phi\hat{\nabla}^{2}_{3}g\right]\\\
&=-c\hat{D}_{3}\phi
g\left(\oint_{S_{3}}\left(\hat{\nabla}_{3}\phi\left(\bm{\eta^{\prime}}_{3},t\right)\right)\cdot\hat{\bm{r}_{3}}dA_{3}^{\prime}\right)-2\hat{D}_{3}\left(\hat{\nabla}_{3}\phi\right)\cdot\left(\hat{\nabla}_{3}g\right)-\hat{D}_{3}\phi\hat{\nabla}^{2}_{3}g,\end{split}$
(65)

where $S_{3}$ is the surface of a sphere of radius $r_{3}=||\bm{\eta}_{3}||$,
$dA_{3}^{\prime}$ is an elemental area on that surface and $\hat{\bm{r}_{3}}$
is the unit outward facing normal vector. The diffusion in the $\bm{\eta}_{3}$
coordinates is isotropic and so $\phi$ is independent of orientation. That is,
$\phi$ only has radial dependence, so we can write,

$\hat{\nabla}_{3}g=\hat{\nabla}_{3}\text{exp}\left(-c\int_{V_{3}}\phi\left(\bm{\eta^{\prime}}_{3},t\right)dV_{3}^{\prime}\right)=-cg\phi\left(\oint_{S_{3}}dA_{3}^{\prime}\right)\hat{\bm{r}_{3}}.$
(66)

By the same reasoning we are able to take the integrand outside of the
integral in Equation Eq. 65 and by substituting in Equation Eq. 66 we obtain,

$\begin{split}\mathcal{L}_{3}\Phi&=\hat{D}_{3}\left[\left(\hat{\nabla}_{3}\phi\right)\cdot\left(\hat{\nabla}_{3}g\right)-\phi\hat{\nabla}^{2}_{3}g\right]\\\
&=-\hat{D}_{3}\hat{\nabla}_{3}\cdot\left(\phi\hat{\nabla}_{3}g\right)\\\
&=\hat{D}_{3}\hat{\nabla}_{3}\cdot\left(cg\phi^{2}\left(\oint_{S_{3}}dA_{3}^{\prime}\right)\hat{\bm{r}_{3}}\right)\\\
&=\hat{D}_{3}\hat{\nabla}_{3}\cdot\left(4\pi
r_{3}^{2}c\phi\Phi\hat{\bm{r}_{3}}\right).\end{split}$ (67)

That is,

$\frac{\partial\Phi\left(\bm{\eta}_{3},t\right)}{\partial
t}=\hat{D}_{3}\hat{\nabla}^{2}_{3}\Phi\left(\bm{\eta}_{3},t\right)+\hat{D}_{3}\hat{\nabla}_{3}\cdot\left(4\pi
r_{3}^{2}c\phi\Phi\left(\bm{\eta}_{3},t\right)\hat{\bm{r}_{3}}\right),$ (68)

as stated in Equation Eq. 35.

## References

  * [1] N. Agmon, Diffusion with back reaction, The Journal of Chemical Physics, 81 (1984), pp. 2811–2817, https://doi.org/10.1063/1.447954.
  * [2] B. Alberts, Molecular biology of the cell, WW Norton & Company, 2017\.
  * [3] S. S. Andrews, Smoldyn publications.  https://www.smoldyn.org/publications.html.
  * [4] S. S. Andrews and D. Bray, Stochastic simulation of chemical reactions with spatial resolution and single molecule detail, Physical Biology, 1 (2004), pp. 137–151, https://doi.org/10.1088/1478-3967/1/3/001.
  * [5] G. E. Briggs and J. B. S. Haldane, A Note on the Kinetics of Enzyme Action, Biochemical Journal, 19 (1925), pp. 338–339, https://doi.org/10.1042/bj0190338.
  * [6] J. Carrera and M. W. Covert, Why build whole-cell models?, Trends in Cell Biology, 25 (2015), pp. 719–722, https://doi.org/10.1016/j.tcb.2015.09.004.
  * [7] W. W. Chen, M. Niepel, and P. K. Sorger, Classic and contemporary approaches to modeling biochemical reactions, Genes & development, 24 (2010), pp. 1861–1875.
  * [8] F. C. Collins and G. E. Kimball, Diffusion-controlled reaction rates, Journal of colloid science, 4 (1949), pp. 425–437.
  * [9] A. Cornish-Bowden, Chapter 1 - basic principles of chemical kinetics, in Fundamentals of Enzyme Kinetics, A. Cornish-Bowden, ed., Butterworth-Heinemann, 1979, pp. 1–15, https://doi.org/10.1016/B978-0-408-10617-7.50006-7.
  * [10] P. Debye, Reaction rates in ionic solutions, Transactions of The Electrochemical Society, 82 (1942), p. 265, https://doi.org/10.1149/1.3071413.
  * [11] M. Doi, Stochastic theory of diffusion-controlled reaction, Journal of Physics A: Mathematical and General, 9 (1976), pp. 1479–1495, https://doi.org/10.1088/0305-4470/9/9/009.
  * [12] 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, https://doi.org/10.1063/1.465983.
  * [13] R. Erban, J. Chapman, and P. Maini, A practical guide to stochastic simulations of reaction-diffusion processes, 2007, https://doi.org/10.48550/ARXIV.0704.1908.
  * [14] R. Erban and S. J. Chapman, Stochastic modelling of reaction–diffusion processes: algorithms for bimolecular reactions, Physical biology, 6 (2009), p. 046001.
  * [15] J. Ëtiles, D. Van Helden, T. Bartol, E. Ëalpeter, and M. Ëalpeter, Miniature end-plate current rise times¡ 100 microseconds from improved dual recordings can be modeled with passive acetylcholine diffusion form a synaptic vesicle, Proc. Natl. Acad. Sci. USA, 93 (1996), pp. 5745–5752.
  * [16] M. Feig and Y. Sugita, Whole-cell models and simulations in molecular detail, Annual review of cell and developmental biology, 35 (2019), p. 191.
  * [17] M. B. Flegg, Smoluchowski reaction kinetics for reactions of any order, SIAM Journal on Applied Mathematics, 76 (2016), pp. 1403–1432.
  * [18] M. A. Hamburg and F. S. Collins, The path to personalized medicine, New England Journal of Medicine, 363 (2010), pp. 301–304.
  * [19] L. H. Hartwell, J. J. Hopfield, S. Leibler, and A. W. Murray, From molecular to modular cell biology, Nature, 402 (1999), pp. C47–C52.
  * [20] S. A. Hollingsworth and R. O. Dror, Molecular dynamics simulation for all, Neuron, 99 (2018), pp. 1129–1143, https://doi.org/10.1016/j.neuron.2018.08.011.
  * [21] E. Honig, G. Roebersen, and P. Wiersema, Effect of hydrodynamic interaction on the coagulation rate of hydrophobic colloids, Journal of Colloid and Interface Science, 36 (1971), pp. 97–109, https://doi.org/10.1016/0021-9797(71)90245-1.
  * [22] S. E. Jackson and J. D. Chester, Personalised cancer medicine, International journal of cancer, 137 (2015), pp. 262–266.
  * [23] J. Karr, J. Sanghvi, D. Macklin, M. Gutschow, J. Jacobs, B. Bolival, N. Assad-Garcia, J. Glass, and M. Covert, A whole-cell computational model predicts phenotype from genotype, Cell, 150 (2012), pp. 389–401, https://doi.org/10.1016/j.cell.2012.05.044.
  * [24] J. R. Karr, J. C. Sanghvi, D. N. Macklin, A. Arora, and M. W. Covert, WholeCellKB: model organism databases for comprehensive whole-cell models, Nucleic Acids Research, 41 (2012), pp. D787–D792, https://doi.org/10.1093/nar/gks1108.
  * [25] J. R. Karr, K. Takahashi, and A. Funahashi, The principles of whole-cell modeling, Current Opinion in Microbiology, 27 (2015), pp. 18–24, https://doi.org/10.1016/j.mib.2015.06.004.
  * [26] H. Kramers, Brownian motion in a field of force and the diffusion model of chemical reactions, Physica, 7 (1940), pp. 284–304, https://doi.org/10.1016/S0031-8914(40)90098-2.
  * [27] J. Lipková, K. C. Zygalakis, S. J. Chapman, and R. Erban, Analysis of brownian dynamics simulations of reversible bimolecular reactions, SIAM Journal on Applied Mathematics, 71 (2011), pp. 714–730.
  * [28] R. Losick and C. Desplan, Stochasticity and cell fate, Science, 320 (2008), pp. 65–68, https://doi.org/10.1126/science.1147888.
  * [29] L. Marucci, M. Barberis, J. Karr, O. Ray, P. R. Race, M. de Souza Andrade, C. Grierson, S. A. Hoffmann, S. Landon, E. Rech, J. Rees-Garbutt, R. Seabrook, W. Shaw, and C. Woods, Computer-aided whole-cell design: Taking a holistic approach by integrating synthetic with systems biology, Frontiers in Bioengineering and Biotechnology, 8 (2020).
  * [30] L. Michaelis, M. L. Menten, et al., Die kinetik der invertinwirkung, Biochem. z, 49 (1913), p. 352.
  * [31] J. D. Murray, Mathematical biology: An introduction, volume 17, 2003\.
  * [32] J. Nielsen, Systems biology of metabolism: a driver for developing personalized and precision medicine, Cell metabolism, 25 (2017), pp. 572–579.
  * [33] O. Purcell, B. Jain, J. R. Karr, M. W. Covert, and T. K. Lu, Towards a whole-cell modeling approach for synthetic biology, Chaos: An Interdisciplinary Journal of Nonlinear Science, 23 (2013), p. 025112, https://doi.org/10.1063/1.4811182.
  * [34] S. A. Rice, Chapter 2 diffusion-controlled reactions in solution, in Diffusion-Limited Reactions, C. Bamford, C. Tipper, and R. Compton, eds., vol. 25 of Comprehensive Chemical Kinetics, Elsevier, 1985, pp. 3–46, https://doi.org/10.1016/S0069-8040(08)70252-8.
  * [35] E. Roberts, A. Magis, J. O. Ortiz, W. Baumeister, and Z. Luthey-Schulten, Noise contributions in an inducible genetic switch: A whole-cell simulation study, PLOS Computational Biology, 7 (2011), pp. 1–21, https://doi.org/10.1371/journal.pcbi.1002010.
  * [36] J. Schöneberg and F. Noé, Readdy-a software for particle-based reaction-diffusion dynamics in crowded cellular environments, PloS one, 8 (2013), p. e74261.
  * [37] L. A. Segel and M. Slemrod, The quasi-steady-state assumption: A case study in perturbation, SIAM Review, 31 (1989), pp. 446–477.
  * [38] M. L. Shuler, S. Leung, and C. C. Dick, A mathematical model for the growth of a single bacterial cell*, Annals of the New York Academy of Sciences, 326 (1979), pp. 35–52, https://doi.org/10.1111/j.1749-6632.1979.tb14150.x.
  * [39] M. v. Smoluchowski, Versuch einer mathematischen theorie der koagulationskinetik kolloider lösungen, Zeitschrift für physikalische Chemie, 92 (1917), pp. 129–168.
  * [40] 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, https://doi.org/10.1063/1.5064867.
  * [41] 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.
  * [42] J. R. Stiles, T. M. Bartol, et al., Monte carlo methods for simulating realistic synaptic microphysiology using mcell, Computational neuroscience: realistic modeling for experimentalists, (2001), pp. 87–127.
  * [43] B. Szigeti, Y. D. Roth, J. A. Sekar, A. P. Goldberg, S. C. Pochiraju, and J. R. Karr, A blueprint for human whole-cell modeling, Current Opinion in Systems Biology, 7 (2018), pp. 8–15.
  * [44] K. Takahashi, S. Tănase-Nicola, and P. R. ten Wolde, Spatio-temporal correlations can drastically change the response of a mapk pathway, Proceedings of the National Academy of Sciences, 107 (2010), pp. 2473–2478, https://doi.org/10.1073/pnas.0906885107.
  * [45] Y. Taniguchi, P. J. Choi, G.-W. Li, H. Chen, M. Babu, J. Hearn, A. Emili, and X. S. Xie, Quantifying e. coli proteome and transcriptome with single-molecule sensitivity in single cells, Science, 329 (2010), pp. 533–538, https://doi.org/10.1126/science.1188308.
  * [46] Z. R. Thornburg, D. M. Bianchi, T. A. Brier, B. R. Gilbert, T. M. Earnest, M. C. Melo, N. Safronova, J. P. Sáenz, A. T. Cook, K. S. Wise, et al., Fundamental behaviors emerge from simulations of a living minimal cell, Cell, 185 (2022), pp. 345–360.
  * [47] M. Tomita, Whole-cell simulation: a grand challenge of the 21st century, Trends in Biotechnology, 19 (2001), pp. 205–210, https://doi.org/10.1016/S0167-7799(01)01636-5.
  * [48] J. S. Van Zon and P. R. Ten Wolde, Green’s-function reaction dynamics: a particle-based approach for simulating biochemical networks in time and space, The Journal of chemical physics, 123 (2005), p. 234910.
  * [49] 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, https://doi.org/10.1103/PhysRevLett.94.128103.
  * [50] A. Vijaykumar, P. G. Bolhuis, and P. R. ten Wolde, Combining molecular dynamics with mesoscopic green’s function reaction dynamics simulations, The Journal of Chemical Physics, 143 (2015), p. 214102, https://doi.org/10.1063/1.4936254.
  * [51] P. G. Wolynes and J. M. Deutch, Slip boundary conditions and the hydrodynamic effect on diffusion controlled reactions, The Journal of Chemical Physics, 65 (1976), pp. 450–454, https://doi.org/10.1063/1.432788.


For the objective part, we need to implement $ n^2 $ unitaries to construct the edge representation through all-to-all strategy. Each unitary requires $ logn$ depth and $log^2 n$ gates to set the flag, $logn$ depth and $logn$ gates to implement the controlled copy, hence overall we need $ n logn $ depth and $n^2 log^2n$ gates. Implementation details can be found in Sec.~\ref{sec:xor} and~\ref{sec:controlled-copy}. The last step which implements the cost requires $ n^2 $ unitaries each with $ logn $ depth and $log^2n$ gates. Note that these unitaries can be implemented in parallel giving $nlogn$ depth. We conclude that the overall depth is $nlogn $, and the number of gates is $ n^2 log^2 n$.

In the LNN architecture, we need to use the procedures described in Theorems~\ref{theorem:reginterlace} and \ref{theorem:all-to-all}. This puts additional cost of swapping the registers which is $logn$, which results in depth $ n log^2 n $.

The number of parameterized gates is 1 for permutation checking, $n^2$ for objective part implementation, and $1$ for Grover Mixer using a single-ancilla decomposition from Lemma 7.11 from~[54], hence $n^2$ in total. For the energy span, observe that
\begin{equation}
\sum_{t=0}^{n-1}\sum_{\substack{i,j=0\\i\neq j}}^{n-1} W_{ij}\delta(b_t, i) \delta(b_{t+1}, j) + A  [\{b_t\}\textrm{ is not a permutation}]  \leq  n\max_{i\neq j}W_{ij} + A
= \order{n}.
\end{equation}

Similarly, as it was in Sec.~\ref{sec:sim-maxkcut}, the cost of the initial state and mixer is negligible.

\section{Superposition of permutations in binary encoding}\label{sec:sup}

Let us recall the original algorithm for generating superposition of permutation in one-hot encoding~[51]. The circuit for generating the superposition of permutations uses $n$ registers, each consisting of $n$ qubits. At each step, one of the registers is fixed so that it complies with the superposition of permutations requirements (the details of the circuit can be found in [51]). The first step fixes the first register, using the second register called the \emph{mask register}. The first step ends with swapping the second and third registers, thus creating space for the \emph{workspace register}.

The next $n-3$ steps fix the consecutive registers as follows. Workspace and mask registers are moved down by $n$ qubits to open up room for the next register to be fixed (so that the new register is between the already fixed registers and the workspace register). Then an appropriate unitary operation is applied so that this new sandwiched register is fixed.

The final step fixes the workspace and mask registers to obtain the desired final state.

In the case of the binary encoding, the following is changed: at the end of each step once a register is fixed, a unary-to-binary encoding circuit is applied as in Fig.~\ref{fig_utob}. As we will show later, encoding change can be accomplished with depth $n$ on LNN architecture. Since the output qubits storing the value in binary encoding are scattered in the register, they are moved to the top of the register leaving $n-⌈logn⌉$ free $|0⟩$ qubits below. Then, the workspace and mask registers need to be moved down only by $⌈logn ⌉$ qubits to create an empty space of $n$ qubits for the next step. Finally, in the end, after the workspace and map registers are fixed, we apply unary-to-binary encoding on both of them and move all the relevant qubits to the top so that the first $n⌈logn⌉$ qubits store the superposition of permutations. Note that since these qubits need to be moved through $n$ qubits, the overall depth is still $ n$. Overall, we need $3n + (n-3) ⌈logn ⌉= nlogn $ qubits, of which the last $3(n-logn)$ qubits are set to $|0⟩$ at the end, and those are the ancilla qubits required for the procedure.

\begin{figure}\centering
\includegraphics{unary_to_binary.pdf}
\caption{\label{fig_utob}Quantum circuit transforming the number in unary encoding into equivalent binary encoding~[44]. $\ket{c_0}$, $\ket{c_1}$, $\ket{c_2}$ are the output bits, the qubits without labels in the output are set to $\ket{0}$. The parts I, II, and III are applied repeatedly for consecutive bits in binary encoding.}
\end{figure}

Let us now show that the procedure presented in Fig.~\ref{fig_utob} can be implemented in $n$ depth even on the LNN network. Note that the procedure consists of several similar steps applied on smaller and smaller circuits, a single step distinguished with a blue block. Each such step consists of three parts: a single CNOT (yellow) which can be implemented by shifting the target qubit down to the one-by-last qubit, applying CNOT, and then undoing the shift. Then there is a sequence of CNOTs represented by the green block, in which the target qubit should be moved upwards after each CNOT and should be moved back to the bottom at the end. Note that in both cases we have $≤CN$ gates, where $C$ is a constant (in the next iterations, it will be $CN/2$, $CN/4$, etc.).

The most complicated procedure is the one distinguished by the red block, where we have a sequence of controlled-SWAPs. In there, we start by interlacing all the qubits. Note that after the interlacing, all qubits on which the swaps need to be applied are neighboring. So we only need to move control qubits through the whole register in a similar fashion, as it was done in the circuit described in Sec.~\ref{sec:controlled-copy}. Then the control qubit should be moved to its original position, and the interlacing should be undone. Note that even after these updates, the circuit can be implemented in $ n$ depth.

Since all the adjustments to the encoding change circuit do not increase the depth at each step, the total depth (and the number of gates) remains the same as in~[51], namely $n^2$ depth and $ n^3$ gates on LNN.

%\section{Numerical analysis -- version with Fuchs-QAOA}\label{sec:appendix-fuchs-numerics}
%Below we present a version of numerical experiment done for MAX-$K$-CUT with Fuchs-QAOA:
%%	\includegraphics{}

%Suppose that instead of starting with the product of $n$ quantum states $\frac{1}{\sqrt n}
%\sum_{i=0}^{n-1} \ket i$, we start with the equal superposition of permutations $\frac{1}{\sqrt n!} \sum_{x \in B_n} \ket{x}$ where each city is encoded in binary. Here, $B_n$ is the set of strings obtained by concatenating $n$ strings each encoding a different city in binary with length $\lceil \log n\rceil$. With this choice of the initial state, the Hamiltonian takes the form
%\sum_{t=0}^{n-1}\sum_{\substack{i,j=0\\i\neq j}}^{n-1} W_{ij}\delta(b_t, i) \delta(b_{t+1}, j).
%Overall, all the cost measures are equal to that of the GM-QAOA, but the number of required qubits is $\order{n \log n}$ instead of $\order{n^2}$. Hence, we can simulate the exact behavior of GM-QAOA using a much less amount of qubit resources. Further improvement of the cost measures depends on more efficient initial state preparation.

\end{document}