Weighted-ensemble Brownian dynamics simulation: Sampling of rare events in non-equilibrium systems

Weighted-ensemble Brownian dynamics simulation: Sampling of rare events in non-equilibrium systems

Justus A. Kromer justuskr@physik.hu-berlin.de Department of Physics, Humboldt-Universität zu Berlin, Newtonstr. 15, 12489 Berlin, Germany    Lutz Schimansky-Geier Department of Physics, Humboldt-Universität zu Berlin, Newtonstr. 15, 12489 Berlin, Germany    Raul Toral IFISC, Instituto de Física Interdisciplinar y Sistemas Complejos, CSIC-UIB, E-07122 Palma de Mallorca, Spain

We provide an algorithm based on weighted-ensemble (WE) methods, to accurately sample systems at steady state. Applying our method to different one- and two-dimensional models, we succeed to calculate steady state probabilities of order and reproduce Arrhenius law for rates of order . Special attention is payed to the simulation of non-potential systems where no detailed balance assumption exists. For this large class of stochastic systems, the stationary probability distribution density is often unknown and cannot be used as preknowledge during the simulation. We compare the algorithms efficiency with standard Brownian dynamics simulations and other WE methods.


Rare events are ubiquitous in many biological, chemical and physical processes Moss and McClintock (1989); G. van Kampen (1992). Whereas the density of states is known in systems at thermal equilibrium, interesting phenomena often occur in non-equilibrium systems Nicolis and Prigogine (1977). Unfortunately, many such problems are inaccessible to analytic methods. Therefore computer simulations are a widely used tool to estimate the density of states or transition rates between them Mannella and Palleschi (1989); Mannella (2000). Since standard Brownian dynamic simulation Burrage (1999); Kloeden and Platen (2011) provides computational costs that are inversely proportional to the state’s probability, specialized methods Bhanot et al. (1987); Giardina et al. (2006); Dickson and Dinner (2010) have to be used to adequately sample rare events, i.e. states with low probability or low transition rates.

In the last decades, flat histogram algorithms Wang and P. Landau (2001) have been developed, allowing one to evenly sample states with highly different probabilities. These algorithms are implementations of the umbrella sampling Torrie and Valleau (1977), where each state is sampled according to a given probability distribution, the so-called umbrella distribution. Within non-equilibrium umbrella sampling (NEUS) Warmflash et al. (2007) the space of interest is divided into different but almost evenly sampled subregions. The interaction between different regions occurs solely due to probability currents between then. Whereby the probability distribution within a region is then calculated by performing Monte Carlo simulations.

In order to calculate low rates between a starting and a final state, forward flux sampling methods can be used (for a review see) Allen et al. (2009). These methods introduce a sequence of surfaces between these states and introduce walkers (copies of the system) to perform weighted trajectories according to the underlying dynamics. If walkers cross one of the surfaces, getting closer to the final state, new walkers with smaller weights are introduced. Finally, many walkers with particular small weights reach the final state. The consideration of the particular weights allows one to calculate very low rates in a finite simulation time. Recently, extensions to these methods have been developed to calculate both, transition rates using umbrella sampling Dickson et al. (2009) and probability distributions using forward flux sampling Valeriani et al. (2007) algorithms.

In this work, we present an algorithm, based on the previously developed weighted-ensemble (WE) Brownian dynamics simulations Huber and Kim (1996); Bhatt et al. (2010); Zhang et al. (2010); Bhatt and Bahar (2012), that allows one to calculate the stationary probability density function (SPDF) as well as transition rates between particular states. Like in WE simulations the space of interest is divided into several subregions and the probability for finding the system in them is calculated by generating equally weighted walkers in each region. By moving to the underlying dynamics, the walkers transport probability between the subregions. Thus, WE methods are usually applied to systems of Brownian particles moving in a potential landscape Bhatt et al. (2010); Bhatt and Zuckerman (2011).

We are interested in an algorithm which allows simulations of stochastic dynamical systems which, apriori do not obey detailed balance for probability fluxes or suppose some special topology of the flow Graham and Haken (1971a, b); Haken (1975); Risken (1984); Schimansky-Geier et al. (1985). Such are given for Brownian particles in conservative force fields under the influence of additive noise. Even canonic dissipative systems possess a vanishing probability flow if transformed to the energy as dynamic variable Haken (1975); Ebeling and Engel-Herbert (1980); Klimontovich (1980). In consequence, both systems allow exact analytic solutions of the SPDF. In contrary, we aim to develop an algorithm which does not assume that neither the deterministic nor the stochastic items (see  (1) below) underly such conditions. Thus, no information on the SPDF can be used for the simulations.

In general the algorithm can be applied to arbitrary dynamical systems of the form:


where is the number of stochastic time-dependent degrees of freedom ; the functions describe the deterministic velocities for the -th direction; represents zero-mean Gaussian noise with delta-like correlation function . The noise intensity along the -th direction is scaled by the functions , which in general depend on the vector . We are interested in high precision sampling of the stationary probability current and the SPDF of finding the system in the -dimensional cube with a finite resolution. We will specify the resolution by the number of evenly spaced supporting points along the -th direction, for which we will determine .

This article is organized as follows: In section I we introduce an algorithm, based on WE methods Huber and Kim (1996), that allows one to calculate low probabilities and rates. Afterwards, we study one- and two dimensional model systems and analyze the algorithms efficiency compared to Brownian dynamics simulation (BDS) and WE techniques.

I The algorithm

First, for sake of clarity, we restrict ourselves here to particle motion in one dimension, , under additive noise. The deterministic part of the dynamics  (1) can be always represented as a conservative force . The noise strength is scaled by the parameter , we put . For such systems, the SPDF is known to be


where is a normalization constant.

We are interested in finding numerically the system’s SPDF, , at a set of evenly spaced supporting points in a finite part of the physical space, given by . The region of interest is divided into subregions of size , the -th subregion is bounded by , , with . Supporting points are given explicitly by , , with , see Fig. 1. Here denotes the largest integer smaller than or equal to .

Figure 1: Scheme of the interpolation points and the supporting points used in the numerical calculations. In this scheme we have used .

Let us introduce the probability for finding a particle in the -th subregion at time , and the corresponding set . Initially, no information on the system is available, thus, each subregion is given an arbitrary amount of probability , simply fulfilling the normalization condition . Naturally, equilibration can be accelerated if one already has information on the SPDF of the system (Eq. (1)). In that case, one can choose close to the set , which optimally approximates the SPDF:


However, in general no such information is required.

i.1 Time evolution

After setting the initial set , the time evolution of the is performed using three different steps.

We start with a redistribution step, in which walkers (copies of the system) are uniformly distributed in each subregion. Besides their individual positions , where denotes the particular subregion and the individual walkers, each walker possesses a given amount of weight . This is nothing but the present probability in the -th subregion distributed to the walkers, which yields


Note that one does not need to introduce walkers in subregions with .

After the redistribution step has been performed, Eq. (1) is integrated for all walkers, using a Brownian dynamic simulation step and an arbitrary integration scheme. This integration step realizes the time evolution . Here walkers transport probability between the subregions. As walkers are independent of each other, it is of importance to note that the particular time evolution of each one of the walkers is due to different sample paths in the stochastic parts of the Langevin equation.

Lastly, an updating step is performed, in which the new probabilities are calculated by summing up the weights of all walkers that are currently located in the particular subregion,


In what follows, we will name the sequence of redistribution, integration and updating step as running step. After an equilibration time , the set reaches a stationary regime, where the ’s fluctuate around their mean values .

i.2 Calculating the stationary probability

The individual are estimated by averaging over a total amount of sets , , taken, after the system has reached the stationary regime, every running steps: . It turns out that the mean probabilities coincide with the stationary probability (see Eq. (3)), for compatibly chosen time step and the size of a subregion (see Sec. I.5).

Finally the SPDF on the supporting points is calculated by adding the adjacent and dividing by the size (in order to have a properly normalized PDF):


i.3 Calculation of the probability current

The stationary probability current at position can be easily calculated by adding up (with the right sign) the weights of all walkers, passing to the right and to the left per unit time. In practice should be the boundary of a subregion. If the current is given by:


and indicate these walkers which cross the boundary moving rightwards, i.e. . Alternatively, assign walkers transporting weight leftwards, .

Averaging over such estimates, taken in the stationary regime, leads to the average current , which converges towards for .

i.4 Implementation of boundary conditions

The implementation of boundary conditions for the probability current or the SPDF is straightforward. Right now, absorbing boundaries are already implemented at and , since walkers that pass these boundaries are not located in any subregion. Therefore, their weights will get lost in the next updating step. Reflecting boundary conditions at can be implemented by setting for all walkers with . Hence, the probability current at will be zero. Reflecting boundaries at can be implemented analogously.

