A stochastic cellular automaton simulation of chemical oscillations in small systems

We propose a stochastic cellular automaton method to simulate chemical reactions in small systems. Unlike the standard Gillespie method, which simulates chemical reactions with a few thousand molecules reacting with each other but without spatial considerations, our systems are divided into independent cells, each containing only a few molecules. Our simulation of the Brusselator produces chemical oscillations that agree extremely well with solutions to deterministic rate equations, and we can see strong oscillations in systems with as few as 10 cells. We are able to study several factors that affect the robustness of these small chemical oscillators: system size, spatial distribution, and correlation of molecules. We have found that non-Poisson particle distributions can greatly suppress chemical oscillations and that chemical reactions can induce correlation between the spatial distributions of particles of different species and create large-scale inhomogeneity in particle concentrations. In addition, incomplete oscillations (misfirings) can appear among strong, regular oscillations when the system size is smaller than a certain threshold, and these misfirings are triggered by random events, with a probability that is related to the system size. Since these effects, resulting from several different physical causes, are difficult to accurately model by adding generic noise factors to deterministic rate equations, as is frequently done in theoretical studies, we argue that our stochastic cellular automaton method is a useful addition to the existing tools for studying small, inhomogeneous, and non-equilibrium reaction-diffusion systems, especially those of biological nature.


I. INTRODUCTION
Oscillatory systems are ubiquitous in the natural world, ranging from simple physical systems, such as pendulums, to complicated biological systems, such as cellular cycles. 1,2 What these oscillators have in common is a restoring mechanism that can pull the system back to a certain configuration-a fixed point-when the system moves away from it. The characteristic frequency of an oscillator depends on system parameters and, in the case of nonlinear oscillators, on initial conditions. Chemical oscillations can occur in systems with numerous chemical species that undergo multistep reactions. In well-stirred systems, the concentrations of the species can change periodically with time, as seen in the Belousov-Zhabotinsky reaction, one of the most famous oscillatory chemical reactions. 3,4 Because the laws of thermodynamics dictate that a closed system eventually reaches a state of equilibrium, where the free energy is a minimum, an oscillatory chemical system is obviously not at equilibrium. In addition, these systems are nonlinear, and the restoring or feedback mechanism is often an autocatalytic reaction, in which a product is also a reactant. 5,6 The scale of biological oscillations varies widely from individual cells to entire ecosystems, with periods ranging from a fraction of a second to days and years. Much progress has been made in understanding the mechanisms responsible for various biological rhythms at the level of enzyme activity, protein synthesis, and gene regulation and expression. 2 Even though one might regard biological oscillators as more complicated versions of chemical oscillators, there is nevertheless a crucial difference in terms of the numbers of molecules involved, due to the difference in the sizes of the molecules. For example, the average number of proteins in an E. coli bacterium, which has a volume of about 1 µm 3 , is about 3 × 10 6 . 7 Consequently, for a proteome size of about 5000, the average number of any particular type of protein in an E. coli bacterium is less than 1000. By contrast, a 100-mM solution of electrolytes contains about 10 8 ions in the same volume. This difference in the numbers of molecules has a major impact on determining the appropriate theoretical tools for studying these systems. Generally speaking, understanding the robustness of biological rhythms requires a good understanding of fluctuations and correlations due to the small numbers of molecules in biological oscillators. 2 The conventional method for studying the kinetics of chemical reactions involves using rate equations, which are ordinary differential equations (ODEs) that describe the rate of change of the concentration of each chemical species based on the chemical reactions involved. The rate equations assume homogeneous concentrations of all chemical species, and the solutions to the rate equations are deterministic for a given set of initial conditions. Applying rate equations to biological systems, however, is problematic, because the fluctuations that naturally arise from small particle numbers are not reflected in these deterministic equations and spatial homogeneity is difficult to achieve in a biological cell crowded with macromolecules. 8 In the 1970s, Gillespie proposed a stochastic algorithm to simulate chemical reactions, and this algorithm has become the de facto standard in the simulation community. 8,9 In the Gillespie algorithm, two random numbers are generated in each simulation step. One of these numbers is used to determine which chemical reaction will occur, with probabilities weighted by the reaction rate constants and particle concentrations. The other number is used to determine the duration of the time step. The simulation terminates when no more reactants remain. Gillespie proved that this algorithm is statistically fully equivalent to rate equations but with the added benefit of having the built-in impact of stochastic events. Similar to deterministic chemical rate equations, the Gillespie algorithm assumes well-mixed, uniformly distributed chemical species. A typical simulation usually involves a few thousand reactants of each species, and averaging a few simulation runs may also be necessary. 8 Despite the enormous success of the Gillespie algorithm and several subsequent modifications that are computationally more efficient, [10][11][12][13][14] the issues of how to treat spatial inhomogeneity and take into account diffusion of molecules remain an ongoing research effort. 15 Methods of simulating reaction-diffusion in heterogeneous systems usually involve dividing systems into many subvolumes and allowing chemical reactions to occur in the subvolumes with diffusion of molecules between subvolumes. 16 For example, in the binomial tau-leap spatial stochastic simulation algorithm, 17 at each time step a chosen subvolume can undergo either a chemical reaction using a Gillespie-based method, or molecules in the subvolume can diffuse into neighboring subvolumes, depending on various parameters specified by the algorithm. In other words, simulation is performed sequentially in the subvolumes, one subvolume for each time step.
In cellular automaton simulations, on the other hand, systems are divided into multiple, independent units called "cells," and events occur simultaneously in discrete time steps in each cell following rules that are applied uniformly to all cells. Adapting Gillespie's method to cellular automata, however, is difficult because time steps in the Gillespie algorithm are determined stochastically and are uneven. This unevenness in time steps does not allow uniform advances in time for all cells, as is customary in a cellular automaton simulation. 18 Nevertheless, researchers have developed several methods with different levels of complexity to study reaction-diffusion systems [18][19][20][21][22] (especially surface reactions [23][24][25][26] ) using cellular automata. Direct comparisons of simulation results and solutions to differential equations, however, are uncommon, 27 except for simple reactions. 22,28 Our method uses standard cellular automaton techniques to divide the "world" into multiple, parallel "cells," each occupied by chemical "agents," and the status of the world is updated in discrete time steps. During each time step, chemical reactions are allowed to take place in each cell independently following rules that are stochastic in nature but that are applied to all cells simultaneously. Chemical agents are also allowed to move between neighboring cells. Our method has elements similar to Gillespie's method and an earlier method by Nakanishi, 29 but with a few major differences in terms of the numbers of molecules involved in each reaction unit and how reactions and time steps are chosen. A major advantage of our method is that it can incorporate spatial inhomogeneity of particles to simulate reaction-diffusion systems when stochastic events are important, and our method is also easy to implement. In the following sections, we will give a rationale for our method and test our algorithm using simple chemical reactions of various orders. We will then extend our method to more complicated chemical oscillators that involve multi-step chemical reactions. We will make direct comparisons of our simulation results with solutions of deterministic rate equations, which we regard as accurate descriptions of the chemical kinetics of macroscopic, homogenous experimental systems. We will also investigate the robustness of these simulated chemical oscillators when the system size is reduced and will report findings that are characteristic of small systems and that are not seen in larger systems. In addition, we will show that chemical oscillations are suppressed or distorted in the non-Poisson particle distributions that we have examined. Even though the scope of the current study is limited principally to homogeneous systems, we hope to further develop this method for application to spatially inhomogeneous reaction-diffusion systems -and particularly to networks of biochemical systems -where fluctuations due to small system size can play critical roles.

II. A STOCHASTIC CELLULAR AUTOMATON SIMULATION OF CHEMICAL REACTIONS
We will introduce our simulation method by examining standard chemical reaction rate equations and giving them statistical interpretations. We will then test our simulation algorithm using specific examples with increasing complexity. Our simulation methods for first-and second-order reactions are similar to those used by other authors, 22,28 but our approach to multi-step reactions is different from the methods of other researchers, to the best of our knowledge.
Since our algorithm simulates chemical reactions on a molecular level, all concentrations used in this paper are number concentrations. We will refer to a reaction that involves a minimal set of reactant molecules as a "stoichiometric reaction." For example, the stoichiometric reaction for the reaction 2A + B → C involves a minimal reactant set of two A molecules and one B molecule. When molecule numbers exceed this minimal set, an important factor, as we will see later, is how many distinct reactant sets, or microstates, these molecules can form. Since the thermal de Broglie wavelength is less than an angstrom at room temperature for all molecules that are of interest here, we will treat the molecules as distinguishable, classical particles when enumerating distinct reactant sets, and we will denote the total number of possible reactant sets that can be formed from a given group of reactant molecules as the "multiplicity factor" (Ω), a concept also used in statistical mechanics.

A. Spatial distribution of molecules
We assume that our chemical system is divided into L cells. For computational simplicity, we assume that the cells lie in one dimension, but our method can easily be generalized to two or three dimensions. If we have N molecules of the same species, then the probability of finding a particular molecule in a cell is 1/L in a homogeneous system. One can easily prove that for large values of N and L, the number of molecules in the cells follows a Poisson distribution: 30 where P(n) is the probability of a cell having n molecules and λ = N/L is the average number of molecules per cell. Note that this distribution depends only on n and λ and does not depend on N and L explicitly. For large values of λ, a Poisson distribution approaches a normal distribution. In this work, however, we are mostly interested in cases where λ is small. To achieve a Poisson distribution in our simulation, we can pick a random cell from L cells and deposit a molecule in it and repeat the process as many times as we have molecules, or we can let molecules perform random walks (RW) multiple times. Our simulation and analysis show that when λ is small, a large number of RW steps are required for a non-Poisson distribution (for example, a flat distribution) to become Poisson. We will experiment with both methods in our work. The first method is computationally less time-consuming, but the second method mimics a diffusion process and is physically more realistic.

Example A first-order chemical reaction
We consider a simple first-order chemical reaction with a forward rate constant k f and a backward rate constant k b : The rate equations for this reaction are where A and B are concentrations of chemical species A and B and are functions of time. (We will use upper-case, italicized letters for concentrations.) To give a statistical interpretation to the above rate equations, we rewrite the differential rate equation for A for the forward reaction (dA/dt = −k f A) as a difference equation for a short time step τ: This equation states that during a time interval τ, the fractional reduction in the number of the molecules is equal to k f τ. Interpreting this fractional change as a probability, we conclude that statistically, the rate equations for the first-order reaction (2) are equivalent to the following statement: During a small time interval τ, any molecule of species A has a probability k f τ of becoming a molecule of species B, and any molecule of species B has a probability k b τ of becoming a molecule of species A (k f τ, k b τ 1; first-order reaction A B).
Since this first-order reaction does not involve any A-B combinations, simulating this reaction is reasonably straightforward using a Monte Carlo method (please see the Appendix for details), and our results confirm previous work by other researchers. 31

Example A second-order chemical reaction (different reactant species)
The following second-order chemical reaction involves two reactants and one product and has a forward rate constant k f and a backward rate constant k b : The rate equations for this reaction are: For the initial conditions A(0) = B(0) = A 0 and C(0) = 0, the analytical solutions to the above rate equations are where and We again rewrite the differential rate equation for A for the forward reaction (dA/dt = −k f AB) as a difference equation for a small time step τ: Comparing the above expression with Eq. (4), we note the appearance of an additional factor of B on the right-hand side here. This indicates that compared to the first-order reaction (2), the probability of a molecule of species A participating in this second-order reaction is increased by a factor of B, the number of molecules of species B in a unit volume. Another way to write this equation is The left-hand side of the above equation is equal to the number of times the stoichiometric reaction will occur in a unit volume (because each reaction consumes one A molecule), and the product AB on the right-hand side is numerically equal to the number of distinct A-B molecular pairs in the same volume. We therefore give the following statistical interpretation to the forward reaction in the second-order reaction (5): In a volume V 1 , the probability that the stoichiometric reaction A + B k f → C will occur during a small time interval τ is equal to k f τ/V 1 , and the maximum number of times that this reaction can occur is equal to the multiplicity factor, which is equal to the number of distinct A-B molecular pairs in the volume:

are the numbers of molecules of species A and B, respectively, in volume
Since rate constants have different units for reactions of different orders, we have included the volume explicitly in the above statement to ensure that the units are correct. In our simulations, however, for convenience, we set V 1 , the volume of each cell, to one unit unless otherwise noted.
Following the steps in the Appendix, we simulate reaction (5), and Fig. 1 shows both the simulation results and the analytical results given by Eq. (7). The agreement is very good.

