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[:]

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 [ ]: