Molecular finite-size effects in stochastic models of equilibrium chemical systems

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 crowded reaction-diffusion master equation (cRDME) 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 cRDME 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 molecular crowding 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; (iv) strongly modulate the Fano factors and coefficients of variation. These crowding-induced 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 [1,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 are 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 [6,7], 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 [11][12][13][14]) or else their solution computed by means of an approximative method such as moment-closure approximation [15][16][17] or the system-size expansion [18][19][20][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][23][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][27][28][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 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 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 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 i denotes species B in voxel i and the notation N e(i) stands for the set of voxels which neighbour voxel i.The parameter k D has units of inverse time and is proportional to the diffusion coefficient D of the species; specifically 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 D .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 the next section.
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 i takes a value of zero if voxel 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 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][35][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 [33]).
The vCME is the non-spatial counterpart of the vRDME.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: 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 the next section.
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 EQUI-LIBRIUM 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 later sections 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) level.We shall start by summarising a result by van Kampen for the CME, which we will subsequently extend to the other three frameworks.
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 ith chemical species.Here k j and k j are the rate constants for the forward and reverse reactions respectively and r ij − s ij is the change in the number of molecules of species X i when reaction 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 i (as given by the conventional non-spatial rate equations).Furthermore we assume that the system has a number, S, of chemical conservation laws of the form: where n = {n 1 , n 2 , ..., 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 wellmixed 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) is the probability that the system is in state n in equilibrium.Hence the equilibrium solution is a product of Poisson distributions constrained by the chemical conservation laws.
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 , ..., n 1 M , ..., n N 1 , ..., n N M }, where n j i 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 − − X 2 has the conservation law n 1 + 2n 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 j 1 + 2n j 2 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 above (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 , i.e., n i = N j=1 n j i .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 the 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 φ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.
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 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., M +1 k=1 n i k = 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 ( φ and φ) 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 Appendix A. 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., Ω φM+1 → N , 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 the vCME and of the CME must be equal, which for system (3) implies: This statement together with Eqs. ( 4) and ( 9) 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 respectively (where X 3 is the empty space species), which implies the relation k It can be shown using the model reduction technique developed in [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 jth reaction in system (3) in volume excluded and dilute conditions are respectively given by: Now by Eq. ( 9) and Eq. ( 14) we then have:

=
E j (fraction of empty space) r M +1,j −s M +1,j . ( 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, Ẽ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

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 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 φ i and φ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.

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 j ) in Eq. ( 6), then the global solution is simply a multivariate Poisson distribution: 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 k ) in Eq. ( 11), then the normalised global probability distribution is given by: which is a Multinomial distribution with parameters ({Ω φ1 /N, ..., Ω φM+1 /N }, N ).Note that Ω φi /N is the fraction of space occupied by particles of species X i and consequently i=1 Ω φi /N = 1.It is well known that the marginal distributions of a multinomial distribution are Binomial [42].For species X i , the marginal distribution is thus Binomial with parameters (N, Ω φi /N ): B. Dilute limit Consider the dilute limit Ω φM+1 → N .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 Ω φi /N → 0, such that Ω φi → Ωφ i for 1 ≤ i ≤ M .Note that the last limit Ω φi → Ωφ 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 [43], that in the dilute limit, the global marginal distribution of the vRDME, Binomial (N, Ω φ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 equations Eqs.(22) and Eq. ( 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,44] (CV i ∝

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- 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 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: (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 ; (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 SPE-CIAL 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  for species X 1 (red) and X 2 (green).From left to right, we plot the Fano factor, coefficient of variation, and skewness as a function of the parameter k 0 , using Eqs.(26) for the upper panel (dilute conditions, CME) and Eqs.(31) for the lower panel (volume exclusion conditions, vCME).
The parameter values are: k 1 = k 2 = k 3 = 0.1, Ω = 10, N = 200.We note that in the dilute limit k 0 ≈ 0, the two systems have the same behaviour.
account.We shall call these Statement 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 L , but that the sum of the number of molecules of species X L+1 , ..., X M is constant at all times.
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 the previous section, that if a chemical species X i is not involved in the chemical conservation law, then the marginal distribution is Poisson 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 parameters ({Ω φ1 and a multinomial with parameters ({Ω φL+1 /k, ..., Ω φ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, M +1 i=1 n i = N , which together imply the combined conservation law 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 − k, Ω φi /(N − k) whereas if it is involved in the chemical conservation law then the marginal distribution is Binomial (k, Ω φ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 − k, Ωφ i /(N − 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 M , the statistics are given by Eq. ( 23) with the parameter N replaced by k and φ 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 the previous section.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 if it is not involved in the conservation laws, is Binomial (s, Ω φi /s) if it is involved in the first conservation law and Binomial (k, Ω φ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).
Explicit solution of Eqs.(37) and Eqs.(39) together with the relevant conservation laws leads to: 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 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 the previous section, (i) the marginal global distribution of species X 1 is Poisson (Ωφ 1 ) according to the RDME and Binomial (N − k, Ω φ1 /(N − 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 × √ 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  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 , 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 the 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 the type (32).(38).While (b) is obtained by varying the ratio k4 / k5 which controls the mean number of inert external crowders n 5 in system (41).The parameters k, k 0 , k 1 , k 2 , k 3 , Ω for both (a) and (b) are the same as in Fig. 5. See text for discussion.

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 + 2n 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 Eq. ( 6) and Eq. ( 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 R 2F 1 to denote Tricomi's confluent hypergeometric function and the regularised hypergeometric function respectively (these are the functions HypergeometricU and Hypergeometric2F1Regularised in Mathematica's notation [45]).Note that for the vCME, we have here considered the volume excluded version of reaction scheme (42), namely X 1 + X 1 k0 − − k1 X 2 + X 3 with X 3 representing the empty space species and the equilibrium constant k0 k1 = k 0 N k 1 Ω (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 + 2n 2 = k.
There are here clearly differences from what we previously found for chemical species involved in chemical conservation laws of the 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 the type (32) (see Fig. 4b), 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 M i=L+1 n i = k is "nested" within the volume exclusion law M +1 i=1 n i = N , while no such nesting is present in the current example where the chemical conservation law is n 1 + 2n 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 = 4k 0 /(k 1 Ω).The parameter L contains information about all the rate constants of the system; it is proportional to the equilibrium The distribution is significantly shifted by excluded volume effects; this is for the case where the maximum number of molecules which can be put inside a compartment of unit volume is N = 100.
The reaction rate constants are k 0 = 0.001 and k 1 = 1, while the conservation law constant is k = 100.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. 7a); (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 lastly a phase in which volume exclusion leads to a change from sub-Poissonian to super-Poissonian statistics (for 7 ≤ L ≤ 11, see Fig. 7b).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. 7c); (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. 7d).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 [46,47] and Smoluchowski models [48,49] for lattice spacings that are neither too small nor too large [50].In the limit of small lattice spacing, the statistics of the RDME do not converge to those of the lattice-free approach [9,51] but it is possible in this case to derive a new convergent RDME called the CRDME which does not suffer from this issue [52].
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 In (c) we compute the difference between the CV in dilute and high volume exclusion conditions, ∆CV i , for both species.In (d) we compute the difference between the mean number of molecules in dilute and high volume exclusion conditions, ∆ n i .See text for discussion.
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 Brownian dynamics (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 dynamics agrees with that given by the bimolecular propensity in the CME (for a derivation see Appendix D of [53]).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 [54].
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 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.

FIG. 1 :
FIG. 1: Schematic illustrating the differences between the RDME and vRDME spatial modeling of the reaction A + A − − B. For the RDME (a), particles can react inside each of the 4 voxels and diffuse between neighbouring voxels.The red and blue circles denote particles of species A and B respectively.The particles are empty to denote that they occupy no volume (point particles) and can pass through each other.For the vRDME (b), the colour coding is the same except that we have also yellow circles denoting the "empty space" E. Each voxel is occupied by at most one particle.Bimolecular reactions occur between neighbouring voxels.Diffusion involves the switching of an empty space molecule and a chemical molecule between two neighbouring sites.On the right of both (a) and (b), is an illustration of what happens when the dimerisation reaction occurs.

FIG. 2 :
FIG.2:The Fano factor (a), coefficient of variation (b), and skewness (c) of a species X i in a chemical system without chemical conservation laws.The red and green lines correspond to the statistics predicted by assuming a finite molecular radius (as given by Eqs.(23)) and by assuming point particles (as given by Eqs.(22)), respectively.Volume exclusion effects become more appreciable as the occupied volume fraction of space tends to unity ( n i /N → 1) which causes the Fano factor to decrease (a), deviations from the inverse square root law for noise strength (b) and the skewness to switch from positive to negative (c).

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.

(
lower panel).The analytical formulae are compared with data from the Stochastic Simulation Algorithm (SSA, open circles) sampling the CME and the vCME, as evidence of their exactness.

FIG. 3 :
FIG. 3: Our theoretical predictions (lines) for the crowded and non-crowded models of the open homodimerisation reaction compared with data from the SSA of the CME and vCME (circles)

2 FIG. 4 :
FIG. 4: Comparisons of vRDME and RDME simulations with our theoretical predictions for the local and global distributions of molecule numbers of species X 1 and X 2 in the heterodimerisation system.Parameter values are k 0 = 20, k 1 = 1, k 2 = 1 and k 3 = 20 and the chemical conservation law is n 2 + n 3 = k = 15.The global compartment volume is Ω = 1 and the total number of voxels for both the RDME and vRDME is N = 49.In all cases there is excellent agreement between simulations and theory.
FIG. 5: Variation of the vRDME global distribution of molecule numbers for species X 1 in the open heterodimerisation reaction, as a function of the occupied volume fraction of space.The degree of volume exclusion is controlled by varying the size of the reactant molecules in (a) and by introducing inert molecules into the system in (b).Specifically (a) is obtained by keeping the compartment size constant and varying the maximum number of molecules N (voxels) which can be accommodated in the compartment for system (38).While (b) is obtained by varying the

FIG. 6 :
FIG.6:The steady-state probability distribution according to the CME (dilute conditions) and to the vCME (volume exclusion is taken into account).The plots are generated using Eqs.(43)-(44).

FIG. 7 :
FIG. 7: Comparison of statistics of intrinsic noise in high volume exclusion and dilute conditions for the closed reversible dimerisation reaction.The statistics are all numerically calculated from Eq. (44); the dilute case is obtained by choosing R = N/k = 1000 and the high volume exclusion case by choosing R = 1.The constant k is fixed to 50 in all cases.The non-dimensional parameter L which is an aggregate of all rate constants and the volume is varied and the statistics plotted as a function of L. In (a) and (b) we show the variation of the Fano factors of the two species.