Webpages about Smoldyn and spatial simulation

Page contents

Meaning of algorithm accuracy

The Smoldyn algorithms are described in several publications and in the Smoldyn User's Manual.

There are two distinct levels of approximation between physical reality and a simulation of it. First, the physical system is simplified to an ideal model system. This is the job of the modeler, who is trying to build a better understanding of the physical system. Secondly, the ideal model system is discretized or otherwise modified to enable its simulation. In general, a good simulator is one that minimizes these modifications, but that also runs sufficiently quickly to be useful. It's important to minimize these modifications because they lead to differences between the simulation results and what the ideal model would do, thus making it hard to understand how the model actually behaves.

The correct way to consider the accuracy of stochastic algorithms that are simulated with discrete time steps is to consider each simulation observation as a snapshot of the simulated system and to then ask whether these snapshots are statistically idential to those that would be observed with the ideal model system. How the simulator actually computes the observable results is irrelevent; what matters is whether those results are statistically correct. Algorithms that produce statistically correct are called exact algorithms.

Most individual Smoldyn algorithms are exact for all length time steps. However, these algorithms are typically not exact when used in combination with others. Importantly though, all Smoldyn algorithms approach exactness as the user reduces the time steps towards zero. As a result, it is always possible to reduce simulation errors and thus evaluate whether the errors that arise with larger time steps are likely to be problematic.

Legitimacy of the unbinding radius

The rational for the unbinding radius, as I put in my 2005 paper, is that we treat the separation between reactants as a bistable system, so that reactants diffuse on some energy surface, but we only consider a bond as having formed when they are closer than a binding radius and we only consider the bond as having been broken when they separate by more than the unbinding radius. This bistability does not violate the physics at all but only simplifies our observation of it. For this reason, in the limit of short time steps, reversible bimolecular reactions are physically exact (within the assumptions of the ideal model system). In other words, they get the physics completely correct, despite the fact that our choice of how we label the molecules between the binding and unbinding radii makes it look questionable. Because the physics is correct, the model is completely correct both inside the binding radius and outside of the unbinding radius, which are places where the molecule labeling is well defined, rather than bistable.

Does Smoldyn obey detailed balance?

Detailed balance means that when a system is at equilibrium, every elementary process (e.g. diffusion, chemical reaction, etc.) is exactly matched by a reverse process. It arises from the principle of microscopic reversibility. It can also be seen as time-reversibility, meaning that detailed balance is obeyed if one cannot tell if a movie of an equilibrium system is being played in the forwards or reverse direction. All physical systems obey detailed balance. For this reason, physicists often find models concerning if they do not obey detailed balance.

Most biological models, independent of simulation method, do not obey detailed balance for the simple reasons that they are not limited to elementary processes, they do not treat the entire system, and they are not at equilibrium. For example, there are many so-called futile cycles in biological reaction networks, of which a simple one would be represented by the irreversible reactions A -> B, B -> C, and C -> A. This could be a legitimate and useful biological model. However, it clearly ignores ATP and other energy sources, it comes to steady-state but not to thermodynamic equilibrium, and these are not elementary reactions (due to the missing ATP or other inputs). At steady-state, molecules will cycle through the different species in a specific sequence, thus violating the time-reversibility of detailed balance. Thus, this model violates detailed balance, whether Smoldyn or another simulator is used, but still might be biologically reasonable.

The next question is whether Smoldyn obeys detailed balance, if the model does. The answer is not quite, although Smoldyn approaches detailed balance in the limit of short time steps. Following are specific algorithm analyses.

Diffusion. Simulated diffusion is exactly reversible, so it obeys detailed balance. This is true both for molecules in solution and bound to surfaces.

Surface interactions. Reflection is exactly reversible, so it obeys detailed balance. Likewise, partially transmitting surfaces are also simulated exactly reversibly, so they also obey detailed balance. However, adsorption and desorption (which have to be included together for equilibrium) do not quite obey detailed balance. This can be seen because they almost match the physics exactly, with any length time step, but not quite. Adsorption differs in that it ignores the so-called Andrews-Bray correction, in which molecules that do not end up across a surface at the end of a time step still might have crossed and recrossed that surface during the time step. Smoldyn corrects for this by doubling the adsorption probability for those molecules that do end up across the surface, thus leading to the correct adsorption rate. However, this is a slight physical inaccuracy, which has the effect of leading to slightly incorrect concentrations very close to the surface (but it reduces to zero in the limit of short time steps). This error violates detailed balance. Another minor physical inaccuracy is that Smoldyn does not allow a molecule to both adsorb and desorb in a single time step. However, this treatment is symmetric, so it does not affect detailed balance.

Zero order reactions. Provided these are balanced with matching unimolecular decay reactions, they are exactly reversible, so they obey detailed balance.

Unimolecular reactions. - Provided that each reaction is made reversible, these are exactly reversible and so obey detailed balance.

Reversible bimolecular reactions. At first glance, these clearly violate detailed balance because reactants bind at the binding radius and unbind at the unbinding radius. Between these two radii there is a net flux of reactants towards each other, clearly violating time-reversibility. However, the situation is not so simple when thinking about the ends of time steps as discrete snapshots of the system (described above), and thus asking whether these snapshots are statistically identical to those from a system that obeys detailed balance. Again, at first glance the answer is no for the region between the binding and unbinding radii. This is true even in the limit of small time steps because there is still a net flux between the binding and unbinding radii. An even more careful treatment accounts for molecule labeling isses between the binding and unbinding radii discussed above. With this, the physics are simulated exactly, in the limit of zero time steps, despite the molecule labeling between the binding and unbinding radii. For this reason, reversible bimolecular reactions do actually obey detailed balance, in the limit of short time steps. With longer time steps, I don't believe that detailed balance is perfect, but I think it is very close.