i.5 Convergence criteria

i.5.1 One-dimensional systems

In order to ensure, that and converge towards the stationary probability distribution and the probability current of Eq. (1) for , the time step and the size of a subregion have to fulfill specific criteria. This is due to the redistribution step, where walkers are uniformly distributed in each subregion. This implicates statistical errors, since they can reach positions in a subregion that, are inaccessible or at least more improbable. Therefore, walkers can more easily escape from potential minimums or reach regions of low probability. This effectively flats the probability distribution, leading to more probability in regions of low probability, for instance, around local maximums of , and less probability in the potentials minimums. In order to overcome this problem, earlier works Huber and Kim (1996); Bhatt et al. (2010) have stored the positions and weights of all walkers. In the next redistribution step, walkers were only spaced on the stored positions according to the weights belonging to them. This requires a lot of computer memory, especially for large and .

However, we found that one does not need to store these information, if the subregions are small enough to ensure that walkers have a non-negligible probability to leave them during one integration step. As a measure of how far a walker can step, due to the fluctuations, in one time step, we use the diffusion length . Thus, the size of a subregion should be small compared to the diffusion length


The distance a walker can pass during an integration step is not only determined by , but also by the deterministic dynamics, leading to a step length in first order. Usually changes very fast at the boundaries of the simulation area, which produces high deterministic velocities and regions of low probability. Walkers can only reach these regions, if the fluctuation are strong enough to balance the deterministic force, i.e. if


for all . This leads to a condition for the time step :


Hence, lower time steps allow one to sample regions, far from the potential extrema and for instance, the tails of the SPDF. If a larger time step is chosen, the will run to zero in subregions with larger deterministic force.

Since it is often difficult to fulfill Eq. (10) in the entire simulation area, one should choose a time step, which allows one to fulfill Eq. (8).

i.5.2 Multidimensional systems

In general, our method can be applied to stochastic dynamical systems in arbitrary dimension (Eq. (1)). For the foundation of the algorithm, we refer to the Appendix  B. However, in order to ensure convergence in a finite region , the criteria for the time step and the size of a subregion along the -th direction should be fulfilled. For noise dominated directions the criteria hold, therefore, the condition for the time step (Eq. (10)) becomes:


and the criteria for the size of a subregion along the -th directions (Eq. (8)) reads:


However, there might be directions without any noise (). In that case the length of a subregion should ensure, that walkers can leave it due to the deterministic term , leading to


otherwise information on the deterministic dynamics gets lost during the redistribution step, since walker that stay into a subregion do not produce any change in the and are again randomly placed in their subregion during the next redistribution step. For equally sized subregions should be the minimum value of Eq. (13) with respect to all x in the simulation area for each of the directions. These criteria can lead to a huge number of subregions, especially in high dimensional spaces.

To appropriately reduce the number of subregions, the should be chosen in order to locally fulfill the criteria. This was implemented by a grouping algorithm, that groups original subregions into ”larger” ones as long as walkers can leave these due to the deterministic term (Eq. 13). Depending on the system, this procedure highly reduces the total number of ”larger” subregions . If the grouping algorithm was used, walkers are randomly placed in each of these ”larger” subregions and a probability of finding a walker in the corresponding area was introduced. We find that such grouping highly reduces the computational costs, since less subregions and therefore less walkers are required. Since walkers jump out of these regions until the next redistribution step starts, the algorithm still approximates the correct SPDF.

i.6 Simulation techniques

Simulations were performed on a Intel®Xeon ®CPU E31245 @ 3.30GHz processor with 16 Gb DDR-3 RAM. The algorithm described above was implemented in a program for one- and two-dimensional systems. Runs of the algorithm are specified by the time step , the size of a subregion , (, in two-dimensional problems), the simulation area, given by and (, ), the number of walkers per subregion , the thermalization time , and the number of running steps between two sets of P denoted as . The numerical integration of the Langevin Eq. (1) was done using a Heun scheme. A resolution of was used in any direction.

To compare the results with other methods, we also perform Brownian dynamics simulation (BDS) using initially uniformly placed particles in the simulation area. Integration was done using the Heun scheme Kloeden and Platen (2011); Gilsing and Shardlow (2007) with integration time step . After a thermalization time the particles positions were recorded after time intervals . The BDS was given a running time (real CPU time) which usually equals the time our algorithm needs to produce its results. After the BDS was stopped and the SPDF was calculated using the recorded particle positions. We set , , to make results comparable.

Ii Model systems and results

ii.1 One dimensional system

In order to demonstrate the implementation of the algorithm, we study overdamped Brownian motion in a bistable potential . Correspondingly, we put and in Eq. (1) which results in a bistable system which is often used to study bistable systems or stochastic resonance therein Risken (1984); Gammaitoni et al. (1998). The two stable states come up to the potentials minimums, located at and , respectively. The corresponding SPDF is given by Eq. (2). For low noise strength, the SPDF attains sharp peaks at the potentials minimums and decreases down to low values at the borders and the local maximum, for instance for .

ii.1.1 Equilibration

At first, we study the equilibration process, performing simulations with different numbers of walkers per subregion . Results are shown in Fig. 2. Analyzing the time dependence of the probability , we find that longest thermalization time occurs at the local maximum of . Note that runs with larger thermalize at lower , but one needs more integration steps.

After thermalization has been achieved, we evaluate the coefficient of variation of the probability, given by


where averages are performed for a fixed number of sets and different and . We find it to scale according to (data not shown).

Figure 2: at obtained from by for the subregion containing for and different numbers of walkers per subregion . Parameters are chosen as in run 1 (see Tab. 1).

ii.1.2 Stationary probability density function

We start to calculate the SPDF in the region for a small noise strength (). The time step and the box size are set according to the criteria (see Tab. 1).

1 -1.4 1.4 0.011 0.00104869 2670
2 -1.75 1.75 0.0015 0.000387297 9037
3 -2.5 2.5 0.0001 0.00010775 46404
Table 1: Time step and box size according to convergence criteria Eq. (10) and (8), where we choose and .

Using the results shown in Fig. 2, we set the thermalization time for a run with . Time averages after thermalization were calculated over an ensemble of sets P.

Figure 3: Estimates for the stationary probability density obtained from the algorithm for run (see Tab. 1) and by using a Brownian dynamics simulation for . Analytic results are obtained from Eq. (2).

Results for the SPDF are shown in Fig. 3. The algorithm calculates the tails of the distribution down to correctly, after a running time . We also calculate the SPDF using Brownian dynamics simulation using , which stops estimating at a level of after the same running time. Further runs were performed (see Tab. 1 (run ) and (run )), approximating the tails down to () and () after a running time of and , respectively (data not shown).

Simulations for different values of indicate that insignificant deviations from the analytic SPDF occur for .

ii.1.3 Probability current

Next, we present that our algorithm can be used to calculate the escape rate to pass the energy barrier at . Such problems are typical for chemical reactions Hänggi et al. (1990) and in the field of neuroscience Tuckwell (1988).

Initially, only particles are assigned at the subregion including , so approximating an initial delta-like probability distribution for . Furthermore, an absorbing boundary right behind the local maximum () is included. To fulfill normalization of the SPDF, walkers that reach are reinjected immediately at . The escape rate to pass the barrier is given by the probability current . For small noise intensities, the probability current on top of the potential barrier can be described using Arrhenius law, namely:


where . Since strong fluctuations are rare, but possible, we will use to approximate . Probability currents were recorded for each time step and averaged over a sequence of running steps, resulting in .

Figure 4: Time dependence of the probability current obtained from our algorithm for , and decreasing noise intensities (from top to bottom). The time step is set and . Note that and vary according to Eq. (10) and (8), respectively, resulting in larger running times for smaller .

Figure 4 shows the time dependence of . After a relaxation regime, where the current decays exponentially, the current reaches its stationary value.

Figure 5: Stationary probability current as a function of the inverse noise strength obtained by time averaging the data partly shown in Fig. 4 in the stationary regime. Error bars show three standard deviations of the stationary data.

The values of , averaged over the stationary regime, are shown in Fig. 5 for different noise intensities. Fulfilling the criteria described above, the algorithm reproduces well Arrhenius law down to corresponding to a current .

ii.1.4 Efficiency compared to weighted-ensemble Brownian dynamics simulation

In order to compare the efficiency of two algorithms, important quantities are the transient time required to first reach the steady state. Once the algorithm reaches the steady state, we quantify the size of the fluctuations by the coefficient of variation Eq. (14) at the potential’s local maximum . The computation time mainly depends on the number of integrations needed, to reach the stationary regime.

Figure 6: Relative fluctuations plotted over the number of integration steps, needed for equilibration, for the algorithm (black) and standard (WE) simulation (red). The number of walkers per subregion is shown for each point. Simulations were done for run and (compare Tab. 1).

In order to compare the efficiency, the size of fluctuations (Eq. (14)) in the local minimum relative to during the stationary regime is plotted over in Fig. 6. The most effective algorithm would be located close to the origin. Comparing the efficiency of standard WE simulations and our algorithm, we find that WE simulations with low equilibrate faster. However, the precision highly depends on the fluctuations during the stationary regime. To produce results of same precision (same ) both algorithms approximately need the same . Usually WE simulations were performed using thousands of walkers per subregions, resulting in a low value of . Here runs of our algorithm producing the same need much less and have memory requirements independent of .

ii.2 Two-dimensional systems

ii.2.1 Poincaré Oscillator

As an example of a two-dimensional system with known SPDF, we consider the Poincaré oscillator Ebeling and Engel-Herbert (1980); Ebeling et al. (1986); Lekkas et al. (1988), represented by the dynamical system:


where represent delta correlated white Gaussian noise with zero mean. Using the energy function , which only depends on the distance to the origin, one can calculate the associated SPDF:


(see appendix  A). Since noise only applies to the -direction, the lengths of a subregion and in - and -direction are calculated by Eq. (13) and (8), respectively. The minimum of is equal to , therefore, we choose , which corresponds to a first order approximation of in the next subregion.

Figure 7: Contour plots of the SPDF ( and ) obtained from the algorithm (black), using the parameters , , , , , , , , , () and the analytic solution (red) Eq. (17). Contour lines are labeled according to represented values of the SPDF and show the rotational symmetry. The SPDF possesses its global maximum at , corresponding to the unit circle for our choice of , and a local minimum in the origin. Running time .

Results for the SPDF are shown in Fig. 7. Interestingly, the algorithm approximates better the SPDF along the direction where no noise was applied. Here the SPDF is sampled down to . We found that the algorithm slightly oversamples the analytic SPDF in the tails. This is due to the statistical errors, made during the redistribution step. By reducing the size of a subregion, this error can be reduced further. Along the -direction, noise is applied. Here the behavior is similar as in the one dimensional example (see above). For runs with larger (results not shown) the algorithm slightly oversamples the SPDF in the origin.

ii.2.2 Bistable system with colored noise

As a further example, we calculate the SPDF of the two dimensional system:


where denotes the time scale separation between and the colored noise . The white Gaussian noise has been already described above. This system has been studied previously in Debnath et al. (1990); Hänggi and Jung (1995). Like in the bistable system we have studied above, the SPDF has maximums at and . However, for some combinations of and , the SPDF possesses a local minimum at .

Figure 8: Contour plots of the SPDF of Eq. (18) (, and ) obtained from the algorithm (black) and BDS (red) of the same running time. Parameters: (top) , , , , , , , , ,(, ), , and (bottom) , , , , , , , , ,(, ), . Running times are (top) and (bottom)

Plots of the SPDF are depicted in Fig. 8. Our algorithm samples the SPDF down to , whereas BDS breaks down at a level of . The minimum, occurring for was clearly found by the algorithm.

ii.2.3 FitzHugh-Nagumo-system

As a last example, we consider the widely used FitzHugh-Nagumo-system Fitzhugh (1961), which is often used in the field of Neuroscience Lindner et al. (2004) or to study synchronization Gunton et al. (2003); Toral et al. (2003) and coherence phenomena Toral et al. (2007), represented by:


Here denotes the timescale separation between the activator variable and the inhibitor variable . , represent independent zero-mean delta-correlated Gaussian white noises. We want to study the stationary probability density in the case of for a time scale separation . We set the parameters according to Ref. Kostur et al. (2003) to , and . Thus, the system is in the excitable regime. Since the deterministic part of the equation for the activator variable increases very fast if is increased, we have to choose a time step , which is small enough, that the walkers’ steps are small compared to , but allows us to fulfill the criteria for the size of the subregions.

Figure 9: Contour plots of the the SPDF of Eq. (19) fod ( (top) and (bottom)) obtained from the algorithm (black) and BDS (red) with the same running time. Parameters: (top) , , , , , , , ,(), , , and (bottom) , , , , , , , , (), , . Running times are (top) and (bottom). Simulation times for the BDS were chosen three times larger.

Contour plots of the SPDF are shown in Fig. 9. Especially regions of low probability are much better sampled, using the algorithm. In the case of low diffusion (, Fig. 9 (top)) the algorithm runs down to , whereas BDS stops at a level of . Especially the minimum is much better sampled by our algorithm. Note that the local maximum located in the surroundings of was not found by BDS. For higher diffusion values (, Fig. 9 (top)), significant differences can be found only in the tails.

Iii Discussion

iii.1 One dimensional system

In the one dimensional system (see section II.1) we succeeded to approximate the SPDF down to . If equally sized subregions are used, a huge number of subregions has to be implemented, in order to fulfill the convergence criteria (see section I.5). However, by evaluating the SPDF according to Eq. (6) the additional computational costs also reduce the fluctuations of the estimated SPDF. If the subregions are to large, the theoretical SPDF is overestimated by the algorithm in potential minimums and underestimated in the potentials maximums.

We also managed to calculate escape rates of size . Although the finite size of a subregion leads to small errors, in the calculated SPDF, Arrhenius law is well reproduced for such small probability currents.

At the stationary regime the fluctuate around their mean value, estimating the SPDF. By either increasing the number of averages or the number of walkers per subregion these fluctuations can be reduced, leading to higher precision. The estimation error scales with and , respectively. However, increasing the number of walkers per subregion also effects the computational costs during the thermalization. Simulations for different show that runs with higher become stationary faster, but this does not compensate for the additional computational costs. We also find, that increasing slightly improves the sampling of the SPDF’s tails. Once the system is in the stationary regime the increase of the computational costs scale linearly with and . An advantage of simulations with small is, that one does not need to perform running steps for subregions with and the number of such regions naturally increases for small . We usually use the smallest possible and scale the estimation error by increasing .

iii.2 Comparison with weighted-ensemble Brownian dynamics

Since the general idea of our method was adapted from prior simulation techniques known as weighted-ensemble (WE) sampling Huber and Kim (1996); Bhatt et al. (2010); Bhatt and Bahar (2012), we want to discuss advantages and disadvantages of our algorithm in comparison with these techniques. The main difference in WE techniques is the redistribution step. Here, using WE techniques, positions and weights of all walkers are stored, and new walkers are introduced on the stored positions considering their particular weights.

In our algorithm, there is no need to store any position or weight, since walkers are randomly placed in each subregion. The resulting statistical errors can be neglected, if the size of the subregions fulfills conditions which, unfortunately, lead to much larger numbers of subregions. By averaging the probabilities, these extra computational costs contribute to an reduction of fluctuations in the stationary regime.

The comparison of the computational costs until equilibration and the achieved precision shows, that both algorithms posses the same efficiency for high precision runs.

iii.3 Two dimensional systems

In the case of two dimensional systems, we found that the algorithm outperforms Brownian dynamics simulations. However, the number of subregions needed according to the criteria can be really high, especially if there are some directions without any noise. Here it is necessary to size the subregions to fulfill the criteria locally. A first step in that direction has already been done by the implementation of the grouping algorithm. We are quite confident that it is possible to reduce the computational costs by optimizing the Mesh. The analyzed examples show, that the running times highly depend on the investigated system. For the bistable system with colored noise and the FitzHugh-Nagumo-system, we obtained good results within minutes, whereas the algorithm needed about hours for the Van-der-Pol oscillator.

Iv Conclusion

We provided and tested an algorithm that allows the calculation of low probabilities and low rates. The algorithm is based on WE Brownian dynamic simulations, but uses a uniform distribution of walkers within each subregion. To our findings, the resulting statistical errors can be neglected if one uses subregions small compared to the diffusion length. In contrast to WE methods, the required memory does not depend on the number of walkers, which leads to less memory requirements for runs with large numbers of walkers.

Special attention was payed to non-equilibrium dynamical systems. Applying the method to one- and two-dimensional model systems, we analyze its efficiency compared to standard Brownian dynamics simulation. Our method outperforms Brownian dynamics simulation by several orders of magnitude and its efficiency is comparable to weighted-ensemble Brownian dynamic simulations in all studied systems and lead to impressive results in regions of low probability and small rates.

This paper was developed within the scope of the IRTG 1740 / TRP 2011/50151-0, funded by the DFG / FAPESP. LSG acknowledges IFISC of the University of Balearic Island for cordial hospitality and support. He also acknowledges Dr. Volkhard Bucholtz from Logos-Verlag Berlin for earlier participation and work on the project. RT acknowledges financial support from MINECO (Spain), Comunitat Autònoma de les Illes Balears, FEDER, and the European Commission under project FIS2007-60327, as well as the warm hospitality at Humboldt University.

Appendix A Stationary probability density of the Poincaré oscillator

Using the energy function , we get from  (16) to a representation as a canonic dissipative system:


The corresponding Fokker Planck equation in the phase space for the SPDF reads Ebeling and Engel-Herbert (1980):


Using the ansatz , the first two items at the r.h.s. cancel. The remaining second line can be integrated once. Assuming an exponentially decaying SPDF at infinitely large energies yield the disappearance of the irreversible probability flux in y-direction Graham and Haken (1971a); Haken (1975). One finds, afterwards :


It leads to Eq. (17).

Appendix B Foundation of the algorithm

In this Appendix we show that the presented algorithm is described by a corresponding Master equation for the probability distribution density for the case of a Markovian hopping process between boxes. We use a single index to label the boxes, but the argument applies to any spatial dimension . We identify the dynamics of the stochastic system which shall be simulated with the discrete stochastic dynamics of the walkers. The latter is defined via the matrices of probabilities per unit time which describe the hopping in the given discretized space. We assume that it shall converge for sufficiently small time scales and box lengths to the outgoing dynamics.

Let assume that we have the probability distribution density given at time in every box with index . It holds


We redistribute the probability in every box to walkers. In result any walker of the same box gets an identical weight


until the next new redistribution.

At later time the walkers may be still inside the box or may have jumped to other boxes. Let denote all possible box-indices which can be reached during a single step from the -box. Then in accordance with the assumption above, the probability is that during a single walker leaves -box and jumps to the box with index . By this hopping the walker carries the weight to the new box, which step is, obviously, connected with a lost in the outgoing box. Therefore, the lost per particle for this specific hopping from can be expressed as


The whole lost of weight will be realized on all possible hopping channels. It is identical for all particles being located in the present box. Hence, the full lost becomes


Alternatively, one can also introduce as the boxes where from walkers can reach the -box. Then, gain of weight transferred by every walker arriving at the -box is expressed, if , by


Again, summing over the different hopping steps and considering that all the walkers inside the box with reach the box, yields the gain


Therefore, the balance of transported weight results in the following shift of the full weight in the -box at time compared to the former one


plus corrections of order corresponding to cases in which two or more walkers jump from one box to another during the time interval . After dividing by the time step and taking the limit we obtain the wanted Master equation for the discretized dynamics in the boxed phase space



  • Moss and McClintock (1989) F. Moss and P. V. McClintock, Noise in nonlinear dynamical systems, Vol. I-III (Cambridge University Press, New York, 1989).
  • G. van Kampen (1992) N. G. van Kampen, Stochastic processes in physics and chemistry (North-Holland, Amsterdam, 1992).
  • Nicolis and Prigogine (1977) G. Nicolis and I. Prigogine, Self-Organization in Nonequilibrium Systems: From Dissipative Structures to Order Through Fluctuations (Wiley-Interscience, New York, 1977).
  • Mannella and Palleschi (1989) R. Mannella and V. Palleschi, Phys. Rev. A 40, 3381 (1989).
  • Mannella (2000) R. Mannella, in Stochastic Processes in Physics, Chemistry, and Biology, Lecture Notes in Physics, Vol. 557, edited by J. A. Freund and T. Pöschel (Springer, Berlin, 2000) p. 353.
  • Burrage (1999) P. Burrage, Runge-Kutta methods for stochastic differential equations, Ph.D. thesis, University of Queensland (1999).
  • Kloeden and Platen (2011) P. Kloeden and E. Platen, Numerical solution of stochastic differential equations, Vol. 23 (Springer, Berlin, 2011).
  • Bhanot et al. (1987) G. Bhanot, R. Salvador, S. Black, P. Carter,  and R. Toral, Phys. Rev. Lett. 59, 803 (1987).
  • Giardina et al. (2006) C. Giardina, J. Kurchan,  and L. Peliti, Phys. Rev. Lett. 96, 120603 (2006).
  • Dickson and Dinner (2010) A. Dickson and A. Dinner, Annu. Rev. Phys. Chem. 61, 441 (2010).
  • Wang and P. Landau (2001) F. Wang and D. P. Landau, Phys. Rev. Lett. 86, 2050 (2001).
  • Torrie and Valleau (1977) G. Torrie and J. Valleau, J. Comp. Phys. 23, 187 (1977).
  • Warmflash et al. (2007) A. Warmflash, P. Bhimalapuram,  and A. Dinner, J. Chem. Phys. 127, 154112 (2007).
  • Allen et al. (2009) R. Allen, C. Valeriani,  and P. ten Wolde, J. Phys.: Condens. Matter 21, 463102 (2009).
  • Dickson et al. (2009) A. Dickson, A. Warmflash,  and A. Dinner, J. Chem. Phys. 131, 154104 (2009).
  • Valeriani et al. (2007) C. Valeriani, R. Allen, M. Morelli, D. Frenkel,  and P. Wolde, J. Chem. Phys. 127, 114109 (2007).
  • Huber and Kim (1996) G. Huber and S. Kim, Biophys. J. 70, 97 (1996).
  • Bhatt et al. (2010) D. Bhatt, B. Zhang,  and D. Zuckerman, J. Chem. Phys. 133, 014110 (2010).
  • Zhang et al. (2010) B. Zhang, D. Jasnow,  and D. Zuckerman, J. Chem. Phys. 132, 054107 (2010).
  • Bhatt and Bahar (2012) D. Bhatt and I. Bahar, J. Chem. Phys. 137, 104101 (2012).
  • Bhatt and Zuckerman (2011) D. Bhatt and D. Zuckerman, J. Chem. Theory Comput. 7, 2520 (2011).
  • Graham and Haken (1971a) R. Graham and H. Haken, Z. Phys. 243, 289 (1971a).
  • Graham and Haken (1971b) R. Graham and H. Haken, Z. Phys. 245, 141 (1971b).
  • Haken (1975) H. Haken, Rev. Mod. Phys. 47, 67 (1975).
  • Risken (1984) H. Risken, The Fokker Planck equation, Vol. 23 (Springer, Berlin, 1984).
  • Schimansky-Geier et al. (1985) L. Schimansky-Geier, A. Tolstopjatenko,  and W. Ebeling, Phys. Lett. 108A, 329 (1985).
  • Ebeling and Engel-Herbert (1980) W. Ebeling and H. Engel-Herbert, Physica A 104, 378 (1980).
  • Klimontovich (1980) Y. Klimontovich, Kinetic theory of elecromagnetic processes (Springer, Berlin, 1983 (in Russian: Nauka, Moscow, 1980)).
  • Gilsing and Shardlow (2007) H. Gilsing and T. Shardlow, J. Comput. Appl. Math. 205, 1002 (2007).
  • Gammaitoni et al. (1998) L. Gammaitoni, P. Hänggi, P. Jung,  and F. Marchesoni, Rev. Mod. Phys. 70, 223 (1998).
  • Hänggi et al. (1990) P. Hänggi, T. P.,  and N. Borkovec, Rev. Mod. Phys. 62, 251 (1990).
  • Tuckwell (1988) H. Tuckwell, Introduction to theoretical neurobiology, Vol. 2 (Cambridge University Press, New York, 1988).
  • Ebeling et al. (1986) W. Ebeling, H. Herzel, W. Richert,  and L. Schimansky-Geier, Zeitschrift f. angew. Math. und Mech. 66, 141 (1986).
  • Lekkas et al. (1988) K. Lekkas, L. Schimansky-Geier,  and H. Engel-Herbert, Z. Phys. B - Condensed Matter 70, 517 (1988).
  • Debnath et al. (1990) G. Debnath, F. Moss, T. Leiber, H. Risken,  and F. Marchesoni, Phys. Rev. A 42, 703 (1990).
  • Hänggi and Jung (1995) P. Hänggi and P. Jung, Adv. Chem. Phys. 89, 239 (1995).
  • Fitzhugh (1961) R. Fitzhugh, Biophys. J. 1, 445 (1961).
  • Lindner et al. (2004) B. Lindner, J. Garcıa-Ojalvo, A. Neiman,  and L. Schimansky-Geier, Phys. Rept. 392, 321 (2004).
  • Gunton et al. (2003) J. Gunton, R. Toral, C. Mirasso, M. Gracheva, et al., Recent. Res. Devel. Applied Phys. 6, 497 (2003).
  • Toral et al. (2003) R. Toral, C. Masoller, C. Mirasso, M. Ciszak,  and O. Calvo, Physica A 325, 192 (2003).
  • Toral et al. (2007) R. Toral, C. Mirasso,  and J. Gunton, Europhys. Lett. 61, 162 (2007).
  • Kostur et al. (2003) M. Kostur, X. Sailer,  and L. Schimansky-Geier, Fluct. Noise Lett. 3, 155 (2003).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description