Tutorial about fluopy

Here we outline the principle use of fluopy.

from pprint import pprint

%matplotlib inline

import matplotlib.pyplot as plt
import numpy as np

import fluopy
import fluopy.analysis as an
import fluopy.blinking as bl
import fluopy.emissions as em
import fluopy.fcs as fcs_p
import fluopy.figure as fi
import fluopy.fluorophores as fl
import fluopy.miscellaneous as mi
import fluopy.prediction as pr
import fluopy.simulation as si
import fluopy.transitions as tr

fluopy.__version__

'0.4.0.dev4+gc3c2cb30d'

For random number generation a Generator is initialized.

rng = np.random.default_rng(seed=1)

Modules

Fluopy contains the following modules. For detailed information look up the API documentation.

print("Fluopy modules:")
print()
pprint([item for item in dir(fluopy) if not item.startswith("_")])

Fluopy modules:

['Analysis',
 'Blinking',
 'Emissions',
 'ExponentialMixtureMarginalModel',
 'ExponentialMixtureModel',
 'FCS',
 'Fluorophore',
 'FluorophoreData',
 'FluorophoreSystem',
 'PairedState',
 'Photoswitching_fingerprint_model',
 'Prediction',
 'Simulation',
 'SingleState',
 'Transition',
 'TransitionSet',
 'analysis',
 'atto643',
 'blinking',
 'calculate_emission_rate',
 'calculate_excitation_rate',
 'calculate_fret_efficiency',
 'calculate_fret_rate',
 'calculate_internal_conversion_rate',
 'calculate_pet_rate',
 'calculate_photon_collection_rate',
 'calculate_photon_flux',
 'calculate_spectral_overlap_integral',
 'check_graph_suitable',
 'construct_state_graphs',
 'construct_transition_graph',
 'convert_wavenumber_wavelength_frequency',
 'cy5_dna',
 'determine_node_order',
 'distributions',
 'draw_networkx_curved_edge_labels',
 'emissions',
 'fcs',
 'figure',
 'fitting',
 'fluo_data',
 'fluorophores',
 'formulas',
 'generate_combinations',
 'get_pis',
 'henderson_hasselbalch_equation',
 'hypoexponential_distribution_cdf',
 'hypoexponential_distribution_pdf',
 'hypoexponential_distribution_pdf_1st_order_derivative',
 'hypoexponential_distribution_pdf_2nd_order_derivative',
 'import_module',
 'integral_kappa_squared',
 'kappa_squared',
 'logger',
 'logging',
 'map_to_lambdas',
 'miscellaneous',
 'module_',
 'network',
 'photoswitching_fingerprint_prepare',
 'plot_graph',
 'prediction',
 'routines',
 'sample_kappa_squared_distribution',
 'simulate_rotational_motion',
 'simulation',
 'simulation_tcspc',
 'submodule',
 'submodules',
 'transitions',
 'universal_figure']

Introduction

Fluopy helps to run Markov chain simulations for a photophysical system of isolated or communicating fluorophores.

To get started the photophysical system has to be defined and instantiated.

This consists of the following steps:

• Define a fluorophore system that contains fluorophores at positions in space.

• Define a set of possible transitions for the fluorophore system.

A simulation provides a series of transition events that occured in a single realization of the Markov chain.

A simulation is run for a given fluorophore system and transition set with a random number generator initialized as seed.

For a fluorophore system and transition set that fullfill certain requirements a prediction for the simulation outcome can be derived based on computational estimation of a steady-state equilibrium.

Finally, the acquired data is analysed and presented in various formats.

Analysis procedure that focus on experimentally observable transitions are available:

• Time traces

• Fluorescence correlation spectroscopy

• ON/Off time distributions

A simulation with pulsed excitation can be carried out. The analysis resembles a TSCPC experiment.

A typical setup, simulation and analysis procedure for selected photophysical systems is presented in the other tutorials.

In the following we will present the relevant code blocks.

Transitions

There are various transition types that connect the photophysical states:

list(tr.TransitionType.__members__)

['EXCITATION',
 'FLUORESCENT_EMISSION',
 'SINGLET_QUENCHING',
 'INTERSYSTEM_CROSSING_ST',
 'INTERSYSTEM_CROSSING_TS',
 'INTERNAL_CONVERSION_S',
 'REVERSE_INTERSYSTEM_CROSSING',
 'PHOTOBLEACHING_1',
 'PHOTOBLEACHING_2',
 'ET_CYCLE_T',
 'ET_CYCLE_S',
 'ADDUCT_T',
 'ADDUCT_S',
 'THERM_ELIMINATION',
 'PHOTO_UNCAGING',
 'RAD_ESCAPE',
 'RAD_RELAX',
 'ISOMERIZATION',
 'PHOTO_BISO',
 'THERM_BISO',
 'FRET',
 'CIS_FRET_1',
 'CIS_FRET_2',
 'OFF_FRET_1',
 'OFF_FRET_2',
 'S_S_ANNIHILATION',
 'S_T_ANNIHILATION',
 'S_T_ANNI_RISC',
 'S_T_ANNI_BLEACH',
 'R_FRET_1',
 'R_FRET_2',
 'H2O_ATTACK_S',
 'H2O_ATTACK_T',
 'BACK_REACTION',
 'S1_S0_TRANSITIONS',
 'CIS_S0_TRANSITIONS',
 'T1_S0_TRANSITIONS',
 'OFF_S0_TRANSITIONS']

