The model
In the present approach we consider proteins to be freely diffusing point-like objects. On the scale of a whole cell which is the scale we are interested in, long range forces between the objects are shielded by the solvent and can therefore safely be ignored. Diffusion is formally modelled by Brownian dynamics, taking intermolecular forces explicitly into account, and integrating over the velocities of the object [23]. In the absence of long-range forces, the Brownian dynamics treatment of diffusion reduces to a random walk process. The random walk process only considers the position of the objects and not their velocities thereby reducing considerably the computational cost; this approach was therefore used in this work. It is also assumed that any differences in reaction kinetics resulting from the different possible orientations of two reacting molecules relative to each other at the time of encounter can safely be integrated into an average reaction kinetic, so that the objects can be treated as point-like. Interactions between these point-like objects are governed by a set of reaction rules (described in detail below) that are designed to emulate the biological system of interest as closely as possible, as illustrated in Figure 2.
Figure 2  Possible applications of the method. An example of the kind of simulation this approach is designed for. This example illustrates a simulation of Schizosaccharomyces pombe yeast cell 10 μm long. Different types of proteins are shown in different colours, each has its own diffusion, reaction and location (nuclear or cytosolic) characteristics.

Reaction rules
We are interested in the reaction probability: the probability that two entities interact during a time step given that they can interact with reaction rate k, their diffusion rates are D1 and D2 and they start a distance d apart at the beginning of the time interval Δt. This probability is illustrated in Figure 1C. The reaction between two diffusing particles can be considered to occur in two steps: firstly the encounter of the two entities through diffusion, followed by the actual chemical reaction. Let us consider two freely diffusing chemical entities, A and B, starting a distance d0 from each other at time t0. At any time t later the rate κAB (t|d0, t0) of the reaction between entities A and B can be expressed as:
κ A B   ( t | d 0  , t 0  ) = p C A B   ( t | d 0  , t 0  ) ⋅ k R A B         ( 2 )    MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaaiiGacqWF6oWAdaahaaWcbeqaaiabdgeabjabdkeacbaakiabcIcaOiabdsha0jabcYha8jabdsgaKnaaBaaaleaacqaIWaamaeqaaOGaeiilaWIaemiDaq3aaSbaaSqaaiabicdaWaqabaGccqGGPaqkcqGH9aqpcqWGWbaCdaqhaaWcbaGaem4qameabaGaemyqaeKaemOqaieaaOGaeiikaGIaemiDaqNaeiiFaWNaemizaq2aaSbaaSqaaiabicdaWaqabaGccqGGSaalcqWG0baDdaWgaaWcbaGaeGimaadabeaakiabcMcaPiabgwSixlabdUgaRnaaDaaaleaacqWGsbGuaeaacqWGbbqqcqWGcbGqaaGccaWLjaGaaCzcamaabmaabaGaeGOmaidacaGLOaGaayzkaaaaaa@567F@
where pCAB MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGWbaCdaqhaaWcbaGaem4qameabaGaemyqaeKaemOqaieaaaaa@3169@ (t|d0, t0) is the probability of the two entities coming into contact at time t and kRAB MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGRbWAdaqhaaWcbaGaemOuaifabaGaemyqaeKaemOqaieaaaaa@317D@ is the rate of the reaction once in contact, averaged over all possible orientations of the two entities relative to each other. Both parts of this equation: pCAB MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGWbaCdaqhaaWcbaGaem4qameabaGaemyqaeKaemOqaieaaaaa@3169@ (t|d0, t0) and kRAB MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGRbWAdaqhaaWcbaGaemOuaifabaGaemyqaeKaemOqaieaaaaa@317D@ can be estimated as described below.
The reaction rate κAB (t|d0, t0) can be integrated over a simulation timestep Δt to provide the probability of at least one reaction taking place in that timestep. We are interested in the probability of at least one event taking place during the time interval Δt, i.e. 1 - P(no event during Δt). The process under consideration is a Poisson process with a time dependent rate of the event taking place. Given the rate κ(t) of an event taking place at a time t, the probability, PAB, of at least one reaction taking place during that time interval takes the general form [24,25]:
PAB (Δt) = 1 - e-I(Δt)     (3)
where
I ( Δ t ) = ∫ 0 Δ t   κ ( t ) d t          ( 4 )    MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGjbqscqGGOaakcqqHuoarcqWG0baDcqGGPaqkcqGH9aqpdaWdXaqaaGGaciab=P7aRjabcIcaOiabdsha0jabcMcaPiabbsgaKjabdsha0bWcbaGaeGimaadabaGaeuiLdqKaemiDaqhaniabgUIiYdGccaWLjaGaaCzcamaabmaabaGaeGinaqdacaGLOaGaayzkaaaaaa@44AD@
Such that the probability of a reaction taking place during timestep δt can be expressed as:
P A B   ( Δ t ) = 1 − e − ∫ 0 Δ t   κ A B   ( t ) d t            ( 5 )    MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGqbaudaahaaWcbeqaaiabdgeabjabdkeacbaakiabcIcaOiabfs5aejabdsha0jabcMcaPiabg2da9iabigdaXiabgkHiTiabdwgaLnaaCaaaleqabaGaeyOeI0Yaa8qmaeaaiiGacqWF6oWAdaahaaadbeqaaiabdgeabjabdkeacbaaliabcIcaOiabdsha0jabcMcaPiabbsgaKjabdsha0badbaGaeGimaadabaGaeuiLdqKaemiDaqhaoiabgUIiYdaaaOGaaCzcaiaaxMaadaqadaqaaiabiwda1aGaayjkaiaawMcaaaaa@4DA9@
where κAB (t) is given in equation (2).

