Skip to content

Commit

Permalink
General update
Browse files Browse the repository at this point in the history
- Headers
- Introduction
- Equilibrium equation.
  • Loading branch information
Luc Blom authored and Luc Blom committed Aug 28, 2017
1 parent 7a13e44 commit ff63ba4
Show file tree
Hide file tree
Showing 10 changed files with 26 additions and 20 deletions.
1 change: 1 addition & 0 deletions _config.yml
Expand Up @@ -18,5 +18,6 @@ github_download_zip: "http://github.com/gfrd/modern_egfrd/zipball/master"
github_download_tar: "http://github.com/gfrd/modern_egfrd/tarball/master"
github_install_instructions: "https://github.com/gfrd/modern_egfrd/blob/master/INSTALL.md"
github_issue_tracker: "https://github.com/gfrd/modern_egfrd/issues"
github_old_repository: "https://github.com/gfrd/egfrd.git"

visualstudio_website: "https://www.visualstudio.com"
6 changes: 3 additions & 3 deletions algorithm.md
Expand Up @@ -3,7 +3,7 @@ title: background
layout: default
---

# Green’s Function Reaction Dynamics
## Green’s Function Reaction Dynamics
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 \[[1,2](#references)\]. The Green’s functions allow GFRD to make large jumps in time and space when the particles are far apart from each other.

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 \[[1,2](#references)\]. Although already up to five orders more efficient than conventional Brownian Dynamics \[[1](#references)\] and also very accurate by its own right, the original GFRD algorithm has three drawbacks: (1) Because of the synchronous nature, the decomposition into one and two-body problems has to happen at every simulation step; (2) All components in the system are propagated according to the smallest tentative reaction time, making the performance sub-optimal; (3) the decomposition into single particles and pairs of particles involves cut-off distances, which makes the algorithm inexact; the systematic error is controlled by a parameter that determines the probability that during a time step a particle travels a distance further than the maximum distance set by the requirement that each particle can interact with at most one other particle; this means that there is a trade-off between performance and error.
Expand Down Expand Up @@ -41,7 +41,7 @@ A Pair is a domain that contains a pair of particles. The reaction-diffusion pro
Fig.3: A Pair, 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 <b>R</b> and one for the inter-particle vector <b>r</b>. The sizes of these domains are determined such that when both <b>R</b> and <b>r</b> would reach their maximum values, |<b>R</b>|<sup>max</sup> and |<b>r</b>|<sup>max</sup>, respectively, the two particles would still be within the original protective domain. The dynamics of <b>R</b> is a diffusion problem of a random walker in a spherical domain with absorbing boundary conditions. The dynamics of <b>r</b> 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 <b>R</b> from its domain, the escape of <b>r</b> from its domain, a bimolecular reaction, or a unimolecular reaction.
</div>

# Performance
## Performance
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 <br/>
δ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/<Δt>, where <Δt> is the magnitude of the average time step in the eGFRD simulations; the size of the average time step scales as <Δ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.

Expand All @@ -51,7 +51,7 @@ Fig.4. CPU time for simulating a system of hard spheres for a fixed amount of re
δ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.
</div>

# References
## References
1. Van Zon JS, Ten Wolde PR (2005) Simulating biochemical networks at the particle level in time and space: Green’s Function Reaction Dynamics. Phys Rev Lett, 94: 128103. ([doi](https://dx.doi.org/10.1103/PhysRevLett.94.128103))
2. Van Zon JS, Ten Wolde PR (2005) Green’s Function Reaction Dynamics: A particle-based approach for simulating biochemical networks in time and space. J Chem Phys, 123: 234910. ([doi](https://dx.doi.org/10.1063/1.2137716), [arXiv](https://arxiv.org/abs/q-bio/0404002))
3. Opplestrup T, Bulatov VV, Gilmer GH, Kalos MH, Sadigh B (2006) First-passage Monte Carlo algorithm: diffusion without all the hops. Phys Rev Lett, 97:230602. ([doi](https://dx.doi.org/10.1103/PhysRevLett.97.230602), [arXiv](https://arxiv.org/abs/0905.3576))
Expand Down
Binary file added includes/images/arrow.jpg
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified includes/images/movie.jpg
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified includes/images/movie_arrow.jpg
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified includes/movies/movie.mp4
Binary file not shown.
Binary file removed includes/packages/WinSupport.zip
Binary file not shown.
10 changes: 7 additions & 3 deletions index.md
Expand Up @@ -11,14 +11,18 @@ GFRD decomposes the many-body reaction-diffusion problem into one- and two-body
<b>Movie 1. eGFRD in action.</b>
</p>

# Applications
The eGFRD algorithm is generic and can be applied to a wide variety of reaction-diffusion problems, including those in population dynamics, evolution, and soft-condensed matter physics. The scheme presented here has been specifically designed to simulate biochemical networks. Recently, we have completely rewritten the code in Modern C++, resulting a very fast simulator.
## Applications
The eGFRD algorithm is generic and can be applied to a wide variety of reaction-diffusion problems, including those in population dynamics, evolution, and soft-condensed matter physics. The scheme presented here has been specifically designed to simulate biochemical networks.

## eGFRD in all dimensions
The original code to simulate reactions and diffusion in 3D (cytoplasm) has been rewritten in Modern C++, resulting in a very fast simulator. A [prototype]({{site.github_old_repository}}) has been developed that can also simulate systems in 2D (membranes) and 1D (filaments) \[[5](#references)\].

## Developers
The eGFRD algorithm was originally developed by the group of Takahashi at the [Riken institute]({{site.riken_website}}) in Japan and the group of Ten Wolde at [AMOLF]({{site.company_website}}) in The Netherlands.

# References
## References
1. Van Zon JS, Ten Wolde PR (2005) Simulating biochemical networks at the particle level in time and space: Green’s Function Reaction Dynamics. Phys Rev Lett, 94: 128103. ([doi](https://dx.doi.org/10.1103/PhysRevLett.94.128103))
2. Van Zon JS, Ten Wolde PR (2005) Green’s Function Reaction Dynamics: A particle-based approach for simulating biochemical networks in time and space. J Chem Phys, 123: 234910. ([doi](https://dx.doi.org/10.1063/1.2137716), [arXiv](https://arxiv.org/abs/q-bio/0404002))
3. Takahashi K, Tanase-Nicola S, Ten Wolde PR (2010) Spatio-temporal correlations can drastically change the response of a MAPK pathway. Proc. Natl Acad Sci USA, 107: 2473 — 2478. ([doi](https://dx.doi.org/10.1073/pnas.0906885107), [arXiv](https://arxiv.org/abs/0907.0514))
4. Opplestrup T, Bulatov VV, Gilmer GH, Kalos MH, Sadigh B (2006) First-passage Monte Carlo algorithm: diffusion without all the hops. Phys Rev Lett, 97:230602. ([doi](https://dx.doi.org/10.1103/PhysRevLett.97.230602), [arXiv](https://arxiv.org/abs/0905.3576))
5. Sokolowski TR, Ten Wolde PR (2017) Spatial-Stochastic Simulation of Reaction-Diffusion Systems. ([arXiv](https://arxiv.org/abs/1705.08669))
23 changes: 12 additions & 11 deletions installation.md
Expand Up @@ -3,27 +3,28 @@ title: installation
layout: default
---

The modern C++ implementation of the eGFRD has been published as an open source package on GitHub. Please follow the [installation instructions]({{site.github_install_instructions}}) on the repository to setup the simulator. Problems can be reported via the eGFRD [issue tracker]({{site.github_issue_tracker}}).

The modern C++ implementation of the eGFRD has been published as open source package on GitHub. Please follow the [installation instructions]({{site.github_install_instructions}}) on the repository to setup the simulator. Problems can be reported on the eGFRD [issue tracker]({{site.github_issue_tracker}}).

# Examples
## Examples
To help you getting started, we have assembled some simple models. These can be used to verify the implementation of the code and to build more complex models.

## Equilibrium
Consider a set of particles A, B and C with the following reaction: A + B ⇾ C with binding rate constant k<sub>a</sub> [m<sup>3</sup>/s] and C ⇾ A + B and unbinding rate k<sub>d</sub> [s<sup>-1</sup>]. In a box with length L [m] and periodic boundary conditions, the reaction will relax towards equilibrium. For this system, we can calculate the number of C particles analytically.
```
<#C> = (NA + NB) + KD*V*sqrt[(NA + NB + KD*V)^2 - 4*NA*NB]/2
```

where NA, NB and NC are the number of A, B, and C respectively in the volume V = L<sup>3</sup>, and KD = k<sub>d</sub>/k<sub>a</sub> is the dissociation constant. The number of C particles thus depends on the dissociation constant, that is the ratio of the association and dissociation rate, and not on the absolute values of the rate constants. This can be used to test the code. To run an eGFRD simulation of this reaction, enter the following command:
<code><#C> = (N<sub>A</sub> + N<sub>B</sub>) + K<sub>D</sub> &bull; V &bull; sqrt[(N<sub>A</sub> + N<sub>B</sub> + K<sub>D</sub> &bull; V)<sup>2</sup> - 4 &bull; N<sub>A</sub> &bull; N<sub>B</sub>]/2</code>

where N<sub>A</sub>, N<sub>B</sub> and N<sub>C</sub> are the number of A, B, and C respectively in the volume V = L<sup>3</sup>, and K<sub>D</sub> = k<sub>d</sub>/k<sub>a</sub> is the dissociation constant. The number of C particles thus depends on the dissociation constant, that is the ratio of the association and dissociation rate, and not on the absolute values of the rate constants. This can be used to test the code. To run an eGFRD simulation of this reaction, enter the following command in the Unix shell after installing the package:
```
./RunGfrd Equilibrium -ka 1e-19 -kd 2e-2 -p 100 -e 200 > data.out
```

The first command line argument sets the simulation type. The value “Equilibrium” represents a build-in simulation of the described association-dissociation reaction. The values of the association rate k<sub>a</sub> and dissociation rate k<sub>d</sub> can also be set as command line argument. The argument `-p` sets the preperation time. This allows the system to relax to an equilibrium. The argument `–e` sets the simulation time during which the statistics are accumulated. To view all arguments for the Equilibrium simulation type --help:
The variable “Equilibrium” describes the simple association-dissociation reaction. The values of the association rate ka and dissociation rate kd are set as command line argument. The parameter -p sets the equilibration time during which the system is allowed to relax to equilibrium, while –e sets the run time during which the statistics is accumulated, both in seconds. All parameter settings can be viewed by running:
```
./RunGfrd Equilibrium --help
```
The result will be stored in the file data.out and shows the simulation time against the instantaneous number of A, B and C particles. Finally, the analytical prediction <#C> is also stored in the file.
The result will be stored in the file data.out and shows the simulation time against the total number of A, B and C molecules. Finally, the analytical prediction <#C> is also stored in the file.

## Power Spectrum
Besides the dissociation constant, it may also prove useful to compute the power spectrum of the simple association-dissociation reaction. The following command prompts GFRD to compute the power spectrum for this reaction, with the same parameters as in the paper of Kaizu et. al. \[[1](#references)\]:
Expand All @@ -45,21 +46,21 @@ The Mitogen-Activated Protein Kinase (MAPK) cascade is one of most studied and b
## Custom
Users can also use the package to set up their own model with ini-file like scripts. A simple example with file name *rebind.gfrd*, is given below. To run the modern eGFRD simulator, copy the script to the same folder as where the RunGfrd executable is located and enter the following command:
```
.\RunGfrd Custom -f rebind.gfrd > data.out
.\RunGfrd Costum -f rebind.gfrd > data.out
```

```
; eGFRD rebinding simulation
[Simulator]
Seed = 0xA2529C06; pseudorandom number generator.
Seed = 0xA2529C06
PrepareTime = 0.5; sec.
EndTime = 15.0; sec.
ShowProgress = False
MaintenanceStep = 10000
MaintenanceFile = sim_dump.out
[World]
[World] ;this is some comment
Matrix = 8
Size = 1e-6; m
Expand Down Expand Up @@ -97,6 +98,6 @@ B = 50
C = 0
```

# References
## References
1. Kaizu K, de Ronde W, Paijmans J, Takahashi K, Tostevin F, ten Wolde PR (2014) The Berg-Purcell Limit Revisited Biophys J, 106:976-985. ([doi](https://dx.doi.org/10.1016/j.bpj.2013.12.030))
2. Takahashia K, Tănase-Nicolad S, ten Wolde PR (2009) Spatio-temporal correlations can drastically change the response of a MAPK pathway ([doi](https://dx.doi.org/10.1073/pnas.0906885107))

0 comments on commit ff63ba4

Please sign in to comment.