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.
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
/home/anto/Documents/ucla/src/brownian
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)
# Simulation time step
t_step = 0.5e-6 # seconds
# Time duration of the simulation
t_max = .2 # seconds
# Particles definition
P = gen_particles(2, box)
# 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 2, 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: 2 Number of time steps: 400000 Emission array size: 3.0 MB (total_emission=True) Emission array size: 6.0 MB (total_emission=False) Position array size: 18.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)
[20982] Simulating particle 0 [20982] Simulating particle 1
In very long simulations, setting delete_pos=True
allows to save RAM by only saving
the emission array.
Here, since we kept the trajectories (delete_pos=False
), we can plot them:
# This requires delete_pos=False in .sim_motion_em()
plot_tracks(S)
Of the two particles only the second (in green) happens to pass through the excitation volume (approximately shown as a black ellipsis). Plotting the emission we see only one "burst".
plot_emission(S)
Generating timestamps from the emission requires two steps.
Firstly we extract the number of emitted photons in each time bin (here 0.5us)
using a Poisson distribution. Photons can be generated by molecules emission or by
a user-defined (constant) background rate. The array of emission events is saved in
S.tt
(stands for timetrace):
S.sim_timetrace(max_em_rate=3e5, bg_rate=2e3)
plot_timetrace(S, scroll_gui=False)
plt.legend(['Emission rate', 'Emitted photons + BG'])
<matplotlib.legend.Legend at 0x7de9750>
Note: in the previous timetrace, photons are rebinned with a 1 ms bin width, much larger than the simulation time step (0.5us).
Once we generated the timetrace (statistically) we compute the timestamps (deterministically)
using for each photon its bin position (so we end up with timestamps with 0.5us resolution).
Timestamps are saved in S.ph_times
:
S.gen_ph_times()
print S.ph_times[:10], '...'
[ 9.00000000e-06 2.54000000e-04 2.59000000e-04 3.69500000e-04 5.87500000e-04 1.99250000e-03 2.12150000e-03 2.80350000e-03 3.45250000e-03 3.72400000e-03] ...
plt.hist(S.ph_times, bins=np.r_[:200]*1e-3, histtype='step');
Let verify that timetrace and timestamps are coherent. We can compute a 1ms-bin-width
histogram from the two arrays and obtain the same result. For the timestamps it's easy,
we just need to use the histogram
function. For the timetrace, that is an already binned
array, we need to rebin with a 2000 times larger bin (2000 = 1ms/0.5us):
# Binning of timestamps
timestamps_hist = np.histogram(S.ph_times, bins=np.r_[:201]*1e-3)
# Rebinning of timetrace
rebin = 2000
rebinned_timetrace = S.tt.reshape(-1, rebin).sum(axis=1)
rebinned_time = (np.arange(rebinned_timetrace.size)+1)*(S.t_step*rebin)
# Comparison plot
plt.plot(timestamps_hist[1][1:], timestamps_hist[0])
plt.plot(rebinned_time, rebinned_timetrace)
# Rigorous check
assert ((timestamps_hist[0] - rebinned_timetrace) == 0).all()
The two plots are identical and the assertion passes demonstrating that the two array are numerically equal.
NOTE on path: The file path_def.py
(loaded by brownian.py
) defines the SIM_DATA_DIR
variable that contains the path where to save the simulation. You can edit path_def.py
to choose a folder and use here SIM_DATA_DIR
. You can also skip that step and paste your
path in the following command (replacing SIM_DATA_DIR):
S.dump(dir_=SIM_DATA_DIR)
Saved to /home/anto/Documents/ucla/src/data/sim/brownian/objects//bromo_sim_D1.2e-11_2P_4pM_step0.5us_t_max0.2s_ID0-0.pickle
Let load the simulation from a file. We know we saved in folder SIM_DATA_DIR
with ID=0, EID=0
and the file name contains 0.2s
(simulation duration):
S2 = load_sim_id(ID=0, EID=0, glob_str='*0.2s*', dir_=SIM_DATA_DIR)
Loaded: /home/anto/Documents/ucla/src/data/sim/brownian/objects/bromo_sim_D1.2e-11_2P_4pM_step0.5us_t_max0.2s_ID0-0.pickle
S2
Box: X 8.0um, Y 8.0um, Z 12.0um D 1.2e-11, #Particles 2, t_step 0.5us, t_max 0.2s EID_ID 0 0
S2.em
array([[ 7.84983849e-07, 7.84919671e-07, 7.92120295e-07, ..., 1.76593796e-05, 1.76628897e-05, 1.76252189e-05]])
The arrays are different but they have the same content:
print id(S2.em), id(S.em)
130640560 45489520
(S.em == S2.em).all()
True
from IPython.core.display import HTML
def css_styling():
styles = open("./styles/custom2.css", "r").read()
return HTML(styles)
css_styling()