{'fraction_non_alphanumeric': 0.02593039564238886, 'fraction_numerical': 0.00929580221142242, 'mean_word_length': 4.75234521575985, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Excluded volume effects in on- and off-lattice reaction-diffusion models

Mathematical models are important tools to study the excluded volume effects on reaction-diffusion systems, which are known to play an important role inside living cells. Detailed microscopic simulations with off-lattice Brownian dynamics become computationally expensive in crowded environments. In this paper we therefore investigate to which extent on-lattice approximations, so called Cellular Automata models, can be used to simulate reactions and diffusion in the presence of crowding molecules. We show that the diffusion is most severely slowed down in the off-lattice model, since randomly distributed obstacles effectively exclude more volume than those ordered on an artificial grid. Crowded reaction rates can be both increased and decreased by the grid structure and it proves important to model the molecules with realistic sizes when excluded volume is taken into account. The grid artifacts increase with increasing crowder density and we conclude that the computationally more efficient on-lattice simulations are accurate approximations only for low crowder densities.


Introduction
Living cells are regulated by complicated signaling pathways that control which genes are expressed and how the cells behave. These reaction networks are spatially organized with important reaction complexes often bound to the cell membrane and with the DNA being confined inside the nucleus in eukaryotes. Moreover, the cytoplasm of cells is highly crowded, meaning that up to 40% of the volume is occupied by macromolecules [18,28], while individual species are present only at very low concentrations. On membranes crowding is even more severe [16], due to attaching actin filaments creating static barriers [15,20]. This macromolecular crowding leads to three excluded volume effects: (i) it forces moving molecules to diffuse around obstacles, or "crowders", slowing down the diffusion; (ii) it either decreases (diffusion limited case) or increases (reaction rate limited case or dimerizations) the reaction rates; and (iii) due to the impeded diffusion, it increases the spatial heterogeneity inside the cells, leading to spatial self organization [5,14]. However, the effect of macromolecular crowding has not yet been fully understood, since single molecule tracking inside cells is still difficult. Consequently, mathematical models are a crucial tool to understand reaction-diffusion processes in the intracellular environment.
One class of mathematical models are so called particle based reactiondiffusion models (PBRD), where we resolve individual molecules by following their diffusive paths and by modeling reactions as random events. These are considered as the microscopic modeling level in reaction-diffusion simulations. One group of PBRD are Brownian dynamics (BD) simulations, where the diffusive motion is modeled as a continuous time continuous space stochastic process. When the particles are modeled as hard spheres these simulations inherently account for the excluded volume effects. But the catch with these methods is, that they become computationally very expensive due to the high number of collision events when applied to dense environments. An approximation to BD simulations are lattice based approaches, such as cellular automata (CA), where we still follow individual trajectories, but the possible particle positions are restricted to an artificial grid. These methods are cheaper to simulate and in this paper we will investigate what influence the artificial lattice has on the excluded volume effects for the reaction and diffusion rates and the stochastic noise. We therefore simulate both models in a two-dimensional plane representing cellular membranes.
This microscopic modeling framework captures the internal end external noise, that often play an important role in computational systems biology. Reaction-rate equations (RREs) are a macroscopic approximation to this microscopic description. These are deterministic ordinary differential equations describing the concentrations of the reacting molecules, and they are applicable when the law of large numbers holds, meaning that the molecules are present at large copy numbers, and when they are well-mixed in a dilute system. They consequently do not capture the stochastic or space dependent effects inside cells, but provide suitable reference solutions in dilute media to compare the effect of macromolecular crowding to.
In the next section we will present the on-and off-lattice microscopic tools in more detail. We will then perform simulations in environments with various crowder densities for static crowders and show that the diffusion is more severely slowed down when simulated in continuous than in discrete space. In Section 4 we will extend the simulations to reactions in a crowded environment, and find that they are non-linearly affected by the artificial lattice. We finally draw conclusions on the agreement between the BD and CA models when excluded volume effects are important in Section 5.

Particle based reaction-diffusion models
Particles in solution move due to their thermal energy, and the collisions with the smaller solvent molecules result in a random walk often modeled by Brownian motion. This causes them to diffuse from high to low concentrations. We will now present how to create sample trajectories of a reaction-diffusion system with on-and off-lattice simulations, where the particles diffuse with diffusion coefficient γ 0 . We demonstrate how to treat reaction events using the example of the association event where k A is the intrinsic reaction rate.

Continuous space
Brownian dynamics (BD) simulations sample trajectories of Brownian motion in continuous space, or off-lattice. To simulate a trajectory one can discretize time with a small enough fixed time step ∆t. We then draw a normally distributed random number ξ ∼ N (0, 1) to compute the new particle position as and equivalently for the other space coordinates, since Brownian motion is independent along the Cartesian axes. This algorithm is implemented among others in the freely available software package Smoldyn [3,2]. Here, a time step dependent binding radius is chosen such that the simulated reaction rate equals k A . Within this binding radius two reaction partners react with probability one and the method has been used for simulations in a crowded environment [9,Ch. 4].
Similarly, in MCell [30] the position of a diffusing molecule after the discrete time step ∆t is sampled from the probability distribution of its position. Another algorithm for simulations in dilute media is Green's functions reaction dynamics (GFRD) [34], where single particles or particle pairs are surrounded by protective domains. An asynchronous time step is chosen in an event driven algorithm to be the time when the first particle leaves its protective domain and in this way one avoids simulating all the relatively uninteresting jumps between collision events. An exact implementation is the so called first-passage kinetic Monte Carlo (FPKMC) algorithm [7,24,31], used in the software ECell [33]. The GFRD algorithms, however, become computationally very expensive [17], when applied to non-dilute systems and we will perform BD simulations with the software Smoldyn to investigate the excluded volume effects in continuous space.
Volume exclusion effects are inherently simulated with these algorithms [32], when the particles are represented as hard spheres rather than point particles. The review article [29], summarizes which of the available software packages allow for this feature.

Discrete space
To gain further efficiency when simulating PBRD we will now restrict the particles' trajectories to an artificial lattice. These so called cellular automata (CA) or lattice gas automata models [6] are widely used to investigate excluded volume effects on both the diffusion rates [8] and the reaction rates [5,12,28]. First, the domain is discretized into a grid or lattice. Each lattice site can then contain at most one molecule and at each time step the molecules are picked in random order to jump to a neighboring site. If the sampled target site is occupied, either a reaction happens with probability p A or the move is rejected, which models the excluded volume effects. Since the mean square displacement (MSD) of a diffusing molecule in d dimension is we choose the time step for a lattice spacing h. The algorithm for the simple reaction (2.1) until final time T then reads as follows.
Algorithm 1 Cellular Automata 1: Initialize system by placing initial numbers of A, B and C molecules randomly on the grid. 2: while t < T do 3: Choose molecules in random order.

4:
for each molecule do

5:
Randomly choose a nearest neighbor site as target.

6:
if target site is empty then 7: Move molecule. if molecule is A(B) and target is occupied by B(A) then 10: Generate a random number ξ.

11:
if ξ < p A then 12: Replace A and B with a C molecule at target site.     Update t := t + ∆t.

22: end while
This algorithm has been used for Cartesian and hexagonal lattices and in [12] it is shown that the algorithm is sensitive to the choice of grid in a crowded environment and in general overestimates the excluded volume effects on the reaction rates as compared to BD simulations. In [8] the algorithm has been extended to model different molecule shapes, which consist of combinations of Cartesian lattices, and the diffusivity of moving molecules among static crowders is computed for a mix of shapes.
In the next sections we will perform both, BD and CA simulations in crowded environments to examine if CA are a good approximation of BD in this setting. The crowders are here represented as static and inert particles, that do not actively participate in the reaction processes.

Diffusion simulations
In all the experiments we choose γ 0 = 1 and perform the on-lattice simulations according to Algorithm 1, referred to as CA, and use the open source software Smoldyn for the off-lattice simulations, referred to as BD. We investigate three different lattices (Cartesian, Cartesian with diagonal jumps and hexagonal) in the CA simulations, see Fig. 3.1, where the coefficients for the diagonal jumps are chosen as in [22]. We discretize the two-dimensional domain 50 × 50 with h = 1, so that all molecules (crowders and diffusing) have size 1 in the CA simulations and are spheres with r = 0.5 in the BD simulations. The resolution of the on-lattice simulations is defined by h and according to (2.4) this leads to the CA time step ∆t CA = 0.25. For the BD simulations we choose a finer resolution with ∆t BD = 0.001, to guarantee that the off-lattice simulations are time step independent. In order to investigate the diffusive motion over long times without encountering boundary effects we implement periodic boundary conditions on all boundaries.  10 0 10 0 more accurate BD simulations with the small time step ∆t BD run ca. 58 times slower than the CA simulations for φ = 0.4. If no crowders are placed in the system (φ = 0) all models agree with the theoretical MSD in dilute medium (2.3) (plot not shown). But, for an increasing φ the on-and off-lattice models increasingly differ and the speed of diffusion decreases, which is expected. The unexpected finding is that the diffusion is more heavily slowed down with BD than with CA simulations. This appears counterintuitive at first, since the artificial grid decreases the available directions of movement and the off-lattice model with infinitely many degrees of freedom (dof) is expected to simulate more mobile particles. But, the lattice also orders the particles, so that they effectively exclude less space, see . This principle can be understood intuitively when considering a parking lot with predefined parking spots, imagine finding a spot or leaving the lot if all cars were parked randomly. The restriction of the degrees of freedom also makes it more probable to choose the possible direction out of a finite number of lattice directions as compared to the probability to sample a jump in the small angular direction ϕ, see Fig. 3.3(c). The increased flexibility on the hexagonal lattice leads to even faster diffusion than on a Cartesian, but allowing for diagonal jumps does not increase diffusion any further, since most of the jumps (82%) are sampled along the Cartesian axes. Hence, the Cartesian grid with only 4 dof has the closest agreement with the off-lattice simulations with infinitely many dof, contrary to the findings in [12] for fractal-like reaction kinetics. For very high and unbiological crowder densities (φ = 0.6), particles in BD and CA simulations on a Cartesian mesh effectively no longer move and only the on-lattice models with more dof allow the particles to find a free passageway. For normal diffusion in a dilute medium in integer dimension d the MSD grows linearly with time (2.3). In the second row in Fig. 3.2 we further observe that the MSD is no longer linear, but that diffusion becomes anomalous with a time varying diffusion constant. This behavior is explained in [16] and illustrated in Fig. 3.4, where the change in the diffusion constant is due to encountering the surrounding obstacles, before converging to an average long time behavior. Note, that we do not depict the small solvent molecules causing the random walk. However, when the crowder density is higher than the percolation threshold, the obstacles form a cluster traversing the whole domain, so that disconnected subdomains form and space has a non-integer dimension smaller than two, consequently the MSD is non-linear for all times [4,35]. The percolation thresholds for lattices are φ = 40.73% (Cartesian) and φ = 50.30% (Hexagonal) [23]. For the off-lattice case with partially overlapping disks for the excluded volume as illustrated in Fig. 3.3(b) it is more difficult to find exact values, but if the free space would consist of fully overlapping spheres 67.63% [25] would have to be freely available which would leave us with an occupancy fraction of 32.37% for the effective excluded volume. In the second row in Fig. 3.2 we see that the simulations reveal sustained anomalous diffusion for the cases when the percolation threshold has been exceeded.
To quantify the excluded volume effects on the effective diffusivity, we evaluate the diffusion constant γ(φ) after the transient phase by taking the last value of the plot x 2 (t) /(4t), for the cases converging to normal slower diffusion, see Fig. 3.5. As expected, the effective diffusion constant γ(φ) decreases with increasing crowder density and we can clearly observe the lattice artifact of underestimating the excluded volume effect on the diffusivity. It also appears that γ(φ) depends linearly on φ for all models, and the decrease in diffusivity for BD agrees with the findings in [21] for equally sized spheres. We now simulate 100 trajectories for one specific crowder distribution and plot the mean together with the 95% confidence interval in Fig. 3.6, in order to investigate the excluded volume effects on the variance of the diffusive motion. As φ increases, less and less space becomes available for the particles' diffusion and hence the variance and diffusion speed decrease. Another interpretation is that the higher number of particles in the system leads to a more deterministic behavior. Similarly, the BD simulations resulting in slower diffusion have a smaller trajectory to trajectory variance, but we do not observe a grid effect for the variance otherwise.

Reaction-diffusion simulations
The excluded volume also has a thermodynamic effect on the reaction rates. For a high diffusion constant the system can be considered well-mixed, and the rate of reaction rate limited reactions is increased since the obstacles decrease the reaction volume. Moreover, dimers effectively exclude less space than two monomers, so that dimerization is also favored by crowding effects [13]. Diffusion-limited reactions on the other hand are impeded due to the slowed down diffusion and the increased time for the reaction partners to meet. There exist many models for reactions in the crowded cell environment, such as fractional dynamics [12,28], a power law approximation of the RREs [27], and fractional [16] and multifractional Brownian motion [19], or scaled particle theory (SPT) to compute the reaction rates using statistical physics [10,11,13,26].
To examine the excluded volume effects on the reaction rates and in particular if CA can capture the same effects as BD, we compare association, dissociation and reversible reactions when simulated with BD and with CA. All CA simulations are here performed on a Cartesian grid, since it shows the best agreement with off-lattice simulations in the pure diffusive case. We model the complex C in two ways: (i) C has the same size as A and B; (ii) C has double the size of A and B, meaning it is a 2 × 1 molecule in CA and a sphere with radius r C = √ 2r A in BD, see Table 4.1. The latter model appears to be more realistic when we simulate a system where the molecules occupy volume, since no mass is lost in a binding reaction and the resulting complex occupies the same volume as the two reacting molecules together. The restriction that molecules have to be composed of Cartesian voxels in CA, however, forces us to choose a rather unrealistic representation of C in that case. To account for the larger size of C we adjust its diffusion constant γ C 0 = γ 0 / √ 2, and we extend the CA model to simulate particles of different sizes by adding an extra time step ∆t CA,C = √ 2∆t CA for the jumps of C.
In this section a(t) denotes the number concentration of A molecules where A(t) is the number of A molecules in the system with volume V at time t.
We focus on diffusion-limited reactions by choosing the reaction probability to be one and compute 10 3 sample trajectories in 10 3 different crowder distributions and plot their mean values in all experiments. The crowders are represented as an additional molecular species, which is static and inert and hence does not actively affect the reactions.

Association events
We first consider the bimolecular binding reaction When volume exclusion is taken into account the reaction radius in Smoldyn is not well calibrated [1]. Hence, we first experimentally find the reaction rate k BD A ( where k A is the inverse of the mean binding time scaled by the system volume V . If we assume that there are initially equally many A and B molecules, a d0 = b d0 , the solution is a d (t) = a d0 a d0 k A t + 1 . (4.5) We introduce static crowders at different occupancies φ to the BD and CA simulations. Note that the produced C molecules rest in the system as moving obstacles and here have the same size as A and B (Model I). In Fig. 4.1 we plot the concentration a(t) and the difference to the dilute solution a d (t) for each φ.
For very short times the reaction speed is increased in Fig. 4.1, because initially close reaction partners are kept in vicinity of each other by the surrounding obstacles. But, generally the diffusion limited reactions simulated here are slower for higher crowder densities. Introducing the artificial lattice, leading to a faster diffusion than in BD simulations, has a non-linear effect on the reaction rates. For early times it slows down the reactions, since close reaction partners have a higher chance of escaping each other. Then, a cross-over between the CA and BD curves occurs where the increased diffusivity of the on-lattice  simulations leads to a faster encounter of initially distant reaction partners and the reaction rate is higher in the CA simulations. The discrepancy between on-and off-lattice simulations increases for higher crowder densities and higher initial concentrations. We call the latter effect self-crowding, meaning we can observe excluded volume effects even without adding explicit obstacles, since the A, B and C molecules themselves act as obstacles to one another. Since C is not actively participating in the reaction, we only simulate Model I, but the self-crowding effect of C would be increased with Model II. In Fig. 4.2 we plot the concentration a(t) simulated with both BD and CA in a highly crowded environment to explain the different phases. We only observe Phase III when the starting concentrations of A and B are 0.2 and φ = 0.4, leading to the unbiological overall occupancy of 60%. In this case the diffusivity in CA is as obstructed as in BD (see Fig. 3.2(f)), but the grid makes it more difficult for molecules to pass each other, such that the moving C molecules can permanently block A and B molecules from reacting, whereas in the off-lattice case there are more chances to pass moving obstacles.

Dissociation events
In this section we examine the dissociation reaction for the two different models for the size of C. We first choose the dissociation probability in CA to be p D = 0.1 and then compute the mean time until a dissociation event happens by From this we can derive the macroscopic reaction rate for the dissociation event We observe a self-crowding effect, since the simulations with φ = 0 differ from the theory for dilute media and this effect increases with increasing initial concentration.
Since an additional lattice site is needed for the dissociation to be successful, excluded volume effects considerably obstruct this reaction, when modeled with Model I, and a(t) decreases for increasing φ.
and use it in the BD simulations. In dilute medium the mean concentration of A then follows the macroscopic RREs, which lead to the dilute concentration When C follows Model I extra space is needed to place one of the reaction products A or B and the dissociation might be rejected if another molecule occupies the position sampled for the extra particle. This effect, however, is not possible to achieve with the software Smoldyn, as it always places the new particle in the system and decides in the next time step if a possible rebinding or a diffusive jump happen. Hence, we only examine CA simulations for different fractions of occupied volume φ and compare them to the analytic solution in the dilute medium in Fig. 4.3. Again, we observe self-crowding effects, as the simulations in the dilute medium for φ = 0 do not follow the dilute concentration (4.9). The A and B molecules here act as moving obstacles that that block the extra lattice site needed for the dissociation of C, and we see that the discrepancy between the simulations and the analytic solution (4.9) is more pronounced for an increased initial concentration c 0 . Introducing obstacles further stabilizes the complex C, so that fewer A and B molecules are created, leading to a lower concentration a(t).
When C is double the size of A and B on the other hand (Model II), the dissociation event is always successful and independent of the occupancy φ, since no additional site is needed. In Fig. 4.4 we show some sample paths of both BD and CA simulations, they both agree with the analytic concentration in (4.9). Thus, the choice of model for C plays an important role for the outcome of the dissociation reaction and macromolecular crowding has no effect on the outcome of a dissociation reaction in discrete and continuous space models if the more realistic Model II is chosen. This model is natural to implement in BD, but leads to the unnatural representation of C in CA, where molecules generally are assumed to all have the same size.

Reversible reactions
We now combine our results to examine the reversible binding reaction where we choose the reaction constants in the same way as in the previous sections. To guarantee the correct rebinding probability in the Smoldyn simulations, we manually set it to 0.25 to agree with the rebinding probability in CA, where a newly produced particle jumps with probability 0.25 into the adjacent reaction partner and they react with p A = 1. The RREs for (4.10) are non-linear and hence not analytically solvable, but we can derive expressions for the steady state concentrationā of Ā a = 1 2 We first investigate Model I for C by simulating the reversible reaction only with CA for different a 0 = b 0 and c 0 = 0. In Fig. 4.5 we observe that increasing crowding stabilizes the complex C and decreases the amount of A and B molecules in the system at steady state as compared to the dilute case (4.11). There are two reasons for this: (i) a particle already occupies the lattice for the newly created molecule and the dissociation event is rejected; (ii) the rebinding time for A and B is decreased, since they escape from each other more rarely. The decreased time for A and B to meet in a crowded environment does here not effect the steady state concentrations. The self crowding effect is of the same order as the crowding effect, where the simulated steady state for φ = 0 is ∼ 0.9ā for a 0 = 0.05 and ∼ 0.79ā for a 0 = 0.2, indicating that the A and B molecules themselves heavily stabilize the complex C. . The self-crowding effect for φ = 0 stabilizes the complex C, leading to a lower concentration of A molecules than predicted by the theory for dilute media. This effect is increased when additional crowding molecules are introduced into the system.
To investigate the grid effect on the reversible reaction, we will now perform the same experiments for Model II, Fig. 4.6. As compared to the CA simulations for Model I we see that the steady state levels are slightly elevated, which is expected since the complex C already occupies the space needed to dissociate into A and B. The steady state level, however, still lies beneath the analytically predicted one for dilute media, indicating that the predominant effect is the decreased rebinding time of A and B in a crowded environment. The on-lattice simulations here underestimate the excluded volume effect as it is easier for A and B to diffuse away from each other. In the case of the reversible reaction, investigated here, the grid artifact on the reaction rates is linear with no crossover between CA and BD, and increases with increasing φ. The computational time between the two models differs considerably, for the case with a 0 = b 0 = 0.2 the CA simulations were ca. 42 faster than the more accurate BD simulations. with Model II, simulated with CA (solid lines) and BD (dashed lines). No additional space is needed to place the new molecule when the complex C decays, but since the steady state concentration of A is lower than the one predicted for dilute media, we can conclude that it is mostly the decreased rebinding time for A and B, that leads to the difference between the theory for dilute media and the simulated steady states in a crowded medium. This effect is stronger when simulated off-lattice, since it is more difficult for newly created A and B molecules to escape from each other.
Last, we investigate the variance in the reversible reaction system and how it depends on the level of crowding. We performed the same experiments as for Fig. 4.6, but with only 100 trajectories in one crowder distribution, but there was no visible difference in the variance of the overall concentration of A. In Fig. 4.7 we depict snapshots of the distributions of A, B and C molecules together with the obstacles when the system has reached steady state. Here, we observe that an increased crowder density leads to spatially more inhomogeneous distributions of A and B, since close reactants are stabilized in the complex C.

Conclusions
To understand complex gene regulatory networks in the crowded cell it is essential to perform realistic and computationally efficient reaction-diffusion simulations capturing the excluded volume effects, due to the high concentration of crowding macromolecules. In this paper we perform rigorous reactiondiffusion simulations in discrete and continuous space to compare how well the on-lattice models approximate the more accurate off-lattice models, when applied to crowded environments. Due to the computational complexity we hereby restrict our study to static crowders, to diffusion-limited reactions with reaction probability one, and to the two-dimensional case, providing insight into reaction-diffusion processes on biological membranes. In a pure diffusion system we first observe that all models result in slower diffusion for an increased fraction of occupied volume, as was expected. This effect, however, is more pronounced with off-lattice simulations, which was not expected as they are more flexible in the number of directions a molecule can move. A possible explanation to this phenomenon is, that the artificial grid orders the particles such that they effectively excluded less space. The consequences of the lattice on the outcome of diffusion limited association reactions is twofold: for short times the bimolecular reaction rates are decreased by the artificial grid, but for long times they are increased due to the faster hitting times between the reactants. For a reversible reaction we find that the onlattice simulations underestimate the crowding effects compared to the more detailed BD model. In all experiments the excluded volume effects increase for higher initial concentrations of the reactants, indicating that they themselves act as crowders, an effect we call self-crowding. When modeling dissociation or reversible reactions we illustrated that it is important to model the molecules with their actual size, a feature usually not considered in the CA models, where all molecules are assumed to be the size of one lattice site. On the other hand, we have given examples of the considerable decrease in computation time for the CA model as compared to the costly off-lattice simulations.
We observe that the overall variance of the concentration of molecules is independent of the presence of crowding molecules, but the spatial variation in a reversible reactive system is increased, since close reactants are more likely to be bound in a complex and distant reactants are more hindered to meet each other.
The main aim of this article is to investigate the accuracy of on-lattice approximations to BD simulations, since they are popular for investigating excluded volume effects. We find that the artificial lattice in the computationally more efficient CA model has significant artifacts especially on the diffusive behavior and the steady-state concentration of a reversible reaction. Hence it can be regarded as an alternative to the more accurate off-lattice simulations only for low crowder distributions.