Example A second-order chemical reaction (same reactant species)
The following second-order chemical reaction involves reactants of the same species and one product: Assuming the forward rate constant is k f and the backward rate constant is k b , the rate equations for this reaction are: For the initial conditions A(0) = A 0 and C(0) = 0, the analytical solutions to the above rate equations are where We now try to give a statistical interpretation for the forward reaction in the rate equation . Suppose that in the ith cell, there are N A,i molecules of species A. If we follow arguments similar to those presented in Example 2, we would then conclude that the multiplicity factor should be N 2 A,i . However, the number of distinguishable pairs in this cell is actually N A,i N A,i − 1 /2. It seems as if two issues need careful consideration: Is enumerating particle pairs the correct statistical interpretation for rate equations when stoichiometric coefficients are larger than 1, and does the order of the molecules in a pair matter? For example, if there are two A molecules in the cell (N A,i = 2), should the multiplicity factor Ω i be 1, 2, or even 4? We will answer this question by running computer experiments with the values of Figure 2 shows simulation results of reaction (9) compared to the analytical results given by Eq. (11). We can see that only Ω i = N A,i N A,i − 1 produces results that agree with analytical solutions. This means that if, for example, we have two A molecules (N A,i = 2) reacting with each other, then the multiplicity factor is Ω i = 2. In other words, the order of the molecules in the pair does matter. In the previous example for reactants of different species, the multiplicity factor is Ω i = 1 for a molecule A reacting with a molecule B (N A,i = 1 and N B,i = 1); the order of the molecules in the pair does not matter in this case. This inconsistency in the ways we count A-B pairs and A-A pairs, however, is simply the result of how the rate equations are written and is not due to some conflict with fundamental physical laws. Instead of using a consistent counting method and introducing new rate constants that are defined differently from the ones in conventional rate equations (as Gillespie did 8,9 ), we choose to use the conventional rate constants but give them different statistical interpretations depending on whether the reactant species are the same or different. (

Example A third-order chemical reaction
For a third-order reaction, we consider Assuming the forward rate constant is k f and the backward rate constant is k b , the rate equations for this reaction are: We will solve the above differential equations numerically. To simulate this reaction, we use the multiplicity factor Ω i = N A,i N A,i − 1 N B,i for the forward reaction in the ith cell (i = 1, 2, . . ., L). In Fig. 3, we plot simulation results and numerical solutions to Eq. (13); we again see excellent agreement between the two. Furthermore, if we set the volume of a cell to V 1 = 2 and double the initial molecule numbers (thus keeping the initial concentrations unchanged), we see no change in Fig. 3, indicating that the volume normalization used in calculating the reaction probability (step 4 in the Appendix) is appropriate.
Another useful result of the rate equations (13) is that at equilibrium (where dA/dt = dB/dt = dC/dt = 0), the following product-to-reactant concentration ratio is a constant: Here K is the equilibrium constant, and this is called the law of mass action in chemistry.
To test this, we run our simulation using randomly generated initial concentrations for each chemical species. Equilibrium concentrations A ∞ , B ∞ , and C ∞ are obtained by averaging the concentrations over many time steps after the system has reached equilibrium. In Fig. 4, we plot simulation results of C ∞ versus A 2 ∞ B ∞ , which are in good agreement with the law of mass action. It is worth noting that the concentrations in the simulations in Fig. 4 are so low that it is much less likely to find two A molecules in a cell than it is to find one. Since it takes a minimum of two A molecules in a cell for the forward reaction in reaction (12) to occur, this means that the forward reaction cannot occur in most cells. The calculation of average concentrations, however, uses molecules in all cells, including these "inactive" cells. Nevertheless, the statistics still work out extremely well for our simulation. Furthermore, rate equations assume that chemical concentrations can make arbitrarily small changes over arbitrarily small time intervals. In contrast, particle numbers  in our simulations are integers, and time steps are discrete. It is reassuring that if we count the particles correctly using appropriate multiplicity factors, our simulations do produce results that agree well with those obtained from continuous equations.

C. Simulation of multiple chemical reactions
In the simple chemical reactions discussed in the previous section, each chemical species can only participate in one reaction, and the reactants and products are of different chemical species. In this section, we extend our simulation to more complicated systems in which a single chemical species can participate in multiple chemical reactions and product molecules can contain the same species as reactant molecules (autocatalytic reactions).
The example we use is the Brusselator, a name given by Tyson 32 to a theoretical system of chemical oscillators developed by Lefever, Nicolis, and Prigogine in Brussels. 33 A Brusselator consists of the following four chemical reactions: We can see that in reaction (15), molecules X are produced by precursors A, but X can also decay to C through reaction (18). In the presence of B, X can produce Y through reaction (16). Reaction (17) is autocatalytic, because X is both a reactant and a product. In this reaction scheme, X is considered an "activator," because it produces, or "activates," Y, and Y is considered an "inhibitor," because it consumes, or "inhibits," X. For the purposes of our discussion here, we ignore backward reactions and assume that some external agents keep the concentrations of reactants A and B constant.
We can write individual rate equations for each reaction in the Brusselator and combine these to obtain the following overall rate equations for X and Y : and These coupled nonlinear differential equations for X and Y can be solved numerically for given rate constants, values of A and B, and initial conditions. Before doing this, however, it is helpful to review some general properties of the Brusselator. First, we note that it is possible for X and Y to remain constant with time (dX/dt = dY /dt = 0). Setting Eqs. (19) and (20) to zero, we find these steady-state values: Whether this steady-state point (also call a fixed point) is stable or not can be analyzed using standard techniques for dynamical systems. 1 Linearizing the system near the fixed point, we find the Jacobian of the system at the fixed point (which involves computing the partial derivatives of dX/dt and dY /dt in Eqs. (19) and (20) with respect to X and Y and evaluating at the fixed point): Since the determinant of the Jacobian here is positive, the fixed point here is unstable when the trace of the Jacobian is positive. In other words, we expect to see chemical oscillations in the Brusselator when For a first simulation of the oscillations of a Brusselator, we choose the parameters A = 3, B = 11, and k 1 = k 2 = k 3 = k 4 = 1. Numerical solutions to the rate equations (19) and (20) using these parameters are shown in Fig. 5(a). We can see that both X(t) and Y (t) have periodic sharp changes over very narrow time durations, which can be interesting -and challenging -to simulate.
We can also see in Fig  concentration of the activators (X), in the meantime, decreases steadily to a very low value and stays there for a while. When the concentration Y reaches a high enough level, however, the autocatalytic reaction (17) starts to occur at a rapid rate, and the system sees a sharp increase in X and an equally sharp decrease in Y within a very narrow time window. The cycle then repeats. We note that an activator X in the Brusselator can participate in several different reactions. In a real chemical reaction, molecules move about randomly in space, and chemical reactions may occur when molecules come into contact with each other. If there are multiple chemical reactions, there is no specific order in which the reactions should occur. After some experimentation, we have found that the following method produces results that compare well with numerical solutions to rate equations (also see Fig. 6): 1. We randomly deposit molecules of species X and Y into L cells to achieve the desired initial concentrations. The numbers of molecules of species A and B, however, are the same in all cells and are kept constant throughout the entire simulation process. 2. During each time step of interval τ, we simulate chemical reactions in each cell independently.
In each cell, we randomly pick a set of four reactions from reactions (15)- (18) and simulate the selected reactions consecutively; some of these reactions may be duplicates by chance. For each of these reactions, the rate constant and the multiplicity factor for the particular reaction in question are used to calculate how many times the stoichiometric reaction should happen, and molecule numbers are updated accordingly. (Please see the Appendix for details.) 3. After all of the selected reactions are simulated in each of the cells, the average concentrations are calculated using where N X ,i and N Y ,i are the numbers of X and Y molecules in the ith cell, respectively. These molecules are then redistributed randomly to achieve a Poisson distribution. Steps 2 and 3 are repeated many times until a desired total simulation time has been reached.
Simulation results using the same parameters as those in Fig. 5(a) are shown in Fig. 5(b); the agreement between simulation and numerical results is excellent. The reason that we are able to use four sequential reactions to simulate four simultaneous reactions is because we choose τ such that the probability of any single reaction occurring (kτ/V n−1 1 ) is so small that the chance of multiple reactions actually taking place in the same cell is almost negligible. In fact, when we choose four fixed reactions in the order (15)-(18), we see no significant difference in the simulation results compared to those presented in Fig. 5(b). Simulating a fixed set of reactions would certainly reduce computation time, but we opted not to do this because of its potential to introduce bias into systems by preferring some reactions to others.   It is worth pointing out that our method of choosing reactions randomly with equal probabilities during each time step works even when the reaction rates are different. In Gillespie's method, on the other hand, reaction rate constants affect which reactions are more likely to be chosen. In addition, a chosen reaction will occur in Gillespie's method in a time interval determined by another random number, whereas in our simulation, whether a chosen reaction will occur -and how many reactant molecules will participate in the reaction -is determined by the Monte Carlo simulation in step 4 in the Appendix. In a sense, our method is less efficient than Gillespie's method, because in many cells, no reaction actually occurs in a given time step. But the advantage of our method is that we can use it to explore spatial parameters by associating different cells with different spatial coordinates, whereas Gillespie's method does not have any "built-in" spatial variables.

A. The effects of particle distributions
We have shown that our stochastic cellular automaton simulation method can produce results that compare well with numerical solutions to deterministic rate equations. Now we want to address some specific aspects of our simulation.
Firstly, how important is it to maintain a Poisson particle distribution when simulating a homogeneous system? We explore this issue by conducting several computer experiments. Naively, we may expect a flat particle distribution to also accurately mimic a well-stirred chemical system (allowing a maximum variation of 1 among particle numbers in different cells, of course, because particle numbers are integers and average concentrations can be non-integers). However, when we change the Poisson distribution in our simulation in Fig. 5(b) to a flat distribution (but leave all other simulation parameters unchanged), we obtain the results shown in Fig. 8(a), where no chemical oscillations but with periods and amplitudes larger than those obtained by solving rate equations numerically (Fig. 5(a)). We do not see such V 1 dependence, however, in simulations using the Poisson distribution.
are observed. The reason is not hard to understand: The autocatalytic reaction (17) requires two X molecules. By forcing a flat distribution, however, no cell can have more than one X molecule after the average concentration X drops below 1, and the autocatalytic reaction (17) can no longer occur. As a result, Y continues to grow due to reaction (16) while X remains constant at a very low value, and chemical oscillations cannot occur. The problem of not having enough X molecules in a cell for the autocatalytic reaction can be solved, of course, by assuming a larger cell volume V 1 , thereby having more molecules in the cells to maintain the same concentrations. When we set V 1 = 4, we indeed see chemical oscillations in our simulation with a flat distribution ( Fig. 8(b)), but the amplitude and period are significantly larger than the expected values from numerical solutions to rate equations ( Fig. 5(a)). When we further increase V 1 to 14 (Fig. 8(c)), the oscillation amplitude and period are closer to the expected values, although they are still a bit larger. If we use a Poisson distribution with these increased values of V 1 , however, we do not see such V 1 dependence and continue to obtain consistent simulation results that all agree with numerical solutions to rate equations. This is an indication that the Poisson distribution is appropriate for capturing the essence of a homogeneous system from a statistical point of view, even when the numbers of molecules in each cell are very low. In other words, the statistics work correctly when the probability of finding a molecule in any cell is the same for all cells (Poisson distribution), but this does not necessarily mean that every cell actually has nearly the same number of molecules (flat distribution). The distinction between these two distributions disappears when large numbers of molecules are present in each cell, which explains why the simulation error using flat distributions decreases when V 1 becomes larger.
We have also tried to redistribute particles after each simulation step using lognormal distributions. A lognormal distribution is asymmetric, with a long tail towards the high-number end. This long tail ensures that the autocatalytic reaction, unlike the simulation using a flat distribution discussed in the previous paragraph, will occur in some cells during the period when X is small. Figure 9 shows the results using lognormal distributions with different variances, and these results deviate significantly from those using a Poisson distribution (Fig. 5(b)), even though all other simulation parameters are the same. We can see that for a variance of 1, stable chemical oscillations are produced with a lognormal distribution, but the period and the amplitude of the oscillations are smaller than those in Fig. 5(b), where the distribution is Poisson. When the variance is set to 2, chemical oscillations are still visible, but the amplitude and period are even smaller and more irregular. A variance of 3, however, almost completely destroys chemical oscillations, although the fluctuations in X and Y seem out of phase. Thus, our computer experiments indicate that of the particle distributions we have tested, only the Poisson distribution produces accurate chemical oscillations in our cellular automaton simulation of small systems. Interestingly enough, however, using a lognormal redistribution after each simulation step does not alter simulation results of reactions that reach equilibrium (reactions in Section II B).

B. Random walks between simulated chemical reactions
There are, of course, other ways to achieve Poisson distributions besides reassigning particle positions after each simulated reaction step. In a real chemical system, molecules undergo diffusion motion, which can be simulated using random walks (RWs). For each RW step of our simulation, every particle has the same 1/3 probability to move to the cell on the left, to move to the cell on the right, or to stay in the same cell; periodic boundary conditions are applied. 34 Figure 10 shows simulation results when particles are allowed to undergo 5 and 500 RW steps at the end of each chemical reaction time step. After completing the RW steps for each time step, we collect a set N X,i of X molecule counts in each of the L cells and a set N Y ,i of Y molecule counts and calculate the Pearson correlation coefficient between the two sets. From Fig. 10, we can see that chemical reactions cause X and Y particle numbers to become correlated, but more RW steps help to reduce this correlation and produce better chemical oscillations. With 5 RW steps, the average X-Y correlation coefficient in Fig. 10(a) is -0.44, and chemical oscillations are greatly suppressed. In Fig. 10(b), we can see that even though 500 RW steps produce simulation results very close to those obtained by randomly reassigning particle positions (Fig. 5(b)), pronounced negative X-Y correlations are still evident during the sharp population reversal periods, with peak values of about -0.8. 500 RW steps after each time step. The correlation coefficient between the numbers of X and Y molecules in each cell is also plotted. Every 10th simulation data point is shown. All simulation parameters are the same as those in Fig. 5(b). Chemical reactions can cause a strong X-Y correlation that can suppress chemical oscillations. This X-Y correlation can be reduced by taking more RW steps for the particle redistribution. The correlation, however, remains strong during the sharp population reversal periods, even with 500 RW steps.
This negative X-Y correlation indicates that in cells with fewer X molecules, there are more Y molecules, and vice versa. After reaction (16), a cell will have less X and more Y, and after reaction (17), a cell will have less Y and more X, both causing negative X-Y correlations. RW steps tend to randomize particle distributions and therefore reduce X-Y correlations. During the sharp population reversal periods, however, this correlation seems difficult to erase, even with 500 RW steps. A closer examination of the distributions of X and Y molecules during this transition period reveals the presence of spatial inhomogeneities with characteristic length scales that can reach a few hundred cells. A few hundred RW steps are obviously insufficient for erasing X-Y correlations under these circumstances.
Our simulation results indicate that a strong correlation between the inhibitor and activator populations plays a role in suppressing chemical oscillations. This, however, is not the only factor. Chemical oscillations are greatly suppressed by lognormal redistributions of particles (Section III A), even though the X-Y correlation is also statistically zero there. Since simulating diffusion by RW is more computationally time-consuming, we argue that reassigning particle positions using a Poisson distribution is a more efficient way to simulate a homogeneous system.

C. The robustness of chemical oscillations in small systems
After gaining some insight into factors that suppress chemical oscillations in the Brusselator in our simulation, we now examine how robust the oscillations are when the size of the system is reduced. Figure 11 shows results of simulations using L = 10, 100, and 1000 cells. Simulation parameters used in Fig. 11(c) are the same as those in Fig. 5(b), but the total simulation time is longer. We can see that with 1000 cells, simulation results are almost identical to solutions to rate equations. It is quite remarkable that the oscillations remain steady and strong even with 100 cells, although some peaks decay prematurely, resulting in incomplete cycles ("misfirings"), as seen in Fig. 11(b). These misfirings, however, do not seem to affect the height or the duration of the subsequent full peaks. With only 10 cells, we still see strong oscillations, but the period and amplitude become more irregular (Fig. 11(a)). Because of computational power limitations, computer simulations are often carried out on small systems, and it is a common practice to average simulation results of small systems to gain better statistics, with the expectation that the averaged results will more closely resemble the behavior of larger systems. Let X and Y be the average values of X and Y, respectively, over many simulations of the same system. Figure 12 shows X and Y from 1000 simulations for (a) a 10-cell system and (b) a 100-cell system. Both X and Y decay to constant values after a few oscillations, and the functional forms of the curves have little resemblance to those in Fig. 11. Other researchers have reported similar damped chemical oscillations in small systems after averaging over several thousand simulations. 35 We argue, however, that the meaning of the mathematical exercise of averaging for chemical oscillators needs careful examination: The damping of the oscillations is caused by destructive superposition due to fluctuations in the periods of the oscillations and the phase shifts that occur after misfirings, and not by the lack of oscillations, which are strong at all times, as seen in Fig. 11.
To further investigate the misfirings, we ran a simulation for L =100 cells for a total time duration of 4000 units and found 198 misfirings. We then divided the total time duration into 160 smaller consecutive intervals of ∆t = 25 units each and counted the number of misfirings in each of these intervals. Figure 13 is a histogram of the results, with the misfiring counts in an individual time interval on the horizontal axis and the number of times these counts occur over the total time duration on the vertical axis. Since we have an average of 198/160 = 1.24 misfirings per 25-unit time interval, we can also use Eq. (1) to calculate the predicted frequency of misfirings according to Poisson statistics. We see in Fig. 13 that our simulation data fit very well to a Poisson distribution, indicating that these misfirings are random events. This also means that we can find the probability of misfirings, but cannot predict precisely when they will actually occur. It is not a coincidence that it takes about 20 units of elapsed time for the averaged oscillations in Fig. 12(b) to die out, because we expect about one misfiring during that amount of time for a system of L = 100 cells. The inset in Fig. 13 shows that the frequency of misfirings decreases when the system size L increases. The linear fit to the data seems to indicate that there should be no misfirings for L greater than ∼ 1400, and indeed we do not see misfirings in our simulations of 2000-cell systems.
Instead of averaging over many simulation runs, a more appropriate method for analyzing the behavior of chemical oscillators in small systems is to compile histograms of amplitudes and periods, as shown in Fig. 14 second, and third rows, respectively. In the Y amplitude histograms, we can see two clearly distinct clusters, with the smaller clusters on the left due to misfirings. 36 The period histograms are plotted in the column on the right, with the periods associated with misfirings excluded. The means and standard deviations of the amplitudes and periods for full cycles are listed in the figures, as well. The amplitude (Y ODE = 13.84) and period (T ODE = 6.10) obtained by numerically solving the Brusselator rate equations (Eqs. (19) and (20)) are indicated in the histograms using dotted vertical lines for comparison. In the column on the left, we have included trajectories of Y -vs.-X in phase space, which also show larger fluctuations in smaller systems. When L = 1000 cells, the mean values for amplitudes (13.82) and periods (6.02) produced by our simulations are almost identical to Y ODE and T ODE . When L = 10 cells, even though the oscillations look irregular in Fig. 11(a), the mean values for amplitudes (13.1) and periods (5.0) are still reasonably close to Y ODE and T ODE . When L = 100 cells, the oscillations are very robust despite the occasional misfirings. It is interesting that the misfirings for all three of the different cell numbers have a characteristic amplitude of about 5. It is also interesting that the period histograms for L = 100 and L = 1000 in the column on the right in Fig. 14 look asymmetric, as if they might be a superposition of two Gaussians of different widths. In the future, we would like to investigate this further to gain a better theoretical understanding of this asymmetry. Another interesting observation is that during the phase of the oscillation when the number of inhibitors (Y) increases steadily and the number of activators (X) is low, the average concentration of X is only about 0.34. According to the Poisson distribution (Eq. (1)), the probability of finding a cell with two X molecules (the minimum number for the autocatalytic reaction to occur) is only about 4%. This means that statistically speaking, only four cells out of 100 have enough X molecules to support the autocatalytic reaction during this phase. As we have discussed earlier, the sharp chemical population reversals and chemical oscillations cannot occur without the autocatalytic reaction. It is quite remarkable that chemical oscillations are still robust in our simulations of such small systems.

IV. CONCLUSION
We have developed a stochastic cellular automaton method to simulate multi-step chemical reactions, and our simulation results agree extremely well with analytical and numerical solutions to deterministic rate equations. We have conducted detailed analyses of the effects of fluctuations and particle spatial distributions on chemical oscillations in small systems, obtaining oscillations in the Brusselator even when average molecule numbers per cell in simulations fall significantly below 1 and when systems contain as few as 100 cells. This is quite remarkable considering that rate equations are usually written for systems with particle concentrations measured in moles per liter.
Of the different particle distribution schemes we have studied, we have found that only the Poisson distribution produces correct chemical oscillations in small systems. In addition, we have seen that chemical reactions cause the numbers of particles of different chemical species in each cell to become correlated, which also suppresses oscillations. Random walks help to reduce these correlations, but this is particularly difficult to do during the sharp activator-inhibitor population reversal period due to the development of large-scale particle inhomogeneity. Since our goal is to simulate homogeneous systems in the present work, we opt to redistribute particles after simulated reactions in each time step to reduce computation time. This approach is valid as long as diffusion of chemicals in the system is fast or the system is well stirred.
We have observed incomplete oscillations (misfirings) in our simulations of small systems, which are triggered by stochastic events but which have no effect on oscillations that reach full height at later times. Since aberrations such as misfirings are common in oscillations in biological systems, we believe our method provides a useful tool to further investigate these aberrations.
Our work has shown that system sizes play an important role in the behavior of chemical oscillations in small systems and that the characteristics of the oscillations are not "scalable," meaning that one cannot produce behaviors of larger systems by averaging over smaller systems. This means that one needs to be very careful when choosing theoretical tools to study biochemical oscillators, where particle numbers are small and particle distributions are inhomogeneous. We further argue that cellular automaton simulations, which can incorporate these kinds of inhomogeneity, have advantages over the more traditional rate equations and Gillespie's method, where spatial variations are not taken into account.
Compared to other simulation methods, our method is conceptually straightforward and has a number of notable differences and advantages. For example, when Gillespie's algorithm is extended to reaction-diffusion systems, uneven time steps for all subvolumes need to be sorted, and only one reaction occurs in one subvolume in each time step. In our method, we use fixed time steps, and as a result, we do not need to sort time steps, and multiple chemical reactions are able to take place simultaneously in all cells at each time step. By taking advantage of this independence of chemical reactions in each cell, we can vectorize computer code to improve computational efficiency, 37 even though code optimization has not been the main focus of the current work. This cellular automaton approach is also very convenient for allowing us to take spatial inhomogeneity into consideration.
In addition, our simulation method offers concrete interpretations of the meanings of chemical rate equations on a molecule-by-molecule basis, which is helpful to students of chemistry, biology, and physics.
We want to point out that despite the fact that our simulation method can produce results that agree well with numerical solutions to chemical rate equations, numerical methods have an advantage in computational efficiency over any Monte Carlo method, ours included. The particular merit of computer simulation is revealed when fluctuations due to small system size are important, because numerical solutions to rate equations are intrinsically deterministic and are therefore not appropriate for stochastic systems.
In the future, we hope to expand our work to study the effects of particle spatial inhomogeneity and diffusion in greater depth. We would also like to explore other reaction-diffusion systems, including systems with more complex oscillators that can be used to model real biological systems. We believe that our method is a useful tool for studying the effects of diffusion, particle distribution, and stochastic events in small systems, especially those of biological nature.

SUPPLEMENTARY MATERIAL
See supplementary material for sample computer code used in our simulations, please see our MATLAB that produces a simulation of the third-order reaction 2A + B C (Fig. 3). The numerical solutions are obtained using the MATLAB command ode45, and it takes about 3.6 seconds to run this simulation on a MacBook Pro using the parameters specified in the caption of Fig. 3.