Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
branch: master
Fetching contributors…

Cannot retrieve contributors at this time

212 lines (188 sloc) 13.976 kb
<!DOCTYPE html>
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8"/>
<title>GFRD algorithm</title>
<link rel="stylesheet" href="/css/reset.css" type="text/css"/>
<link rel="stylesheet" href="/css/style.css" type="text/css"/>
<link rel="stylesheet" href="/css/syntax.css" type="text/css"/>
<link rel="stylesheet" href="/css/menu.css" type="text/css"/>
<link rel="stylesheet" href="/css/footer.css" type="text/css"/>
<script src="/html5media/html5media.min.js"></script>
<!--[if lt IE 9]>
<script src="/ie7/IE9.js"></script>
<![endif]-->
</head>
<body>
<div class="algorithm">
<div id="globalheader">
<!--googleoff: all-->
<ul id="globalnav">
<li id="gn-home"><a href="/">Home</a></li>
<li id="gn-algorithm"><a href="/algorithm">Algorithm</a></li>
<li id="gn-background"><a href="/background">Background</a></li>
<li id="gn-downloads"><a href="/downloads">Downloads</a></li>
<li id="gn-documentation"><a href="/documentation">Documentation</a></li>
<li id="gn-support"><a href="/support">Support</a></li>
<li class="gn-empty"><a></a></li>
</ul>
<!--googleon: all-->
<div id="lastbutton"> </div>
</div>
</div>
<h1>Algorithm</h1>
<div id="container">
<div class="grid2col"><blockquote><h2>Green’s Function Reaction Dynamics</h2></blockquote>
<div class="column first"><p>A reaction-diffusion system is a many-body problem that cannot be solved
analytically. The key idea of GFRD is to decompose the many-body problem into
one- and two-body problems that can be solved analytically via Green’s
functions, and use these Green’s functions to set up an event-driven
algorithm<sup><a href="#c1">1</a>,<a href="#c2">2</a></sup>. The Green’s functions allow GFRD to make
large jumps in time and space when the particles are far apart from each
other.</p>
<p>In the original version of the algorithm, the many-body problem was solved
by determining at each step of the simulation a maximum time step such that
each particle could interact with at most one other particle during that time
step<sup><a href="#c1">1</a>,<a href="#c2">2</a></sup>. This scheme was a synchronous event-driven
algorithm, because at each time step all the particles were propagated
simultaneously. Moreover, the scheme was not exact, because the decomposition
into single particles and pairs of particles involved cut-off distances,
introducing a trade-off between speed and error.</p>
<p>In the new version of the algorithm, we have implemented ideas of Opplestrup
and coworkers<sup><a href="#c3">3</a></sup>. In this new scheme, protective domains are put
around single particles and pairs of particles, as shown in Fig. 1. For each
of the domains, the reaction-diffusion problem is solved analytically using
Green’s functions. This yields for each domain an <em>event type</em> and an <em>event
time</em>; as described below, the set of possible <em>event types</em> depends on
whether the domain is a <em>Single</em>, meaning it contains a single particle, or a
<em>Pair</em>, meaning it contains a pair of particles. The <em>event times</em> are put in
a chronologically ordered <em>event list</em>, and the events are then executed in
chronological order. When an event is executed, the particles of the
corresponding domain are propagated, for the propagated particles new domains,
with new <em>event types</em> and new <em>event times</em>, are determined, and the new
events are put back in the <em>event list</em>. The new version of GFRD, called
eGFRD, is thus an exact, event-driven, asynchronous algorithm. Its
asynchronous nature makes eGFRD similar in spirit to event-driven MD
simulations of hard spheres and the Gibson-Bruck scheme, which is an exact,
event-driven asynchronous algorithm for simulating the <em>zero-dimensional</em>
chemical master equation<sup><a href="#c4">4</a></sup>. A movie of eGFRD in action is shown in
Movie 1.</p>
<h3><em>Single</em></h3>
<p>A <em>Single</em> is a domain that contains only a single particle. The next <em>event
type</em> is either a unimolecular reaction of the type A ⇾ B + C + … or A ⇾ Ø, or
an escape from the domain, as sketched in Fig. 2. To determine which of the
two possible events will occur, a tentative <em>event time</em> is determined for
each of them; the event with the smallest tentative <em>event time</em> is the one
that will occur for that domain; it is the one that is put in the <em>event
list</em>.</p>
<h3><em>Pair</em></h3>
<p>A <em>Pair</em> is a domain that contains a pair of particles. The reaction-diffusion
problem of two particles that can react with each other and diffuse in a
spherical domain with absorbing boundary conditions can, to our knowledge, not
be solved analytically. We therefore apply a coordinate transformation leading
to a diffusion problem for the center-of-mass coordinate <strong>R</strong> and a
reaction-diffusion problem for the inter-particle vector <strong>r</strong>, as sketched in
Fig. 3. The <em>event types</em> that can occur are <strong>R</strong> leaving its domain, <strong>r</strong>
leaving its domain, a unimolecular reaction or a bimolecular reaction. For
each of these events a tentative <em>event time</em> is determined, and the event
with the smallest tentative <em>event time</em> is the one that will happen, and put
in the <em>event list</em>.</p>
<h2>Performance</h2>
<p>Fig.4 shows the power of eGFRD. Plotted is the CPU time for simulating a
system consisting of hard spheres for a fixed amount of real time as a
function of the number of particles N, for two scenarios: one in which the
concentration C is kept constant (solid line) and one in which the volume V is
kept constant (dashed lines). For comparison, two curves (in blue) obtained
with Brownian Dynamics (BD) are also shown: one corresponding to a BD
simulation with a conventional, small time step of δt = 10<sup>-6</sup> σ<sup>2</sup> / D, and one
corresponding to a (relaxed) BD simulation using a relatively large time step
of δt = 10<sup>-3</sup> σ<sup>2</sup> / D. It is seen that the CPU time of eGFRD scales linearly
with N when the concentration is constant, as in BD; this is because a cell
list is used for determining which particles or domains interact with a given
particle or domain. However, the CPU time of eGFRD scales with N<sup>5/3</sup> when the
volume is kept constant. This can be understood by noting that the CPU time
scales with N / &lt;Δt>, where &lt;Δt> is the magnitude of the average time step in
the eGFRD simulations; the size of the average time step scales as &lt;Δt> ~ 1 /
l<sup>2</sup> ~ N<sup>2/3</sup>, where l is the average distance between the particles. Indeed,
when N is increased while V is kept constant, the concentration C increases
and the average distance between the particles decreases. This decreases the
jumps in time and space that can be made during the eGFRD simulations. Since
the CPU time of BD scales with N wile that of eGFRD scales with N<sup>5/3</sup> when V
is constant, there is a concentration above which BD becomes more efficient
than eGFRD. This is demonstrated in the inset of Fig.4. It is seen that the
cross-over occurs at mM concentration. However, the concentrations of
components in signal transduction pathways are typically in the µM range,
while those in gene regulation networks are typically in the nM range. Under
these biologically relevant conditions, eGFRD is up to 6 orders of magnitude
more efficient than BD.</p></div>
<div class="column last"><div class="media"><p><embed width="441px" height="281px" src="/images/gfrd.svg"
title="Fig.1. GFRD" type="image/svg+xml" pluginspage="http://www.adobe.com/svg/viewer/install/"/><br/>
Fig.1: <em>Overview of GFRD. The principal idea of GFRD is to decompose the
many-body reaction-diffusion problem into one- and two-body problems that can
be solved analytically using Green’s Functions. To this end, protective
domains are put around single particles and pairs of particles, leading to
so-called </em>Singles<em> and </em>Pairs<em>, respectively. For each domain, a next </em>event
type<em> and a next </em>event time<em> is determined. The </em>event times<em> are put in a
chronologically ordered </em>event list<em>, and the events are then executed in
chronological order. When an event is executed, the particles of the
corresponding domain are propagated, new domains with new events are
determined, and the events are put back in the event list.</em></p></div>
<div class="media"><p><video class="video" width="440px" height="260px" controls preload="none" title="Movie 1. GFRD in 2D">
<source src="/movies/planar-surface_6-particles_200-steps/planar-surface_6-particles_200-steps.ogv" type='video/ogg; codecs="theora, vorbis"'/>
<source src="/movies/planar-surface_6-particles_200-steps/planar-surface_6-particles_200-steps.mp4" type='video/mp4; codecs="avc1.42E01E, mp4a.40.2"'/>
</video><br/>
Movie 1 (<a href="/movies/planar-surface_6-particles_200-steps/planar-surface_6-particles_200-steps.ogv">ogv</a>,
<a href="/movies/planar-surface_6-particles_200-steps/planar-surface_6-particles_200-steps.mp4">mp4</a>).
<em>A movie of eGFRD in action in 2D. The cylinders are protective domains for single particles
(light blue) or pairs of particles (dark blue). With the color red the
protective domain that will be updated next is highlighted.</em></p></div>
<div class="media"><p><embed width="441px" height="155px" src="/images/single.svg"
title="Fig.2. Single" type="image/svg+xml" pluginspage="http://www.adobe.com/svg/viewer/install/"/><br/>
Fig.2: <em>A </em>Single<em>, a protective domain that contains only a single particle.
The possible </em>event types<em> are either a unimolecular reaction or the escape of
the single particle from the protective domain.</em></p></div>
<div class="media"><p><embed width="441px" height="155px" src="/images/pair.svg"
title="Fig.3. Pair" type="image/svg+xml" pluginspage="http://www.adobe.com/svg/viewer/install/"/><br/>
Fig.3: <em>A </em>Pair<em>, a protective domain that contains a pair of particles. A
coordinate transformation is applied, and one protective domain is defined for
the center-of-mass coordinate <strong>R</strong> and one for the inter-particle vector <em>r</em>.
The sizes of these domains are determined such that when both <strong>R</strong> and <strong>r</strong>
would reach their maximum values, </em>|<em><strong>R</strong></em>|<em><sup><strong>max</strong></sup> and </em>|<em><strong>r</strong></em>|<em><sup><strong>max</strong></sup>,
respectively, the two particles would still be within the original protective
domain. The dynamics of <strong>R</strong> is a diffusion problem of a random walker in a
spherical domain with absorbing boundary conditions. The dynamics of <strong>r</strong>
takes into account that the two particles not only diffuse, but can also react
with each other; this is the problem of a random walker between two concentric
spherical surfaces, with a radiation boundary condition at the inner surface,
describing the bimolecular reaction, and an adsorbing boundary condition at
the outer surface, describing the escape from the domain. The next event is
thus either the escape of <strong>R</strong> from its domain, the escape of <strong>r</strong> from its
domain, a bimolecular reaction, or a unimolecular reaction.</em></p></div></div></div>
<p><img src="/images/performance.png" title="Fig.4. Performance"><br/>
Fig.4. <em>CPU time for simulating a system of hard spheres for a fixed amount of
real time as a function of the number of particles N for two scenarios: one in
which the concentration C is kept constant (solid line), and one in which the
volume V is kept constant (dashed line). For comparison, we also show the
results for BD simulations (blue curves), with one curve corresponding to a
simulation with a conventional small time step of δt = 10<sup>-6</sup> σ<sup>2</sup> / D, and one
corresponding to a (relaxed) BD simulation using a relatively large time step
of δt = 10<sup>-3</sup> σ<sup>2</sup> / D. The inset shows the CPU time as a function of the
concentration C, for N=300 and N=3000 (the volume is thus varied). It is seen
that above mM concentrations, BD is more efficient than eGFRD. However, in the
biologically relevant concentration range of nM to µM eGFRD can be up to 6
orders of magnitude more efficient than conventional BD.</em></p>
<div class="references"><h3>References</h3>
<p><a href="#r1"><span id="c1">1</span></a> van Zon JS, tenWolde PR (2005) Simulating biochemical networks at the particle level and in time and space: Green’s function reaction dynamics. <em>Phys Rev Lett</em>, 94:128103. (<a href="http://dx.doi.org/10.1103/PhysRevLett.94.128103">doi</a>)<br/>
<a href="#r2"><span id="c2">2</span></a> van Zon JS, ten Wolde PR (2005) Green’s-function reaction dynamics: A particle-based approach for simulating biochemical networks in time and space. <em>J Chem Phys</em>, 123:234910. (<a href="http://dx.doi.org/10.1063/1.2137716">doi</a>, <a href="http://arxiv.org/abs/q-bio/0404002">arXiv</a>)<br/>
<a href="#r3"><span id="c3">3</span></a> Opplestrup T, Bulatov VV, Gilmer GH, Kalos MH, Sadigh B (2006) First-passage Monte Carlo algorithm: diffusion without all the hops. <em>Phys Rev Lett</em>, 97:230602. (<a href="http://dx.doi.org/10.1103/PhysRevLett.97.230602">doi</a>, <a href="http://arxiv.org/abs/0905.3576">arXiv</a>)<br/>
<a href="#r4"><span id="c4">4</span></a> Gibson MA, Bruck J (2000) Efficient exact stochastic simulation of chemical systems with many species and many channels. <em>J Phys Chem A</em>, 104: 1876 — 1888. (<a href="http://dx.doi.org/10.1021/jp993732q">doi</a>)</p></div>
</div>
<div id="globalfooter">
<p>Copyright &copy; 2010 FOM institute AMOLF</p>
<ul class="piped">
<li><a href="/about">About</a></li>
<li><a href="/credits">Credits</a></li>
</ul>
</div>
</body>
</html>
Jump to Line
Something went wrong with that request. Please try again.