The contact probability: pCAB MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGWbaCdaqhaaWcbaGaem4qameabaGaemyqaeKaemOqaieaaaaa@3169@ (t|d0, t0)
The probability of contact, pCAB MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGWbaCdaqhaaWcbaGaem4qameabaGaemyqaeKaemOqaieaaaaa@3169@ (t|d0, t0), is determined by the diffusion process of the two entities A and B during the time interval Δt = t - t0. The interacting bodies follow the laws of diffusion such that the probability of finding a given entity in an infinitesimal volume element dV, a distance d away from its starting position a time Δt later, is given by the well known Gaussian distribution: (4πDΔt)−3/2e−d24DΔtdV MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqGGOaakcqaI0aaniiGacqWFapaCcqWGebarcqqHuoarcqWG0baDcqGGPaqkdaahaaWcbeqaaiabgkHiTiabiodaZiabc+caViabikdaYaaakiabdwgaLnaaCaaaleqabaGaeyOeI0YaaSaaaeaacqWGKbazdaahaaadbeqaaiabikdaYaaaaSqaaiabisda0iabdseaejabfs5aejabdsha0baaaaGccqWGKbazcqWGwbGvaaa@4557@, where D is the diffusion constant of the entity [26].
The present approach is illustrated in Figure 3. The two entities diffuse freely starting a distance d0 apart. The probability of them coming into contact increases with time and reaches a maximum. Subsequently the two entities diffuse further and the probability of them coming into contact decreases with time. A mathematically equivalent description is given if A is considered to be diffusing with diffusion constant D+ = DA + DB while B remains stationary. It follows that the probability of contact, pCAB MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGWbaCdaqhaaWcbaGaem4qameabaGaemyqaeKaemOqaieaaaaa@3169@ (t|r0, t0), is given by:
Figure 3  Probability density functions. The diffusion of the two entities A and B (red and blue) gives rise to a certain probability of them coming into contact (green). The two entities diffuse freely starting a distance d apart. As time goes by the probability of them coming into contact increases and reaches a maximum. Subsequently the two entities diffuse further and the probability of them coming into contact decreases with time (note: this is a 1D projection of the probability density profile, in 3D the integral over the probability density is correctly normalised to 1). p C A B   ( t | d 0  , t 0  ) = ( 4 π D +  Δ t )  − 3 / 2   e − d 0 2   4 D +  Δ t     δ V C        ( 6 )    MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGWbaCdaqhaaWcbaGaem4qameabaGaemyqaeKaemOqaieaaOGaeiikaGIaemiDaqNaeiiFaWNaemizaq2aaSbaaSqaaiabicdaWaqabaGccqGGSaalcqWG0baDdaWgaaWcbaGaeGimaadabeaakiabcMcaPiabg2da9iabcIcaOiabisda0GGaciab=b8aWjabdseaenaaBaaaleaacqGHRaWkaeqaaOGaeuiLdqKaemiDaqNaeiykaKYaaWbaaSqabeaacqGHsislcqaIZaWmcqGGVaWlcqaIYaGmaaGccqWGLbqzdaahaaWcbeqaaiabgkHiTmaalaaabaGaemizaq2aa0baaWqaaiabicdaWaqaaiabikdaYaaaaSqaaiabisda0iabdseaenaaBaaameaacqGHRaWkaeqaaSGaeuiLdqKaemiDaqhaaaaakiab=r7aKjabdAfawnaaBaaaleaacqWGdbWqaeqaaOGaaCzcaiaaxMaadaqadaqaaiabiAda2aGaayjkaiaawMcaaaaa@5E2C@
where D+ = DA + DB, Δt = t - t0 and δVC is a small contact volume defined such that if two entities are found to be within this volume, they are considered to be in contact. This small contact volume, δVC, will be considered further below.

The reaction rate: kRAB MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGRbWAdaqhaaWcbaGaemOuaifabaGaemyqaeKaemOqaieaaaaa@317D@
In order to get a good first approximation for pCAB MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGWbaCdaqhaaWcbaGaem4qameabaGaemyqaeKaemOqaieaaaaa@3169@ (t|d0, t0), we initially consider the well-mixed limit. In this limit the distribution of the two entities A and B is uniform over space. This approximation is thus only valid in the long timestep limit. For short timesteps, the approximations break down and a correction to the reaction rate, kRAB MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGRbWAdaqhaaWcbaGaemOuaifabaGaemyqaeKaemOqaieaaaaa@317D@, has to be introduced. An exact analytical solution has recently been presented that is also correct for short time steps [27]. However, the mathematics involved are complex and difficult to implement; the approach used here, although approximate, is simpler and is not expected to alter the findings significantly.