Photophysical states are named according to the Enum class tr.SingleState and tr.PairedState, which link state IDs to photophysical state names.

A single fluorophore can occupy the following photophysical states:

list(tr.SingleState.__members__)

['S0', 'S1', 'S2', 'T1', 'T2', 'B', 'cis', 'OFF', 'OFF2', 'R']

Two interacting fluorophores can occupy the following (paired) photophysical states:

list(tr.PairedState.__members__)

['S1_S0',
 'S0_S1',
 'S1_T1',
 'S1_Cis',
 'S0_Cis',
 'S1_OFF',
 'S0_S0',
 'S0_T2',
 'S1_S1',
 'S0_T1',
 'S0_OFF2',
 'S0_OFF',
 'S0_B',
 'S1_R',
 'S0_R']

A specific transition of a certain type from state 1 to state 2 is defined as:

transition = tr.Transition(
    transition_type=tr.TransitionType.EXCITATION, rate=1e6, fluorophore_ids=[1]
)
pprint(transition)

Transition(identity=None,
           transition_type=<TransitionType.EXCITATION: TransitionAttributes(abbreviation='EXC', initial_state=<SingleState.S0: 0>, final_state=<SingleState.S1: 1>, photon=False)>,
           abbreviation='EXC',
           initial_state=<SingleState.S0: 0>,
           final_state=<SingleState.S1: 1>,
           rate=1000000.0,
           photon=False,
           fluorophore_ids=[1])

A photophysical system is made up of one or more fluorophores or fluorophore pairs with multiple transitions:

transition_1 = tr.Transition(
    tr.TransitionType.EXCITATION, rate=1e6, fluorophore_ids=[0]
)
transition_2 = tr.Transition(
    tr.TransitionType.FLUORESCENT_EMISSION, rate=1e9, fluorophore_ids=[0]
)
transitions = {
    "cy5_dna": [transition_1, transition_2],
}
pprint(transitions)

{'cy5_dna': [Transition(identity=None,
                        transition_type=<TransitionType.EXCITATION: TransitionAttributes(abbreviation='EXC', initial_state=<SingleState.S0: 0>, final_state=<SingleState.S1: 1>, photon=False)>,
                        abbreviation='EXC',
                        initial_state=<SingleState.S0: 0>,
                        final_state=<SingleState.S1: 1>,
                        rate=1000000.0,
                        photon=False,
                        fluorophore_ids=[0]),
             Transition(identity=None,
                        transition_type=<TransitionType.FLUORESCENT_EMISSION: TransitionAttributes(abbreviation='FLU', initial_state=<SingleState.S1: 1>, final_state=<SingleState.S0: 0>, photon=True)>,
                        abbreviation='FLU',
                        initial_state=<SingleState.S1: 1>,
                        final_state=<SingleState.S0: 0>,
                        rate=1000000000.0,
                        photon=True,
                        fluorophore_ids=[0])]}

The fluorophore system

A photophysical system is defined by one or more fluorophores that occupy electronic states and go through transitions from one state to another.

A number of input parameters are kept in fluopy.fluo_data as FluorophoreData instances. These are automatically inserted if the name matches.

fluorophore = fl.Fluorophore(name="some_dye", position=[0, 0])
pprint(fluorophore)

There is no FluorophoreData for Fluorophore some_dye in fluopy.fluo_data. Parameters have to be defined manually.

Fluorophore(identity=None,
            name='some_dye',
            position=array([0, 0]),
            constants=None)

fluorophore = fl.Fluorophore(name="cy5_dna", position=[0, 0])
pprint(fluorophore)

Fluorophore(identity=None,
            name='cy5_dna',
            position=array([0, 0]),
            constants=FluorophoreData(data_files='cy5_data',
                                      QUANTUM_YIELD=0.27,
                                      FLUORESCENCE_LIFETIME=1.7e-09,
                                      S1_QUENCH_RATE=0,
                                      ISC_ST_RATE=830000.0,
                                      ISC_TS_RATE=5000.0,
                                      RISC_RATE=0,
                                      STA_EFFICIENCY=0,
                                      PHOTOBLEACH_T1_RATE=10.0,
                                      PHOTOBLEACH_T2_RATE=0,
                                      DSTORM_PET_T_RATE_MOL=100000000.0,
                                      DSTORM_PET_S_RATE_MOL=1000000000.0,
                                      DSTORM_PET_SUCCESS_RATE=0.001,
                                      DSTORM_TH_EL_RATE_1=0.01,
                                      DSTORM_TH_EL_RATE_2=0,
                                      DSTORM_P_EL_CROSS_SECTION=6e-24,
                                      RAD_ESCAPE_EFFICIENCY=0.01,
                                      RAD_RELAX_RATE=1300.0,
                                      OFRET_EFFICIENCY=0.001,
                                      ISO_RATE=4000000.0,
                                      BISO_CROSS_SECTION=6e-18,
                                      BISO_THERMAL_RATE=5000.0,
                                      BISO_EFFICIENCY=0.04,
                                      H2O_ATTACK_S=0,
                                      H2O_ATTACK_T=0,
                                      BACK_REACTION=0))

A fluophore system contains the fluorophore identities, their spatial arrangement and transition rate constants.

fluorophores = [fluorophore]
fluorophore_system = fl.FluorophoreSystem(fluorophores=fluorophores)

fluorophore_system.plot(scale=1, quadratic=True);

Fluorophore data

The transition rate constants are specific for a kind of fluorophore and the chemical environment.

Standard datasets with all constant photophysical attributes of the fluorophore are defined for various fluorophores.

pprint(fluopy.fluo_data.atto643)

FluorophoreData(data_files='atto643_data',
                QUANTUM_YIELD=0.6,
                FLUORESCENCE_LIFETIME=3e-09,
                S1_QUENCH_RATE=0,
                ISC_ST_RATE=1000000.0,
                ISC_TS_RATE=100000.0,
                RISC_RATE=0,
                STA_EFFICIENCY=0,
                PHOTOBLEACH_T1_RATE=1,
                PHOTOBLEACH_T2_RATE=0,
                DSTORM_PET_T_RATE_MOL=0,
                DSTORM_PET_S_RATE_MOL=0,
                DSTORM_PET_SUCCESS_RATE=0,
                DSTORM_TH_EL_RATE_1=0,
                DSTORM_TH_EL_RATE_2=0,
                DSTORM_P_EL_CROSS_SECTION=0,
                RAD_ESCAPE_EFFICIENCY=0,
                RAD_RELAX_RATE=0,
                OFRET_EFFICIENCY=0,
                ISO_RATE=0,
                BISO_CROSS_SECTION=0,
                BISO_THERMAL_RATE=0,
                BISO_EFFICIENCY=0,
                H2O_ATTACK_S=30000.0,
                H2O_ATTACK_T=0,
                BACK_REACTION=0.1)

pprint(fluopy.fluo_data.cy5_dna)

FluorophoreData(data_files='cy5_data',
                QUANTUM_YIELD=0.27,
                FLUORESCENCE_LIFETIME=1.7e-09,
                S1_QUENCH_RATE=0,
                ISC_ST_RATE=830000.0,
                ISC_TS_RATE=5000.0,
                RISC_RATE=0,
                STA_EFFICIENCY=0,
                PHOTOBLEACH_T1_RATE=10.0,
                PHOTOBLEACH_T2_RATE=0,
                DSTORM_PET_T_RATE_MOL=100000000.0,
                DSTORM_PET_S_RATE_MOL=1000000000.0,
                DSTORM_PET_SUCCESS_RATE=0.001,
                DSTORM_TH_EL_RATE_1=0.01,
                DSTORM_TH_EL_RATE_2=0,
                DSTORM_P_EL_CROSS_SECTION=6e-24,
                RAD_ESCAPE_EFFICIENCY=0.01,
                RAD_RELAX_RATE=1300.0,
                OFRET_EFFICIENCY=0.001,
                ISO_RATE=4000000.0,
                BISO_CROSS_SECTION=6e-18,
                BISO_THERMAL_RATE=5000.0,
                BISO_EFFICIENCY=0.04,
                H2O_ATTACK_S=0,
                H2O_ATTACK_T=0,
                BACK_REACTION=0)

The data is collected in a dataclass specifying possible constants fluopy.fluo_data.FluorophoreData.

fluorophore_data = fluopy.fluo_data.FluorophoreData()
pprint(fluorophore_data)

FluorophoreData(data_files=None,
                QUANTUM_YIELD=0,
                FLUORESCENCE_LIFETIME=0,
                S1_QUENCH_RATE=0,
                ISC_ST_RATE=0,
                ISC_TS_RATE=0,
                RISC_RATE=0,
                STA_EFFICIENCY=0,
                PHOTOBLEACH_T1_RATE=0,
                PHOTOBLEACH_T2_RATE=0,
                DSTORM_PET_T_RATE_MOL=0,
                DSTORM_PET_S_RATE_MOL=0,
                DSTORM_PET_SUCCESS_RATE=0,
                DSTORM_TH_EL_RATE_1=0,
                DSTORM_TH_EL_RATE_2=0,
                DSTORM_P_EL_CROSS_SECTION=0,
                RAD_ESCAPE_EFFICIENCY=0,
                RAD_RELAX_RATE=0,
                OFRET_EFFICIENCY=0,
                ISO_RATE=0,
                BISO_CROSS_SECTION=0,
                BISO_THERMAL_RATE=0,
                BISO_EFFICIENCY=0,
                H2O_ATTACK_S=0,
                H2O_ATTACK_T=0,
                BACK_REACTION=0)

Transition set

A set of possible transitions can be generated from a given fluorophore system by helper functions:

transitions = fluorophore_system.load_transitions(
    summarize=False,
    irradiance=2.5,
    wavelength=640,
    bleaching=False,
    energy_transfer=False,
    dstorm=False,
)
pprint(transitions)

{'cy5_dna': [Transition(identity=None,
                        transition_type=<TransitionType.EXCITATION: TransitionAttributes(abbreviation='EXC', initial_state=<SingleState.S0: 0>, final_state=<SingleState.S1: 1>, photon=False)>,
                        abbreviation='EXC',
                        initial_state=<SingleState.S0: 0>,
                        final_state=<SingleState.S1: 1>,
                        rate=np.float64(7269625.549132028),
                        photon=False,
                        fluorophore_ids=[0]),
             Transition(identity=None,
                        transition_type=<TransitionType.FLUORESCENT_EMISSION: TransitionAttributes(abbreviation='FLU', initial_state=<SingleState.S1: 1>, final_state=<SingleState.S0: 0>, photon=True)>,
                        abbreviation='FLU',
                        initial_state=<SingleState.S1: 1>,
                        final_state=<SingleState.S0: 0>,
                        rate=np.float64(158823529.4117647),
                        photon=True,
                        fluorophore_ids=[0]),
             Transition(identity=None,
                        transition_type=<TransitionType.INTERSYSTEM_CROSSING_ST: TransitionAttributes(abbreviation='ISC_ST', initial_state=<SingleState.S1: 1>, final_state=<SingleState.T1: 3>, photon=False)>,
                        abbreviation='ISC_ST',
                        initial_state=<SingleState.S1: 1>,
                        final_state=<SingleState.T1: 3>,
                        rate=830000.0,
                        photon=False,
                        fluorophore_ids=[0]),
             Transition(identity=None,
                        transition_type=<TransitionType.INTERSYSTEM_CROSSING_TS: TransitionAttributes(abbreviation='ISC_TS', initial_state=<SingleState.T1: 3>, final_state=<SingleState.S0: 0>, photon=False)>,
                        abbreviation='ISC_TS',
                        initial_state=<SingleState.T1: 3>,
                        final_state=<SingleState.S0: 0>,
                        rate=5000.0,
                        photon=False,
                        fluorophore_ids=[0]),
             Transition(identity=None,
                        transition_type=<TransitionType.ISOMERIZATION: TransitionAttributes(abbreviation='ISO', initial_state=<SingleState.S1: 1>, final_state=<SingleState.cis: 6>, photon=False)>,
                        abbreviation='ISO',
                        initial_state=<SingleState.S1: 1>,
                        final_state=<SingleState.cis: 6>,
                        rate=4000000.0,
                        photon=False,
                        fluorophore_ids=[0]),
             Transition(identity=None,
                        transition_type=<TransitionType.PHOTO_BISO: TransitionAttributes(abbreviation='PBISO', initial_state=<SingleState.cis: 6>, final_state=<SingleState.S0: 0>, photon=False)>,
                        abbreviation='PBISO',
                        initial_state=<SingleState.cis: 6>,
                        final_state=<SingleState.S0: 0>,
                        rate=np.float64(48327.51904841),
                        photon=False,
                        fluorophore_ids=[0]),
             Transition(identity=None,
                        transition_type=<TransitionType.THERM_BISO: TransitionAttributes(abbreviation='TBISO', initial_state=<SingleState.cis: 6>, final_state=<SingleState.S0: 0>, photon=False)>,
                        abbreviation='TBISO',
                        initial_state=<SingleState.cis: 6>,
                        final_state=<SingleState.S0: 0>,
                        rate=5000.0,
                        photon=False,
                        fluorophore_ids=[0]),
             Transition(identity=None,
                        transition_type=<TransitionType.INTERNAL_CONVERSION_S: TransitionAttributes(abbreviation='IC', initial_state=<SingleState.S1: 1>, final_state=<SingleState.S0: 0>, photon=False)>,
                        abbreviation='IC',
                        initial_state=<SingleState.S1: 1>,
                        final_state=<SingleState.S0: 0>,
                        rate=np.float64(424581764.7058823),
                        photon=False,
                        fluorophore_ids=[0])]}

A transition set has to be defined from the individual transitions for a given fluorophore system:

transition_set = tr.TransitionSet(transitions, fluorophore_system)

TransitionSet provides helper methods so that a transition_set can be easily modified with the following methods:

transition_set_no_bleaching = transition_set.remove_absorbing_states()
transition_set_no_ets = transition_set.remove_energy_transfers()
transition_set_filtered = transition_set.filter_by_identity(remove_list=[6])

transition_set_adj = transition_set.adjust_rates(
    change_dict={1: 1e3}, keep_zero_rates=False
)

If for no other reasons you should remove transitions with rate constants of zero (this has possibly already been done above through keep_zero_rates=False):

transtition_set = transition_set.remove_zero_rates()

A transition set has to be finalized before a simulation is started:

transition_set.finalize()
list(vars(transition_set).keys())

['transitions',
 'fluorophore_system',
 'transition_df',
 'single_states',
 '_combined_state_transitions_df',
 '_row_sums',
 '_transition_matrix']

A graphical representation of all possible transitions is given:

transition_set.plot(graph_type="shell", colors=None, scale=1);

All transitions are listed in the transitionsattribute or in a dataframe:

transition_set.transition_df

		transition_type	abbreviation	initial_state	final_state	rate	photon	fluorophore_ids	absorbing
Fluorophore	identity
cy5_dna	0	TransitionType.EXCITATION	EXC	SingleState.S0	SingleState.S1	7.269626e+06	False	[0]	False
	1	TransitionType.FLUORESCENT_EMISSION	FLU	SingleState.S1	SingleState.S0	1.588235e+08	True	[0]	False
	2	TransitionType.INTERSYSTEM_CROSSING_ST	ISC_ST	SingleState.S1	SingleState.T1	8.300000e+05	False	[0]	False
	3	TransitionType.INTERSYSTEM_CROSSING_TS	ISC_TS	SingleState.T1	SingleState.S0	5.000000e+03	False	[0]	False
	4	TransitionType.ISOMERIZATION	ISO	SingleState.S1	SingleState.cis	4.000000e+06	False	[0]	False
	5	TransitionType.PHOTO_BISO	PBISO	SingleState.cis	SingleState.S0	4.832752e+04	False	[0]	False
	6	TransitionType.THERM_BISO	TBISO	SingleState.cis	SingleState.S0	5.000000e+03	False	[0]	False
	7	TransitionType.INTERNAL_CONVERSION_S	IC	SingleState.S1	SingleState.S0	4.245818e+08	False	[0]	False

Note: In transition_set.single_states, only the ‘relevant’ states are mentioned - the ones that actually occur in transition_set.plot().

Note: Transitions that have a rate of 0 or are not usable by the fluorophore do not occur transition_set.single_states; however, transitions that may not be reachable still do.

Single dye simulation setup

For further demonstration we set up a simple fluorophore system that contains a single Cy5 fluorophore.

Define the fluorophore system

fluorophore = fl.Fluorophore(name="cy5_dna", position=[0, 0])
fluorophore_system = fl.FluorophoreSystem(fluorophores=[fluorophore])

Define the transition set

transitions = fluorophore_system.load_transitions(
    summarize=True,
    irradiance=2,
    wavelength=640,
    bleaching=False,
    energy_transfer=False,
    dstorm=False,
)
transition_set = tr.TransitionSet(transitions, fluorophore_system)
transition_set = transition_set.remove_energy_transfers()
transition_set.finalize()

<fluopy.transitions.TransitionSet at 0x7827f0f58050>

Make a prediction

For a fluorophore system and transition set that fullfill certain requirements a prediction for the simulation outcome can be derived based on computational estimation of a steady-state equilibrium.

prediction = pr.Prediction(transition_set)
list(vars(prediction).keys())

['energy_transfer',
 'absorbing_chain',
 'transition_set',
 'frequency_transitions',
 'frequency_states',
 'transition_time_distributions',
 'lifetime_distributions',
 'mean_transition_times',
 'mean_lifetimes',
 'state_occupations']

The results are accessible through methods of the prediction instance.

prediction.plot_frequency_transitions(scale=0.5)
prediction.plot_frequency_states(scale=0.5)
prediction.plot_mean_lifetimes(scale=0.5)
prediction.plot_mean_transition_times(scale=0.5)
prediction.plot_state_occupations(scale=0.5)
prediction.plot_transition_time_distributions(
    fluorophore="cy5_dna", transition_id=0, scale=0.5
);

Note for multi-fluorophore systems:

• if energy transfers occur, pairs have separate color in plot but are normalized with other transitions of their donor

• absorbing states: the prediction works using the fundamental matrix. For Q, the transitions leading to an OVERALL absorbed state are dropped. This can be done since bleaching is usually by far the least occurring transition, meaning each fluorophore will reach its limiting distribution (as if it was no absorbing Markov Chain) the frequencies and state occupations of absorbing transitions and states are set to 0, their lifetimes to inf (this is different to analysis, where - if absorbing transitions happened - only the lifetime and the state occupation is 0)

Run a simulation

To carry out the complete Markov chain simulation, a Simulation object has to be instantiated from the transition set.

The simulation can then run for a given number of transition steps, or until it reaches a defined end time.

%%time
simulation = si.Simulation(transition_set)
list(vars(simulation).keys())

CPU times: user 25 μs, sys: 3 μs, total: 28 μs
Wall time: 31.9 μs

['transition_set',
 'time_series',
 'transition_series',
 'state_series',
 'memmap_path']

%%time
# simulate until it reaches given end_time
simulation.run(start_at=None, size=1e6, end_time=1, seed=rng, use_memmap=None)
mi.print_class(simulation)

Floating point precision error warning:
 The smallest safe increment is 2.22e-16.
 Everything drawn below this number might be rounded to zero
 when approaching the time limit of this simulation.
 Using the highest possible rate which occurs for example in state combination [1]
 gives a probability of 1.31e-07 for a smaller increment to be drawn.

Attributes of <fluopy.simulation.Simulation object at 0x782822240d70>:
.................................................................
transition_set = <fluopy.transitions.TransitionSet object at 0x7827f0f58050>
_________________________________________________________________
time_series = array([0.00000000e+00, 1.15167401e-07, 1.1846001....99893691e-01, 1.00000000e+00], shape=(3302282,))
_________________________________________________________________
transition_series = array([0, 4, 0, ..., 4, 0, 2], shape=(3302280,), dtype=uint32)
_________________________________________________________________
state_series = array([[0, 1, 0, ..., 0, 1, 3]], shape=(1, 3302281), dtype=int8)
_________________________________________________________________
memmap_path = None
_________________________________________________________________


CPU times: user 9.83 s, sys: 37 ms, total: 9.87 s
Wall time: 9.87 s

If the simulation goes through many transtition events, the arrays may become too large for the available RAM. In that case, use memory maps to store the data on the local drive.

Run an approximate simulation

Under certain requirements (similar to prediction) a simulation approximation can can be carried out, saving time.

%%time
approximation = si.Simulation(transition_set=transition_set)
approximation.approximate(prediction=prediction, size=1e6, seed=rng)
list(vars(approximation).keys())

Floating point precision error warning:
 The higher limit of smallest increment with a probability of 1.00e-03 is 1.70e-12.
 This was estimated using the highest possible rate which occurs for example in state combination [1].
 Everything drawn below this number will be rounded to zero starting somewhere between 1.00e+03 - 1.00e+04.

CPU times: user 89.8 ms, sys: 12.9 ms, total: 103 ms
Wall time: 103 ms

['transition_set',
 'time_series',
 'transition_series',
 'state_series',
 'memmap_path']

Analyze the simulation

Analyze the data from a given simulation for state probabilities and lifetimes, among others.

analysis = an.Analysis(simulation=simulation)
list(vars(analysis).keys())

['simulation',
 'frequency_transitions',
 'frequency_states',
 'transition_time_distributions',
 'lifetime_distributions',
 'mean_transition_times',
 'mean_lifetimes',
 'state_occupations']

mi.print_class(analysis)

Attributes of <fluopy.analysis.Analysis object at 0x7827f3183e00>:
.................................................................
simulation = <fluopy.simulation.Simulation object at 0x782822240d70>
_________________________________________________________________
frequency_transitions = array([0.4979499 , 0.13429449, 0.00072193, 0.00337857, 0.35955491,
       0.00072162, 0.00337857])
_________________________________________________________________
frequency_states = {'cy5_dna': array([0.49794975, 0.49794975, 0.00072193, 0.00337857])}
_________________________________________________________________
transition_time_distributions = [array([1.15167401e-07, 2.00370096e-07, 1.0292415....05830117e-08, 1.19170847e-08], shape=(1644370,)), array([2.70692572e-09, 1.04365849e-09, 1.1245765...1.39646961e-10, 7.48977991e-10], shape=(443478,)), array([8.13131906e-11, 1.12316807e-09, 3.2367626..., 5.24545962e-10, 1.79783488e-09], shape=(2384,)), array([2.78520369e-11, 4.96192071e-11, 5.2462151... 5.82091275e-09, 1.47267176e-09], shape=(11157,)), array([3.29261809e-09, 3.21472338e-10, 4.8111506....36261496e-10, 5.08187048e-09], shape=(1187351,)), array([2.09049325e-04, 4.70720765e-05, 3.1177581..., 2.92710445e-05, 7.05170142e-06], shape=(2383,)), ...]
_________________________________________________________________
lifetime_distributions = {'cy5_dna': [array([1.15167401e-07, 2.00370096e-07, 1.0292415....05830117e-08, 1.19170847e-08], shape=(1644370,)), array([3.29261809e-09, 3.21472338e-10, 4.8111506....08187048e-09, 1.79783488e-09], shape=(1644370,)), array([2.09049325e-04, 4.70720765e-05, 3.1177581..., 2.92710445e-05, 7.05170142e-06], shape=(2383,)), array([2.30303056e-05, 1.90388195e-05, 9.5060310... 1.93452378e-05, 2.26451732e-06], shape=(11157,))]}
_________________________________________________________________
mean_transition_times = array([1.72132328e-07, 1.70004097e-09, 1.6551900... 1.69931154e-09, 1.93393009e-04, 2.26937551e-05])
_________________________________________________________________
mean_lifetimes = {'cy5_dna': array([1.72132328e-07, 1.69954958e-09, 1.93393009e-04, 2.26937551e-05])}
_________________________________________________________________
state_occupations = {'cy5_dna': array([0.28302459, 0.00279444, 0.46100879, 0.25317218])}
_________________________________________________________________

In the following analysis plots the simulation results are overlaid with prediction outcomes (crosses) if a prediction is provided.

analysis.get_fluorescence_lifetimes(fluorophore="cy5_dna")
analysis.get_emitting_transition_lifetimes(fluorophore="cy5_dna")

analysis.plot_frequency_transitions(scale=0.5, prediction=prediction)
analysis.plot_frequency_states(scale=0.5, prediction=prediction)
analysis.plot_mean_transition_times(scale=0.5, prediction=prediction)
analysis.plot_mean_lifetimes(scale=0.5, prediction=prediction)
analysis.plot_state_occupations(scale=0.5, prediction=prediction)
analysis.plot_lifetime_distributions(
    scale=0.5, prediction=prediction, fluorophore="cy5_dna", state_identity=1
)
analysis.plot_transition_time_distributions(
    scale=0.5, prediction=prediction, fluorophore="cy5_dna", transition_id=0
);

Note for multi-fluorophore systems:

• The analysis follows these rules:

  • time to transitions of energy transfers are collected only from the donor’s point of view.

  • time to transition of e.g., fluorescence, does not differentiate whether energy transfer was also an option or not one energy transfer transition occurrence resembles the transition of the donor and the acceptor state lifetimes and occurrences do not differentiate whether energy transfer was involved or not.

• Energy transfer has an impact on true fluorescence lifetime only if the acceptors are truly available. If the potential acceptor is in a non-accepting state, the energy transfer is not available. If multiple fluorophores are potential acceptors, all of their respective rates are taken into account. Depending on the photophysical system, energy transfers may only have a marginal impact on the mean fluorescence lifetime.

Simulation of experimentally observable (photons per frames) only

If you are not interested in all hidden state transitions but only in the observable photon emission events, you can limit the simulation to the relevant transitions.

This approach requires less memory than a full simulation.

However, the complete statistical analysis is not available. A photon time series, FCS and blinking analysis is available.

The observation of emitted photons is further influenced by spectral filters, optical elements and detector noise contributions.

For multicolor fluorophore-systems, instantiate two Emissions objects with different bandpasses.

Extract photon emission events from simulation

%%time
emissions = em.Emissions(frame_time="5ms", seed=rng, bandpass=(600, 800))
emissions.extract(simulation=simulation)
emissions

CPU times: user 142 ms, sys: 6 μs, total: 142 ms
Wall time: 142 ms

<fluopy.emissions.Emissions at 0x7827f3182900>

Simulate photon emission events

Helper routines are available to further modify the emission time series and correct for detection efficiency and noise contributions:

emissions.add_photon_collection_objective(p=0.1, seed=rng)  # 1.
emissions.add_quantum_efficiency(p=0.9, seed=rng)  # 3.1.
emissions.add_poisson_noise(
    rate=0.05, seed=rng
)  # 3.2. (dark noise), note the frame time
emissions.add_emccd_gain(emccd_gain=10, seed=rng)  # 4.
emissions.add_gaussian_noise(mean=10, std=1, seed=rng)  # 5. (readout noise)

emissions.plot_time_series();

emissions.plot_histogram();

emissions.plot_cumulative_events();

Simulation of pulsed excitation

A simulation of observable photon emissions with pulsed laser excitation can be carried out. The analysis resembles a TSCPC experiment.

In this simulation, the returned lifetimes are the time differences of photon emission to the last laser pulse. I.e., if the time between laser pulses is very short, it may shorten the observed lifetimes.

The simulation does not discriminate between multiple identical fluorophores (i.e., if one fluorophore gets excited, it transfers the energy to another fluorophore, and this fluorophore emits a photon).

However, if the fluorophores are different, only one of them may be directly excitable by the laser pulse (due to the wavelength) and, depending on the provided bandpass, only the fluorescence of the other may be detected.

If details is True, a complete simulation object will be returned with transition_series, state_series and time_series (requiring more memory).

%%time
emissions_tcspc = em.Emissions(frame_time="10us", seed=rng, bandpass=(600, 800))
lifetimes_DA, lifetimes_D, lifetimes_all, simulation_object = emissions_tcspc.tcspc(
    transition_set=transition_set,
    number_pulses=1e5,
    pulse_duration=1e-11,
    time_between_pulses=1e-7,
    excitation_rates={"cy5_dna": 1e11},
    size=1e5,
    store_time_points=True,
    # details = True
)

the last frame (of index 0.01) has 1.00e+00 times the pulses of other frames.

CPU times: user 3.94 s, sys: 71.4 ms, total: 4.01 s
Wall time: 3.94 s

emissions_tcspc.plot_time_series()
fi.universal_figure(
    type_="hist", data=lifetimes_all, ylabel="PD", density=True, xlabel="Lifetime (s)"
);

Note:

Since the TCSPC simulation assumes a laser pulse to be instantaneous, a laser pulse that excited multiple fluorophores will lead to duplicates in the time series. This is handled by introducing minimal amount of spaces between these duplicates. The excitations are added in the correct order (fluorophore-specific) such that these series exactly represent the transition chain, even with multiple different distances between energy transfer partners. Note that when monitoring the fluorescence lifetime via the simulation object it will be the true one (as opposed to the generally returned lifetimes, which are the observed ones).

Fluorescence correlation spectroscopy

Observed fluorescence emission events can be analyzed by a correlation analysis.

fcs = fcs_p.FCS(emissions)
list(vars(fcs).keys())

['emissions', 'autocorrelation', 'tau']

Autocorrelation of time points

The autocorrelation can be computed on the exact time points (being limited by the smallest available time difference).

%%time
fcs.autocorrelate_time_points(
    exp_min=-16, exp_max=0, points_per_base=4, base=4, normalize=True
)

The exp_max 0 yields a base to the power of exp_max 1 that is larger than the last time point 0.9998935849738816. Therefore, exp_max is adjusted to -1.

CPU times: user 3.5 s, sys: 362 ms, total: 3.86 s
Wall time: 3.86 s

<fluopy.fcs.FCS at 0x7827f3182e40>

mi.print_class(fcs)
fcs.plot(normalize_to=None, unit="s", scale=1);

Attributes of <fluopy.fcs.FCS object at 0x7827f3182e40>:
.................................................................
emissions = <fluopy.emissions.Emissions object at 0x7827f3182900>
_________________________________________________________________
autocorrelation = array([0.68531317, 0.48458959, 0.68531317, 1.230... 1.0166399 , 1.02127411, 1.03142128, 1.05049859])
_________________________________________________________________
tau = array([2.81051449e-10, 3.97466771e-10, 5.6210289... 1.06694174e-01, 1.50888348e-01, 2.13388348e-01])
_________________________________________________________________

Autocorrelation of time series

The autocorrelation can be computed on the binned time series (being limited by the time bin).

There is not much to be seen in this example, since the bin time is larger than all relevant fluctuation time scales.

fcs.autocorrelate_time_series(log=True, m=4, normalize=True)

<fluopy.fcs.FCS at 0x7827f3182e40>

mi.print_class(fcs)
fcs.plot(normalize_to=None, unit="s", scale=1);

Attributes of <fluopy.fcs.FCS object at 0x7827f3182e40>:
.................................................................
emissions = <fluopy.emissions.Emissions object at 0x7827f3182900>
_________________________________________________________________
autocorrelation = array([1.00211007, 1.00070264, 1.0035714 , 0.997... 0.99998845, 0.99884571, 0.99879731, 1.00040408])
_________________________________________________________________
tau = array([0.005, 0.01 , 0.015, 0.02 , 0.03 , 0.04 ,...0.12 ,
       0.16 , 0.24 , 0.32 , 0.48 , 0.64 ])
_________________________________________________________________

Some fcs fits are available:

# fcs_predict = fcs_p.fit_dark(tau, dark_lifetime, dark_occupation)
# fcs_predict = fcs_p.fit_antibunching(tau, excitation_rate, s1_lifetime)
# fcs_predict = fcs_p.fit_triplet_cis(tau, k_isc, k_T, k_01, k_10, k_iso, k_biso_eff)

Antibunching

Alternatively, you can focus on fast time scales in a linear scale and observe photon bunching or antibunching.

# sensible to tau_max and bin_width, see coincidence notebook
hist, bins = fcs_p.coincidence(
    emissions.event_time_points[: int(2e5)], tau_max=1e-7, bin_width=1e-9, seed=rng
)
fi.universal_figure(
    type_="line",
    data=[bins, hist],
    xlabel=r"$\tau$ (s)",
    ylabel=r"$g^{(2)}(\tau)$",
    scale=1,
);

Blinking

For analyzing a fluorescence on/off process, the blinking time scales can be estimated by thresholding an emission time series.

An ON-period is a number of consecutive frames where each frame contains a minimum amount of emissions (> threshold).

An OFF-period is a number of consecutive frames where each frame contains a maximum amount of emissions (≤ threshold).

Each ON-period is followed by an OFF-period and vice versa.

Emissions from a short simulation

We limit the dataset to 2000 frames for illustration purposes.

%%time
emissions = em.Emissions(frame_time="10us", seed=rng, bandpass=None)
emissions.simulate(transition_set=transition_set, store_time_points=False, frames=2000)
emissions

CPU times: user 863 ms, sys: 932 μs, total: 864 ms
Wall time: 864 ms

<fluopy.emissions.Emissions at 0x7827f0e6dd10>

threshold: int = 5

emissions.plot_time_series(scale=1)
plt.hlines(threshold, xmin=0, xmax=0.02)

<matplotlib.collections.LineCollection at 0x7827eeb96900>

blinks = bl.Blinking(emissions, threshold=threshold)
blinks

<fluopy.blinking.Blinking at 0x7827eeb96a50>

mi.print_class(blinks)

Attributes of <fluopy.blinking.Blinking object at 0x7827eeb96a50>:
.................................................................
emissions = <fluopy.emissions.Emissions object at 0x7827f0e6dd10>
_________________________________________________________________
on_periods = array([ 2,  2,  2,  1,  5,  1,  3,  1,  2,  3,  ...,  1,  2,  3,
        6,  5,  2,  7,  6, 10,  2])
_________________________________________________________________
off_periods = array([ 2,  3,  2,  8,  2, 12, 14, 18,  1, 10,  ...,  5,  1,  1,  1,
        3,  1,  3,  7,  3,  1])
_________________________________________________________________
on_periods_frames = array([   1,    5,   10,   14,   23,   30,   43,...1944, 1953, 1959, 1964, 1978, 1987,
       1998])
_________________________________________________________________
off_periods_frames = array([   3,    7,   12,   15,   28,   31,   46,... 1939, 1943, 1950, 1958, 1961, 1971, 1984, 1997])
_________________________________________________________________

# plot a histogram of OFF times
blinks.plot(
    mode="off_histogram", density=True, display_mean=True, as_time="ms", scale=0.5
)

# plot a histogram of ON times
blinks.plot(
    mode="on_histogram", density=True, display_mean=True, as_time="ms", scale=0.5
)

# plot a time series of OFF times
blinks.plot(mode="off_frame_series", scale=0.5)

# plot a time series of ON times
blinks.plot(mode="on_frame_series", scale=0.5);

You can extract the analytical OFF statistics from the emission time series without differentiating between fluorophores.

on_off_times_analytic, on_off_values_analytic = bl.get_analytical_off_statistics(
    off_frames=blinks.off_periods_frames,
    off_periods=blinks.off_periods,
    on_frames=blinks.on_periods_frames,
    frame_time=blinks.emissions.parameters["frame_time"],
)

bl.plot_off_statistics(
    on_off_times_analytic, on_off_values_analytic, scale=1, title="analytical OFF"
);

Emissions from the long simulation

Get more detailed information from a complete simulation:

%%time
emissions = em.Emissions(frame_time="200ns", seed=rng, bandpass=None)
emissions.extract(simulation=simulation)
emissions

CPU times: user 1.13 s, sys: 106 ms, total: 1.23 s
Wall time: 1.23 s

<fluopy.emissions.Emissions at 0x7827f30995b0>

blinks = bl.Blinking(emissions, threshold=threshold)
blinks

<fluopy.blinking.Blinking at 0x7827ec678e10>

mi.print_class(blinks)

Attributes of <fluopy.blinking.Blinking object at 0x7827ec678e10>:
.................................................................
emissions = <fluopy.emissions.Emissions object at 0x7827f30995b0>
_________________________________________________________________
on_periods = array([1, 1])
_________________________________________________________________
off_periods = array([1920216])
_________________________________________________________________
on_periods_frames = array([1462038, 3382255])
_________________________________________________________________
off_periods_frames = array([1462039])
_________________________________________________________________

# plot a histogram of OFF times
blinks.plot(
    mode="off_histogram", density=True, display_mean=True, as_time="ms", scale=0.5
)

# plot a histogram of ON times
blinks.plot(
    mode="on_histogram", density=True, display_mean=True, as_time="ms", scale=0.5
)

# plot a time series of OFF times
blinks.plot(mode="off_frame_series", scale=0.5)

# plot a time series of ON times
blinks.plot(mode="on_frame_series", scale=0.5)

array([[<Axes: xlabel='identity', ylabel='consecutive ON frames'>]],
      dtype=object)