Hybrid stochastic simulation

Hybrid stochastic simulations are a sub-class of stochastic simulations. These simulations combine existing stochastic simulations with other stochastic simulations or algorithms. Generally they are used for physics and physics-related research. The goal of a hybrid stochastic simulation varies based on context, however they typically aim to either improve accuracy or reduce computational complexity. The first hybrid stochastic simulation was developed in 1985.

History

The first hybrid stochastic simulation was developed by Simon Duane at the University of Illinois at Urbana-Champaign in 1985. It combined the Langevin equation with microcanonical ensembles. Duane's hybrid stochastic simulation was based upon the idea that the two algorithms complemented each other. The Langevin equation excelled at simulating long-time properties, but the addition of noise into the system created inefficient exploration of short-time properties. The microcanonical ensemble approach meanwhile excelled at exploring short-time properties, but became less reliable for long-time properties. By combining the two methods, the weakness of each could be mitigated by the strength of the other. Duane's initial results using this hybrid stochastic simulation were positive when the model correctly supported the idea of an abrupt finite-temperature transition in quantum chromodynamics, which was an controversial subject at the time.

Since then many hybrid stochastic simulations have been developed, aiming to overcome deficiencies in the stochastic simulations that they were based upon.

Methods

The Dobramysl and Holcman Method

The Dobramysl and Holcman mixed analytical-stochastic simulation model was published in 2018 by Ulrich Dobramysl and David Holcman, from the University of Cambridge and University of Oxford respectively. It simulates parts of Brownian trajectories, instead of simulating the entire trajectory. This approach is particularly relevant when a Brownian particle evolves in an infinite space. Trajectories are then simulated only in the neighborhood of small targets. Otherwise, explicit analytical expressions are used to map the initial point to a distribution located on an imaginary surface around the targets. This method has many possible applications, including generating gradient cues in an open space and simulating the diffusion of molecules that have to bind to cell receptors.

Principle of the algorithm

The algorithm avoids the explicit simulation long trajectories with large excursions and thus it circumvents the need for an arbitrary cutoff distance for the infinite domain. The algorithm consists of mapping the source position to a half-sphere containing the absorbing windows. Inside the sphere, classical Brownian simulations are run until the particle is absorbed or exits through the sphere surface. The algorithm consists of the following steps:

1. The source releases a particle at position [math]\displaystyle{ x(t=0)=x_0 }[/math].
2. If [math]\displaystyle{ |x(t)|\gt R' }[/math], map the particle's position to the surface of the sphere S(R), using the distribution of exit point [math]\displaystyle{ P_{map} }[/math]. In three dimensions, there is a finite probability for a Brownian particle to escape to infinity upon which the trajectory is terminated.
3. In the first time step, [math]\displaystyle{ t=\Delta t, |x^t|\gt R', }[/math] use the mapping, [math]\displaystyle{ P_{map} }[/math], to map the particle's position to the sphere S(R). This leads to a sequence of mapped position [math]\displaystyle{ x_1,..x_n }[/math] until the particle is absorbed. Note that for the mapping, there is again a finite probability that the particle escapes to infinity, in which case, the trajectory is terminated.
4. The Euler-Maruyama scheme can be used to perform a Brownian step: [math]\displaystyle{ x(t)=x(t-\Delta t)+\sqrt{2D\Delta t}\mathbf{R}, }[/math] where [math]\displaystyle{ \mathbf{R} }[/math] is a vector of standard normal random variables.
5. When either [math]\displaystyle{ x(t)\cdot \mathbf{e}_z \le 0 }[/math] (in the case of half-space) or[math]\displaystyle{ |x(t)|\leq 1 }[/math] (in the case of the sphere), and [math]\displaystyle{ |x(t)-x_i|\lt \epsilon }[/math] for any value i, one considers that the particle is being absorbed by window i and terminate the trajectory.
6. If the particle crosses any reflective boundary, go back to step 3 to generate a new position. Otherwise return to step 2.

Mapping the source for a ball in 3D

One can map the source for a ball in 3D to get the first passage probability for hitting a ball before escaping to infinity. The mapping is as follows:

[math]\displaystyle{ P_{map}(x|y)=\frac{|y|}{R}\frac{1}{4\pi}\frac{\beta^2-1}{(1+\beta^2-2\beta\cos\gamma)^{3/2}} }[/math], with [math]\displaystyle{ \int_{\partial B(R)}P(x,y)dS_{x}=\beta^{-1}=\frac{R}{|y|}, }[/math] and [math]\displaystyle{ |x||y|\cos\gamma=x\cdot y. }[/math]

The probability distribution of hitting is obtained by normalizing the integral of the flux.

Remarks

The choice of the radius R is arbitrary as long as the sphere S(R) encloses all windows with a buffer of at least size [math]\displaystyle{ \epsilon }[/math]. The radius R' should be chosen such that frequent re-crossings are avoided, e.g. [math]\displaystyle{ R'\leq R+10\sqrt{2D\Delta t}. }[/math] This algorithm can be used to simulate trajectories of Brownian particles at steady-state close to a region of interest. Note that there is no approximation involved.

Two-regime method

The Two-Regime Method for reaction–diffusion simulations was created by Mark Flegg, Jonathan Chapman and Radek Erban at the University of Oxford. It combines molecular-based algorithms with compartment-based approaches at ideal points during calculations to reduce computational cost. The molecular-based algorithms are great at giving highly accurate detail on localized regions of interest. Compartment-based models excel at efficient simulations of large regions. The main use for this model is to increase both the speed and accuracy of reaction–diffusion simulations, and provide more control to the simulator over methods to characterize regions of interest.

Principle of the algorithm

The Two-Regime Method works by having two regimes of interest. One region is event-based and primarily uses compartment-based approaches, while the other region is time-based and relies on molecular-based regimes. The steps of the algorithm are as follows:

1. Divide the computational domain [math]\displaystyle{ \Omega }[/math] into two parts. The parts should be non-overlapping.
2. Determine which part of the domain would be better fit for a compartment-based approach, label that [math]\displaystyle{ \Omega_C }[/math]. The other domain will be [math]\displaystyle{ \Omega_M }[/math]. Ensure the following is true: [math]\displaystyle{ \Omega = \Omega_M \cup \Omega_C }[/math].
3. Treat molecules in [math]\displaystyle{ \Omega_M }[/math] as free molecules in a continuous space. These molecules will diffuse and react through molecular-based approaches.

Molecules jump between compartments while in region [math]\displaystyle{ \Omega_C }[/math] with the chance of jumping into [math]\displaystyle{ \Omega_M }[/math] where movement will then be simulated using Brownian motion. Many possibilities exist to couple these regions, which can vary based on the purpose of the simulation.

Remarks

This algorithm and ones built upon it are used to study the conversion of species. They can also be coupled with the Fokker-Planck equation to simulate population and single trajectories using Brownian simulations.

Applications

Hybrid stochastic simulations have been used to:

• Predict the impact of HIV-prophylaxis on patients, aiding in the development of exposure prophylaxis medications.
• Model tumor suppression mechanisms in cancer research.
• Simulate train trajectories, which helps in the development of railway traffic schedules.