The well-mixed limit
Let us consider the simple reaction: A+B→kY MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqqGbbqqcqGHRaWkcqqGcbGqdaGdKaWcbaGaem4AaSgabeGccaGLsgcacqqGzbqwaaa@33D2@ occurring in a finite volume V. To first order, the rate of change in the number of molecules Y is:
d n Y   d t   = k N A    n A  n B   V        ( 7 )    MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaadaWcaaqaaiabbsgaKjabd6gaUnaaBaaaleaacqqGzbqwaeqaaaGcbaGaeeizaqMaemiDaqhaaiabg2da9maalaaabaGaem4AaSgabaGaemOta40aaSbaaSqaaiabdgeabbqabaaaaOWaaSaaaeaacqWGUbGBdaWgaaWcbaGaeeyqaeeabeaakiabd6gaUnaaBaaaleaacqqGcbGqaeqaaaGcbaGaemOvayfaaiaaxMaacaWLjaWaaeWaaeaacqaI3aWnaiaawIcacaGLPaaaaaa@42D2@
where nA, nB and nY are the number of molecules of A, B and Y (respectively), NA is Avogadro's constant and k is rate of the reaction. The rate of change in nY can also be expressed as:
d n Y   d t   = p ^  C A B   ⋅ k R A B         ( 8 )    MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaadaWcaaqaaiabbsgaKjabd6gaUnaaBaaaleaacqqGzbqwaeqaaaGcbaGaeeizaqMaemiDaqhaaiabg2da9iqbdchaWzaajaWaa0baaSqaaiabdoeadbqaaiabdgeabjabdkeacbaakiabgwSixlabdUgaRnaaDaaaleaacqWGsbGuaeaacqWGbbqqcqWGcbGqaaGccaWLjaGaaCzcamaabmaabaGaeGioaGdacaGLOaGaayzkaaaaaa@446C@
where p^CAB MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacuWGWbaCgaqcamaaDaaaleaacqWGdbWqaeaacqWGbbqqcqWGcbGqaaaaaa@3179@ is here the ensemble average probability of any two entities A and B being in contact in volume V and kRAB MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGRbWAdaqhaaWcbaGaemOuaifabaGaemyqaeKaemOqaieaaaaa@317D@ the rate of reaction if they are in contact.
The objects are considered to be uniformly distributed over the volume V such that the probability of A and B occupying the same contact volume δVC, p^CAB MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacuWGWbaCgaqcamaaDaaaleaacqWGdbWqaeaacqWGbbqqcqWGcbGqaaaaaa@3179@, can be expressed as:
p ^  C A B   = n A  n B   V  δ V C        ( 9 )    MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacuWGWbaCgaqcamaaDaaaleaacqWGdbWqaeaacqWGbbqqcqWGcbGqaaGccqGH9aqpdaWcaaqaaiabd6gaUnaaBaaaleaacqqGbbqqaeqaaOGaemOBa42aaSbaaSqaaiabbkeacbqabaaakeaacqWGwbGvaaacciGae8hTdqMaemOvay1aaSbaaSqaaiabdoeadbqabaGccaWLjaGaaCzcamaabmaabaGaeGyoaKdacaGLOaGaayzkaaaaaa@410B@
where δVC is the same infinitesimal volume as in equation (6). By combining equations (7), (8), and (9), we can extract the rate of the reaction if A and B are in contact:
k R , mixed  A B   = k N A  δ V C          ( 10  )    MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGRbWAdaqhaaWcbaGaemOuaiLaeiilaWIaeeyBa0MaeeyAaKMaeeiEaGNaeeyzauMaeeizaqgabaGaemyqaeKaemOqaieaaOGaeyypa0ZaaSaaaeaacqWGRbWAaeaacqWGobGtdaWgaaWcbaGaemyqaeeabeaaiiGakiab=r7aKjabdAfawnaaBaaaleaacqWGdbWqaeqaaaaakiaaxMaacaWLjaWaaeWaaeaacqaIXaqmcqaIWaamaiaawIcacaGLPaaaaaa@46E4@
Combining equation (10) and pCAB MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGWbaCdaqhaaWcbaGaem4qameabaGaemyqaeKaemOqaieaaaaa@3169@ from equation (6) into equation (2), the rate of the reaction between entities A and B, starting a distance d0 apart at a time t later, is given by:
κ mixed  A B   ( t | d 0  , t 0  ) = ( 4 π D +  Δ t )  − 3 / 2   e − d 0 2   4 D +  Δ t     ⋅ k N A          ( 11  )    MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaaiiGacqWF6oWAdaqhaaWcbaGaeeyBa0MaeeyAaKMaeeiEaGNaeeyzauMaeeizaqgabaGaemyqaeKaemOqaieaaOGaeiikaGIaemiDaqNaeiiFaWNaemizaq2aaSbaaSqaaiabicdaWaqabaGccqGGSaalcqWG0baDdaWgaaWcbaGaeGimaadabeaakiabcMcaPiabg2da9iabcIcaOiabisda0iab=b8aWjabdseaenaaBaaaleaacqGHRaWkaeqaaOGaeuiLdqKaemiDaqNaeiykaKYaaWbaaSqabeaacqGHsislcqaIZaWmcqGGVaWlcqaIYaGmaaGccqWGLbqzdaahaaWcbeqaaiabgkHiTmaalaaabaGaemizaq2aa0baaWqaaiabicdaWaqaaiabikdaYaaaaSqaaiabisda0iabdseaenaaBaaameaacqGHRaWkaeqaaSGaeuiLdqKaemiDaqhaaaaakiabgwSixpaalaaabaGaem4AaSgabaGaemOta40aaSbaaSqaaiabdgeabbqabaaaaOGaaCzcaiaaxMaadaqadaqaaiabigdaXiabigdaXaGaayjkaiaawMcaaaaa@671D@
In doing this, notice that the contact volume δVC cancels out of the equations. This effectively removes any information about the size of the particles from subsequent considerations.
Inserting kmixedAB MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGRbWAdaqhaaWcbaGaeeyBa0MaeeyAaKMaeeiEaGNaeeyzauMaeeizaqgabaGaemyqaeKaemOqaieaaaaa@3721@ from equation (11) into equation (4), I(Δt) can be solved analytically using the standard integral:
∫ 0 Δ t   t − 3 2       e − 1 t    d t = π   Erfc ( 1 Δ t      )        ( 12  )    MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaadaWdXaqaaiabdsha0naaCaaaleqabaGaeyOeI0YaaSaaaeaacqaIZaWmaeaacqaIYaGmaaaaaaqaaiabicdaWaqaaiabfs5aejabdsha0bqdcqGHRiI8aOGaemyzau2aaWbaaSqabeaacqGHsisldaWcaaqaaiabigdaXaqaaiabdsha0baaaaGccqqGKbazcqWG0baDcqGH9aqpdaGcaaqaaGGaciab=b8aWbWcbeaakiabbccaGiabbweafjabbkhaYjabbAgaMjabbogaJnaabmaabaWaaSaaaeaacqaIXaqmaeaadaGcaaqaaiabfs5aejabdsha0bWcbeaaaaaakiaawIcacaGLPaaacaWLjaGaaCzcamaabmaabaGaeGymaeJaeGOmaidacaGLOaGaayzkaaaaaa@51C8@
where Erfc is the complimentary error function defined by:
Erfc ( x ) = 2 π  ∫ x ∞  e − z 2       d z       ( 13  )    MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqqGfbqrcqqGYbGCcqqGMbGzcqqGJbWycqGGOaakcqWG4baEcqGGPaqkcqGH9aqpdaWcaaqaaiabikdaYaqaaGGaciab=b8aWbaadaWdXaqaaiabdwgaLnaaCaaaleqabaGaeyOeI0IaemOEaO3aaWbaaWqabeaacqaIYaGmaaaaaaWcbaGaemiEaGhabaGaeyOhIukaniabgUIiYdGccqqGKbazcqWG6bGEcaWLjaGaaCzcamaabmaabaGaeGymaeJaeG4mamdacaGLOaGaayzkaaaaaa@4A63@
such that I(Δt) has the analytical form:
I ( Δ t ) = k N A    ⋅ 1 4 d π D +    Erfc ( d 4 Δ t D +       )        ( 14  )    MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGjbqscqGGOaakcqqHuoarcqWG0baDcqGGPaqkcqGH9aqpdaWcaaqaaiabdUgaRbqaaiabd6eaonaaBaaaleaacqWGbbqqaeqaaaaakiabgwSixpaalaaabaGaeGymaedabaGaeGinaqJaemizaqgcciGae8hWdaNaemiraq0aaSbaaSqaaiabgUcaRaqabaaaaOGaeeyrauKaeeOCaiNaeeOzayMaee4yam2aaeWaaeaadaWcaaqaaiabdsgaKbqaamaakaaabaGaeGinaqJaeuiLdqKaemiDaqNaemiraq0aaSbaaSqaaiabgUcaRaqabaaabeaaaaaakiaawIcacaGLPaaacaWLjaGaaCzcamaabmaabaGaeGymaeJaeGinaqdacaGLOaGaayzkaaaaaa@5368@
Finally we can express the probability PmixedAB MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGqbaudaqhaaWcbaGaeeyBa0MaeeyAaKMaeeiEaGNaeeyzauMaeeizaqgabaGaemyqaeKaemOqaieaaaaa@36EB@ of a reaction taking place between entities A and B, starting a distance d0 apart, during the time interval Δt as:
P mixed  A B   ( t | d 0  , t 0  ) = 1 − e − k N A    ⋅ 1 4 d 0  π D +    Erfc ( d 0   4 Δ t D +       )          ( 15  )    MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGqbaudaqhaaWcbaGaeeyBa0MaeeyAaKMaeeiEaGNaeeyzauMaeeizaqgabaGaemyqaeKaemOqaieaaOGaeiikaGIaemiDaqNaeiiFaWNaemizaq2aaSbaaSqaaiabicdaWaqabaGccqGGSaalcqWG0baDdaWgaaWcbaGaeGimaadabeaakiabcMcaPiabg2da9iabigdaXiabgkHiTiabdwgaLnaaCaaaleqabaGaeyOeI0YaaSaaaeaacqWGRbWAaeaacqWGobGtdaWgaaadbaGaemyqaeeabeaaaaWccqGHflY1daWcaaqaaiabigdaXaqaaiabisda0iabdsgaKnaaBaaameaacqaIWaamaeqaaGGacSGae8hWdaNaemiraq0aaSbaaWqaaiabgUcaRaqabaaaaSGaeeyrauKaeeOCaiNaeeOzayMaee4yam2aaeWaaeaadaWcaaqaaiabdsgaKnaaBaaameaacqaIWaamaeqaaaWcbaWaaOaaaeaacqaI0aancqqHuoarcqWG0baDcqWGebardaWgaaadbaGaey4kaScabeaaaeqaaaaaaSGaayjkaiaawMcaaaaakiaaxMaacaWLjaWaaeWaaeaacqaIXaqmcqaI1aqnaiaawIcacaGLPaaaaaa@6942@
The probability of two chemical entities interacting in a timestep Δt is thus expressed in equation (15) in terms of the reaction rate k, the sum of the diffusion constants of the two entities D+ and the time interval Δt. This approach provides a good description of the interactions in terms of the underlying diffusion process and reactivity of the two entities.

