Fast simulation of Brownian dynamics in a crowded environment

Brownian dynamics simulations are an increasingly popular tool for understanding spatially-distributed biochemical reaction systems. Recent improvements in our understanding of the cellular environment show that volume exclusion effects are fundamental to reaction networks inside cells. These systems are frequently studied by incorporating inert hard spheres (crowders) into three-dimensional Brownian dynamics simulations, however these methods are extremely slow owing to the sheer number of possible collisions between particles. Here we propose a rigorous"crowder-free"method to dramatically increase simulation speed for crowded biochemical reaction systems by eliminating the need to explicitly simulate the crowders. We consider both the case where the reactive particles are point particles, and where they themselves occupy a volume. We use simulations of simple chemical reaction networks to confirm that our simplification is just as accurate as the original algorithm, and that it corresponds to a large speed increase.


I. INTRODUCTION
The fact that living cells constitute crowded cytoplasmic and nuclear environments has been appreciated for several decades. 1,2 However, the significance of excluded volume effects to specific biochemical processes has recently been highlighted by a multitude of experimental and theoretical observations. It is now established that crowding by large inert molecules can place limits on the total number of transcription factors in a cell, 3 can cause DNA to change its shape, 4 can encourage protein structure self-assembly, 5 and can both enhance and diminish transcription factor binding rates. 6 Correspondingly, several authors have recently proposed a variety of mathematical descriptions of crowding effects. Many of these are modifications of the compartment-based reaction-diffusion master equation, [7][8][9] which divides space into a lattice and models diffusion as particles hopping between neighbouring lattice sites. Lattice-based models have, however, been shown to underestimate the effects of crowding compared to more detailed descriptions. 10,11 Some authors have proposed introducing crowding effects directly into nonspatial descriptions such as the chemical master equation 12 or the deterministic reaction rate equations. 13,14 Alternatively, a popular lattice-free spatial technique involves Brownian dynamics (BD) simulations. [15][16][17] BD simulations explicitly track the positions of particles and model diffusion as a Brownian random walk in continuous space. Several popular modern BD simulators do not model crowding explicitly, since they assume particles to be point-particles with no physical volume. 18,19 Highly detailed molecular dynamics simulators are also popular, incorporating particle shapes, charge distributions, and hydrodynamic interactions, 17,[20][21][22] but their increased accuracy comes at the cost of considerably longer simulation times. However, designing algorithms to accurately study the behaviour of hard sphere colloids (uniform suspensions of insoluble particles) without hydrodynamic interactions were a popular problem in chemical physics long before the biochemical implications of volume exclusion were fully appreciated. [23][24][25] One such algorithm was proposed by Cichocki and Hinsen. 26 The idea behind the Cichocki-Hinsen algorithm is simple to state: only one particle is moved at a time, and if the attempted move results in a collision the particle is simply placed back in its previous position, thereby crudely modelling a steric repulsion. Despite its relative simplicity, the Cichocki-Hinsen algorithm has been proved to converge to the Smoluchowski equation in the limit of short simulation time steps. 26 Furthermore, it has been shown to agree perfectly with far more detailed algorithms which incorporate particle velocity and momentum. 27 It is therefore commonly used to simulate Brownian diffusion of hard spheres in the study of both physical chemistry [28][29][30] and cell biology. 15,31,32 However because of its fine-grained detail each simulation is computationally expensive, and many independent simulations are required to get good statistical samples.
In this article, we propose a modification to the Cichocki-Hinsen algorithm for reaction-diffusion systems. Our simplification arises from distinguishing between reactive particles (which may either be point particles or have a finite volume) and hard sphere crowders. Assuming that the crowder particles are uniformly distributed, we rigorously derive the probability that a reactive particle will collide with a crowder in a single time step, and use this to write a modified Cichocki-Hinsen algorithm which does not explicitly simulate crowders: we call this the crowder-free algorithm. We show that the crowder-free algorithm results in a dramatic speed increase over the original Cichocki-Hinsen algorithm of up to three orders of magnitude. Perhaps more surprisingly, the output data of the two algorithms are near-indistinguishable in terms of short-time diffusion coefficients, long-time diffusion coefficients, and reaction dynamics for each example that we test.
In Section II we propose the crowder-free algorithm for a system of reactive point particles in a sea of hard sphere crowders. We first outline the Cichocki-Hinsen algorithm for a point particle reaction-diffusion system. We then derive the probability that a small diffusive jump by a reactive point particle results in a collision with a crowder. Using this expression, we outline the crowder-free algorithm. We subsequently test our algorithm's speed and accuracy in modelling pure diffusion, zero, first, and second-order reactions, and the reactiondiffusion system A + B C in the presence of crowders. In Section III we analogously propose the crowder-free algorithm for a system of finite-size reactive particles in a sea of hard sphere crowders. We then derive the probability that a small diffusive jump by a finite-size reactive particle results in a collision with a crowder: this is shown to be very similar to the point particle expression. We again test our algorithm's speed and accuracy in modelling pure diffusion, zero, first, and second-order reactions, and the reaction-diffusion system ∅ − → X, X + X − → ∅ in the presence of crowders. We conclude with a discussion in Section IV.

