Enzyme kinetics simulation at the scale of individual particles

TL;DR

This paper develops a new proximity-based reaction condition for particle simulations that accurately models enzyme kinetics involving multiple timescales without explicitly simulating fast reactions.

Contribution

It introduces a novel reaction condition that captures short-timescale enzyme reactions in particle-based models, improving their accuracy.

Findings

Successfully reproduces non-linear enzyme reaction rates

Validates the new reaction condition through particle-based simulations

Enhances modeling accuracy for enzyme kinetics in reaction-diffusion systems

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 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 accurate modelling of enzyme kinetics by many particle-based simulations of reaction-diffusion systems as 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 and demonstrate that non-linear reaction rates typical of enzyme kinetics can be reproduced without needing to explicitly simulate reactions on the short timescale.

Equations

$$K=\frac{k}{V}=\frac{4\pi \sigma {\hat{D}}_2}{V}.$$

$${\mathbf{\eta }}_1$$

$${\mathbf{\eta }}_i$$

$${\bar{\mathbf{x}}}_i$$

$${\hat{D}}_i=\{\begin{array}{l}{\bar{D}}_N\text{ when }i=1,\\ D_i+{\bar{D}}_{i-1}\text{ when }i>1,\end{array}$$

$${\bar{D}}_j=\frac{1}{{\sum }_{i=1}^j{D}_i^{-1}}.$$

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

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

$$P\left(\mathbf{\eta }\in \partial {\Omega }_R,t\right)=0.$$

$$\underset{\mathbf{\eta }\to \infty }{\lim }P\left(\mathbf{\eta },t\right)=P_{\infty }=\frac{1}{V^{N-1}}.$$

$${\Omega }_R=\left\{\mathbf{\eta }={\mathbf{\eta }}_2:r_2\le \sigma \right\},\text{where }{r}_2=||{\mathbf{\eta }}_2||.$$

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

$$K=\frac{4\pi \sigma {\hat{D}}_2}{V}\left(1+\frac{\sigma }{\sqrt{4{\hat{D}}_{2}t}}\right).$$

$$\underset{t\to \infty }{\lim }P\left(r_2,t\right)=\frac{1}{V}\left(1-\frac{\sigma }{r_2}\right).$$

$$\text{ce}A+B<=>[k_1][k_{-1}]X\text{and}\text{ce}X+C->[k_2]P+C,$$

$$\frac{da}{dt}$$

$$\frac{dx}{dt}$$

$$k_{-1}=\frac{{\bar{k}}_{-1}}{\epsilon }\text{and}{c}_0=\frac{{\bar{c}}_0}{\epsilon },$$

$$\bar{a}=\frac{a}{a_0},\bar{b}=\frac{b}{b_0},\bar{c}=\frac{c}{c_0},\bar{x}=\frac{x}{a_0},\text{and}T=k_1{b}_{0}t,$$

$$\epsilon \frac{d\bar{a}}{dT}$$

$$\epsilon \frac{d\bar{x}}{dT}$$

$$\mu =\frac{{\bar{k}}_{-1}}{k_1{b}_0},\text{and}\nu =\frac{k_2{\bar{c}}_0}{k_1{b}_0}.$$

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

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

$$\Gamma =\frac{k_{-1}}{k_2},$$

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

$$k_3\left(c\right)=\frac{k_1}{\Gamma +c}.$$

$$\left[\frac{\partial }{\partial t}-{\hat{D}}_2{\hat{\nabla }}_2^2\right]P_2\left({\mathbf{\eta }}_2,t\right)$$

$$\left[\frac{\partial }{\partial t}-{\hat{D}}_3{\hat{\nabla }}_3^2\right]P_3\left({\mathbf{\eta }}_3,t\right)$$

$$P\left(,{\mathcal{H}}_3,t\right)=\frac{1}{N_C}.$$

$$P\left(,{\mathcal{H}}_3,t|{\mathbf{\eta }}_3\right)={\left[1-{\int }_{V_3}P\left(,{{\mathbf{\eta }}^{\prime }}_3,t\right)d{V}_3^{\prime }\right]}^{N_C-1},$$

$$P\left(,{\mathbf{\eta }}_3,t|{\mathcal{H}}_3\right)=N_{C}P\left(,{\mathbf{\eta }}_3,t\right){\left[1-{\int }_{V_3}P\left(,{{\mathbf{\eta }}^{\prime }}_3,t\right)d{V}_3^{\prime }\right]}^{N_C-1}.$$

$$\Phi \left({\mathbf{\eta }}_3,t\right)=\varphi \left({\mathbf{\eta }}_3,t\right){\left[1-\frac{c}{N_C}{\int }_{V_3}\varphi \left({{\mathbf{\eta }}^{\prime }}_3,t\right)d{V}_3^{\prime }\right]}^{N_C-1}.$$

Full text

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

1 Introduction

Whole cell models promise to revolutionise systems biology. The capability to simulate the integrated function of every gene and molecule in a cell would assist clinicians in individualising 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 neighbourhood 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 [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 idealised 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 Section 2.

Smoluchowski’s original reaction condition has been criticised 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 reaction condition 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 summarises 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 reaction condition.

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

[TABLE]

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 (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 initialised 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 reaction condition 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 initialisation. 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 characterised 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 model of bimolecular reactions. 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 would have to wait a prohibitively long time before the substrate concentration had 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 modeller. 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.

Motivated by these issues, we propose a modification to Smoluchowski’s original reaction condition that allows us to reproduce non-linear reaction rates in trimolecular systems that are reminiscent of those in observed in Michaelis-Menten kinetics. Our new reaction condition is informed by the singular perturbation analysis of a trimolecular enzymatic system and incorporates the influence of the complex formation without requiring events on $t_{c}$ to be explicitly modelled. This completely circumvents all the issues posed 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 Generalised Smoluchowski theory

Smoluchowski only considered bimolecular reactions, but his work can be generalised to incorporate reactions of any order, as demonstrated 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}$ 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

[TABLE]

is the ‘centre of diffusion’ of the first $i$ molecules (analogous to the centre 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 centre of diffusion of the $N$ molecules and $\bm{\eta}_{i}$ for $i>1$ describes the separation between $\bm{x}_{i}$ and the centre 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

[TABLE]

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

[TABLE]

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

[TABLE]

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 normalisation,

[TABLE]

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,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

Before the system reaches this steady state the reaction rate can be arbitrarily high, diverging to infinity as we approach $t=0$. This transient behaviour is an artefact of the initial condition, where two reactants have a non-zero probability of being initialised such that $r_{2}\leq\sigma$ (or even an elevated chance of being initialised with $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 (1).

3 Generalised reaction conditions

To illustrate the aforementioned issues with Smoluchowski’s reaction condition 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 summarised succinctly in chemical shorthand as

[TABLE]

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 (15) yields a system of five ODEs that can be reduced to just two equations,

[TABLE]

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 Equations (16) and (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,

[TABLE]

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

To proceed we define the dimensionless variables

[TABLE]

where $a_{0}$ and $b_{0}$ are the initial concentrations of $A$ and $B$ respectively. Applying our change of variables to Equations (16) and (17) gives

[TABLE]

where we have defined the dimensionless parameters

[TABLE]

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

[TABLE]

Substituting Equation (23) into Equation (20) and redimensionalising gives to leading order in $\varepsilon$,

[TABLE]

where

[TABLE]

is analogous to the Michaelis constant for Reaction (15).

The preceding analysis shows that on the long timescale Reaction (15) can be reduced to a single trimolecular reaction

[TABLE]

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

[TABLE]

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.

We have shown that under the pseudo-steady-state approximation Reaction (15) is equivalent to Reaction (26) on the long timescale. The disparate timescales required to make this approximation valid, also preclude accurate particle-based simulation of Reaction (15) by any methods that rely on Smoluchowski’s bimolecular reaction condition given in Equation (11). For instance, the fast bimolecular reaction between $X$ and $C$ cannot be accurately simulated on the long timescale using this condition as shown in Fig. 1(a). To avoid this issue, we seek instead construct a particle-based simulation of Reaction (26), which requires the development of a new trimolecular reaction condition that reproduces the non-linear reaction rate in Equation (27) as shown in Fig. 1(b).

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

[TABLE]

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 (where each state consists 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 (28). On the manifold, states diffuse independently in the $\bm{\eta}_{3}$ space in accordance with Equation (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 the $\bm{\eta}_{3}$ coordinates 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 it 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

[TABLE]

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}||$ centred on the origin, i.e.

[TABLE]

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

[TABLE]

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

[TABLE]

which when substituted into Equation (32) give,

[TABLE]

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

[TABLE]

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

[TABLE]

where $\hat{\bm{r}_{3}}$ is the unit outward facing normal vector of a sphere of radius $r_{3}$; a derivation of this result can be found 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 the $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)$,

[TABLE]

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

[TABLE]

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 Equation (36) becomes negligible as was the case for the trimolecular generalisation of Smoluchowski 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 non-linear reaction kinetics in the concentration of $C$ we propose a long thin boundary $\partial\Omega_{NL}$ where

[TABLE]

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 the normal to the boundary is $\mathbf{\hat{n}}=\bm{\eta}_{2}$ and the steady-state solution to Equation (37) with $P\left(,\bm{\eta},t|\bm{\mathcal{H}}_{3}\right)=0$ on $\partial\Omega_{NL}$ where $\Omega_{NL}$ is given by Equation (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

[TABLE]

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

[TABLE]

Finding the flux in the $\mathbf{\hat{n}}$ direction over $\partial\Omega_{NL}$ 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,

[TABLE]

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$,

[TABLE]

where $\mathcal{L}$ is the Laplace transform and $F(4\pi r_{3}^{3}/3)=f(r_{3})$. For Reaction (26) the reaction rate is, $K(c)=ck_{3}(c)$ where $k_{3}(c)$ is given by Equation (27). In this case 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),

[TABLE]

Finally, substituting this result into Equation (38) we obtain a reaction boundary $\partial\Omega_{NL}$,

[TABLE]

that reproduces the non-linear kinetics of Reaction (26) and thus can be used to construct a particle-based simulation of Reaction (15).

We will focus on investigating the kinetics that result from the reaction boundary in Equation (44) as to our knowledge this constitutes the first example of a proximity-based reaction condition capable of directly reproducing non-linear kinetics. Before proceeding however, it is worth considering what kinds of kinetics are attainable via our method. In the current presentation, our method requires that the original reaction network can be reduced to an equivalent trimolecular reaction in the same form as Reaction (26). Treating bimolecular reactions with non-linear reaction rates - such as the original Michaelis-Menten system [30] - is more difficult, since removing the third reactant would also reduce the spatial degrees of freedom that can be utilised when designing the reaction boundary. In principle, the method could be generalised to higher order reactions, i.e. those involving four or more reactants, and since the additional reactants introduce new spatial degrees of freedom it is conceivable that this might allow more exotic non-linear rates to be simulated. An additional restriction can be inferred directly from Equation (42), namely that the inverse Laplace transform of the reaction rate must exist. Moreover, the resulting function $f$ must be non-negative for all values of $r_{3}$ since it denotes a distance.

4 Simulation of non-linear kinetics

To construct a particle-based simulation of trimolecular reactions like Reaction (26) we can make use of the methods typically employed by simulations of Smoluchowski’s framework. These simulations are often divided into time-driven (TD), and event-driven (ED) approaches. 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 Smoldyn software package, which has been used widely in literature [3]. Smoldyn 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]

[TABLE]

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 (3) and the reaction boundary in Equation (44) is tested to determine if any reaction events took place during the last timestep. Equation (45) 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 avoided by making the timestep 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 the size of the reaction boundary in Equation (44) is a function of $r_{3}$. As a consequence, we need to compute a series of corrections so that the appropriate value may be retrieved for any value of $r_{3}$.

Reaction (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

[TABLE]

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

[TABLE]

A molecule of $AB$ is initialised at the centre of diffusion of the reacting $A$ and $B$ molecules, which is defined by Equation (4) as

[TABLE]

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 (3),

[TABLE]

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

Reaction (47) is deceptively similar to Reaction (15) with the complexes $AB$ and $X$ appearing interchangeable. However, $AB$ is not a chemical species. Instead a molecule of $AB$ is constructed within the simulation purely as a convenient way to track any pair of $A$ and $B$ molecules that are close enough to react with a molecule of $C$. Pairs of $A$ and $B$ molecules that are not close enough to form one of these fictitious 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 $A$ and $B$ molecules are reintroduced at the positions

[TABLE]

respectively. When viewed in this way, Reaction (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 $AB$ complex to the closest molecule of $C$. The position of the complex is the centre of diffusion of the $A$ and $B$ molecules allowing for simple calculation of the separation to each molecule of $C$ using Equation (3),

[TABLE]

where the subscript $3$ is associated with a particular molecule of $C$. Once the minimum value of $\bm{\eta}_{3}$ has been found, the reaction boundary in Equation (44) can be checked to determine if $A$ and $B$ are close enough for a reaction to occur, as shown in Fig. 2. 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(r_{3})$ can be retrieved efficiently.

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 Fig. 3. 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 a molecule in TD mode may approach a protective domain before the domain’s escape time, it is possible for a reaction to occur between the TD molecule and the molecule within the domain. 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 the fictitious $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 ‘centre 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 and due to this we did not implement protective domains for $AB$ molecules during the simulation.

5 Numerical results

In this section, we present the results of a series of particle-based simulations of Reaction (26) conducted using an implementation of the simulation framework discussed in Section 4. The first test considers a well-mixed trimolecular system and is designed to validate our implementation. While the second test explores whether the reaction boundary in Equation (44) reproduces the non-linear behaviour described in Equation (24) as $c$ is varied.

5.1 Well-mixed trimolecular test system

To verify our simulation method we consider the trimolecular system

[TABLE]

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 (27). The population of $A$ molecules in the system is governed by the Law of Mass Action and when expressed in terms of molecular populations rather than concentrations may be written

[TABLE]

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

[TABLE]

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 the dimensionless variables

[TABLE]

Equation (53) becomes

[TABLE]

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],

[TABLE]

We choose the scaled rate constants $k_{0}^{\prime}=1$ and $k_{1}^{\prime}=0.2$ such that $\bar{N}_{A}=0.2N_{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 $3\times 10^{4}$ times and the population of $A$ molecules was sampled at the conclusion of each simulation.

In Fig. 4, the sampled distribution of $N_{A}$ (blue bars) is compared with the expected Poisson distribution (red dots). Neither the mean, $6.00\pm 0.01$, nor the variance, $5.98\pm 0.05$, of the simulated distribution deviate significantly from the initial Poisson distribution. This indicates that the kinetics observed in our particle-based simulation closely resemble the theoretical kinetics for the choices of reaction parameters considered. However, if the reaction parameters are not chosen carefully, one may expect to observe some errors arising from the fact that the boundary in Equation (44) only approximates the actual reaction boundary required to reproduce the kinetics of Reaction (26). Moreover, the reaction boundary is further altered during the simulation in an attempt to correct for the errors introduced by the use of finite timesteps.

5.2 Non-linear kinetics

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

[TABLE]

To demonstrate that our framework reproduces this expected non-linear behaviour, we use the same methodology described in 5.1, to simulate System (52) with $\Gamma=2$, and $k_{1}=0.3$, while $k_{0}$ and $V$ are chosen so that $\bar{N}_{A}=0.3N_{A}$. These reaction parameters were selected since they yield convenient dimensional reaction rates for several of the $\bar{N}_{C}$ values considered. Although, any set that did not violate the assumptions used to derive the reaction boundary in Equation (44) - namely that the boundary is thin in $r_{2}$ and long in $r_{3}$ - could be considered.

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 reaction rate is calculated by taking the inverse of the mean of the steady-state distribution obtained for $\bar{N}_{A}$ for each value of $\bar{N}_{C}$. These simulated rates are shown by the blue points in Fig. 5, where they are compared against a plot of Equation (58) shown by the red dashed line. The simulated reaction rate agrees with the theoretical rate for all the values of $\bar{N}_{C}$ considered, and can be seen to reproduce the characteristic non-linear behaviour predicted by Equation (58). However, there is some indication that the simulated reaction rate has a lower horizontal asymptote than predicted, and it is possible that the rates would eventually diverge if $\bar{N}_{C}$ was increased further. 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. In addition, it should be noted that initially we considered $k_{0}=1$ and $V=1$ for all of our simulations, but for these parameters we found that the simulated rate significantly exceeded the theoretical rate for $\bar{N}_{C}<5$. 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 the reaction boundary in Equation (44) is only valid in the limit that $V\to\infty$ so that an effectively infinite number of independent $C$ molecules are present regardless of the concentration of $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. Therefore, if for example $N_{C}=2$, then the simulation would only contain two independent molecules of $C$. This can be corrected by increasing $N_{C}$, but necessitates an increase in $\Gamma^{\prime}$ if $\bar{N}_{C}$ is to be kept constant. Since $\Gamma$ controls the shape of the reaction boundary, we leave it unchanged at $\Gamma=2$ and increase $\Gamma^{\prime}$ by expanding the simulation volume instead. Similarly, we retain $k_{1}=0.3$ and alter $k_{0}$ - which does not impact the shape of the reaction boundary - so that $\bar{N}_{A}=0.3N_{A}$ for all values of $\bar{N}_{C}$ considered. In this way we were able to obtain the simulated reaction rates for $\bar{N}_{C}\in\left\{1,2,3,4\right\}$ that are shown in Fig. 5. Fig. 6 demonstrates the corrections obtained via this method for $\bar{N}_{C}=1$. Here we plot the simulated reaction rate as a function of $\Gamma^{\prime}$, which is altered by changing $V$ and keeping $\Gamma$ fixed. We then set $N_{C}=\Gamma^{\prime}$ so that $\bar{N}_{C}$ remains $1$. The value of $\bar{K}$ is calculated from $5\times 10^{3}$ simulations for each value in $\Gamma^{\prime}\in\left\{1,2,4,8,16\right\}$ and is shown in blue, while the theoretical rate given by Equation (58) is shown by the red dotted line. We can see that as $\Gamma^{\prime}$ - and hence the number of independent $C$ molecules - decreases, the simulated reaction rate increases and exceeds the theoretical reaction rate for $\Gamma^{\prime}\leq 4$ ($N_{C}\leq 4$). For $\Gamma^{\prime}=8$ and $\Gamma^{\prime}=16$ the simulated reaction rate agrees with the theoretical rate, indicating that the simulation contains sufficiently many independent molecules of $C$ to be a good approximation of the effectively infinite population considered in Section 3. Finally, we note that the simulated reaction rate does appear to decrease slightly when $\Gamma^{\prime}$ - and hence $N_{C}$ - is increased from $8$ to $16$ and this is likely a consequence of the crowding observed earlier.

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 bimolecular reaction condition is unable to correctly incorporate the fast reactions associated with the formation of enzyme-substrate complexes, the reaction boundary in Equation (44) allows us to reproduce the non-linear kinetics of Reaction (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 the size of the boundary is determined by the relative proximity of the reactants. Although Reaction (15) is trimolecular, we have shown that it can be thought of as a bimolecular reaction between $A$ and $B$ where the size of the reaction boundary is determined by the proximity of the third molecule $C$. This view allows trimolecular reactions 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 implemented in the presence of our generalised reaction boundaries. Leveraging these ideas we have conducted proof of concept simulations that demonstrate our reaction boundary reproduces the expected non-linear kinetics. While the theory presented here is limited to systems that can be reduced to a single trimolecular reaction similar in form to Reaction (26), 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 density function, $\Phi\left(\bm{\eta}_{3},t\right)$, defined originally in Equation (LABEL:eq:Phi_definition),

[TABLE]

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}\to\infty$, so that $\Phi\left(\bm{\eta}_{3},t\right)$ is given by Equation (34),

[TABLE]

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 (3),

[TABLE]

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

[TABLE]

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

[TABLE]

where we have used the $t$ subscript to denote differentiation with respect to time and defined

[TABLE]

for notational convenience. We can also take Laplacian of $\Phi$ with respect to the coordinates of $\bm{\eta}_{3}$

[TABLE]

which when combined with Equation (A.5) yields

[TABLE]

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

[TABLE]

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 have

[TABLE]

By the same reasoning we are able to take the integrand outside of the integral in Equation (A.9) and by substituting in Equation (A.10) we obtain

[TABLE]

That is,

[TABLE]

as stated in Equation (35).