The reaction-diffusion master equation (RDME) is a standard modelling approach for understanding stochastic and spatial chemical kinetics. An inherent assumption is that molecules are point-like. Here, we introduce the excluded volume reaction-diffusion master equation (vRDME) which takes into account volume exclusion effects on stochastic kinetics due to a finite molecular radius. We obtain an exact closed form solution of the RDME and of the vRDME for a general chemical system in equilibrium conditions. The difference between the two solutions increases with the ratio of molecular diameter to the compartment length scale. We show that an increase in the fraction of excluded space can (i) lead to deviations from the classical inverse square root law for the noise-strength, (ii) flip the skewness of the probability distribution from right to left-skewed, (iii) shift the equilibrium of bimolecular reactions so that more product molecules are formed, and (iv) strongly modulate the Fano factors and coefficients of variation. These volume exclusion effects are found to be particularly pronounced for chemical species not involved in chemical conservation laws. Finally, we show that statistics obtained using the vRDME are in good agreement with those obtained from Brownian dynamics with excluded volume interactions.

## I. INTRODUCTION

A large number of studies have investigated the properties of noisy chemical dynamics (for recent reviews see, for example, Refs. 1 and 2). The importance of the topic stems from an increasing interest in understanding the dynamics of chemical systems with small numbers of molecules for one or more species, for which stochastic effects are important. A natural example of such chemical systems is biochemical pathways inside cells;^{3} artificial examples include reactions occurring inside nano-spaces such as nano-reactors^{4} and carbon nanotubes.^{5}

The approaches at the heart of these studies include Brownian dynamics (BD)^{6,7} and the reaction-diffusion master equation (RDME)^{8,9} and its non-spatial counterpart, the chemical master equation (CME).^{1} Brownian dynamics typically models point or hard spherical molecules which diffuse and interact with each other via chemical and steric interactions. The RDME provides an approximate spatially discretised version of Brownian dynamics, whereby space is divided into small volume elements (voxels), reactions occur between point molecules inside each voxel, and diffusion of molecules is simulated by “hopping” of molecules between neighbouring voxels. The CME is a non-spatial approximation of the RDME, valid in the limit of well-mixed dynamics throughout the whole compartment. While Brownian dynamics is clearly the most realistic, the RDME and CME are far superior in terms of computational efficiency and have enabled the simulation of complex biochemical systems via the stochastic simulation algorithm (SSA) and its variants.^{1,10} Another advantage of master equations is that in many cases, they either can be solved exactly (see, for example, Refs. 11–14) or else their solution computed by means of an approximative method such as moment-closure approximation^{15–17} or the system-size expansion^{18–21} leading to insight which cannot be easily obtained by tediously long simulations using Brownian dynamics.

Nevertheless, a convincing argument can be made that the assumption of point molecules by the RDME and CME is highly unrealistic, given that several experimental studies^{22–24} have suggested that volume exclusion effects due to molecular crowding strongly modulate intracellular chemical equilibria and even play an important role in the regulation of gene expression rates.^{25} Brownian dynamics does not necessarily ignore such volume exclusion effects but is not an ideal simulation tool due to its heavy computational demand, not to mention its analytical impenetrability. A considerable number of studies have ignored chemical reaction kinetics and focused on understanding the diffusion of a tracer molecule in a sea of inert hard sphere molecules.^{26–29} A few studies have, in contrast, sought to understand the effect of crowding on the stochastic chemical properties of very simple chemical systems in the reaction-limited regime, by renormalising the propensities of the CME to account for volume exclusion effects.^{30,31} However to-date no general conclusions have been made, to our knowledge, about the impact of volume exclusion on the statistics of intrinsic noise in chemical systems. In other words, we would like to obtain insight into how the predictions of the RDME and CME for the distributions of molecule numbers of a general chemical system are modified, if interacting molecules are modelled as hard particles with a finite radius.

In this paper, we take a step in this direction. We assume that all the molecules in a general chemical system are roughly of equal molecular size and devise a version of the RDME (the excluded volume reaction-difusion master equation (vRDME)) which models reactions between such particles. Of course the assumption of a population of molecules with equal sizes is rough, however as we shall see, it enables us to carry analytical calculations and to get a general idea of the impact of volume exclusion on the statistics of intrinsic noise. The paper is organised as follows. In Section II, we discuss in detail the RDME and the vRDME, and their non-spatial counterparts, the CME and vCME, pointing out their crucial differences. In Section III, we use these master equations to derive exact closed-form expressions for the local and global distributions of molecule numbers in the presence and absence of volume exclusion effects. The relationship between the rate constants of the volume excluded and dilute approaches is discussed in Section IV. Next, we use the results of Sections III and IV to explore the stochastic properties of chemical systems with no chemical conservation laws (Section V), with chemical conservation laws of a special type (Section VI) and with chemical conservation laws of a more general type (Section VII). The validity of the vRDME as an accurate approximation to a spatially continuous microscopic description is explored in Section VIII. We finally conclude by a summary and discussion in Section IX.

## II. THE CME, RDME, vRDME, AND vCME

In this section, we concisely describe the four mathematical frameworks used in this article: the CME, the RDME, and modified versions of these two, which take into account volume exclusion effects, and which we call the vCME and vRDME, respectively. To clarify the differences between the four mathematical frameworks, we use the example of a simple reversible dimerisation whereby two molecules of a monomer (species A) diffuse and eventually bind to form a single molecule of the dimer (species B) and which at a later time dissociates back into the constituent monomers.

**The CME** describes the stochastic time evolution of the molecule numbers of each chemical species in a well-mixed compartment. A major simplifying assumption is that the molecules are point particles. For the dimerisation reaction, the CME models the chemical process $A+A \u21cc k 1 k 0 B$, where *k*_{0} and *k*_{1} are the rate constants for the forward and backward reactions.

**The RDME** is the spatial counterpart of the CME. The compartment is divided into *N* subvolumes called voxels, each well-mixed (well-mixing is not assumed throughout the whole volume, only locally). The RDME describes the stochastic time evolution of the molecule numbers of each chemical species in each voxel, with the assumption that the particles are point-like. For the dimerisation reaction, the RDME models the processes,

where *A _{i}* denotes species A in voxel

*i*,

*B*denotes species B in voxel

_{i}*i*, and the notation

*Ne*(

*i*) stands for the set of voxels which neighbour voxel

*i*. The parameter

*k*has units of inverse time and is proportional to the diffusion coefficient

_{D}*D*of the species; specifically

*k*=

_{D}*D*/Δ

*x*

^{2}where Δ

*x*is the side length of a voxel. The first reaction corresponds to the dimerisation reaction taking place in voxel

*i*, while the second and third reactions model the diffusion of the monomer and the dimer between neighbouring boxes

*i*and

*j*with rate

*k*. The RDME model with 4 voxels is schematically represented in Fig. 1(a). The particles are empty to underline that they occupy no volume, and the grid corresponds to the voxels. The relationship between the RDME and CME will be clarified further in Sec. III.

_{D}**The vRDME** is a modified version of the RDME, which we introduce in this paper as a means to take into account volume exclusion effects. In the vRDME, molecules are assumed to have roughly the same diameter and the voxel size is fixed to this length scale (unlike the RDME where the voxel size is arbitrary). A volume exclusion rule is implemented such that each voxel can accommodate at most one chemical molecule. An “empty space species” is defined whose molecule numbers reflect whether a voxel is empty or occupied. The volume exclusion rule is then implemented via the interaction of the empty space species and a chemical species. Bimolecular reactions involve the interaction of two chemical molecules in neighbouring voxels. For the dimerisation reaction, the vRDME models the processes,