II. POINT PARTICLES IN A CROWDED ENVIRONMENT
We first describe the Cichocki-Hinsen algorithm as applied to a system of reactive point particles in a sea of inert spherical crowders of radius R. The boundaries of the reaction volume can be of any type (periodic, reflective, etc.) as long as the number of crowder particles remains constant in time (i.e., no absorbing boundaries). Since the original Cichocki-Hinsen algorithm was written for purely diffusive systems, we have added some steps for reactive systems. The reactive method we use is the Doi model, 33,34 which assigns each reaction j a rate λ j .
If reaction j is a bimolecular reaction, it is also assigned a reaction distance r j . Bimolecular reaction j occurs with rate λ j when two reactive particles of the relevant type come within a distance r j of each other. Particles created by bimolecular reactions are typically placed midway between the two parent particles (this is the method we employ in our examples), though different placements may be appropriate for different examples. Unbinding reactions are assigned a rate λ j and an unbinding distance σ j . These reactions occur with rate λ j and normally the daughter particles are placed diametrically opposite each other on a sphere of diameter σ j centered around the parent particle, at a uniformly distributed angle (this again is the method we employ in our examples). Other standards exist for unbinding reactions (including those with more than two daughter particles) and the choice of which to implement is up to the user. Other monomolecular and zero-order reactions are simply assigned a rate λ j . Note that reaction distances and unbinding distances are not physical radii and do not exclude any volume.
Cichocki-Hinsen algorithm with reactive point particles 1. Uniformly distribute the reactive particles and the crowders in the volume, such that no crowders are intersecting each other and no reactive particles lie inside a crowder. Let N be the total number of particles (reactive and crowders), and randomly assign each particle a unique index 1, . . . , N. 2. For each i = 1, . . . , N, propose a new position for particle i at a random Normal(0, √ 2D i ∆t) displacement in each spatial dimension, where D i is the diffusion coefficient of particle i and ∆t is the simulation time step. If this new position causes an intersection between any particle (reactive and crowder), place particle i back in its original position. If not, place particle i in the new position. 3. For each reactive particle involved in a bimolecular reaction j, check if any reactive particle of the appropriate types lies inside a sphere of radius r j around the particle. For each appropriate reactive particle inside this sphere, propose a reaction with probability λ j ∆t. If successful, check if any daughter particle would intersect a crowder. If so, skip the reaction; if not, allow the reaction to proceed. 4. For each reactive particle of a type involved in a unimolecular reaction j, propose a reaction with probability λ j ∆t. If successful, check if any daughter particle would intersect a crowder. If so, skip the reaction; if not, allow the reaction to proceed. 5. For each zero-order reaction, propose a reaction with probability λ j ∆t. If successful, check if any of the new particles would intersect a crowder. If so, skip the reaction; if not, allow the reaction to proceed. 6. Advance time by ∆t. Let N be the new total number of particles and randomly reassign each particle a unique index 1, . . . , N. Return to (2) and repeat until a target time has elapsed.
The overwhelmingly time-consuming step of this algorithm is step (2), in which potential particle overlaps must be checked N times. The reaction steps (3)-(5) also involve potential overlaps, but as ∆t should typically be taken small enough that at most one reaction could plausibly happen per time step, these should not be particularly time-consuming. Our aim in Subsection II A is therefore to reduce the time taken by step (2). Note that step (1) can also be particularly timeconsuming: although our simplification does not particularly aim to fix that problem, it happens that by increasing the speed of step (2) we also dramatically shorten step (1).

