This notebook is part of PyBroMo a python-based single-molecule Brownian motion diffusion simulator that simulates confocal smFRET experiments. You can find the full list of notebooks in Usage Examples.
In this notebook we test that the diffusion coefficient of the simulated data matches with the value used as input. See also [Theory - On Browniam motion and Diffusion cofficient](Theory - On Browniam motion and Diffusion cofficient.ipynb)
Here we will see a simple example usage of the PyBroMo software.
The main file is brownian.py
: it contains the simulator code,
helper functions and the default simulation parameters.
For convenience, here we use an initialization script that will enter into the PyBroMo folder,
load brownian.py
and also import some useful IPython functions:
%run -i load_bromo.py
C:\Data\Antonio\software\src\pybromo
The default simulation parameters are defined in brownian.py
.
Let take a look at them:
# Simulation time step
t_step = 0.5e-6 # seconds
# Diffusion coefficient
Du = 12.0 # um^2 / s
D = Du*(1e-6)**2
# Time duration of the simulation
t_max = 0.1 # seconds
N_samples = int(t_max/t_step)
# PSF definition
psf = NumericPSF()
# Box definition
box = Box(x1=-4.e-6,x2=4.e-6,y1=-4.e-6,y2=4.e-6,z1=-6e-6,z2=6e-6)
# Particles definition
P = gen_particles(1, box)
# Particle simulation definition
S = ParticlesSimulation(D=D, t_step=t_step, particles=P, box=box, psf=psf)
After running brownian.py
we can define some custom parameters:
%run brownian.py
# Initialize the random seed to a fixed value (for notebook reproducibility)
np.random.seed(6)
# Diffusion coefficient
Du = 12.0 # um^2 / s
D = Du*(1e-6)**2
# Simulation time step
t_step = 0.5e-6 # seconds
# Time duration of the simulation
t_max = .2 # seconds
# All the particle starts in the origin for this test
num_part = 20
P = [Particle(x0=0, y0=0, z0=0) for _ in xrange(num_part)]
# Particle simulation definition
S = ParticlesSimulation(D=D, t_step=t_step, t_max=t_max, particles=P, box=box, psf=psf)
The most important line is the last line which creates an object S
that contains all the simulation parameters (it also contains methods to run
the simulation). You can print S
and check the current parameters:
S
Box: X 8.0um, Y 8.0um, Z 12.0um D 1.2e-11, #Particles 20, t_step 0.5us, t_max 0.2s EID_ID 0 0
or check the required RAM for the current parameters:
S.print_RAM()
Number of particles: 20 Number of time steps: 400000 Emission array size: 3.0 MB (total_emission=True) Emission array size: 60.0 MB (total_emission=False) Position array size: 180.0 MB
Now we launch the simulation with the previous parameters (2 particles, 0.2s total time).
We simulate the brownian motion trajectories and the emission at the same time.
The trajectories are saved in the array S.pos
, the emission in the array S.em
.
The emission represents the emission rate in each time step (0.5us).
S.sim_motion_em(total_emission=True, delete_pos=False, wrap=False)
[3968] Simulating particle 0 [3968] Simulating particle 1 [3968] Simulating particle 2 [3968] Simulating particle 3 [3968] Simulating particle 4 [3968] Simulating particle 5 [3968] Simulating particle 6 [3968] Simulating particle 7 [3968] Simulating particle 8 [3968] Simulating particle 9 [3968] Simulating particle 10 [3968] Simulating particle 11 [3968] Simulating particle 12 [3968] Simulating particle 13 [3968] Simulating particle 14 [3968] Simulating particle 15 [3968] Simulating particle 16 [3968] Simulating particle 17 [3968] Simulating particle 18 [3968] Simulating particle 19
x, y, z = S.pos[:, 0, :], S.pos[:, 1, :], S.pos[:, 2, :]
dx, dy, dz = np.diff(x), np.diff(y), np.diff(z)
r_squared = x**2 + y**2 + z**2
dr_squared = dx**2 + dy**2 + dz**2
dr_squared.mean()/(2*3*t_step), D
(1.1995038121250671e-11, 1.2e-11)
time = S.time()
time[1] - time[0]
4.9999999999999998e-07
t_step
5e-07
plt.plot(time, r_squared.T, alpha=0.1);
plt.plot(time, r_squared.mean(0), 'r', lw=3, alpha=0.8);
plt.plot(time, (2*3*D)*time, '--k', lw=3, alpha=0.5)
[<matplotlib.lines.Line2D at 0x1334e908>]