where *E _{i}* denotes an “empty space molecule” in voxel

*i*(the molecule number of species

*E*takes a value of zero if voxel

_{i}*i*is occupied and one if it is empty). The first process models the chemical reaction between two

*A*particles in neighbouring voxels and the other two processes model the diffusion of molecules between neighbouring voxels. Note that because we can interchange the indices

*i*and

*j*, the chemical reaction between two

*A*particles in neighbouring voxels

*i*and

*j*leads to either a

*B*molecule in voxel

*i*or a

*B*molecule in voxel

*j*. The reaction rates have a tilde to denote that these quantities are different than the rates used for the RDME (see later for the relationship between the rate constants of the RDME and the vRDME). The vRDME model with

*N*= 36 voxels is illustrated schematically in Fig. 1(b).

We note that the microscopic stochastic processes modelled by the vRDME have been previously simulated by means of Monte Carlo simulations on a two dimensional lattice, specifically applied to understanding diffusion-limited kinetics in crowded environments.^{34–36} As well, the vRDME is a special case of a class of stochastic population models based on “patch dynamics,” a framework developed by McKane and Newman in the context of ecological systems.^{32} Specifically the vRDME corresponds to one of two types of spatial patch models, the case called “micro-patch” where each patch (each voxel in our terminology) can hold at most one individual. The bulk of studies to date have however focused on the “mesoscopic-patch” approach whereby each patch can hold at most a number *N* of individuals where *N* is typically a number much greater than one, and in which one assumes well-mixing and reactions occurring inside each patch, rather than between neighbouring patches (see, for example, Ref. 33).

**The vCME** is the non-spatial counterpart of the vRDME. The vCME is to the vRDME, what the CME is to the RDME. Hence, the vCME is basically the CME but with two additional properties: (i) besides tracking the total number of molecules of each chemical species in the compartment, it also tracks the total number of empty space molecules in the compartment (this new non-chemical species is denoted as E); (ii) a global exclusion volume rule is imposed, namely, that the total number of molecules of all species (chemical species and the empty space species) adds to a time-independent constant *N* (which corresponds to the number of voxels in the vRDME). For the dimerisation reaction, the vCME models the processes: $A+A \u21cc k \u0303 1 k \u0303 0 B+E$. Note the difference between the CME and vCME; the rate constants are also not the same, hence the tildes.The relationship between the vRDME and vCME will be clarified further in Sec. III.

We have in this section introduced the various mathematical frameworks by means of a simple chemical reaction system but these are applicable to more general systems of chemical reactions.

## III. EXACT SOLUTION OF THE CME, RDME, vRDME, AND vCME IN EQUILIBRIUM CONDITIONS

We will here focus on reversible chemical systems in equilibrium conditions, i.e., those in detailed balance.^{37} The reason for this restriction is that as we shall see, it enables us to write an explicit exact solution of the CME, RDME, vRDME, and vCME, which will be crucial in Secs. III A–III D to understand the differences between them, i.e, the impact of molecular crowding on the stochastic dynamics of biochemical systems at the local (voxel) and global (compartment) levels. We shall start by summarising a result by van Kampen for the CME, which we will subsequently extend to the other three frameworks.

### A. Global distribution of molecule numbers assuming point particle interactions

Consider a well-mixed reversible system of *M* chemical species interacting via *R* chemical reactions, where the *j*^{th} reaction has the form

where *X _{i}* denotes the

*i*th chemical species. Here,

*k*and $ k j \u2032 $ are the rate constants for the forward and reverse reactions, respectively, and

_{j}*r*−

_{ij}*s*is the change in the number of molecules of species

_{ij}*X*when reaction

_{i}*j*occurs. We consider the system to be confined in a compartment of volume Ω. The set of deterministic equilibrium constants

^{39}characterising this mass-action system are

where *ϕ _{i}* is the deterministic concentration of species

*X*(as given by the conventional non-spatial rate equations). Furthermore, we assume that the system has a number,

_{i}*S*, of chemical conservation laws of the form

where $ n \u2192 = { n 1 , n 2 , \u2026 , n M } $ describes the number of molecules of each chemical species in the compartment of volume Ω and the *K _{j}*’s are time-independent constants set by the initial conditions and stoichiometry of the reaction system. Now the time-evolution of the global (whole compartment) distribution of molecule numbers assuming point particles and well-mixed conditions is given by the CME. Assuming mass-action kinetics, van Kampen showed that the exact equilibrium solution of the CME for system (3) is given by

^{40}

where *C* is a normalisation constant, *δ*(., .) is a Kronecker delta, and $P ( n \u2192 ) $ is the probability that the system is in state $ n \u2192 $ in equilibrium. Hence, the equilibrium solution is a product of Poisson distributions constrained by the chemical conservation laws.

### B. Local distribution of molecule numbers assuming point particle interactions

The result is also easily extensible to the RDME. The latter is a master equation describing the time-evolution of the probability that the system is in state $ { n 1 1 , \u2026 , n M 1 , \u2026 , n 1 N , \u2026 , n M N } $, where $ n i j $ is the number of molecules of species *X _{i}* in voxel

*j*and

*N*is the total number of voxels. This is a local description since it describes what happens at each point in space inside the compartment. Now at the voxel level, no chemical conservation laws hold; such laws are only global. For example, the reaction $ X 1 + X 1 \u21cc X 2 $ has the conservation law

*n*

_{1}+ 2

*n*

_{2}= constant, which is defined on the total number of molecules of

*X*

_{1}and

*X*

_{2}in the compartment, but locally in voxel

*j*, $ n 1 j +2 n 2 j $ is not a constant due to the diffusive crosstalk with neighbouring voxels. It also follows that since we are considering a system in equilibrium, the deterministic concentration of a species in each voxel is the same as the deterministic concentration of the species in the whole compartment (that is equal to the solution of the non-spatial deterministic rate equations). Hence, given that there are only global conservation laws, that the local deterministic concentration is the same as the global deterministic concentration and that the voxel volume is Ω/

*N*, by analogy to the CME result (Eq. (6)), it follows immediately that the equilibrium probability distribution solution of the RDME is given by

where *n _{i}* is the global concentration of species

*X*, i.e.,$ n i = \u2211 j = 1 N n i j $.

_{i}### C. Global distribution of molecule numbers for finite size particle interactions

The result of van Kampen can also be straightforwardly extended to the vCME. We will assume that the degree of molecular crowding is not so high that it prevents well-mixing in the limit of long times; this is, the case if all molecules are mobile. The reactions here are modified than those for the CME because of the interaction of the chemical and empty space species. Hence, chemical system (3) is now modified to

where *X _{i}*,

*i*= 1, …,

*M*are the chemical species, and

*X*

_{M+1}is the empty space species. The deterministic equilibrium constants are then given by

where $ \varphi \u0303 i $ is the deterministic concentration of species *X _{i}* according to the deterministic rate equations one would write for the reaction scheme (8). Another difference from the CME is that besides the

*S*chemical conservation laws given by Eq. (5), we now also have an additional global conservation law stemming from volume exclusion, namely,