A. Derivation
We first make two observations which form the basis of our method of reducing the time taken by the Cichocki-Hinsen algorithm. First, the crowders are inert and contribute little to the actual reactive behaviour of the system; their only function is to occasionally prevent a reactive particle from moving or the reaction from happening. Second, the crowders are uniformly distributed in space: this implies that each proposed reactive particle movement has roughly the same chance of being impeded by a crowder.
One common method of modelling diffusion in a crowded environment, based on the crowder uniformity assumption, is to simply replace the diffusion coefficient D with D(1 − φ), where φ is the proportion of the total volume occupied by crowders. 35 The idea is that if a particle attempts to move to a new location, there is a 1 − φ probability of that location not being occupied by a crowder. This is a valid assumption if the random particle displacement at a time step δx R, that is, if the particle moves by a distance much greater than the crowder size, such that its new location can be roughly considered a uniform random variable. However, it makes little sense to permit δx R, because that would allow particles to pass through crowders with a single jump.
On the other hand, permitting δx R makes physical sense, because the tiny perturbations which make up Brownian motion are much smaller than any particle radius. Furthermore, this is precisely the limit in which Cichocki and Hinsen proved their algorithm to be exact. 26 In that limit, however, we cannot use the 1− φ assumption. To understand why not, consider that the particle is already in a permitted location: this implies that with probability 1 there is a small sphere with radius > 0 around the particle which does not intersect any crowder. This local effect implies that the particle's new position cannot be treated as uniformly distributed: if δx is small enough (δx < ), the particle's new position is guaranteed to not intersect any crowder. In summary, if we require that δx R, then the probability that the particle's new position is illegal (intersects a crowder) is not given by 1 − φ but by some function of δx. We now attempt to derive that function.
Consider what happens when a point-particle proposes to move by a displacement δx. This is illustrated in Fig. 1. The particle's proposed new position will be on the surface of a sphere of radius δx around its current position. There will be no crowders with their centres in a sphere of radius R around the particle (otherwise the point particle could not be where it is currently), however there is a non-zero probability that there are crowders with their centres inside the spherical shell between the sphere of radius R + δx and the sphere of radius R (the grey region in Fig. 1). If there are crowders in this region, then there is some probability that the point particle's proposed FIG. 1. Diagram of a point particle attempting to move near a crowder of radius R. The particle attempts to displace itself a distance δx, such that its future position is on the surface of a sphere of radius δx around its current position. There may be crowders with their centres in the spherical shell of radius R + δx (grey region), which could prevent the particle displacement. The proposed position will be illegal if it is on the dotted segment of the sphere of radius δx. new position is illegal: this is precisely the probability that the proposed position intersects the crowder (the dotted segment in Fig. 1). Now, suppose that there are N C crowders of radius R inside a volume V. Assuming a uniform crowder distribution, the probability that a given crowder could collide with the point particle in a single time step is simply the ratio of the volume of the grey region to the total volume, (1) The probability of finding n crowders in the grey region is then given by the Binomial distribution, Of course, Eq.
(2) is only valid for small n, because there is a physical limit to how many crowders can fit in the relevant region. However, this is of little concern, since we are only concerned with the probabilities up to o δx R , which turns out to correspond only to n = 0 and n = 1, We now consider the probability that the proposed new point particle position intersects the crowder. This is given by the surface area of the spherical cap of the sphere of radius δx which lies inside the sphere of radius R around the crowder (the dotted segment in Fig. 1) divided by the total surface area of the sphere of radius δx. This is given by where d is the separation between the centres of the point particle and the crowder. 36 The expected value of d is simply R + δx 2 , so inserting this into Eq. (5) gives Combining Eq. (4) with Eq. (6) gives the probability that the proposed move is illegal, (7) Writing this in terms of the proportion of occupied volume, φ = 4 3 πN C R 3 V , leads to the simplified expression We can therefore write a much faster version of Cichocki-Hinsen algorithm which does not include any crowders. Only point particles need to be modelled explicitly in our algorithm, while the effect of crowders is incorporated by denying a point particle's proposed movement with probability P(illegal). For obvious reasons, we call this a crowder-free algorithm. This idea is shown in Fig. 2. The left panel shows the Cichocki- Hinsen algorithm with crowders (red) and point particles (green, purple). The points are not allowed to intersect the crowders, but the reaction radii are. The right panel shows the crowder-free algorithm, which looks identical to Cichocki-Hinsen without crowders. It is clear that the crowder-free algorithm will be easier to simulate.
Since none of the remaining particles in the crowder-free algorithm occupy any volume, we can move all particles simultaneously. The algorithm therefore essentially reduces to the classical Doi algorithm, with an extra clause for preventing particle movement. Some minor changes must also be made to the reaction parts of the algorithm (steps (3)-(5)), which originally prevented a reaction if a newly created particle would intersect a crowder. Since we no longer explicitly model crowders, we must modify this step. If the reaction is either bimolecular or unbinding, the new particle will be placed at a small displacement σ from a previous particle location (In the case of bimolecular reactions, σ should be the smaller of the two possible σ's). If σ R 1, then we can simply modify the diffusion formula to become P(illegal) = 3φσ 4R . Can we assume that σ R 1? In some cases, such as monomolecular conversion reaction of the type A → B, we will have σ = 0, and it would be absurd to prevent such reactions due to crowding. However, some reactions may have quite a large unbinding distance, and the diffusion formula may prove to be invalid. At each such reaction, we therefore check if σ R < 0.1. If this condition is true, we use the formula P(illegal) = 3φσ 4R , otherwise we use the formula , which is the probability that a uniformly distributed point particle would intersect a crowder. The choice of 0.1 is essentially arbitrary and can obviously be made smaller if required; we find that it gives good results, however. For zero-order reactions, we always use the formula , since particles created by these reactions have no parent particles.
Crowder-free algorithm with reactive point particles 1. Uniformly distribute the reactive particles in the volume. 2. Propose new positions for all particles at a random Normal(0, √ 2D i ∆t) displacement in each spatial dimension, where D i is the diffusion coefficient of particle i and ∆t is the simulation time step. Calculate δx, the length of the displacement, for each particle. With probability 3φδx 4R reject the proposed move, otherwise accept it. 3. For each particle of a type involved in a bimolecular reaction j, check if any particle of the appropriate types lies inside a sphere of radius r around the particle, where r is the reaction radius for the relevant reaction. For each appropriate particle inside this sphere, propose the reaction with probability λ j ∆t, where λ j is the corresponding reaction rate. For each daughter particle, calculate σ, the length of the displacement from the nearest parent particle. If σ R < 0.1, with probability 3φσ 4R reject the proposed reaction, otherwise accept it. Otherwise if σ R ≥ 0.1, with probability (where N C is the number of crowders) reject the proposed reaction, otherwise accept it. 4. For each reactive particle of a type involved in a unimolecular reaction, propose a reaction with probability λ j ∆t, where λ j is the reaction rate. For each daughter particle, calculate σ, the length of the displacement from the parent particle. If σ reject the proposed reaction, otherwise accept it. 5. For each zero-order reaction, propose a reaction with probability λ j ∆t, where λ j is the reaction rate. With prob- reject the proposed reaction, otherwise accept it. 6. Advance time by ∆t. Return to (2) and repeat until a target time has elapsed.
In Sec. II B, we confirm that the crowder-free algorithm is orders of magnitude faster than Cichocki-Hinsen, while retaining its accuracy.

B. Comparative tests
In our first test of the crowder-free algorithm, we consider a single point particle diffusing in space, surrounded by a uniform distribution of crowders. In this scenario, the crowderfree algorithm should show a dramatic improvement over the original Cichocki-Hinsen algorithm in terms of computation time.
Indeed, as shown in Fig. 3, we find that the crowder-free algorithm is at least an order of magnitude faster than the standard algorithm when there are only 10 crowders, this increases to three orders of magnitude when there are 500 crowders. A significant advantage is that the crowder-free algorithm does not scale with the number of crowders, making it particularly useful for studying high levels of crowding.
Note that the computation speed for the crowder-free algorithm is slightly faster for higher crowding levels. The reason for this is that as the number of crowders increases, the probability that the diffusing particle does not move on a given time step increases. It takes marginally less computational time to not move a particle than it does to move one (since we do not need to update the particle position).
Of course, fast simulation is of little use if the results of the algorithm are inaccurate. In our second test, we therefore use sample paths from both algorithms to compute the effective short-time diffusion coefficient D * of a single point particle in crowded space. 38 This is done by performing a simulation with input diffusion coefficient D, computing the squared displacement of the particle at each time step, and taking the mean of that value over the entire simulation. This value is equated to 6D * ∆t to find an estimate for the effective short-time diffusion coefficient D * .
The non-dimensional parameter D * D is the effective reduction in short-time diffusion coefficient due to crowding. For no crowding, we expect D * D = 1, and the value should decrease as crowding increases. This is because large jumps are more likely to result in a collision with a crowder than small jumps, so the effective diffusion coefficient appears to be reduced. In Fig. 4 we plot D * D as a function of the proportion of occupied volume φ. As expected, both algorithms show a reduction in the effective short-time diffusion coefficient as crowding increases, and both algorithms give very similar results, with their error bars always intersecting. Each data point is an average of 10 simulations, each simulation ran until the point FIG. 4. Relative reduction in the short-time diffusion coefficient for both the Cichocki-Hinsen algorithm (blue) and the crowder-free algorithm (red), for a single point particle diffusing in space, as a function of the proportion of occupied volume φ. All data points are an average of 10 simulations, error bars are 1 standard deviation. Parameter values are V = 1, R = 0.05, ∆t = 10 −5 , D = 0.1 for the point particle, D = 0.01 for the crowders. particle, initially located at ( 1 2 , 1 2 , 1 2 ), left the unit cube with corners at (0,0,0) and (1,1,1).
It has been shown, however, that short-time diffusion coefficients can differ strongly from long-time diffusion coefficients measured over a whole trajectory. 38 In Fig. 5, therefore, we also plot the long-time diffusion coefficients estimated using the unbiased diffusion estimator developed in Ref. 37. Each point in Fig. 5 is the mean of 3 independent simulations of length 10 3 time steps, and each error bar corresponds to the standard deviation.
We have confirmed that the crowder-free algorithm simulates diffusion as accurately as the original Cichocki-Hinsen algorithm, but we have not tested whether it accurately simulates reactions. In our next test, we check that some basic reactions happen with the same frequency for both algorithms. In Fig. 6 we show the results of these tests. In Fig. 6(a) we compare the Cichocki-Hinsen algorithm with the crowder-free algorithm for a zero-order reaction under a variety of crowding levels. The quantity Λ on the y-axis is the ratio between the actual frequency of the reaction and the rate specified in the algorithm, Λ = 1 therefore corresponds to no effect of crowding. We observe that both algorithms show the same linear reduction in the effective rate as crowding increases. In Fig. 6(b) we perform the same comparison for an unbinding reaction of the form X → Y + Z, for two different unbinding distances σ = 0.001 and σ = 0.05. The σ = 0.001 case corresponds to a short distance (as defined in the algorithm) and so we would use Eq. (8) in the crowder-free algorithm, whereas σ = 0.05 corresponds to a large distance, and so we would use φ, the volume occupied, in the crowder-free algorithm. For both parameter sets, the crowder-free algorithm agrees well with Cichocki-Hinsen. The unbinding reaction with a larger unbinding distance is naturally more affected by crowding than the reaction with a smaller unbinding distance, since large jumps are more likely to be impeded by a crowder. In Fig. 6(c) we perform the same test for a binding reaction of the form X + Y → Z, again both algorithms   6. (a) Λ for a zero-order reaction ∅ → X. (b) Λ for an unbinding (first order) reaction X → Y + Z, for different unbinding distances σ = 0.001 and σ = 0.05. (c) Λ for a binding (second order) reaction X + Y → Z. For all plots, R = 0.05, D = 1, ∆t = 10 −4 , V = 1. All data points are an average of 100 independent simulations, error bars are 1 standard deviation. The quantity Λ is the ratio between the actual frequency of the reaction and the rate specified in the algorithm; thus Λ = 1 corresponds to crowding having no effect on the reaction frequency. agree well, though since these are point-particles there does not appear to be a significant effect of crowding on the binding rate. This differs from the finite-size particle case shown in Fig. 11.
To confirm the above test, we finally use both algorithms to compute the equilibrium distribution of the reaction A + B C in the presence of low and high levels of crowding. We expect the typical number of C molecules to be higher for high crowding, because the unbinding reaction will occur less frequently.
For each algorithm, we simulated two long trajectories of a system initially consisting of 30 uniformly distributed A molecules and 30 uniformly distributed B molecules, in a sea of 10 (low crowding) and 700 (high crowding) crowders. The simulation time was much longer than the time for the system to reach equilibrium. In Fig. 7 we show the equilibrium distribution for the number of C molecules. The mean number of C molecules shifts from around 6 with low crowding to around 11 with high crowding. The crowder-free algorithm agrees almost perfectly with the Cichocki-Hinsen algorithm for both examples, thus confirming that the crowder-free algorithm accurately imitates the Cichocki-Hinsen algorithm, but with a dramatic reduction in computation time.

C. A note on more complex systems
The crowder-free algorithm proposed above specifically concerns a uniform distribution of crowders with the same radius; however the results can equally be applied to more complex systems.
For sets of crowders with different radii, say N (i) C crowders of radius R i for i = 1, . . . , k, we can simply use the formula which will give the probability of a move δx resulting in a collision. Of course, this formula relies on the assumption that δx R i for all i = 1, . . . , k. For systems with a non-uniform distribution of crowders of radius R, the algorithm can still be used if the crowder distribution is locally uniform. In that case, we can divide the volume up into k subvolumes V i with N (i) C crowders for i = 1, . . . , k, where V 1 + · · · + V k = V and N (1) C + · · · + N (k) C = N C . Then we can apply the formula for a point particle in the ith subvolume. However, this method will only work if the crowder distribution remains roughly constant in time. If the crowders are diffusing fast enough that the overall distribution flattens on the time scale of the simulation, then subvolume i will not always contain N (i) C crowders. Since we do not know how N (i) C will change a priori, we cannot use the crowder-free algorithm for such examples.

III. FINITE-SIZE PARTICLES IN A CROWDED ENVIRONMENT
Studying the behaviour of reactive point particles in the presence of crowders provides useful information about real biochemical systems in which the reactive particles are much smaller than the crowders they encounter. This is an accurate description of, for example, small proteins or amino acids diffusing in the vicinity of ribosomes or large enzymes. However, biochemical particles also encounter crowders with a similar size to themselves. In order to study these examples effectively, we must also be able to simulate reactive particles which occupy a non-zero volume. A version of the Cichocki-Hinsen algorithm for which the reactive particles occupy a non-zero volume is given below. Since reactive particles now have a physical radius, we no longer need to define a reaction distance for bimolecular reactions: particles react with a rate λ j if they physically intersect. This is known as partial-absorption Smoluchowski binding. 39 Cichocki-Hinsen algorithm with finite-size reactive particles 1. Uniformly distribute the reactive particles and the crowders in the volume, such that no particles (reactive or crowder) are intersecting each other. Let N be the total number of particles, and randomly assign each particle a unique index 1, . . . , N. 2. Uniformly sample an integer i from 1, . . . , N. Propose a new position for particle i at a random Normal(0, √ 2D i ∆t) displacement in each spatial dimension, where D i is the diffusion coefficient of particle i and ∆t is the simulation time step. If particle i is a crowder, check if this new position causes an intersection between any particle. If so, place particle i back in its original position, if not, place particle i in the new position. Otherwise if particle i is a reactive particle, check if this new position causes an intersection between i and exactly one other reactive particle and no crowders. If so, and if that particle can react with i, proceed to (3). Otherwise, if the new position causes any other type of intersection, place the particle back in its original position; if not, place the particle in its new position. Proceed to (4).

Propose a bimolecular reaction j with probability λ j ∆t,
where λ j is the corresponding reaction rate. If unsuccessful, place particle i back in its previous position and proceed to (4). Otherwise if successful, check if any daughter particle would intersect another particle. If so, skip the reaction, place particle i back in its original position; if not, allow the reaction to proceed. 4. For each reactive particle of a type involved in a unimolecular reaction j, propose a reaction with probability λ j ∆t/N, where λ j is the reaction rate. If successful, check if any daughter particle would intersect any other particle. If so, skip the reaction; if not, allow the reaction to proceed. 5. For each zero-order reaction j, propose a reaction with probability λ j ∆t/N, where λ j is the reaction rate. If successful, check if any of the new particles would intersect another particle. If so, skip the reaction; if not, allow the reaction to proceed. 6. Advance time by ∆t/N. Let N be the new total number of particles and randomly reassign each particle a unique index 1, . . . , N. Return to (2) and repeat until a target time has elapsed.
Note that this algorithm is distinct from the Cichocki-Hinsen algorithm in Section II in several ways, mainly because in this algorithm time is advanced by ∆t N at each time step. This is because here step (3) is nested inside step (2). The reason for this is that bimolecular reactions occur in this algorithm when two reactive particles physically intersect. This is an illegal move, and if the particles do not react then they must not be allowed to remain in that position, but rather revert to the previous position, hence bimolecular reactions and diffusion are closely coupled in this algorithm. It follows that N can change during steps (2)-(3), and so it does not make sense to place step (2) inside a for-loop over i = 1, . . . , N.
Again, step (2) is the overwhelmingly time consuming step for this algorithm, so as before we will attempt to find an expression giving the probability that a given jump causes an intersection with a crowder. However, we will not be able to get substantial speed gains on the same scale that we obtained with point-particles because now even a crowder-free algorithm will contain finite-size reactive particles. Our speed increase will arise from removing a subset of the volume-occupying particles (the crowders) rather than all of them, as before. Obviously, our method will work best if there are many more crowders than reactive particles, though it will always be faster than the standard algorithm.

A. Derivation
To derive an analogous formula to Eq. (8) for the finitevolume case, consider a reactive particle with radius r > 0 attempting to move a distance δx in a sea of N C uniformly distributed crowders of radius R. In Section II, we observed that, to first order in δx R , the probability of a reactive particle performing an illegal move depends only on its behaviour in the vicinity of a single crowder. However, a particle of radius r moving near a single crowder of radius R is identical to a point-particle moving near a crowder of radius R + r: in both cases, the two particle centres are forbidden from being nearer than R + r from each other. It follows that Eq. (7) can be easily adapted for use in this section, but with R replaced by R + r. In other words, we can simply write Observe that we do not need to consider the probability of intersecting reactive particles here. This is because the reactive particles will all be simulated explicitly, so a collision between reactive particles in the crowder-free algorithm will be simulated identically to the original algorithm. As before, we will also need to moderately adapt the reaction part of our algorithm. Again, if a daughter particle is created at a small distance σ from a parent particle, and σ R + r, then we can use the formula P(illegal) = πN C (R + r) 2 σ V . Note, however, that this is much less likely to occur with finitesize particles, since σ will typically be a similar order of magnitude to r, which is in turn typically a similar order of magnitude to R. We could alternatively use the probability that a uniformly distributed point in space can accommodate a particle of radius r. Users may wish to simulate reactions where it makes more sense to use P(illegal) = πN C (R + r) 2 σ V as the probability that a new particle is obstructed by a crowder (for instance, an enzyme releasing a much smaller substrate). Since these matters are system-specific, we leave it up to the user to decide which formula is more appropriate. However, for almost all reactions we will use the probability that a uniformly distributed point in space can accommodate a particle of radius r.
This probability is not the simple expression used in Section II, rather it derives from the scaled particle theory (SPT). The reason for this is that there are unoccupied points in space which are inaccessible to the particle of radius r. These are the points which do not lie inside a crowder but do lie within a distance R + r from a crowder's centre. SPT has been used to obtain analytical expressions for the effect of crowding on intrinsic noise in two-dimensional systems and was observed to give very accurate results. 12 In three dimensions, it offers an expression for the probability that a uniformly distributed point in space of volume V can accommodate a particle of radius r, given that the space contains N C crowders of radius R, 40 log The crowder-free algorithm for finite-size reactive particles is then as follows: Crowder-free algorithm with finite-size reactive particles 1. Uniformly distribute the reactive particles in the volume, such that no particles are intersecting each other. Let N be the total number of particles, and randomly assign each particle a unique index 1, . . . , N. 2. Uniformly sample an integer i from 1, . . . , N. Propose a new position for particle i at a random Normal(0, √ 2D i ∆t) displacement in each spatial dimension, where D i is the diffusion coefficient of particle i and ∆t is the simulation time step. With probability , where r is the radius of particle i, put the particle back in its original position. Otherwise, check if this new position causes an intersection between i and exactly one other particle. If so, and if that particle can react with i, proceed to (3). Otherwise, if the new position causes any other type of intersection, place the particle back in its original position; if not, place the particle in its new position. Proceed to (4). 3. Propose a bimolecular reaction j with probability λ j ∆t, where λ j is the corresponding reaction rate. If unsuccessful, place particle i back in its previous position and proceed to (4). Otherwise if successful, evaluate P(legal) according to Eq. (12) for each daughter particle. Let p be the product of each P(legal). With probability 1 p, skip the reaction and place particle i back in its original position. Otherwise check if any daughter particle would intersect another particle. If so, skip the reaction, place particle i back in its original position; if not, allow the reaction to proceed. 4. For each reactive particle of a type involved in a unimolecular reaction j, propose a reaction with probability λ j ∆t/N, where λ j is the reaction rate. If the reaction is of the type A − → B and the radius of B is less than or equal to that of A, allow the reaction to proceed. Otherwise, evaluate P(legal) according to Eq. (12) for each daughter particle. Let p be the product of each P(legal). With probability 1 p, skip the reaction. Otherwise, check if any daughter particle would intersect any other particle. If so, skip the reaction; if not, allow the reaction to proceed. 5. For each zero-order reaction j, propose a reaction with probability λ j ∆t/N, where λ j is the reaction rate. If successful, evaluate P(legal) according to Eq. (12). With probability 1 − P(legal), skip the reaction. Otherwise check if any of the new particles would intersect another particle. If so, skip the reaction; if not, allow the reaction to proceed. 6. Advance time by ∆t/N. Let N be the new total number of particles and randomly reassign each particle a unique index 1, . . . , N. Return to (2) and repeat until a target time has elapsed.
There is one significant case for which our crowder-free algorithm will not give accurate results, namely, if the crowders are stationary and the level of crowding is high. Simulating such systems with Cichocki-Hinsen reveals that reactive particles can get trapped in regions surrounded by stationary crowders and simply stay there for the entirety of the simulation without reacting or moving significantly. Obviously, these cases cannot be covered by the crowder-free algorithm because all reactive particles (of the same radius) have the same probability of diffusing at any time. We therefore recommend not using the crowder-free algorithm for systems with stationary crowders unless the level of crowding is sufficiently low that no trapping regions could exist. Note that this is not a problem if the reactive particles are point-particles, because they occupy no volume and will always be able to escape from a trapping region eventually.

B. Comparative tests
In this section, we perform similar tests on the crowderfree algorithm for finite-size particles to those we performed in Section II B. We initially test the time taken for both methods to simulate pure diffusion in the presence of an increasing number of crowders. To ensure that the results are different from those in Section II B, we now simulate 50 diffusing "reactive" particles (so-called even though they do not react in this example) in a sea of crowders. Of course, we do not expect to get anywhere near the 1000-fold speed increase that we achieved for the point-particle case: even with no crowders, we have to simulate 50 volume-occupying molecules, constantly ensuring that they do not intersect.
The results of this test are plotted in Fig. 8. With 10 crowders, the crowder-free algorithm takes half the time of the Cichocki-Hinsen algorithm, while with 400 crowders, the crowder-free algorithm has a speed increase of over 20 times. Even for finite-size particles, therefore, the crowder-free algorithm offers a considerable speed increase, and its lack of dependence on crowder number makes it especially useful for studying high levels of crowding. The next test we perform compares estimates of short-time diffusion coefficients from the two algorithms. In both cases, we simulate 20 finite-size particles diffusing in a sea of crowders. Because of this, a single simulation gives 20 different estimates of the diffusion coefficient. In Fig. 9 we plot the mean (points) and standard deviation (error bars) of this sample of 20, for a variety of levels of crowding. Since the "reactive" particles themselves occupy a volume, we incorporate this into our calculation of the proportion of occupied volume φ. As in Fig. 4, the two algorithms agree, with error bars intersecting for each data point. Note that, compared to Fig. 4, the diffusion coefficient is reduced more for the same level of crowding. This confirms the intuitive hypothesis that finite-size particles are more influenced by crowding than point particles.
As in the point-particle case, we also study the long-time diffusion coefficients of finite-sized particles. In Fig. 10 we plot the long-time diffusion coefficients estimated using the unbiased diffusion estimator developed in Ref. 37. We calculated the diffusion coefficient for three different particle sizes, ranging from the same size as the crowders to several orders of magnitude smaller than the crowders. Both algorithms agree well, and we observed that, as we might expect, the diffusion  of smaller particles is less affected by crowding than larger particles. Each point in Fig. 10 is the mean of 10 independent simulations of length 10 3 time steps, and each error bar corresponds to the standard deviation.
In our next test, we check that zero, first, and secondorder reactions happen with the same frequency for both the Cichocki-Hinsen algorithm and our crowder-free algorithm for a variety of particle sizes. In Fig. 11 we show the results of these tests. In Fig. 11(a) we compare the Cichocki-Hinsen algorithm with the crowder-free algorithm for a zero-order reaction under a variety of crowding levels and for two particle sizes. As before, the quantity Λ on the y-axis is the ratio between the actual frequency of the reaction and the rate specified in the algorithm, Λ = 1 therefore corresponds to no effect of crowding. We observe that both algorithms show the same linear reduction in the effective rate as crowding increases, and the rate is reduced more for large particles than for small ones. In Fig. 11(b) we perform the same comparison for an unbinding reaction of the form X → Y + Z, for two different parameter sets. First for X unbinding into equal sized Y and Z with radius 0.05 and unbinding distance 0.1, and second for X unbinding into different sized Y and Z, with radii 0.01 and 0.005, respectively, and unbinding distance 0.015. For both parameter sets, the crowder-free algorithm agrees well with Cichocki-Hinsen. In Fig. 11(c) we perform the same test for a binding reaction of the form X + Y → Z, for two different parameter sets. First for X and Y with radius 0.05, and second for X and Y with radius 0.02. Again both algorithms agree well for both parameter sets, and we note that higher crowding reduces the binding rate, because it takes longer for particles to find each other.
Finally, we compare the algorithms' performance at estimating an equilibrium distribution of a chemical reaction. This time we simulate the reaction ∅ − → X, X + X − → ∅, in which particles are created at uniformly distributed points in space and react with a fixed rate when they collide. This system has previously been studied spatially as an example of protein synthesis and degradation. 41 We expect that, contrary to the example in Fig. 7, crowding will reduce the mean number  11. (a) Λ for a zero-order reaction ∅ → X, for different particle radii r = 0.05 and r = 0.005. (b) Λ for an unbinding (first-order) reaction X → Y + Z, for different Y radii r = 0.05 and r = 0.005. (c) Λ for a binding (second-order) reaction X + Y → Z, for different X and Y radii, r = 0.05 and r = 0.02. For all plots, R = 0.05, D = 1, ∆t = 10 −4 , V = 1. All data points are an average of 100 independent simulations, error bars are 1 standard deviation. The quantity Λ is the ratio between the actual frequency of the reaction and the rate specified in the algorithm; thus Λ = 1 corresponds to crowding having no effect on the reaction frequency.
of X, since the creation of X will be less likely in crowded conditions.
In Fig. 12 we plot the equilibrium distribution of the number of X molecules for both algorithms in both low and high crowding conditions. Each distribution is calculated as a time average over a single long trajectory between 10 5 and 10 7 iterations. The crowder-free algorithm clearly requires fewer iterations than Cichocki-Hinsen because each iteration of both algorithms advances time by ∆t N , where N is the total number of particles, and Cichocki-Hinsen generally has many more particles to simulate. As predicted, the mean of the distribution is much lower in the high crowding example than the low crowding example. As with all previous tests, the crowderfree algorithm agrees almost perfectly with the Cichocki-Hinsen algorithm, confirming that our algorithm suffers little apparent loss of accuracy compared to the Cichocki-Hinsen algorithm, despite its considerable speed increases. Note that the value of φ is not given for these examples because the number of reactive particles fluctuates over time, and therefore so does φ.

IV. DISCUSSION
In this paper, we have proposed a modification to the commonly used Cichocki-Hinsen Brownian dynamics algorithm for simulating reaction-diffusion systems in a crowded environment. We call our modified algorithm a crowder-free algorithm because we do not simulate crowders explicitly. Instead, we rigorously derive the probability that a small displacement of size δx would result in a collision with a crowder. This implies that, instead of simulating crowders, we can simply reject each attempted particle displacement with precisely that probability.
We tested our algorithm in terms of both speed and accuracy, both for cases with reactive point particles and with finitesize reactive particles. The crowder-free algorithm always provides a speed increase over the underlying Cichocki-Hinsen algorithm: this speed increase varied from 2 to over 1000 for the set of examples studied in this paper. Furthermore, the crowder-free algorithm provides data which are nearindistinguishable from the data extracted from the corresponding Cichocki-Hinsen algorithm: this was shown to be true for both diffusive and reactive information. The crowder-free algorithm therefore shows no apparent loss of accuracy compared to the Cichocki-Hinsen algorithm, which, coupled with the clear speed increases, makes it a very attractive algorithm for simulating chemical reactions in a crowded environment.
There are two main cases where the crowder-free algorithm is not more effective than the Cichocki-Hinsen algorithm. The first case is if the initial crowder distribution is not uniform and spreads out over time. In that case, our algorithm is inadequate because we do not know a priori how fast the crowders will diffuse. Note, however, that non-uniform crowder distributions are not a problem in themselves: we can simply subdivide the volume into regions where the distribution is locally uniform, and derive separate values of P(illegal) in each region. The second case involves stationary crowders. If the level of crowding is high and the crowders do not diffuse, then some regions of space may be entirely segregated from others. Since the crowder-free algorithm allows all reactive particles to diffuse anywhere in space, it cannot accurately imitate the Cichocki-Hinsen algorithm in this case.
Finally, we note that further speed increases in both the crowder-free algorithm and the Cichocki-Hinsen algorithm may be possible by more efficient methods of measuring the