|
| 1 | +""" |
| 2 | +Flow Cytometry Simulation: Full System Example with Workflow |
| 3 | +============================================================ |
| 4 | +
|
| 5 | +This tutorial demonstrates a complete flow cytometry simulation using the FlowCyPy library. |
| 6 | +It models fluidics, optics, signal processing, and classification of multiple particle populations. |
| 7 | +
|
| 8 | +Steps Covered: |
| 9 | +-------------- |
| 10 | +1. Configure simulation parameters and noise models |
| 11 | +2. Define laser source, flow cell geometry, and fluidics |
| 12 | +3. Add synthetic particle populations |
| 13 | +4. Set up detectors, amplifier, and digitizer |
| 14 | +5. Simulate analog and digital signal acquisition |
| 15 | +6. Apply triggering and peak detection |
| 16 | +7. Classify particle events based on peak features |
| 17 | +""" |
| 18 | + |
| 19 | +from FlowCyPy.workflow import ( |
| 20 | + ureg, |
| 21 | + SimulationSettings, |
| 22 | + Workflow, |
| 23 | + Detector, |
| 24 | + circuits, |
| 25 | + peak_locator, |
| 26 | + triggering_system, |
| 27 | + distributions, |
| 28 | + population, |
| 29 | + GammaModel, |
| 30 | + classifiers, |
| 31 | +) |
| 32 | + |
| 33 | +SimulationSettings.include_noises = False |
| 34 | +SimulationSettings.include_shot_noise = True |
| 35 | +SimulationSettings.include_dark_current_noise = True |
| 36 | +SimulationSettings.include_source_noise = True |
| 37 | +SimulationSettings.include_amplifier_noise = True |
| 38 | +SimulationSettings.assume_perfect_hydrodynamic_focusing = True |
| 39 | +SimulationSettings.population_cutoff_bypass = False |
| 40 | + |
| 41 | +population_0 = population.Sphere( |
| 42 | + name="Pop 0", |
| 43 | + medium_refractive_index=distributions.Delta(1.33 * ureg.RIU), |
| 44 | + concentration=5e10 * ureg.particle / ureg.milliliter, |
| 45 | + diameter=distributions.RosinRammler( |
| 46 | + shape=150 * ureg.nanometer, |
| 47 | + scale=50 * ureg.nanometer, |
| 48 | + low_cutoff=50.0 * ureg.nanometer, |
| 49 | + ), |
| 50 | + refractive_index=distributions.Normal( |
| 51 | + mean=1.44 * ureg.RIU, |
| 52 | + standard_deviation=0.002 * ureg.RIU, |
| 53 | + low_cutoff=1.33 * ureg.RIU, |
| 54 | + ), |
| 55 | +) |
| 56 | + |
| 57 | +population_1 = population.Sphere( |
| 58 | + name="Pop 1", |
| 59 | + medium_refractive_index=distributions.Delta(1.33 * ureg.RIU), |
| 60 | + concentration=5e17 * ureg.particle / ureg.milliliter, |
| 61 | + diameter=distributions.RosinRammler( |
| 62 | + shape=50 * ureg.nanometer, |
| 63 | + scale=50 * ureg.nanometer, |
| 64 | + ), |
| 65 | + refractive_index=distributions.Normal( |
| 66 | + mean=1.44 * ureg.RIU, |
| 67 | + standard_deviation=0.002 * ureg.RIU, |
| 68 | + low_cutoff=1.33 * ureg.RIU, |
| 69 | + ), |
| 70 | + sampling_method=GammaModel(mc_samples=10_000), |
| 71 | +) |
| 72 | + |
| 73 | + |
| 74 | +detector_0 = Detector( |
| 75 | + name="side", |
| 76 | + phi_angle=90 * ureg.degree, |
| 77 | + numerical_aperture=0.3 * ureg.AU, |
| 78 | + responsivity=1 * ureg.ampere / ureg.watt, |
| 79 | +) |
| 80 | + |
| 81 | +detector_1 = Detector( |
| 82 | + name="forward", |
| 83 | + phi_angle=0 * ureg.degree, |
| 84 | + numerical_aperture=0.3 * ureg.AU, |
| 85 | + responsivity=1 * ureg.ampere / ureg.watt, |
| 86 | +) |
| 87 | + |
| 88 | +discriminator = triggering_system.DynamicWindow( |
| 89 | + trigger_detector_name="forward", |
| 90 | + threshold="2sigma", |
| 91 | + pre_buffer=20, |
| 92 | + post_buffer=20, |
| 93 | +) |
| 94 | + |
| 95 | +peak_locator = peak_locator.GlobalPeakLocator(compute_width=False) |
| 96 | + |
| 97 | +analog_processing = [ |
| 98 | + circuits.BaselineRestorator(window_size=10 * ureg.microsecond), |
| 99 | + circuits.BesselLowPass(cutoff=2 * ureg.megahertz, order=4, gain=2), |
| 100 | +] |
| 101 | + |
| 102 | +workflow = Workflow( |
| 103 | + wavelength=405 * ureg.nanometer, |
| 104 | + source_numerical_aperture=0.1 * ureg.AU, |
| 105 | + optical_power=200 * ureg.milliwatt, |
| 106 | + sample_volume_flow=80 * ureg.microliter / ureg.minute, |
| 107 | + sheath_volume_flow=1 * ureg.milliliter / ureg.minute, |
| 108 | + width=200 * ureg.micrometer, |
| 109 | + height=100 * ureg.micrometer, |
| 110 | + populations=[population_0, population_1], |
| 111 | + gain=10 * ureg.volt / ureg.ampere, |
| 112 | + bandwidth=10 * ureg.megahertz, |
| 113 | + bit_depth="14bit", |
| 114 | + sampling_rate=60 * ureg.megahertz, |
| 115 | + saturation_levels="auto", |
| 116 | + background_power=0.001 * ureg.milliwatt, |
| 117 | + detectors=[detector_0, detector_1], |
| 118 | + analog_processing=analog_processing, |
| 119 | + trigger=discriminator, |
| 120 | + peak_locator=peak_locator, |
| 121 | +) |
| 122 | + |
| 123 | +workflow.initialize() |
| 124 | + |
| 125 | +run_record = workflow.run(run_time=1 * ureg.millisecond) |
| 126 | + |
| 127 | +_ = run_record.event_collection.plot(x="Diameter") |
| 128 | + |
| 129 | +# %% |
| 130 | +# Step 5: Plot Events and Raw Analog Signals |
| 131 | +# ------------------------------------------ |
| 132 | +_ = run_record.event_collection.plot(x="forward") |
| 133 | + |
| 134 | + |
| 135 | +# %% |
| 136 | +# Plot raw analog signals |
| 137 | +# ----------------------- |
| 138 | +_ = run_record.plot_analog(figure_size=(12, 8)) |
| 139 | + |
| 140 | + |
| 141 | +# %% |
| 142 | +# Step 6: Plot Triggered Analog Segments |
| 143 | +# -------------------------------------- |
| 144 | +_ = run_record.plot_digital(figure_size=(12, 8)) |
| 145 | + |
| 146 | + |
| 147 | +# %% |
| 148 | +# Step 7: Plot Peak Features |
| 149 | +# -------------------------- |
| 150 | +_ = run_record.peaks.plot(x=("forward", "Height")) |
| 151 | + |
| 152 | + |
| 153 | +# %% |
| 154 | +# Step 8: Classify Events from Peak Features |
| 155 | +# ------------------------------------------ |
| 156 | +classifier = classifiers.KmeansClassifier(number_of_clusters=2) |
| 157 | + |
| 158 | +classified = classifier.run( |
| 159 | + dataframe=run_record.peaks.unstack("Detector"), |
| 160 | + features=["Height"], |
| 161 | + detectors=["side", "forward"], |
| 162 | +) |
| 163 | + |
| 164 | +_ = classified.plot(x=("side", "Height"), y=("forward", "Height")) |
0 commit comments