where *N* is the total number of molecules which can be maximally fit in the compartment. Given this information, by analogy with the CME result above (Eq. (6)), it follows immediately that the equilibrium probability distribution solution of the vCME is given by

Note that the global distribution is explicitly a function of *N*; this is in contrast to the global distribution of the CME which has no such information.

### D. Local distribution of molecule numbers for finite size particle interactions

The result of van Kampen can also be extended to the vRDME. We will assume that molecular crowding does not prevent diffusion between any two voxels in the compartment; this is the case if all molecules are mobile. This requirement is needed to satisfy detailed balance. Since the system is in equilibrium, the deterministic concentrations in each voxel are the same as the deterministic concentrations in the whole compartment according to the vCME. The state vector is $ { n 1 1 , \u2026 , n M + 1 1 , \u2026 , n 1 N , \u2026 , n M + 1 N } $, where $ n i j $ is the number of molecules of species *X _{i}* in voxel

*j*(1 ≤

*i*≤

*M*), $ n M + 1 j $ is the number of empty space molecules in voxel

*j*, and

*N*is the total number of voxels. A crucial difference from the RDME is that in addition to global conservation laws, now we also have a conservation law in each voxel, namely, there can be at most one molecule in each voxel, i.e., $ \u2211 k = 1 M + 1 n k i =1$, for

*i*= 1, …,

*N*. Given this information and taking into account the fact that the voxel volume is Ω/

*N*, by analogy with the CME result above (Eq. (6)), it follows immediately that the equilibrium probability distribution solution of the vRDME is given by

Note that because of the constraints due to conservation laws (chemical or volume exclusion), generally the mean concentration vector calculated using the exact equilibrium solutions of the CME, RDME, vCME, and vRDME are not equal to their deterministic concentration vector ($ \varphi \u2192 $ and $ \varphi \u0303 \u2192 $), respectively, except in the macroscopic limit of large volumes.

Note also that the local distribution solutions are independent of the underlying connectivity of the lattice of the RDME and vRDME. This is because in equilibrium, the solution of a master equation is generally a product of Poissonians constrained by local and global conservation laws,^{40} and these laws are not in any way influenced by the lattice connectivity. Of course as previously mentioned, the condition of detailed balance (equilibrium) is compatible only with those lattices such that there exists a path connecting any two lattice points. Thus, this is the only requirement on a lattice, for our results to hold.

It is also a fact that in detailed balance (equilibrium) conditions, the global probability distribution calculated starting from the local distribution solution of the RDME exactly matches the global distribution solution of the CME, independent of the diffusion coefficients. The same applies to the vRDME and the vCME. This maybe intuitive to some readers but for the sake of completeness we provide a proof in the Appendix. Thus although we initially presented the vCME in Section II via an intuitive approach, the macroscopic solution of the vCME stands on a solid basis since it can be obtained from the microscopic approach of the vRDME.

The rest of this article is devoted to obtaining insight into the effect of volume exclusion on the global distribution of molecule numbers and to a much lesser extent on the local distribution of molecule numbers. Due to the equivalence of the vRDME and vCME in equilibrium conditions and at the global level of description, we shall use both interchangeably when referring to any conclusions made assuming a finite molecular radius. Similarly, we shall use RDME and CME interchangeably when discussing conclusions made assuming point particles.

## IV. RELATIONSHIP OF RATE CONSTANTS IN THE CME AND vCME

Previously, we have denoted the rate constants in the vCME formalism by tildes to clarify that they are different to those in the CME. Here, we show the connection between the two.

We start by noting that the dilute limit of infinitesimally small molecules corresponds to the limit of infinitely large number of voxels *N* (in the vRDME and vCME) at constant compartment volume Ω. Equivalently, this corresponds to the limit, i.e., $\Omega \varphi \u0303 M + 1 \u2192N$, where practically all of space is empty (species *X*_{M+1} is the empty space species). In this limit, the deterministic (global) concentrations of the vCME and of the CME must be equal, which for system (3) implies

This equation encapsulates the relationship between the rate constants of the volume excluded and dilute probabilistic descriptions. For example, for the reversible dimerisation reaction previously considered, the CME and vCME formulations model the reactions $ X 1 + X 1 \u21cc k 1 \u2032 k 1 X 2 $ and $ X 1 + X 1 \u21cc k \u0303 1 \u2032 k \u0303 1 X 2 + X 3 $, respectively (where *X*_{3} is the empty space species), which implies the relation $ k \u0303 1 \u2032 / k \u0303 1 = k 1 \u2032 \Omega / k 1 N$.

It can be shown using the model reduction technique developed in Ref. 38 that in the limit of abundant empty space species (the dilute limit), the global distribution over the molecule numbers of the chemical species as given by the vCME, Eq. (11), tends to the global distribution of the CME, Eq. (6).

Using the relationship between the rate constants derived above, we can also understand how the effective equilibrium constant changes as a function of the strength of volume exclusion effects. According to the standard definition in physical chemistry and thermodynamics,^{39} the effective equilibrium constant of the j-th reaction in system (3) in volume excluded and dilute conditions are, respectively, given by

This implies that the effective equilibrium constant of unimolecular reactions is unaffected by crowding since in this case *r*_{M+1,j} = *s*_{M+1,j} = 0 because no space is involved. The effective equilibrium constant for bimolecular reactions is however increased relative to the one for non-crowded conditions, $ E \u0303 j > E j $, since in this case a single molecule of empty space is produced when two molecules bind to form a single molecule (*r*_{M+1,j} = 1, *s*_{M+1,j} = 0).

This result for bimolecular reactions can indeed be deduced without any calculation but with the application of Le Chatelier’s principle in physical chemistry^{39} to the vCME formalism. This principle states that a system in equilibrium will counteract the effect of an applied change by adjusting to a new equilibrium. Now if we consider the reaction $ X 1 + X 1 \u21cc X 2 $, this is modelled in the vCME formalism as $ X 1 + X 1 \u21cc X 2 + X 3 $ and hence by Le Chatelier’s principle, an increase in volume exclusion, i.e., a decrease in *X*_{3} (the empty space species) will induce the system to shift its equilibrium to the right to counteract this decrease, in the process causing an increase in the amount of species *X*_{2} and a decrease in the amount of species *X*_{1} which amounts to an increase in the effective equilibrium constant.

These results for unimolecular and bimolecular reactions encapsulate the essence of the thermodynamic theory of crowding developed by Minton and co-workers.^{24} However it is the first time, to our knowledge, that they have been obtained using the deterministic limit of a master equation description.

## V. STOCHASTIC DESCRIPTION OF CHEMICAL SYSTEMS WITHOUT CHEMICAL CONSERVATION LAWS

In this section, we use the results of Sections III and IV to show that *if there are no chemical conservation laws, then the marginal distribution of the global molecule numbers of each chemical species X_{i} is Poisson (Ωϕ_{i}) if crowding is ignored and binomial ($N,\Omega \varphi \u0303 i /N$) if crowding is taken into account*. Here, Ω is the compartment volume,

*N*is the maximum number of particles which can be placed in the compartment if volume exclusion is taken into account, and

*ϕ*and $ \varphi \u0303 i $ are the deterministic concentrations of the CME and vCME, respectively. We shall call this Statement 1. In what follows, we discuss the physical implications of this statement, as well as confirm our results using stochastic simulations of the CME and the vCME applied to an open homodimerisation reaction.

_{i}### A. Derivation of Statement 1

As shown in Section III, the global probability distribution assuming point particles, for a system with *M* chemical species, is generally given by the solution of the CME, namely, Eq. (6). Now if there are no chemical conservation laws, i.e., there is no factor *δ*(*f _{j}*(

*n*

_{1},

*n*

_{2}, …,

*n*),

_{M}*K*) in Eq. (6), then the global solution is simply a multivariate Poisson distribution,

_{j} and hence it follows that the marginal distribution of each chemical species *X _{i}* when volume exclusion is ignored is a Poisson with mean Ω

*ϕ*.

_{i}We also showed that the global probability distribution for molecules with a finite radius and assuming *N* of them can at most be packed in the compartment, for a system with *M* chemical species (and an additional species *X*_{M+1} representing free space), is generally given by the solution of the vCME, namely, Eq. (11). Now if there are no chemical conservation laws, i.e., there are no factors of the type *δ*(*f _{k}*(

*n*

_{1},

*n*

_{2}, …,

*n*),

_{M}*K*) in Eq. (11), then the normalised global probability distribution is given by

_{k} which is a multinomial distribution with parameters $ ( { \Omega \varphi \u0303 1 / N , \u2026 , \Omega \varphi \u0303 M + 1 / N} , N ) $. Note that $\Omega \varphi \u0303 i /N$ is the fraction of space occupied by particles of species *X _{i}* and consequently $ \u2211 i = 1 M + 1 \Omega \varphi \u0303 i /N=1$. It is well known that the marginal distributions of a multinomial distribution are binomial.

^{41}For species

*X*, the marginal distribution is thus binomial with parameters $ ( N , \Omega \varphi \u0303 i / N ) $,

_{i}### B. Dilute limit

Consider the dilute limit $\Omega \varphi \u0303 M + 1 \u2192N$. This can equivalently be seen as the limit of large numbers of voxels (at constant compartment volume Ω) in the vRDME such that the occupied volume fractions of all chemical species (except the empty space species) tend to zero and the deterministic solution of the vRDME (vCME) approaches that of the RDME (CME), i.e., *N* → ∞ and $\Omega \varphi \u0303 i /N\u21920$, such that $\Omega \varphi \u0303 i \u2192 \Omega \varphi i $ for 1 ≤ *i* ≤ *M*. Note that the last limit $\Omega \varphi \u0303 i \u2192 \Omega \varphi i $ for 1 ≤ *i* ≤ *M* follows by the specific relationship between the rate constants of the vRDME and of the RDME enforced in Section IV. Note also that the dilute limit implies infinitesimally small molecules, since the volume of each molecule is roughly that of a voxel Ω/*N*. It then follows by the Poisson limit theorem,^{42} that in the dilute limit, the global marginal distribution of the vRDME, binomial $ ( N , \Omega \varphi \u0303 i / N ) $, tends to the global marginal distribution of the RDME, Poisson (Ω*ϕ _{i}*).

### C. Statistical measures and physical implications

The Fano factor (*F*) is defined as the ratio of the variance and the mean and is a measure of how much a distribution differs from a Poisson distribution; the coefficient of variation (*CV*) is defined as the ratio of the standard deviation and the mean and is a measure of how “noisy” a system is; and the skewness (*SK*) of a distribution is the third standardised moment of the distribution and is a measure of how asymmetrical it is. These measures are well known for the Poisson and binomial distributions and hence we can state that assuming point particles, the statistics of chemical species *X _{i}* are given by

while for those modelled assuming a finite molecular radius, the statistics of chemical species *X _{i}* are given by

The differences between Eqs. (22) and (23) are illustrated in Fig. 2 where we plot the qualitative behaviour of the Fano factor, the coefficient of variation squared, and the skewness for a system, in which volume exclusion effects are neglected due to the assumption of point particles (green lines) and when they are taken into account (red lines).

The physical implication of these results is as follows. *As the fraction of occupied space increases, i.e., as 〈 n_{i}〉/N → 1, the fluctuations change from Poissonian (F = 1) to sub-Poissonian (F < 1), the well known classical noise-strength power law^{37,43} (CV_{i} ∝ 〈n_{i}〉^{−1/2}) becomes invalid, and the distribution of fluctuations changes from being skewed to the right (positive skewness) to being skewed to the left (negative skewness). The latter occurs when the occupied volume fraction 〈n_{i}〉/N exceeds 1/2*.

Another interpretation of the results is that if one ignores volume exclusion effects, i.e., employs the CME/RDME to model chemical systems without chemical conservation laws, then the dependence of the Fano factor, coefficient of variation, and the skewness on the mean molecule numbers is qualitatively wrong for high molecule numbers. It also follows from the properties of multivariate Poisson and multinomial distributions that ignoring volume exclusion implies zero covariance between the molecule numbers of different species while taking it into account implies a negative covariance.

### D. Application: Open homodimerisation reaction

We consider the dilute (point particle) chemical system,

whereby a species *X*_{1} is produced and subsequently two molecules of this species reversibly bind to form another molecule of type *X*_{2}. This system follows no chemical conservation laws and hence is of the type discussed above. The Fano factor, coefficient of variation, and skewness for the fluctuations in both species are given by Eqs. (22) together with the deterministic equilibrium constants,

This procedure leads to the following equations:

The volume exclusion version (assuming a finite molecular radius) of the chemical system (24) is given by

where species *X*_{3} is the empty space species. The Fano factor, coefficient of variation, and skewness for the fluctuations in both species are given by Eqs. (23) together with the deterministic equilibrium constants,

and the conservation law due to volume exclusion,

where *N* is the total number of molecules which can be contained in the compartment. Furthermore, we know that in the dilute limit, *ϕ*_{3} → *N*/Ω, the effective equilibrium constants of the crowded system must equal the equilibrium constants of the non-crowded system (as previously discussed at length in Section IV). Thus using Eqs. (25) and (28), we have the following relationship between the rate constants of the crowded system and of the non-crowded system:

The overall procedure described above leads to the following equations:

Comparing the statistical quantities Eqs. (26) and (31), one notices the stark difference in the parametric dependence of these quantities if volume exclusion effects are taken into account. These differences are illustrated in Fig. 3 where the solid lines show the analytical predictions for the Fano factor, coefficient of variation, and skewness of both species as a function of the parameter *k*_{0}, for dilute (upper panel) and volume exclusion conditions (lower panel). The analytical formulae are compared with data from the SSA (open circles) sampling the CME and the vCME, as evidence of their exactness.

As one can appreciate from these plots, the dependence on *k*_{0} is strongly affected by volume exclusion, except of course in the limit of small *k*_{0} where there are few particles in the compartment. Some qualitative differences which are particularly noticeable and interesting are the following: (i) volume exclusion has very little impact on the Fano factor of species *X*_{1} but a strong impact on the Fano factor of species *X*_{2} (a change from constant to strongly monotonic decreasing as a function of *k*_{0}); (ii) in contrast, volume exclusion has a strong impact on the coefficient of variation of species *X*_{1} (a change from a monotonic decreasing function to a parabolic function of *k*_{0}) but little impact on the coefficient of variation of *X*_{2}; and (iii) the skewness of species *X*_{2} becomes negative as *k*_{0} increases beyond a certain threshold, for volume excluded conditions, but remains positive in dilute conditions. These stark differences suggest that the parameter dependencies of various statistical quantities that one obtains using conventional dilute approaches (the RDME and CME) may not always reflect the actual parameter dependencies *in vivo*.

## VI. STOCHASTIC DESCRIPTION OF CHEMICAL SYSTEMS WITH A SPECIAL TYPE OF CHEMICAL CONSERVATION LAWS

In this section, we study systems with a chemical conservation law implying that the sum of the molecule numbers of some of the species is a constant *k*. For these systems, we show that *(i) the marginal distribution of the global molecule numbers of a chemical species X_{i} not involved in the conservation law is Poisson (Ωϕ_{i}) if volume exclusion is ignored and binomial ($N\u2212k,\Omega \varphi \u0303 i / ( N \u2212 k ) $) if it is taken into account. (ii) The marginal distribution of the global molecule numbers of a chemical species X_{i} involved in the conservation law is binomial (k, Ωϕ_{i}/k) if volume exclusion is ignored and binomial ($k,\Omega \varphi \u0303 i /k$) if it is taken into account*. We shall call these Statements 2 and 3, respectively. We also discuss the physical implications of these statements, as well as confirm our results using stochastic simulations of the RDME and the vRDME applied to an open heterodimerisation reaction.

### A. Derivation of Statements 2 and 3 and the dilute limit

We consider a chemical system with *M* chemical species and a chemical conservation law of the form

where *n _{i}* is the number of molecules of species

*i*, and 1 ≤

*L*≤

*M*− 2. This is a special case of the general global conservation law considered in Section III. It implies that there are no restrictions on the number of molecules of species

*X*

_{1}, …,

*X*, but that the sum of the number of molecules of species

_{L}*X*

_{L+1}, …,

*X*is constant at all times.

_{M}The global probability distribution for *M* chemical species, assuming point particles, Eq. (6), is then given by

which is the product of a multivariate Poisson distribution with parameters ({Ω*ϕ*_{1}, …, Ω*ϕ _{L}*}) and a multinomial distribution with parameters ({Ω

*ϕ*

_{L+1}/

*k*, …, Ω

*ϕ*/

_{M}*k*},

*k*). The multinomial originates from the constraint placed by the chemical conservation law Eq. (32). Hence it follows, by the same arguments as in Sec. V, that if a chemical species

*X*is not involved in the chemical conservation law, then the marginal distribution is Poisson (Ω

_{i}*ϕ*) whereas if it is involved in the chemical conservation law then the marginal distribution is binomial (

_{i}*k*, Ω

*ϕ*/

_{i}*k*).

The global probability distribution for *M* chemical species, assuming a finite molecular radius, Eq. (11), specialised to the conservation law Eq. (32), is given by

