PyBroMo - B.1 Disk-single-core - Generate photon timestamps

This notebook is part of PyBroMo a python-based single-molecule Brownian motion diffusion simulator that simulates confocal smFRET experiments.

Overview

In this notebook we show how to generated timestamps of emitted-photons from saved diffusion traces.



In [ ]:

    
%matplotlib inline
import numpy as np
import tables
import matplotlib.pyplot as plt
import seaborn as sns
import pybromo as pbm
print('Numpy version:', np.__version__)
print('PyTables version:', tables.__version__)
print('PyBroMo version:', pbm.__version__)

Timestamps simulation

As a memo, let's write some formulas related to the FRET efficiency:

$$ k = \frac{F_a}{F_d} \quad,\qquad E = \frac{k}{k+1} \qquad\Rightarrow\qquad k = \frac{E}{1-E}$$



In [ ]:

    
S = pbm.ParticlesSimulation.from_datafile('016', mode='w')



In [ ]:

    
def em_rates_from_E(em_rate_tot, E_values):
    E_values = np.asarray(E_values)
    em_rates_a = E_values * em_rate_tot
    em_rates_d = em_rate_tot - em_rates_a

    k_values = E_values/(1 - E_values)
    assert np.allclose((em_rates_a/em_rates_d), k_values)

    em_rates = np.hstack([em_rates_a, em_rates_d])
    return em_rates



In [ ]:

    
em_rate_tot = 200e3
E_list = np.array([0.01, 0.02, 0.05, 0.1, 0.2, 0.4])

em_rate_list = em_rates_from_E(em_rate_tot, E_list)
em_rate_list



In [ ]:

    
# Get the random state at the end of the diffusion simulation
saved_rs_state = S.traj_group._v_attrs['last_random_state']
pbm.hash_(saved_rs_state)

Simulation of the series of emission rates

Here we perform the timetsamps simulation for the list of emission rates previously computed.

At this point we also choose the Poisson background level to add to the simulation.



In [ ]:

    
em_rate_list

Simulate timestamps for background = 1kcps:



In [ ]:

    
rs = np.random.RandomState()
rs.set_state(saved_rs_state)



In [ ]:

    
%%timeit -n1 -r1
for em_rate in em_rate_list:
    print(' Emission rate: ', em_rate, flush=True)
    S.simulate_timestamps_mix(max_rates=(em_rate,), populations=(slice(0, 20),), 
                              bg_rate=1e3, rs=rs)

Simulate timestamps for background = 4kcps:



In [ ]:

    
%%timeit -n1 -r1
for em_rate in em_rate_list:
    print(' Emission rate: ', em_rate, flush=True)
    S.simulate_timestamps_mix(max_rates=(em_rate,), populations=(slice(0, 20),), 
                              bg_rate=4e3, rs=rs)



In [ ]:

    
for k in S.ts_store.h5file.root.timestamps._v_children.keys():
    if not k.endswith('_par'):
        print(k)



In [ ]:

    
ts, ts_par = S.get_timestamps_part('Pop1_P20_Pstart0_max_rate198000cps_BG4000cps_t_1s_rs_8798a6')



In [ ]:

    
ts[:]

array([ 8650, 19320, 22290, ..., 19990800, 19994140, 19996650])



In [ ]:

    
bins = np.arange(0, 1, 1e-3)
plt.hist(ts*ts.attrs['clk_p'], bins=bins, histtype='step');

Verify the simulation

Check that the new arrays show up in the data file:



In [ ]:

    
group = '/timestamps'

print('Nodes in in %s:\n' % group)

print(S.ts_store.h5file.get_node(group))
for node in S.ts_store.h5file.get_node(group)._f_list_nodes():
    print('\t%s' % node.name)
    #print('\t    %s' % node.title)



In [ ]:

    
[t for t in S.timestamp_names if 'BG4000cps' in t]



In [ ]:

    
S.ts_store.close()



In [ ]:



In [ ]: