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

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;³ artificial examples include reactions occurring inside nano-spaces such as nano-reactors⁴ and carbon nanotubes.⁵ 

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).¹ 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.²⁵ 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.

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 ⇌ k 1 k 0 B$⁠, where k₀ and k₁ 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 Ne(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 k_D = D/Δx² 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 Sec. III.

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 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.³² 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 ⇌ k ̃ 1 k ̃ 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.

We will here focus on reversible chemical systems in equilibrium conditions, i.e., those in detailed balance.³⁷ 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.

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³⁹ 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 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⁴⁰ 

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.

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 M 1 , … , n 1 N , … , 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 ⇌ X 2$ has the conservation law n₁ + 2n₂ = constant, which is defined on the total number of molecules of X₁ and X₂ 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, i.e.,$n i = ∑ j = 1 N n i j$⁠.

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 $ϕ ̃ 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.

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 , … , n M + 1 1 , … , n 1 N , … , 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., $∑ 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 (⁠$ϕ →$ 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,⁴⁰ 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.

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 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 $X 1 + X 1 ⇌ k 1 ′ k 1 X 2$ and $X 1 + X 1 ⇌ k ̃ 1 ′ k ̃ 1 X 2 + X 3$⁠, respectively (where X₃ is the empty space species), which implies the relation $k ̃ 1 ′ / k ̃ 1 = k 1 ′ Ω/ 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,³⁹ the effective equilibrium constant of the j-th reaction in system (3) in volume excluded and dilute conditions are, respectively, given by

Now by Eqs. (9) and (14), we then have

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 ̃ 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³⁹ 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 ⇌ X 2$⁠, this is modelled in the vCME formalism as $X 1 + X 1 ⇌ X 2 + X 3$ and hence by Le Chatelier’s principle, an increase in volume exclusion, i.e., a decrease in X₃ (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₂ and a decrease in the amount of species X₁ 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.²⁴ However it is the first time, to our knowledge, that they have been obtained using the deterministic limit of a master equation description.

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,Ω ϕ ̃ 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 ϕ_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.

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₁, n₂, …, 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₁, n₂, …, 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 M + 1 Ω ϕ ̃ i /N=1$⁠. It is well known that the marginal distributions of a multinomial distribution are binomial.⁴¹ For species X_i, the marginal distribution is thus binomial with parameters $( N , Ω ϕ ̃ i / N )$⁠,

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,⁴² 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).

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.

We consider the dilute (point particle) chemical system,

whereby a species X₁ is produced and subsequently two molecules of this species reversibly bind to form another molecule of type X₂. 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₃ 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, ϕ₃ → 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₀, 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₀ is strongly affected by volume exclusion, except of course in the limit of small k₀ 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₁ but a strong impact on the Fano factor of species X₂ (a change from constant to strongly monotonic decreasing as a function of k₀); (ii) in contrast, volume exclusion has a strong impact on the coefficient of variation of species X₁ (a change from a monotonic decreasing function to a parabolic function of k₀) but little impact on the coefficient of variation of X₂; and (iii) the skewness of species X₂ becomes negative as k₀ 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.

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−k,Ω ϕ ̃ i / ( N − 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,Ω ϕ ̃ 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.

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₁, …, 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 ({Ωϕ₁, …, Ωϕ_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_i 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 (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 $( { Ω ϕ ̃ 1 / ( N − k ) , … , Ω ϕ ̃ L / ( N − k ) , Ω ϕ ̃ M + 1 / ( N − k ) } , N − k )$ 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, $∑ i = 1 M + 1 n i =N$⁠, which together implies the combined conservation law $∑ i = 1 L n i + n M + 1 =N−k$⁠. 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.

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 , Ω ϕ ̃ / ( 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 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 − k − s , Ω ϕ ̃ i / ( N − k − s ) )$ if it is not involved in the conservation laws, 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.

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

whereby a species X₁ is produced and subsequently molecules of this species and that of X₂ reversibly bind to form molecules of X₃, a heterodimer. The system has the implicit chemical conservation law n₂ + n₃ = 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₄ is the empty space species and now we have two conservation laws: the chemical law n₂ + n₃ = k and the volume exclusion law n₁ + n₂ + n₃ + n₄ = 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 (39) together with the relevant conservation laws leads to

This implies that the concentrations of species X₂ and X₃ (the species involved in a chemical conservation law) are insensitive to volume exclusion effects but the concentration of species X₁ is found to decrease when crowding is taken into account. Intuitively, this is because species X₂ and X₃ are involved in a chemical conservation law and hence the impact of the second conservation law due to volume exclusion is nullified; species X₁ 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₁ is Poisson (Ωϕ₁) 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₂ is binomial (k, Ωϕ₂/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₁ 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₁, i.e., binomial $( N − k , Ω ϕ ̃ 1 / ( N − 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₁.

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.²³ To this end, we consider a modified version of reaction scheme (38),

where X₄ is the empty space species and X₅ 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₁, i.e., binomial $( N − k , Ω ϕ ̃ 1 / ( N − k ) )$⁠, as a function of the mean number of inert external crowder molecules 〈n₅〉. Note that the effect of increasing molecular crowding by adding more molecules of X₅ is to induce a shift of the probability distribution to the left such that there are fewer molecules, on average, of X₁ 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₂ is unchanged by the degree of volume exclusion since as we noted earlier both the RDME and vRDME give the same result.

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).

Consider the point particle model of the reaction system,

whereby two molecules of X₁ reversibly bind to form a dimer X₂. This system has the global chemical conservation law n₁ + 2n₂ = 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₂ 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⁴⁴). Note that for the vCME, we have here considered the volume excluded version of reaction scheme (42), namely, $X 1 + X 1 ⇌ k ̃ 1 k ̃ 0 X 2 + X 3$ with X₃ representing the empty space species and the equilibrium constant $k ̃ 0 k ̃ 1 = k 0 N k 1 Ω$ (as elucidated in Section IV). All the statistics of the molecule numbers of species X₁ can be deduced from those of X₂ given the conservation law n₁ + 2n₂ = 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 $∑ i = L + 1 M n i =k$ is “nested” within the volume exclusion law $∑ i = 1 M + 1 n i =N$⁠, while no such nesting is present in the current example where the chemical conservation law is n₁ + 2n₂ = k.