which is the product of a multinomial distribution with para-meters $ ( { \Omega \varphi \u0303 1 / ( N \u2212 k ) , \u2026 , \Omega \varphi \u0303 L / ( N \u2212 k ) , \Omega \varphi \u0303 M + 1 / ( N \u2212 k ) } , N \u2212 k ) $ and a multinomial with parameters $ ( { \Omega \varphi \u0303 L + 1 / k , \u2026 , \Omega \varphi \u0303 M / k} , k ) $. The latter multinomial originates from the chemical conservation law Eq. (32). The former multinomial originates from the combination of the chemical conservation law Eq. (32) and the volume exclusion law in the vRDME, $ \u2211 i = 1 M + 1 n i =N$, which together implies the combined conservation law $ \u2211 i = 1 L n i + n M + 1 =N\u2212k$. Hence, it follows that if a chemical species *X _{i}* is not involved in the chemical conservation law, then the marginal distribution is binomial $ ( N \u2212 k , \Omega \varphi \u0303 i / ( N \u2212 k ) $ whereas if it is involved in the chemical conservation law then the marginal distribution is binomial $ ( k , \Omega \varphi \u0303 i / k ) $. It is straightforward to verify using the Poisson limit theorem that in the dilute limit, the global binomial solution of the vRDME approaches the Poisson solution of the RDME.

### B. Statistical measures and Physical implications

For those species not involved in the chemical conservation law, the marginal is Poisson (Ω*ϕ _{i}*) and hence the statistical quantities are given by Eq. (22) if one assumes dilute conditions. If volume exclusion effects are considered then the marginal distribution is binomial $ ( N \u2212 k , \Omega \varphi \u0303 / ( N \u2212 k ) ) $ and hence the statistics are given by Eq. (23) with the parameter

*N*replaced by

*N*−

*k*. Similarly it can be reasoned that for those species involved in the chemical conservation law, i.e. species

*X*

_{L+1}, …,

*X*, the statistics are given by Eq. (23) with the parameter

_{M}*N*replaced by

*k*and $ \varphi \u0303 $ replaced by

*ϕ*if dilute conditions are assumed and by Eq. (23) with the parameter

*N*replaced by

*k*if volume exclusion is taken into account.

The physical implication of these results is as follows. For both species which are involved and not involved in the chemical conservation law, taking into account volume exclusion implies that the marginal distribution is binomial and hence we can make the same statement as for chemical systems without any chemical conservation laws. Namely for chemical systems with a chemical conservation law of the type Eq. (32), *as the extent of volume exclusion increases, i.e, as 〈 n_{i}〉/N → 1, the fluctuations become more sub-Poissonian, deviations from the classical noise-strength power law become more apparent, and the distribution of fluctuations changes from being skewed to the right (positive skewness) to being skewed to the left (negative skewness). The latter occurs when the 〈n_{i}〉/(N − k) exceeds 1/2 for species not involved in the chemical conservation law and when 〈n_{i}〉/k exceeds 1/2 for species involved in the chemical conservation law*.

However, there are also some differences between the results here and those of Sec. VI A. The RDME predicts the wrong qualitative dependence of the Fano factor, coefficient of variation, and the skewness on the mean molecule numbers for all species in the system without any chemical conservation law. For systems with a chemical conservation law, this is still true for those species not involved in a chemical conservation law. However, the RDME does predict the right qualitative dependence for those species involved in the chemical conservation law (since it predicts a binomial marginal distribution, same as the vRDME, though with different parametrisation).

The results here can also be generalised to a system with a number of chemical conservation laws of the sum type. For example, for a system with two conservation laws of the type,

where 1 ≤ *z* ≤ *L* − 2, by a similar reasoning as above, we find, assuming a finite molecular radius, that the marginal distributions of species *X _{i}* is binomial $ ( N \u2212 k \u2212 s , \Omega \varphi \u0303 i / ( N \u2212 k \u2212 s ) ) $ if it is not involved in the conservation laws, binomial $ ( s , \Omega \varphi \u0303 i / s ) $ if it is involved in the first conservation law, and binomial $ ( k , \Omega \varphi \u0303 i / k ) $ if it is involved in the second conservation law.

### C. Application: Open heterodimerisation reaction

We now consider the dilute (point particle) model of the chemical system,

whereby a species *X*_{1} is produced and subsequently molecules of this species and that of *X*_{2} reversibly bind to form molecules of *X*_{3}, a heterodimer. The system has the implicit chemical conservation law *n*_{2} + *n*_{3} = *k* where *k* is a time-independent constant determined by the initial conditions, and hence it is of the type studied above. The deterministic equilibrium constants are

The volume excluded version (assuming finite molecular size) of reaction scheme (36) is

where *X*_{4} is the empty space species and now we have two conservation laws: the chemical law *n*_{2} + *n*_{3} = *k* and the volume exclusion law *n*_{1} + *n*_{2} + *n*_{3} + *n*_{4} = *N* where *N* is the maximum number of molecules which the compartment can accommodate. The deterministic equilibrium constants are

where we used the relationship between the rate constants of the volume excluded and dilute systems (as in the previous example and as elucidated in Section IV).

This implies that the concentrations of species *X*_{2} and *X*_{3} (the species involved in a chemical conservation law) are insensitive to volume exclusion effects but the concentration of species *X*_{1} is found to decrease when crowding is taken into account. Intuitively, this is because species *X*_{2} and *X*_{3} are involved in a chemical conservation law and hence the impact of the second conservation law due to volume exclusion is nullified; species *X*_{1} in contrast is not involved in any chemical conservation law and is hence strongly affected by the conservation law due to volume exclusion.

According to our theory in Sec. VI B, (i) the marginal global distribution of species *X*_{1} is Poisson (Ω*ϕ*_{1}) according to the RDME and binomial $ ( N \u2212 k , \Omega \varphi \u0303 1 / ( N \u2212 k ) ) $ according to the vRDME. The mean of the latter is less than that of the former. This is verified via stochastic simulations of the RDME and vRDME using the SSA—see Fig. 4(a); (ii) the marginal global distribution of species *X*_{2} is binomial (*k*, Ω*ϕ*_{2}/*k*) for both the RDME and vRDME. This is also verified via stochastic simulations using the SSA—see Fig. 4(b). In Fig. 4(c) and 4(d), we also show that stochastic simulations using the SSA agree with the theoretical expressions obtained by marginalising the local (voxel) distributions given by Eqs. (7) and (12) in Section III. Note that for the purpose of stochastic simulations using the RDME and vRDME, we need to specify a lattice type. We choose the RDME and vRDME lattices to be periodic, square, and in two dimensions, with the neighbourhood of a voxel being the four cells orthogonally surrounding it. The compartment volume Ω will be fixed to one, meaning that for *N* voxels, the lattice consists of $ N \xd7 N $ voxels with lattice spacing $1/ N $. We shall use this lattice for all stochastic simulations in this article.

Of interest is how the vRDME probability distribution of the global number of molecules of species *X*_{1} changes as the ratio of molecular diameter to compartment length scale is varied. The ratio of the compartment side length to the molecular diameter (the lattice spacing) is given by $ N $. In Fig. 5(a), we plot the global marginal probability distribution solution of the vRDME for species *X*_{1}, i.e., binomial $ ( N \u2212 k , \Omega \varphi \u0303 1 / ( N \u2212 k ) ) $, as a function of the total number of voxels *N* while keeping the compartment volume constant. Good agreement of the vRDME and RDME solutions is obtained when *N* = 1600, i.e, when the compartment side length is about forty times larger than the molecular diameter; here the molecules are small enough that the system is dilute. In contrast, clear differences exist between the vRDME and RDME predictions when *N* = 100 (and smaller values) which corresponds to the case of molecules whose diameter is at least 1/10 of the compartment side length; for these cases, the RDME overpredicts the true global concentration of *X*_{1}.

It is also interesting to understand the effects of increasing the degree of volume exclusion by adding inert molecules to the chemical reaction system. This is of particular relevance to understanding intracellular reaction systems which typically operate in such conditions, i.e., molecules of other intracellular pathways which are inert with respect to the reaction system of interest exert influence on the latter via volume exclusion effects.^{23} To this end, we consider a modified version of reaction scheme (38),

where *X*_{4} is the empty space species and *X*_{5} is a chemical species which does not chemically interact with the rest of the molecules (an inert external crowder). In Fig. 5(b), we plot the global marginal probability distribution solution of the vRDME for species *X*_{1}, i.e., binomial $ ( N \u2212 k , \Omega \varphi \u0303 1 / ( N \u2212 k ) ) $, as a function of the mean number of inert external crowder molecules 〈*n*_{5}〉. Note that the effect of increasing molecular crowding by adding more molecules of *X*_{5} is to induce a shift of the probability distribution to the left such that there are fewer molecules, on average, of *X*_{1} in the system. This is qualitatively similar to the effect seen in Fig. 5(a). This similarity arises because an increase in the fraction of occupied space can either be induced by increasing the size of the reactant molecules while keeping the compartment size fixed (the case of Fig. 5(a)) or else by introducing inert molecules into the system (the case of Fig. 5(b)). Note that in both cases, the marginal distribution of *X*_{2} is unchanged by the degree of volume exclusion since as we noted earlier both the RDME and vRDME give the same result.

## VII. STOCHASTIC DESCRIPTION OF CHEMICAL SYSTEMS WITH MORE GENERAL CHEMICAL CONSERVATION LAWS

Previously, we have considered chemical conservation laws of type (32). Though common, these are not the only chemical conservation laws in nature. The general theory presented in Section III also applies to these other conservation laws. In what follows, we use the latter results to study an example of a chemical system constrained by a chemical conservation law which is not of the sum type. In particular, we will show that in this case, the global marginal distribution of the number of molecules for a species involved in the conservation law is not binomial, unlike the case of a species involved in a chemical conservation law of type (32).

### A. Closed dimerisation reaction

Consider the point particle model of the reaction system,

whereby two molecules of *X*_{1} reversibly bind to form a dimer *X*_{2}. This system has the global chemical conservation law *n*_{1} + 2*n*_{2} = *k* where *k* is a time-independent constant fixed by the initial conditions and hence it is not of the same type as the chemical conservation laws (32) considered earlier. According to Eqs. (6) and (11), the (normalised) marginal probability distribution solution for species *X*_{2} according to the CME and the vCME is given by

respectively. Here, we have used the notation *H _{U}* and $ H 2 F 1 R $ to denote Tricomi’s confluent hypergeometric function and the regularised hypergeometric function, respectively (these are the functions HypergeometricU and Hypergeometric2F1Regularised in Mathematica’s notation

^{44}). Note that for the vCME, we have here considered the volume excluded version of reaction scheme (42), namely, $ X 1 + X 1 \u21cc k \u0303 1 k \u0303 0 X 2 + X 3 $ with

*X*

_{3}representing the empty space species and the equilibrium constant $ k \u0303 0 k \u0303 1 = k 0 N k 1 \Omega $ (as elucidated in Section IV). All the statistics of the molecule numbers of species

*X*

_{1}can be deduced from those of

*X*

_{2}given the conservation law

*n*

_{1}+ 2

*n*

_{2}=

*k*.

There are here clearly differences from what we previously found for chemical species involved in chemical conservation laws of type (32). While for the latter, the global marginal distributions where binomial independent of whether volume exclusion is taken into account or not (see Section VI), in the example presented in this section, the global marginal distributions are not binomial. This difference can be traced to the fact that a binomial originates as the marginal of a multinomial distribution and the latter is effectively a product of Poissons constrained by a rule stating that the sum of the fluctuating variables is a constant; this rule is naturally obeyed by systems in which the chemical conservation law is of the type (32).

In Fig. 6, we plot the steady-state probability distribution of global molecule numbers according to the CME and vCME for the case when *N* = *k*, i.e., the minimum number of voxels required to accommodate the maximum number of molecules allowed by the dimerisation reaction. We note that while the chemical conservation law shielded the effects of volume exclusion law for those species involved in laws of type (32) (see Fig. 4(b)), no such shielding occurred in the example here, as can be appreciated from the large difference between the two distributions in Fig. 6. Likely, the implicit mathematical reason for these differences is that for systems in Section VI, the chemical conservation law $ \u2211 i = L + 1 M n i =k$ is “nested” within the volume exclusion law $ \u2211 i = 1 M + 1 n i =N$, while no such nesting is present in the current example where the chemical conservation law is *n*_{1} + 2*n*_{2} = *k*.

We now study the effect of volume exclusion on various statistics. We first note that the first and second moments of the vCME solution can be conveniently written in terms of three non-dimensional parameters: *k*, *R* = *N*/*k* and *L* = 4*k*_{0}/(*k*_{1}Ω). The parameter *L* contains information about all the rate constants of the system; it is proportional to the equilibrium constant *k*_{0}/*k*_{1} of the reaction in the absence of volume exclusion. The parameter *R* is an inverse measure of volume exclusion. This is since as *N* increases at constant compartment volume Ω, molecules “become smaller” and hence the system becomes more dilute. The maximum degree of volume exclusion occurs when *N* = *k*, i.e, *R* = 1 and the dilute limit occurs when *R* → ∞. In Fig. 7, we fix *k* = 50 and use Eq. (44) to calculate the statistics in very dilute conditions (*R* = 1000) and in highly crowded conditions (*R* = 1) as a function of the parameter *L*. The dilute statistics agree very well with those which can be calculated directly from the CME using Eq. (43).

In particular we find that (i) the Fano factor of species *X*_{2} is always less than one and hence the distribution is sub-Poissonian in both volume excluded and dilute conditions (see Fig. 7(a)); (ii) the Fano factor of species *X*_{1} can be greater than or less than 1 leading to three distinctive phases: sub-Poisson statistics in both volume excluded and dilute conditions (for *L* < 7), super-Poissonian in both conditions (for *L* > 11) and finally a phase in which volume exclusion leads to a change from sub-Poissonian to super-Poissonian statistics (for 7 ≤ *L* ≤ 11, see Fig. 7(b)). This is in contradistinction to the results in Section VI where we found that a species involved in chemical conservation laws of the type (32) has sub-Poissonian fluctuations in both volume excluded and dilute conditions; (iii) volume exclusion leads to a decrease in the coefficient of variation of species *X*_{2} and to an increase in the coefficient of variation of species *X*_{1} (see Fig. 7(c)); (iv) volume exclusion leads to an increase in the mean number of molecules of species *X*_{2} and to a decrease in the number of molecules of species *X*_{1} (see Fig. 7(d)). Thus, the inclusion of volume exclusion causes a shift of the equilibrium to the right, namely, it leads to the production of more *X*_{2} molecules and of less *X*_{1} molecules. This is in agreement with the standard thermodynamic theory by Minton and co-workers.^{24} We have numerically verified that these results hold for even *k*.

As we saw in this example, the general properties of systems with chemical conservation laws of a general type are not typically open to analytical investigation because of the complicated form of the exact steady-state probability distribution solution of the CME and vCME. Nevertheless, these expressions can be easily investigated numerically which is advantageous compared to lengthy stochastic simulations.

## VIII. COMPARISON OF THE vRDME WITH BROWNIAN DYNAMICS

The vRDME has at least one major disadvantage—it is based on an artificial spatial lattice. Ideally, one would like to derive the vRDME rigorously from a lattice-free approach or at least to show that it is a reliable approximation of a lattice-free description under some conditions.

A derivation of this type has been previously attempted for the dilute case. In particular, it has been shown that for systems with bimolecular reactions, the RDME provides a good approximation to the lattice-free descriptions offered by the Doi^{45,46} and Smoluchowski models^{47,48} for lattice spacings that are neither too small nor too large.^{49} In the limit of small lattice spacing, the statistics of the RDME do not converge to those of the lattice-free approach^{9,50} but it is possible in this case to derive a new convergent RDME called the CRDME which does not suffer from this issue.^{51}

The question of agreement between a lattice and lattice-free approach in the case of volume excluded interactions is relatively simpler than for the dilute case because there is one less parameter: unlike the RDME, the lattice spacing of the vRDME is fixed to equal the molecule diameter. A formal derivation of the vRDME from the volume excluded versions of the spatially continuous Doi or Smoluchowski models is beyond the scope of this paper; here we shall be content with comparing the statistics of the vRDME with those obtained from microscopic BD simulations for a simple example.

In particular, we test the validity of the RDME and vRDME by comparing their global distribution solutions for the closed dimerisation system (42) given by Eqs. (43) and (44), respectively, with the distributions calculated from ensemble averaging BD simulations of the same chemical system. The BD simulations consider particles to be two-dimensional hard disks which move randomly in space according to standard Brownian motion. With parameters defined as in Eqs. (43) and (44), diffusion coefficient *D* and time-step Δ*t*, the BD algorithm we use is as follows:

Place

*k*particles of type*X*_{1}with radius*r*at uniformly distributed points in [0, 1] × [0, 1] such that they do not overlap. Set time*t*= 0. Proceed to (2).Propose a normal random number with mean 0 and standard deviation $ 2 D \Delta t $ to add to each particle coordinate. If no pairs of particles will overlap, accept the new coordinates and proceed to (4). If exactly one pair of particles will overlap and they are both type

*X*_{1}, proceed to (3). Else reject the new coordinates and attempt (2) again.Choose a uniform random number

*rand*between 0 and 1. If*rand*≥*p*Δ*t*reject the new coordinates from (2) and attempt (2) again. Else if*rand*<*p*Δ*t*, remove the overlapping*X*_{1}particles. Place a*X*_{2}particle with radius*r*midway between the centres of the removed particles. Choose an exponential random number*exprand*with mean 1/*k*_{1}. Assign a number*τ*=*t*+*exprand*to the new*X*_{2}particle. Proceed to (4).For each

*X*_{2}particle, check if*t*>*τ*. If not, proceed to (5). Else for each*X*_{2}particle with*t*>*τ*, remove it and place two*X*_{1}particles just touching at a random angle such that their midpoint is the former centre of the*X*_{2}particle. If any of the new*X*_{1}particles overlap other particles, remove them, replace the*X*_{2}particle, and set*τ*=*t*+*exprand*. Proceed to (5).Advance time by setting

*t*=*t*+ Δ*t*. Store the total number of*X*_{1}and*X*_{2}particles in memory. Return to (2) and repeat until a given time has elapsed.

Note that, in the above algorithm, $p= k 0 2 \pi r 2 $ which is the probability per unit time that a given pair of *X*_{1} particles reacts. This choice guarantees that in the limit of well-mixed and dilute conditions, the rate at which dimerisation occurs in the Brownian dynamics agrees with that given by the bimolecular propensity in the CME (for a derivation see Appendix D of Ref. 52). Note also that the precise choice of the distance at which one places the two particles of type *X*_{1} when a dimer *X*_{2} dissociates has little effect on the statistics collected, as long as we have reaction-limited kinetics (probability of the association of two particles of type *X*_{1} is very small). The above algorithm can be considered a volume-excluded version of standard BD algorithms used for dilute reversible systems.^{53}

For accuracy, Δ*t* should be chosen small enough such that at most one reaction normally happens in each time step. To compare BD and vRDME, we choose the particles to have a diameter equal to the width of a vRDME voxel. This implies that the proportion of volume occupied by a BD particle is slightly less than the proportion of volume occupied by a vRDME voxel; however, it is the most natural way of assigning a diameter, and it ensures that BD can feasibly reach the levels of volume exclusion that we want to model with the vRDME.

In Fig. 8, we compare BD simulations with the exact global distributions of the RDME and the vRDME as given by Eqs. (43) and (44), respectively. In panel 8(a), we show the equilibrium global probability distribution of *X*_{2} computed with BD (blue histogram), vRDME (yellow histogram), and RDME (grey dashed line), in dilute conditions. In this case, in BD, the particle diameters were $ 1 20 $ and there were 24 *X*_{1} particles initially; equivalently, in the vRDME, the number of voxels is *N* = 400. It follows that the percentage of occupied volume in this case varies from 3 − 6%, where 3% is reached when all *X*_{1} particles are bound in dimers *X*_{2}. Since this corresponds to fairly dilute conditions, it is unsurprising that BD, the vRDME and the RDME essentially agree. In panel 8(b), we show the same plot in high volume exclusion conditions. In this case, in BD, the particle diameters were $ 1 6 $ and there were 24 *X*_{1} particles initially; equivalently, in the vRDME, the number of voxels is *N* = 36. Therefore the percentage of occupied volume in this case varies from 33 − 67%. Thus, this corresponds to considerably high volume exclusion; the vRDME here agrees with BD but the RDME strongly disagrees with both.

Hence our analysis confirms that for the dimerisation reaction, the vRDME gives global statistics that are in very good agreement with those obtained from a microscopic lattice-free approach, for a parameter set in both low and high volume exclusion scenarios. This is likely mainly due to the fact that the vRDME is a description at the natural length scale of the system (the molecular diameter). Further research is however necessary to clarify whether the agreement between the vRDME and BD holds for a broad range of parameter values and for general chemical systems.

## IX. SUMMARY AND CONCLUSION

In this paper, we have elucidated some of the effects which volume exclusion can have on intrinsic noise in chemical systems which are in equilibrium. In particular, the novelty of our study is that we can make precise statements on the relationship between the probability distribution solution of the master equation and the extent of volume exclusion. This was possible because we obtained an exact solution of the local and global probability distribution of the RDME and of its excluded volume version, the vRDME, in equilibrium (detailed balance) conditions.

A summary of our findings is as follows. We found that the type of the global marginal distributions of the RDME and vRDME varies according to the type of chemical conservation law: (i) for those systems with no chemical conservation law, the global marginal distribution of the RDME and vRDME for all species is Poisson and binomial, respectively; (ii) for those systems with a chemical conservation law of the sum type, Eq. (32), the global marginal distribution of the RDME and vRDME for a species not involved in the chemical conservation law is also Poisson and binomial, respectively; (iii) for those systems with a chemical conservation law of the sum type, Eq. (32), the global marginal distribution of the RDME and vRDME for a species involved in the chemical conservation law is binomial. Taking into account volume exclusion has very little or no impact on the fluctuations in this case; (iv) for those systems with a chemical conservation law of a more general type, nothing can be directly deduced about the type of marginal distributions because of the complexity of the exact normalised probability distributions. However for a specific system of this type, we found that the global fluctuations were neither Poisson nor binomial for species involved in the chemical conservation law and that volume exclusion did have a strong impact on the fluctuations, in contrast to systems with a chemical conservation law of the sum type.

Given points (i)-(iii), we can clearly state that the largest impact of volume exclusion is likely to be on the intrinsic noise statistics of those species not involved in chemical conservation laws; the fact that the RDME solution is Poisson while the vRDME solution is binomial implies that as the extent of molecular crowding increases, the fluctuations become increasingly sub-Poissonian, deviations from the classical inverse square root law for the noise-strength become more apparent, and the marginal distribution of molecule number fluctuations changes from being skewed to the right (positive skewness) to being skewed to the left (negative skewness).

We note that the vRDME used in our study is based on an inherent assumption that the size of all molecules, reactant and inert, is roughly the same and equal to the size of a voxel. This is, of course, a gross simplification of reality, nevertheless the major benefits of this formulation is that (i) the vRDME is exactly solvable in equilibrium conditions and (ii) it appears to be an accurate approximation of microscopic spatially continuous stochastic simulations. Hence a comparison of the exact solution of the vRDME with the exact solution of the RDME (which assumes point particles) allows us to obtain a rough picture of the effects of volume exclusion on intrinsic noise, results which are difficult to obtain if we had to resort to computationally expensive stochastic simulations.

Open questions which remain to be addressed involve understanding the impact of volume exclusion on non-equilibrium steady-states and on the time evolution of moments; these are challenging questions given that exact solutions of master equations are highly unlikely to be found in such conditions. Finally, we expect the extension of the vRDME framework to allow the modelling of chemical reactions involving hard molecules of various sizes to be of paramount importance for the accurate prediction of the effect of volume exclusion on real chemical systems.

## Acknowledgments

This work was supported by the BBSRC EASTBIO Ph.D. studentship to S.S. and by a Leverhulme grant award to R.G. (No. RPG-2013-171). C.C. thanks Philipp Thomas for useful discussions.

### APPENDIX: DERIVATION OF THE GLOBAL DISTRIBUTION OF MOLECULE NUMBERS OF THE vRDME AND RDME

As discussed in Section III in the main text, it is straightforward to show that the solution of the vRDME in equilibrium conditions is

where *n _{i}* is the global concentration of species

*X*, i.e., $ n i = \u2211 j = 1 N n i j $.

_{i}We now use Eq. (A1) to calculate the probability over the vector of the global number of molecules $ n \u2192 = { n 1 , \u2026 , n M } $. We start by noting that the definition of the global concentration of species *X _{i}*, i.e., $ n i = \u2211 j = 1 N n i j $ together with the conservation law factor $ \u220f k = 1 N \delta ( \u2211 i = 1 M + 1 n i k , 1 ) $ is equivalent to the factor $\delta ( \u2211 i = 1 M + 1 n i , N ) $. Thus, we have

The passage from Eq. (A2) to Eq. (A4) can be explained as follows. The sum in Eq. (A2) is over the local molecule numbers only and hence the delta function over the global molecule numbers $\delta ( \u2211 i = 1 M + 1 n i , N ) \delta ( f m ( n 1 , n 2 , \u2026 , n M ) , K m ) $ is unaffected by this sum and can be left outside, which leads to Eq. (A3). Now the term in square brackets in the latter equation is a product of independent Poissonians (the correlation between Poissonians is induced by the delta functions outside of the square brackets). Due to the delta function $\delta ( \u2211 r = 1 N n i r , n i ) $, the sum in the square brackets amounts to calculating the probability distribution of a sum of independent Poisson random variables, which leads to the final result Eq. (A4). Note that Eq. (A4) is the same as the equilibrium solution of the vCME, Eq. (11), which establishes the equivalence of the vRDME and vCME at the global level in equilibrium conditions. By an analogous approach, one can also show the equivalence of the RDME and CME at the global level in equilibrium conditions.