Short timestep correction
The equations above hold for situations where the system can be considered to be well-mixed. However, this assumption breaks down for small timesteps: as chemical entities react over time, there tend to be fewer potentially interacting partners close to each other so that the distribution of the two entities is no longer uniform. In the long timescale limit this is not a problem as the system is well-mixed by each diffusion step and the approximations hold. At each timestep, the reaction process creates 'dips' in the probability distribution of the entities, the spatial extent of these 'dips' is comparable to the spatial extent of the probability of reaction. In order to remain well-mixed the distance covered by one step of diffusion must be greater than the spatial extent of the 'dips' created by the reaction process. For diffusion constants typical of biomolecules, the spatial width at half-maximum of PmixedAB MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGqbaudaqhaaWcbaGaeeyBa0MaeeyAaKMaeeiEaGNaeeyzauMaeeizaqgabaGaemyqaeKaemOqaieaaaaa@36EB@ (t|d0, t0) goes to ~0.1 μm for Δt of the order of seconds. Considering this distance as being covered by diffusion, this gives us a typical timescale of Δt ≥ 0.01 s. The system can therefore be assumed to be well-mixed for timesteps of Δt ≥ 0.01 s. For shorter timesteps, this effect can be corrected for, as will be shown below.
Due to the reaction process, the average concentration around a chemical entity is less than predicted by the uniform distribution. The desired rate can be derived by correcting kR,mixedAB MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGRbWAdaqhaaWcbaGaemOuaiLaeiilaWIaeeyBa0MaeeyAaKMaeeiEaGNaeeyzauMaeeizaqgabaGaemyqaeKaemOqaieaaaaa@392E@ by a scaling factor as follows:
k R A B   = C ⋅ k R , mixed  A B         ( 16  )    MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGRbWAdaqhaaWcbaGaemOuaifabaGaemyqaeKaemOqaieaaOGaeyypa0Jaem4qamKaeyyXICTaem4AaS2aa0baaSqaaiabdkfasjabcYcaSiabb2gaTjabbMgaPjabbIha4jabbwgaLjabbsgaKbqaaiabdgeabjabdkeacbaakiaaxMaacaWLjaWaaeWaaeaacqaIXaqmcqaI2aGnaiaawIcacaGLPaaaaaa@4729@
The procedure we used for doing this is very similar to that used by Andrews and Bray [17] to correct for the same effect in the Smoluchowski approach and is outlined below. More elaborate mathematical considerations of this process can be found in the recent paper by Zon and Wolde 2005 [27].
Using the rate of reaction upon encounter from equation (16), the probability of the reaction taking place after each diffusion step is now given by:
P A B   ( t 0  + Δ t | d 0  , t 0  ) = 1 − e − C ⋅ k N A    ⋅ 1 4 d 0  π D +    Erfc ( d 0   2 s    )          ( 17  )    MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGqbaudaahaaWcbeqaaiabdgeabjabdkeacbaakiabcIcaOiabdsha0naaBaaaleaacqaIWaamaeqaaOGaey4kaSIaeuiLdqKaemiDaqNaeiiFaWNaemizaq2aaSbaaSqaaiabicdaWaqabaGccqGGSaalcqWG0baDdaWgaaWcbaGaeGimaadabeaakiabcMcaPiabg2da9iabigdaXiabgkHiTiabdwgaLnaaCaaaleqabaGaeyOeI0Iaem4qamKaeyyXIC9aaSaaaeaacqWGRbWAaeaacqWGobGtdaWgaaadbaGaemyqaeeabeaaaaWccqGHflY1daWcaaqaaiabigdaXaqaaiabisda0iabdsgaKnaaBaaameaacqaIWaamaeqaaGGacSGae8hWdaNaemiraq0aaSbaaWqaaiabgUcaRaqabaaaaSGaeeyrauKaeeOCaiNaeeOzayMaee4yam2aaeWaaeaadaWcaaqaaiabdsgaKnaaBaaameaacqaIWaamaeqaaaWcbaGaeGOmaiJaem4CamhaaaGaayjkaiaawMcaaaaakiaaxMaacaWLjaWaaeWaaeaacqaIXaqmcqaI3aWnaiaawIcacaGLPaaaaaa@6704@
where Δt is the timestep of the simulation and we use the substitution s=2D+Δt MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGZbWCcqGH9aqpdaGcaaqaaiabikdaYiabdseaenaaBaaaleaacqGHRaWkaeqaaOGaeuiLdqKaemiDaqhaleqaaaaa@352E@.
For the purposes of the correction, we are interested in average effects, so from now on we consider the average concentrations of entities A and B and not the positions of entities A and B. Let us consider the radial concentration ρB(r, t) of entity B around entity A, with A considered to be static at r = 0, while entity B has diffusion constant D+ = DA + DB.
The radial concentration ρB(r, t) of entity B around entity A, at time t is propagated for a simulation timestep, Δt, to give ρB(r, t + Δt) using a Green's function [28]:
ρ B  ( r , t + Δ t ) = ∫ 0 ∞  ρ B     ( r ′  , t ) G s  ( r , r ′  ) 4 π r ′  2  d r ′        ( 18  )    MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaaiiGacqWFbpGCdaWgaaWcbaGaemOqaieabeaakiabcIcaOiabdkhaYjabcYcaSiabdsha0jabgUcaRiabfs5aejabdsha0jabcMcaPiabg2da9maapedabaGae8xWdi3aaSbaaSqaaiabdkeacbqabaaabaGaeGimaadabaGaeyOhIukaniabgUIiYdGccqGGOaakcuWGYbGCgaqbaiabcYcaSiabdsha0jabcMcaPiabdEeahnaaBaaaleaacqWGZbWCaeqaaOGaeiikaGIaemOCaiNaeiilaWIafmOCaiNbauaacqGGPaqkcqaI0aancqWFapaCcuWGYbGCgaqbamaaCaaaleqabaGaeGOmaidaaOGaeeizaqMafmOCaiNbauaacaWLjaGaaCzcamaabmaabaGaeGymaeJaeGioaGdacaGLOaGaayzkaaaaaa@5BE3@
where the Green's function Gs(r, r') is given by:
G s  ( r , r ′  ) = 1 4 π r r ′  s 2 π     ( e − ( r − r ′  )  2   2 s 2      − e − ( r + r ′  )  2   2 s 2      )       ( 19  )    MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGhbWrdaWgaaWcbaGaem4CamhabeaakiabcIcaOiabdkhaYjabcYcaSiqbdkhaYzaafaGaeiykaKIaeyypa0ZaaSaaaeaacqaIXaqmaeaacqaI0aaniiGacqWFapaCcqWGYbGCcuWGYbGCgaqbaiabdohaZnaakaaabaGaeGOmaiJae8hWdahaleqaaaaakiabcIcaOiabdwgaLnaaCaaaleqabaGaeyOeI0YaaSaaaeaacqGGOaakcqWGYbGCcqGHsislcuWGYbGCgaqbaiabcMcaPmaaCaaameqabaGaeGOmaidaaaWcbaGaeGOmaiJaem4Cam3aaWbaaWqabeaacqaIYaGmaaaaaaaakiabgkHiTiabdwgaLnaaCaaaleqabaGaeyOeI0YaaSaaaeaacqGGOaakcqWGYbGCcqGHRaWkcuWGYbGCgaqbaiabcMcaPmaaCaaameqabaGaeGOmaidaaaWcbaGaeGOmaiJaem4Cam3aaWbaaWqabeaacqaIYaGmaaaaaaaakiabcMcaPiaaxMaacaWLjaWaaeWaaeaacqaIXaqmcqaI5aqoaiaawIcacaGLPaaaaaa@6185@
The entity A is then allowed to interact with entity B such that the new concentration of B is given by:
ρ ′  B  ( r , t + Δ t ) = ( 1 − P A B   ( r ) ) ⋅ ρ B  ( r , t + Δ t )       ( 20  )    MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaaiiGacuWFbpGCgaqbamaaBaaaleaacqWGcbGqaeqaaOGaeiikaGIaemOCaiNaeiilaWIaemiDaqNaey4kaSIaeuiLdqKaemiDaqNaeiykaKIaeyypa0JaeiikaGIaeGymaeJaeyOeI0Iaemiuaa1aaWbaaSqabeaacqWGbbqqcqWGcbGqaaGccqGGOaakcqWGYbGCcqGGPaqkcqGGPaqkcqGHflY1cqWFbpGCdaWgaaWcbaGaemOqaieabeaakiabcIcaOiabdkhaYjabcYcaSiabdsha0jabgUcaRiabfs5aejabdsha0jabcMcaPiaaxMaacaWLjaWaaeWaaeaacqaIYaGmcqaIWaamaiaawIcacaGLPaaaaaa@5735@
where PAB (r) is the probability of A and B interacting in the following timestep from equation (17).
The reaction step acts as a sink for the concentration of B, while the concentration of B is assumed to be constant at r = ∞. The long-distance equilibrium solution for ρB(r) is known to be of the form [26,29]:
ρ B  ( r ) = 1 + a r        ( 21  )    MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaaiiGacqWFbpGCdaWgaaWcbaGaemOqaieabeaakiabcIcaOiabdkhaYjabcMcaPiabg2da9iabigdaXiabgUcaRmaalaaabaGaemyyaegabaGaemOCaihaaiaaxMaacaWLjaWaaeWaaeaacqaIYaGmcqaIXaqmaiaawIcacaGLPaaaaaa@3D24@
This allows us to solve numerically for ρB(r, t + Δt) around r = 0 while using an analytical extension for long distances (long distance was defined as r > rP, rP such that PAB (r) = 10-6). Equation (18) is then split into a numerical and an analytical part:
ρ B  ( r , t + Δ t ) = ∫ 0 r P    ρ B     ( r ′  , t ) G s  ( r , r ′  ) 4 π r ′  2  d r ′  + ∫ r P   ∞  ρ B ana      ( r ′  , t ) G s  ( r , r ′  ) 4 π r ′  2  d r ′        ( 22  )    MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaaiiGacqWFbpGCdaWgaaWcbaGaemOqaieabeaakiabcIcaOiabdkhaYjabcYcaSiabdsha0jabgUcaRiabfs5aejabdsha0jabcMcaPiabg2da9maapedabaGae8xWdi3aaSbaaSqaaiabdkeacbqabaaabaGaeGimaadabaGaemOCai3aaSbaaWqaaiabdcfaqbqabaaaniabgUIiYdGccqGGOaakcuWGYbGCgaqbaiabcYcaSiabdsha0jabcMcaPiabdEeahnaaBaaaleaacqWGZbWCaeqaaOGaeiikaGIaemOCaiNaeiilaWIafmOCaiNbauaacqGGPaqkcqaI0aancqWFapaCcuWGYbGCgaqbamaaCaaaleqabaGaeGOmaidaaOGaeeizaqMafmOCaiNbauaacqGHRaWkdaWdXaqaaiab=f8aYnaaDaaaleaacqWGcbGqaeaacqqGHbqycqqGUbGBcqqGHbqyaaaabaGaemOCai3aaSbaaWqaaiabdcfaqbqabaaaleaacqGHEisPa0Gaey4kIipakiabcIcaOiqbdkhaYzaafaGaeiilaWIaemiDaqNaeiykaKIaem4raC0aaSbaaSqaaiabdohaZbqabaGccqGGOaakcqWGYbGCcqGGSaalcuWGYbGCgaqbaiabcMcaPiabisda0iab=b8aWjqbdkhaYzaafaWaaWbaaSqabeaacqaIYaGmaaGccqqGKbazcuWGYbGCgaqbaiaaxMaacaWLjaWaaeWaaeaacqaIYaGmcqaIYaGmaiaawIcacaGLPaaaaaa@8126@
ρ B  ( r , t + Δ t ) = ∫ 0 r P    ρ B     ( r ′  , t ) G s  ( r , r ′  ) 4 π r ′  2  d r ′  + ρ B ana   ( r , t + Δ t )       ( 23  )    MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaaiiGacqWFbpGCdaWgaaWcbaGaemOqaieabeaakiabcIcaOiabdkhaYjabcYcaSiabdsha0jabgUcaRiabfs5aejabdsha0jabcMcaPiabg2da9maapedabaGae8xWdi3aaSbaaSqaaiabdkeacbqabaaabaGaeGimaadabaGaemOCai3aaSbaaWqaaiabdcfaqbqabaaaniabgUIiYdGccqGGOaakcuWGYbGCgaqbaiabcYcaSiabdsha0jabcMcaPiabdEeahnaaBaaaleaacqWGZbWCaeqaaOGaeiikaGIaemOCaiNaeiilaWIafmOCaiNbauaacqGGPaqkcqaI0aancqWFapaCcuWGYbGCgaqbamaaCaaaleqabaGaeGOmaidaaOGaeeizaqMafmOCaiNbauaacqGHRaWkcqWFbpGCdaqhaaWcbaGaemOqaieabaGaeeyyaeMaeeOBa4MaeeyyaegaaOGaeiikaGIaemOCaiNaeiilaWIaemiDaqNaey4kaSIaeuiLdqKaemiDaqNaeiykaKIaaCzcaiaaxMaadaqadaqaaiabikdaYiabiodaZaGaayjkaiaawMcaaaaa@6E2C@
By inserting (21), the analytical extension, ρBana MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaaiiGacqWFbpGCdaqhaaWcbaGaemOqaieabaGaeeyyaeMaeeOBa4Maeeyyaegaaaaa@33A2@ (r, t + Δt), can be derived and is given by:
ρ B ana   ( r , t + Δ t ) = s r 2 π     G s  ( r , r P  ) + 1 2  ( E +  + E −  ) + a 2 r   ( E −  − E +  )       ( 24  )    MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaaiiGacqWFbpGCdaqhaaWcbaGaemOqaieabaGaeeyyaeMaeeOBa4MaeeyyaegaaOGaeiikaGIaemOCaiNaeiilaWIaemiDaqNaey4kaSIaeuiLdqKaemiDaqNaeiykaKIaeyypa0ZaaSaaaeaacqWGZbWCaeaacqWGYbGCdaGcaaqaaiabikdaYiab=b8aWbWcbeaaaaGccqWGhbWrdaWgaaWcbaGaem4CamhabeaakiabcIcaOiabdkhaYjabcYcaSiabdkhaYnaaBaaaleaacqWGqbauaeqaaOGaeiykaKIaey4kaSYaaSaaaeaacqaIXaqmaeaacqaIYaGmaaGaeiikaGIaemyrau0aaSbaaSqaaiabgUcaRaqabaGccqGHRaWkcqWGfbqrdaWgaaWcbaGaeyOeI0cabeaakiabcMcaPiabgUcaRmaalaaabaGaemyyaegabaGaeGOmaiJaemOCaihaaiabcIcaOiabdweafnaaBaaaleaacqGHsislaeqaaOGaeyOeI0Iaemyrau0aaSbaaSqaaiabgUcaRaqabaGccqGGPaqkcaWLjaGaaCzcamaabmaabaGaeGOmaiJaeGinaqdacaGLOaGaayzkaaaaaa@6737@
where: E±=Erfc(rP±rs2) MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGfbqrdaWgaaWcbaGaeyySaelabeaakiabg2da9iabbweafjabbkhaYjabbAgaMjabbogaJnaabmaabaWaaSaaaeaacqWGYbGCdaWgaaWcbaGaemiuaafabeaakiabgglaXkabdkhaYbqaaiabdohaZnaakaaabaGaeGOmaidaleqaaaaaaOGaayjkaiaawMcaaaaa@404B@
A diffusion step is performed using numerical integration for the remainder of equation (23). The value of the constant a in ρBana MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaaiiGacqWFbpGCdaqhaaWcbaGaemOqaieabaGaeeyyaeMaeeOBa4Maeeyyaegaaaaa@33A2@ (r, t) is found by fitting the last 10% of ρB(r, t) relative to rP after each diffusion step. After each diffusion step the two entities are allowed to interact following the probability given in equation (17). In order to achieve a steady state, entity A is left unaffected by the reaction process. The process of diffusion/reaction is repeated until the radial concentration ρB(r, t) reaches a steady state: ρB(r, ∞).
The effective rate of the reaction is given by:
k eff   = ∫ 0 ∞  P A B      ( r ) ρ B  ( r , ∞ ) 4 π r 2  d r       ( 25  )    MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGRbWAdaWgaaWcbaGaeeyzauMaeeOzayMaeeOzaygabeaakiabg2da9maapedabaGaemiuaa1aaWbaaSqabeaacqWGbbqqcqWGcbGqaaaabaGaeGimaadabaGaeyOhIukaniabgUIiYdGccqGGOaakcqWGYbGCcqGGPaqkiiGacqWFbpGCdaWgaaWcbaGaemOqaieabeaakiabcIcaOiabdkhaYjabcYcaSiabg6HiLkabcMcaPiabisda0iab=b8aWjabdkhaYnaaCaaaleqabaGaeGOmaidaaOGaeeizaqMaemOCaiNaaCzcaiaaxMaadaqadaqaaiabikdaYiabiwda1aGaayjkaiaawMcaaaaa@5383@
The values of keff were determined for a range of values of the correction factor, C, and variable s. Using this array of data it is then possible to find the correct value of the correction factor, C, for any pair of values of the desired reaction rate and simulation timestep: (keff, Δt).
Figure 4 shows examples of the probability distribution and corresponding reaction radius, σb, using the Smoluchowski approach for different timesteps. As can be seen, at large timesteps, the distribution of B before the reaction is close to uniform and the correction factor correspondingly small. On the other hand at small timesteps, the probability distribution of B, before reaction, is far from uniform and the correction factor, C, is large. For k = 106 M-1s-1 and DA = DB = 1 μm2s-1 the correction factors, C, are 1.12, 1.46, 4.48 and 4.44 103 for timesteps of 10-1s, 10-2s, 10-3s and 10-4s, respectively. Figure 4 also illustrates the fact that as the timesteps diminish the corrected reaction probability converges with the Smoluchowski cutoff at σb. For short timesteps the two approaches appear to be equivalent.
Figure 4  Radial probability densities and the correction factor C. The radial probability distribution, ρB, before interaction with A (green), and after (blue). The blue distribution evolves to the green distribution after one step of diffusion. The probability of interaction is shown in black. The probability of interaction as given by the corrected Smoluchowski approach of Andrews and Bray is also shown in red (P(r) = 1 for r <σb). Four timesteps are shown: Δt = 0.1s, 0.01s, 0.001s, 0.0001s. The correction factors and corrected Smoluchowski binding radii, σb, corresponding to the different timesteps are also shown, k = 106 M-1s-1, and DA = DB = 1 μm2s-1. Figure 4 also shows the radial probability distribution of the reactants for the different timesteps. These distributions are continuous throughout the reaction region. In contrast, those reported by Andrews and Bray show a strong discontinuity due to the sharply defined reaction radius [17].

The challenge of long timesteps
Relation (15) holds only for a pair of interacting molecules. In the more general context of an actual chemical reaction, the number of potential reactive partners at any one time can be much higher. The assumption is that during a time step Δt the probability of an entity interacting with its closest neighbour, Pclosest, is much higher than the probability of it interacting with its next nearest neighbour, Pnext nearest: Pclosest ≫ Pnext nearest.
This condition will clearly be fulfilled if the average distance, dtravel, travelled by an entity during the time interval Δt is less than the average distance between particles, d. The average travel distance is: dtravel = dtravel=6DΔt MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGKbazdaWgaaWcbaGaeeiDaqNaeeOCaiNaeeyyaeMaeeODayNaeeyzauMaeeiBaWgabeaakiabg2da9maakaaabaGaeGOnayJaemiraqKaeuiLdqKaemiDaqhaleqaaaaa@3C7C@, where D is the diffusion constant. The average interparticle distance given by: d=(VN)1/3 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGKbazcqGH9aqpcqGGOaakdaWcaaqaaiabdAfawbqaaiabd6eaobaacqGGPaqkdaahaaWcbeqaaiabigdaXiabc+caViabiodaZaaaaaa@3616@, where V is the volume and N is the number of particle, can also be expressed in terms of the concentration as: d ∝ C-1/3. where C is the concentration. For typical biological situations, C ≃ μM ml-1 and D ≃ μm2s-1, such that the condition reduces to Δt <~0.002 s.
When this condition is not fulfilled, a given entity can in principle interact with several other entities at any given timestep. This will affect the probability of the reaction with any given particle as reactions are considered mutually exclusive. In principle these probabilities can be calculated at each timestep during the simulation so that the correct statistics are reproduced. However, it was decided that the resulting extension to the code would introduce extra computational overhead, without significant benefit, and was not implemented. On the other hand, simply assuming that entities can interact with at most one other entity in the following timestep, when in fact they can interact with several, leads to the simulated reaction taking place at a rate greater than the expected rate.
Another problem which emerges for long timesteps concerns boundaries: for reactions, the algorithm assumes free diffusion in the space around the entities. This assumption is mostly correct when the timestep is small but at large timestep, the chance of encountering a boundary during that timestep become significant. At that point the free diffusion assumption assumed in the previous equations breaks down leading the reaction happening too fast in the simulation. Simply put, entities close to boundaries have less volume in which to diffuse and, therefore, a higher chance of encounter than entities far from any boundary. We can expect this effect to become important when the scale of the system becomes comparable to the typical distance travelled during a timestep. For biological systems on the μm scale, and chemical entities diffusing with diffusion constants ~μms-1, this sets an absolute upper limit on the timesteps at ~0.1 s. Properly taking into account these boundary effects is beyond the scope of the present work. However, boundary effects are not expected to play an important role as the timescale for these effects is ~0.1 s, which is a much longer timescale than the limit previously set by single particle interaction at ~0.002 s.