In [1]:
%matplotlib inline
import numpy
from ecell4 import *
from ecell4.extra.ensemble import ensemble_simulations
from ecell4_base.core import GSLRandomNumberGenerator, Integer3
Parameters are given as follows. D, radius, N_A, U, and ka_factor mean a diffusion constant, a radius of molecules, an initial number of molecules of A and B, a ratio of dissociated form of A at the steady state, and a ratio between an intrinsic association rate and collision rate defined as ka andkD below, respectively. Dimensions of length and time are assumed to be micro-meter and second.
In [2]:
D = 1
radius = 0.005
N_A = 60
U = 0.5
ka_factor = 0.1 # 0.1 is for reaction-limited
In [3]:
N = 20 # a number of samples
rng = GSLRandomNumberGenerator()
rng.seed(0)
Calculating optimal reaction rates. ka and kd are intrinsic, kon and koff are effective reaction rates.
In [4]:
kD = 4 * numpy.pi * (radius * 2) * (D * 2)
ka = kD * ka_factor
kd = ka * N_A * U * U / (1 - U)
kon = ka * kD / (ka + kD)
koff = kd * kon / ka
Start with no C molecules, and simulate 3 seconds.
In [5]:
y0 = {'A': N_A, 'B': N_A}
duration = 3
T = numpy.linspace(0, duration, 21)
Make a model with effective rates. This model is for macroscopic simulation algorithms.
In [6]:
with species_attributes():
A | B | C | {'radius': radius, 'D': D}
with reaction_rules():
A + B == C | (kon, koff)
m = get_model()
Save a result with ode as obs, and plot it:
In [7]:
obs = run_simulation(numpy.linspace(0, duration, 101), y0, model=m, return_type='observer',
solver='ode')
viz.plot_number_observer(obs)
Simulating with gillespie (Bars represent standard error of the mean):
In [8]:
ensemble_simulations(T, y0, model=m, opt_args=('o', obs, '-'),
solver='gillespie', n=N)
Simulating with meso:
In [9]:
ensemble_simulations(T, y0, model=m, opt_args=('o', obs, '-'),
solver=('meso', Integer3(1, 1, 1), 0.25), n=N)
Make a model with intrinsic rates. This model is for microscopic (particle) simulation algorithms.
In [10]:
with species_attributes():
A | B | C | {'radius': radius, 'D': D}
with reaction_rules():
A + B == C | (ka, kd)
m = get_model()
Simulating with spatiocyte. voxel_radius is given as radius:
In [11]:
ensemble_simulations(T, y0, model=m, opt_args=('o', obs, '-'),
solver=('spatiocyte', radius), n=N)
Simulating with egfrd:
In [12]:
ensemble_simulations(T, y0, model=m, opt_args=('o', obs, '-'),
solver=('egfrd', Integer3(4, 4, 4)), n=N)