Can adaptive mesh refinement regimes be improved with a Voronoi partition?

In this post Robert Ross, a first-year D.Phil. student at the WCMB, discusses a result from a 10-week research project on adaptive mesh refinement regimes that he undertook this summer at the Systems Biology Doctoral Training Centre.

Why are logicians, theoretical physicists and computational biologists all going to be sad on Christmas Day? They all asked Santa for something that does not exist (probably). A complete formal system, a global inertial reference frame, and unlimited computational resources are things apparently even Santa's cleverest elves cannot make.

A woeful and ill-posed (haha! get it?) joke aside, why does it matter to a computational biologist that their computational resources such as RAM and CPU performance are finite? In an ideal mathematical/computational biology research environment we would have the resources, both computational and experimental, along with the collective intelligence to model every aspect of a biological phenomenon. By ignoring nothing, our model would exactly reproduce the behaviour of the biological phenomenon, and we could claim to understand it.

However, because of the aforementioned limitations, we have to create simpler models instead. These simpler models aim to capture the essential behaviour of our system. Yet despite our best efforts, it is often not possible to make a biological system simpler without losing some essential feature. In this case we can turn to "brute-force" computational methods. These allow us to simulate complicated systems, and those not amenable to a deterministic approach. In this scenario using our finite computational resources effectively becomes critical. Methods to reduce computational load while maintaining the accuracy of stochastic simulations have been the focus of much recent research at the WCMB, including but not limited to Erban et al. (2013), Flegg et al. (2013), Franz et al. (2013) and Yates et al. (2012).

What is an adaptive mesh refinement regime?

Adaptive mesh refinement (AMR) regimes concern lattice-based stochastic simulations, in which space is quantised into a series of "compartments". Lattice-based methods are governed by probabilistic rules, which in spatial models, are encapsulated by the reaction diffusion master equation (RDME). Computationally, this is realised using a position-jump model, where species may react or ''jump" between adjacent compartments. The accuracy of a RDME is directly linked to the size of the compartments used (Isaacson, 2013).  

AMR regimes allow the resolution of a simulation to dynamically respond to the problem being modelled, such that computational resources are engaged where they are most required. For instance, a part of the domain where two reactive species meet may be of more interest than another part of the domain where no species reside. Therefore, smaller compartments could be used in this part of the domain to give greater accuracy, and larger compartments could be used elsewhere to reduce computational-load. Importantly, this area of interest may move in the domain during a simulation, meaning that the size of the compartments of which the domain is composed must be able to change. This is particularly relevant to biological modelling, with many biological processes relying on "dynamic" spatio-temporal interpretations to be effectively understood.

Given such an AMR regime, we have considered the following question: How do we best partition the domain to derive transition rates between compartments? This is the question we will address in the remainder of this blog.

Two methods for calculating transition rates

AMR regimes naturally give rise to non-standard domains, i.e. domains with compartments of different sizes.  Much like standard domains, a transition rate between two compartments will be directly linked to the size of the compartments involved in the transition.

We compared two methods, a standard finite-volume associated flux-method with interval-centred partitioning (Bayati et al., 2011), and a previously untested Voronoi domain partition method (Yates et al., 2013)).  Importantly, recent research carried out at the WCMB has suggested that to derive the correct transition rates between compartments on a non-standard domain a Voronoi domain partition must be used. For further details please see Yates et al. (2013).

We thought these two methods would give different results for two reasons:

(1) The finite-volume flux-method on a non-standard domain uses interval-centred compartment-centring. Conversely, a Voronoi partition defines cell intervals to be equidistant from the designated cell centres (i.e. Voronoi, and dependent on the chosen metric). These two different schemes can be seen in Figures 1 and 2.

Figure 1: Interval-centred domain partition. Distance between predefined compartment centres is represented by h. Here a standard interval-centred compartment is shown. Crosses indicate the compartment centre and vertical lines represent compartment boundaries.

Figure 2: Compartment length (l) correction for Voronoi partition. The compartment length and boundaries are corrected according to a Voronoi partition. Crosses indicate the compartment centre, and block vertical lines compartmental boundaries for interval-centred compartment centres. The dashed vertical lines represent Voronoi domain partition compartmental boundaries.

(2) The finite-volume flux-method only takes the distances between the centre of the compartments in which a transition event occurs, as opposed to the compartments both sides of the compartment from which a molecule is jumping.  Whereas the Voronoi domain partition method uses the compartments both sides of the compartment in question to define any transition rate.

Importantly, these two methods give different transition rates for the same non-standard domain.

We compared the performance of these two methods on a non-standard domain with a simple simulation of one-dimensional diffusion. The initial conditions for all simulations presented here are 1000 molecules in compartment 1. The results of this simulation can be seen in Figure 3.

Figure 3: Comparison of Voronoi domain partition and finite-volume flux methods for calculating transition rates on a non-standard domain. The blue circles are the stochastic simulation using the finite-volume flux-method to derive the transition rates, the red circles are the stochastic simulation using the Voronoi domain partition method for deriving transition rates and the green line is the PDE solution.  Inset: This is a magnified region of the same figure covering the non-standard area of the domain. All results presented here are averaged from 1000 repeats of the relevant simulations.

That these two methods gave similar results is surprising. We had postulated that improvements to the AMR regime could be made through the use of the Voronoi domain partition method. Transition rates in compartmental models can be calculated from first passage process arguments, or more generally, by methods which condition a transition probability by other transitions possible for a molecule. In our one-dimensional scenario, this would be a molecule leaving through the left interval as opposed to the right. The finite-volume flux-method does not do this, and yet still provides accurate results compared to the Voronoi domain partition method.  

Voronoi domain partitioning of a non-standard domain results in a compartment centre losing its intuitive "Euclidean" sense. We therefore thought it was relevant to see how the finite-volume flux-method performed on a non-standard domain where the compartment ''centres" have been arbitrarily defined. Figure 4 shows a concentration profile of this simulation with the corresponding numerical solution of the PDE overlaid.

Figure 4: A comparison of arbitrarily-centred finite-volume flux-method with the PDE solution.  The blue circles are the stochastic simulation and the green line is the PDE solution.  This is a magnified region of the same figure covering the non-standard area of the domain.

Interestingly, in this simulation the finite-volume flux-method performs only slightly worse than the interval-centred compartment centre definition (see Figure 3). Still, that it performs more poorly suggests that derivation of transition rates with this the finite-volume flux-method requires the centre of the compartment to be used. If this is the case, it carries an implicit assumption about where the other interval is positioned, i.e. equidistant from the compartment centre. Whether this is what allows it not to consider the other interval explicitly is not known, and is currently under examination.