# MUFASA: A Continuous-Time Stochastic Framework for Realistic Fluorescence Microscopy Simulation

Wessim Omezzine, Sébastien Schaub, Laure Blanc-Féraud, Luca Calatroni

► To cite this version:

Wessim Omezzine, Sébastien Schaub, Laure Blanc-Féraud, Luca Calatroni. MUFASA: A Continuous-Time Stochastic Framework for Realistic Fluorescence Microscopy Simulation. 2025. hal-05112069v2

HAL Id: hal-05112069

https://hal.science/hal-05112069v2

Preprint submitted on 30 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# MUFASA: A Continuous-Time Stochastic Framework for Realistic Fluorescence Microscopy Simulation

Wessim Omezzine¹, Sébastien Schaub², Laure Blanc-Féraud¹, and Luca Calatroni³

¹Université Côte d'Azur, CNRS, I3S

² Sorbonne University, CNRS, LBDV

³MaLGA Centre, DIBRIS, University of Genoa and MMS, Italian Institute of Technology, Genoa, Italy

*Correspondence: omezzine@i3s.unice.fr, sebastien.schaub@imev-mer.fr, blancf@i3s.unice.fr, luca.calatroni@unige.it

## Abstract

We present MUFASA (Multi-Protocol Unified Fluorescence-based Advanced Simulation Algorithm), a physically grounded, continuous-time simulator for super-resolution fluorescence microscopy. By modeling fluorophore dynamics using continuous-time Markov chains, MUFASA’s simulation features yield realistic photon emission behavior across both Single Molecule Localization Microscopy (SMLM) and fluorescence fluctuation-based (FF-SRM) protocols—independently of frame duration and sampling. The framework supports both individual emitters and structure-level simulations, incorporating photophysical transitions, photobleaching, and camera properties.

To quantitatively validate simulations with real data, we introduce a novel validation metric based on the 1-Wasserstein distance between simulated and experimental photon-count distributions. In addition to simulation, another functionality estimates key photophysical parameters (e.g., molar extinction coefficient) and to suggest optimal light-source power ranges from fluctuation data. An intuitive Python-based graphical interface enables real-time parameter tuning, visualization, and TIFF export. Designed for biologists, physicists, microscopists, and numerical imaging engineers, MUFASA offers a practical platform for microscopy experiment design, hypothesis testing and the generation of realistic training data for data-driven microscopy methods across modalities.

## Introduction

Modern Super-resolution Microscopy (SRM) methods aim to overcome the diffraction limit by exploiting the photokinetics of the molecule and their intricate behaviors often represented in terms of electronic state transitions in Jablonski diagram (Figure 1).

As outlined in [1], SRM can be grouped into three main families. Early advancements were made with STimulated-Emission-Depletion (STED) [2], which improves resolution by selectively deactivating fluorophores at the periphery of the excitation spot using a depletion laser, thereby reducing the effective point spread function. Structured Illumination Microscopy (SIM) [3] enhances resolution by pro-
jecting known light patterns and computationally reconstructing high-frequency information beyond the diffraction limit. Subsequently, Single Molecule Localization Microscopy (SMLM) [4] emerged as a prominent family of methods, enabling super-resolution imaging through the precise localization of photoactivatable *blinking* fluorophores that stochastically switch between fluorescent (on) and non-fluorescent (off) regimes—a behavior commonly referred to as *fluorescence blinking*. By temporally isolating the emission of individual molecules, SMLM approaches such as STochastic Optical Reconstruc-tion Microscopy (StORM) [5] and PhotoActivated Localization Microscopy (PALM) [6] can localize each emitter and reconstruct high-resolution images from these localizations. These methods typically operate under high light power conditions (up to several kW cm⁻²) to drive fluorophores into long-lived dark states or induce controlled photobleaching, thereby temporally separating emitters.

Building on this perspective, Fluorescence Fluctuation-based SRM (FF-SRM) techniques [7] have been developed, which exploit the temporal intensity fluctuations of fluorophores to achieve resolution enhancement. Unlike SMLM, which isolates single emitters in time, FF-SRM analyzes the collective temporal behavior of multiple emitters within each diffraction-limited pixel. These methods usually operate under low laser power conditions (e.g., 0.1 W cm⁻²), which preserve fluorophore activity and lead to mean continuous photon emission. These fluctuations are often modeled using Poisson statistics, under the assumption that photon emission is a memoryless process occurring with constant rate over short intervals [8]. Higher-order statistical moments (e.g., cumulants of order $N \ge 2$) are then computed to suppress noise and enhance spatial detail. By leveraging these temporal dependencies, FF-SRM reconstructs super-resolved images without requiring explicit localization of individual molecules.

However, SRM still faces challenges in both re-producibility and interpretability. Imaging results are highly sensitive to experimental parameters and the fluorophores photophysical properties, which are often poorly characterized. For instance, small changes in laser intensity can drastically alter blinking behavior, leading to inconsistent reconstructions [9]. Misinterpretation of photophysical artifacts—such as fluorophore reactivation or triplet state transitions—can result in false-positive localizations or incorrect molecule counts [10, 11]. Physically grounded simulation frameworks would be an excellent tool not only to guide experimental design, enable protocol comparison, support method validation, and generate synthetic training data for data-driven methods. Physically grounded simulators can produce realistic fluorescence signals with known ground truth, estimating key photophysical effects such as blinking kinetics, bleaching, and noise. This enables training models that generalize better across imaging conditions while being benchmarked rigorously using synthetic validation datasets.

Multiple efforts have been made to simulate photophysics. Early tools such as SuReSim [12] and FluoSim [13] relied on simplified photoswitching models, primarily focusing on “on” and “off” regimes to simulate StORM and PALM dynamics, implying a constant photon flux during “on” regimes. The switching behavior is governed by *transition rates*—that is, the probabilities per unit time of transitioning between states. In these simulators, these transition rates were user-defined and not derived from first-principles or measured experimental parameters, limiting their physical interpretability. More recent tools like VirtualSMLM [14] and SMIS [15] have better physical consistency by linking transition rates to experimental parameters. However, these simulators use fixed discrete time grids, which are intrinsically unable to describe rapid transitions, i.e occurring at sub-frame timescales.

The SOFI Simulation tool [16] addresses this lim-
<figure><img src="image_2.png" /><figcaption>**Figure 1:** Fluorophore state transitions in fluorescence microscopy. The Jablonski diagram represents quantum states involved in photon absorption, emission, intersystem crossing, and photobleaching. Vertical solid arrows indicate radiative and non-radiative transitions between the ground singlet state ($S_0$), the first excited singlet state ($S_1$), and the triplet state ($T_1$). Wavy and dotted lines represent non-radiative transitions. Additional transitions from a non-activated state (NA) and to a bleached state (B) are included to reflect specific behaviors observed in methods like PALM. This description of the state space forms the basis for the continuous-time Markov chain modeling used in MUFASA.</figcaption></figure>

itation for fluctuation-based models by employing
Continuous-Time Markov Chains (CTMCs, rather
than modeling fluorophore dynamics as a sequence
of updates on a fixed temporal grid. This enables
a more faithful representation of blinking dynamics,
which occur randomly but with well-defined rates,
independent of any external clock.

However, even this approach remains limited in
scope. The SOFI Simulation tool applies the CTMC
only to the switching between observable "on" and
"off" regimes. Within the "on" regime itself, photon
emission is still modeled as a constant-rate process,
thereby ignoring the rapid quantum jumps between
excited $S_1$ and ground $S_0$ singlet states that govern
the actual photon emission events. Capturing these
sub-frame dynamics requires extending the CTMC
framework to model the full set of quantum state
transitions involved in fluorescence.

To address this gap and encompass different
fluorescence dynamics, we introduce the Multi-

Protocol **Unified Fluorescence-based Advanced** 118
**Simulation Architecture (MUFASA)**, a compre- 119
hensive, physics-inspired simulator that integrates 120
both SMLM and FF-SRM protocols within a single 121
continuous-time framework. Unlike previous tools 122
that model either discrete switching events or sim- 123
plified emission processes, MUFASA simulates quan- 124
tum state transitions directly using CTMCs. This 125
enables the accurate reproduction of both intermit- 126
tent and continuous emission behaviors. MUFASA 127
comprises two core features: 128

* **MUFASA-Sim**, which performs protocol-aware simulations of photon emission based on quantum photophysics, including realistic transition times and photobleaching effects. 129 130 131 132

* **MUFASA-Design**, which extracts relevant photophysical parameters—such as the molar extinction coefficient—from experimental fluctuation data, and provides optimal illumination 133 134 135 136
power range for imaging.

MUFASA incorporates the detailed photophysical dynamics described by the Jablonski diagram, intp a comprehensive simulator, enabling the modeling of rapid transitions between ground and excited states. In addition to its simulation capabilities, MUFASA includes a graphical user interface (GUI) for interactive exploration and a quantitative validation framework for assessing the realism of generated data. Specifically, we introduce the use of the 1-Wasserstein distance—a statistical metric that quantifies the difference between probability distributions based on the minimal transport cost needed to transform one distribution into another. This enables rigorous comparison between simulated photon emission statistics and experimental single-molecule traces.

Designed for both flexibility and realism, MUFASA serves as a valuable tool for a large community of researchers involved in the design of SRM experiments, from experimental sample preparation to numerical simulation, such as hypothesis testing, or generating physically grounded training data and realistic simulations for data-driven SRM methods.

## Results

An overview of the MUFASA simulation framework and its core components is illustrated in Figure 2, which provides visual context for the analyses presented in the following subsections.

Under the SMLM-type dynamics, we simulate isolated fluorophores to examine protocol-specific emission statistics at high temporal resolution (Figure 2b). We qualitatively analyze the distribution of emitted photons over time, allowing us to study how varying experimental conditions influence molecular dynamics. These traces reveal protocol-dependent behaviors: continuous emission for FF, zero-inflated blinking in StORM, and sparse activation in PALM.

We extended the framework to model the stochastic

fluorescence dynamics of multiple independent fluorophores within a biological structure (Figure 2e). This extension builds on the single-molecule simulation pipeline (Figure 2a), but now requires additional input: the spatial coordinates of all fluorophores in the structure. Each fluorophore is then treated as an independent emitter, governed by the same photophysical parameters—such as state transition rates and excitation conditions—as in the single-molecule case. After simulating time dynamics, MUFASA models the optical system and the camera detection process (Figure 2c) using a point spread function (PSF) convolution followed by spatial downsampling and noise modeling. The PSF accounts for diffraction-limited optical blur, while the camera model incorporates detector-specific characteristics such as pixel size, shot noise and gain fluctuations, producing realistic low-resolution image frames. When applied to simulations of microtubule structures labeled with 6,498 fluorophores, MUFASA successfully reproduces contrast patterns and spatial density variations characteristic of each SRM modality, for instance closely mirroring observations from experimental SMLM acquisitions.

The MUFASA-Design feature complements MUFASA-Sim by estimating key photophysical parameters (e.g., molar extinction coefficients) from experimental data and simulating the optimal illumination power range, thereby streamlining experimental planning and optimization. Using enhanced green fluorescent protein (EGFP) and Alexa Fluor 647 as examples (see Supplementary Note 7), we show that the estimated parameters deviate by less than 10% from reference values.

To summarize, these results are validated across four SRM modalities that exploit temporal dynamics—FF, Blinking, StORM, and PALM—using both qualitative and quantitative benchmarks against experimental data. Further details are provided below.
<figure><img src="image_3.png" /><figcaption>**Figure 2:** The MUFASA simulation pipeline for photophysics-aware fluorescence microscopy imaging. (a) A simulated biological structure is populated with discrete fluorophore labels to serve as the input to the MUFASA framework. (b) Fluorophore positions are discretized onto a regular 2D grid of *N* sampling sites. (c) Depending on the protocol, MUFASA generates *N* high-temporal-resolution photon-emission time traces by simulating its photophysical state transitions (excitation, emission, dark states) using a continuous-time Markov chain (CTMC) model and records the resulting per-molecule photon-count time traces. (d) Photon counts from all fluorophores within each grid cell are summed over all discrete time intervals to yield the total fluorescence contribution per pixel—i.e. the raw, per-frame emission signal. (e) Each raw frame is processed through a simulation of the microscope's optical response: first, the high-resolution photon map is convolved with the point spread function (PSF) to model diffraction blur. Then, spatial undersampling is applied by summing photon counts within blocks of high-resolution pixels to match the coarser detector grid. Finally, realistic camera noise is added (e.g., shot noise, EMCCD amplification), yielding a noisy, realistic resolution grayscale image stack. (f) The final output is a time-lapse series of low-resolution images that faithfully conveys both spatial heterogeneity and the temporal photophysics inherent to specific modalities. (g) A zoomed-in camera frame shows two selected fluorophores within a defined region of interest (ROI), marked by red and violet circles. For each fluorophore, the central pixel's intensity is tracked over time under the simulated acquisition sequence, revealing temporal fluorescence traces (middle: red; right: violet). These traces reflect differences in photoswitching behavior and emission variability captured by MUFASA.</figcaption></figure>

## Simulating and approximating single-molecule traces for different SRM protocols

Using MUFASA-Sim, we analyzed the continuous-time behavior of an isolated Alexa Fluor 647 molecule over 200 frames (10 ms per frame). For each frame, the number of emitted photons corresponds to the number of absorption–emission events, occurring within 1 – 2 nanoseconds as shown in Figure 3. This emission pattern directly results from

the underlying quantum state dynamics illustrated in the same figure and formalized through transition rates in the CTMC model, as described in the Methods section on continuous-time dynamics.

We simulated Alexa Fluor 647 fluorescence dynamics under four protocols: FF, Blinking (200 frames 10ms per frame), StORM and PALM (1000 frames, 10 ms per frame), using photophysical rates from the Thermo Fisher Scientific datasheet [17]. We included all states depicted in the Jablonski diagram (figure 1). Transitions between quantum states

225

226

227

228

229

230

231

232

233

234

235
<figure><img src="image_4.png" /><figcaption>Figure 3: Comparative simulation of photophysical state transitions, single-molecule photon emission traces, and photon-count distributions across four fluorescence microscopy techniques: Fluorescence Fluctuation (FF), Blinking, Stochastic Optical Reconstruction Microscopy (StORM), and Photoactivated Localization Microscopy (PALM). The sampling frequency is fixed at 10 ms per frame for all protocols.</figcaption></figure>

are governed by rate constants $k_{ij}$, representing the probability per second of transitioning from state *i* to state *j*. To simplify the model without sacrificing accuracy, transitions with rates $k_{ij} < 10^{-4} s^{-1}$ were pruned. Protocol differences derive solely from variations in experimental conditions like illumination

power (0.1 vs $5kW/cm^2$), intersystem crossing (ISC) rates, which modulate the probability of transitioning from the excited singlet state to the long-lived triplet state.

MUFASA-Sim is built on continuous-time Markov chains rather than discrete updates; hence, it faith-

242

243

244

245

246

247
<figure><img src="image_5.png" /><figcaption>**Figure 4:** Sub-frame photophysical transitions revealed by MUFASA's continuous-time CTMC model. Zoom in on single-molecule state trajectories under FF (left) and StORM (right) conditions, overlaid with the 10 ms frame boundaries (dashed line). MUFASA generates exponential holding time for states capturing rapid transitions, independently of frame duration.</figcaption></figure>

fully captures rapid sub-frame transitions, where the lifetime of each state is randomly drawn from an exponential distribution, as shown in Figure 4. These transitions would not be captured by a frame-based update scheme.

Under the FF protocol, where the light power is relatively low (0.1 W/cm²), the molecule exhibits continuous emission, rapidly alternating between the ground and excited singlet states. For Alexa Fluor 647, the mean lifetime in the excited state ($\tau_1$ = 1 ns) is much shorter than the typical holding time in the ground state, which can reach up to 1 ms. The per-frame photon emission trace fluctuates around a mean of approximately 300 photons. A statistical analysis of the distribution of emitted photon counts reveals that it can be well-approximated by a Poisson distribution. Theoretical results derived in Section 2 show that, under certain conditions, the emission process indeed follows a Poisson distribution. More precisely, using Equation 1, one can compute $k_{01}$ and find that the theoretical Poisson rate $\mu_{theoretical}$ matches the observed mean photon count.

The blinking and StORM protocol share similar photophysical mechanisms: both increase the probability of intersystem crossing (ISC) and extend the time spent in the triplet state. However, they differ significantly in terms of excitation power and the

likelihood of photobleaching. StORM typically uses a higher laser intensity (5 kW/cm²), leading to a greater number of emitted photons during each activation event, compared to blinking (1 kW/cm²). In both cases, the molecule alternates between dark and fluorescent regimes, producing emission traces with prolonged zero-photon intervals during the off state.

PALM, in contrast, relies on the random activation of molecules from a non-activated (NA) state. Activation from this NA state occurs with a probability $p_{act}$ per second. Once activated, Alexa Fluor 647 undergoes multiple absorption–emission cycles, typically within one or two frames (~15 ms) before bleaching. Over 1000 frames (10 ms per frame), the zero-photon count fraction is 99.2%, with a heavy tail consisting of rare 1–2 frame bursts of high photon emission. In large experimental datasets, this pattern persists over tens or even hundreds of thousands of frames, with fluorescent spots generally lasting no more than 1–3 frames.

Across Blinking, StORM, and PALM protocols, the emission traces reveal features such as activation events, long triplet-state dwell times, and sparse burst-like behavior. These characteristics lead to deviations from the standard Poisson assumption, often producing zero-inflated photon count distributions. More appropriate models, such as Poisson–Gamma (negative binomial) mixtures, may better capture the increased variance introduced by blinking and state switching. However, rather than relying on protocol-specific statistical models, MUFASA addresses this variability through a unified framework based on continuous-time Markov chains (CTMCs), enabling realistic simulation of photophysical state transitions across all SRM protocols.

## MUFASA-Sim: Biological structure simulation

We extend the use of MUFASA-Sim beyond single-molecule modeling to biological structures for a given
<figure><img src="image_6.png" /><figcaption>Figure 5: Protocol-specific simulated acquisitions of microtubule structures using MUFASA-Sim. Representative output frames from simulations of a microtubule bundle under four distinct SRM protocols: PALM, StORM, FF, and Blinking. All pixel intensities are expressed in analog-to-digital units (ADU), as output by the simulated camera. Although the same fluorophore model is used across all conditions, protocol-specific behaviors are recapitulated through changes in excitation regime, transition dynamics, and detection noise. PALM and StORM produce sparse, high-intensity single-molecule events due to activation and rapid photobleaching. FF imaging yields dense, continuous emission with spatially structured signal, while Blinking exhibits temporally modulated sparsity reflecting stochastic transitions between emissive and dark states.</figcaption></figure>

spatial pattern.

As a test, we used a synthetic microtubule derived from fluorophore coordinate data available in [18] to generate protocol-specific image sequences. The geometry, illustrated in Supplementary Figure 4, comprises a dense distribution of 6,498 fluorophores arranged over a 2048 × 2048 grid. Local molecular densities vary between 0 and 4 emitters per pixel, mimicking experimentally observed heterogeneity.

To capture realistic temporal behavior, high-

resolution photon emission traces were generated by simulating the continuous-time dynamics of fluorophores using protocol-specific CTMC models. These emission maps were then blurred by a simplified Gaussian point spread function (PSF)—noting that more advanced PSF models can be incorporated, as described in [19], downsampled by a factor $L = 4$, and passed through a noise-aware camera model to produce digitally quantized output. The noise model includes Poisson-distributed shot

325

326

327

328

329

330

331

332

333

334
noise, additive Gaussian readout noise, and EM gain fluctuations according to models established in [20–22]. Acquisition parameters, including optical and camera properties, are summarized in Supplementary Table 1, with protocol-specific experimental and molecular details provided in Supplementary Tables 2 and 3, respectively.

Protocol-specific simulations are shown in Figure 5 and Supplementary Figures 5–9, each demonstrating distinct behavior. All simulations used the same fluorophore model; environmental parameters such as excitation intensity and non-radiative decay rates were adjusted to reflect the operational principles of each SRM modality. The parameters used for the camera models are detailed in Supplementary Note 5.

These qualitative visual results confirm that MUFASA-Sim faithfully reproduces the protocol-specific behaviors of different SRM techniques. Quantitative validation against experimental data is presented in the Validation Section.

## MUFASA-Design: Estimation of molar extinction coefficient $\varepsilon$

The molar extinction coefficient $\varepsilon$ quantifies how efficiently a molecule absorbs light at a given wavelength. This parameter is important for biologists, as it accurately quantifies concentration [23]. For physicists, it plays a central role in modeling light-matter interactions in nanophotonics and quantum optics [24–26].

However, determining $\varepsilon$ is challenging because it is influenced by the fluorophore quantum yield—the probability that an absorbed photon results in the emission of a new photon—which itself depends on environmental factors such as solvent polarity, pH, and temperature. Variations in quantum yield can lead to significant deviations in the observed absorbance, complicating accurate estimation of $\varepsilon$ from experimental data [23].

MUFASA-Design enables the estimation of $\varepsilon$ from

fluorescence fluctuation data by inverting a physically grounded expression for the excitation rate $k_{01}$ given by the Equation 1:

$$ \varepsilon = \frac{\hat{k}_{01} \mathcal{N}_A h c}{P \lambda \cdot 10^{-6} \ln(10)} $$

Here, $\hat{k}_{01}$ is the excitation rate estimated from fluorescence time traces after denoising the experimental signal (as detailed in Supplementary Note 3), $\mathcal{N}_A$ is Avogadro’s number, $h$ is Planck’s constant, $c$ is the speed of light, $P$ is the excitation power density, and $\lambda$ is the excitation wavelength.

<table><caption>**Table 1:** Comparison of true and MUFASA-Design estimated molar extinction coefficients.</caption><thead><tr><td>Fluorophore</td><td>True ε</td><td>Estimated ε</td></tr></thead><tbody><tr><td>EGFP</td><td>55,900 [27]</td><td>52.529</td></tr><tr><td>Alexa Fluor 647</td><td>239,000 [17]</td><td>231,260</td></tr></tbody></table>

We validated this approach by estimating its value for two widely used fluorophores (Table 1). In both cases, MUFASA-Design yielded estimates within 10% of the literature values, despite operating on noisy simulated emission data.

## MUFASA-Design: Illumination power range calibration

Proper calibration of laser light is a critical step in fluorescence microscopy, particularly in super-resolution applications where photobleaching can severely limit data quality and experiment duration. In practice, the selection of illumination intensity often relies on trial and error, leading to potential overexposure or inadequate signal [28].

To address this, MUFASA-Design enables users to simulate the photokinetics of fluorophores across a wide range of light intensities. We applied the simulator to model EGFP behavior for excitation intensities ranging from 0.1 to 10 kW/cm² under the FF protocol. For each setting, MUFASA-Design computes the expected photon emission per molecule.

These simulations offer a practical tool for calibrating illumination intensity in SRM acquisitions
by enabling users to predict fluorophore response prior to experimentation. For example, MUFASA-Design simulations of EGFP photodynamics (Figure 6) reveal that at intermediate light powers (e.g., between 0.5 and 1 kW/cm²), one can achieve stable emission exceeding 3000 photons per frame for over 500 frames before photobleaching becomes relevant. This allows experimentalists to select illumination settings that balance brightness and longevity, reducing photodamage and avoiding suboptimal conditions. Rather than relying purely on experimental trial and error, MUFASA thus provides a principled approach to explore and constrain viable illumination power regimes before real sample acquisition.

<figure><img src="image_7.png" /><figcaption>**Figure 6:** Simulation of EGFP photodynamics across illumination powers using MUFASA-Design. Left: Mean cumulative photons that a single EGFP molecule can emit before reaching the bleaching criterion when illuminated continuously with 488-nm light between 0.1 – 2kW/cm². The initial linear rise ≤ 0.3kW/cm² reflects the excitation-limited regime, whereas the shoulder above (≈ 0.7kW/cm²) marks the cycle-limited behavior. Right: Expected frame index at which the bleaching event occurs for the same simulations. Because bleaching is modelled as a fixed maximum number of excitation cycles, increasing light power shortens the time required to reach this limit. This curve allows users to identify illumination regimes that ensure stable emission while minimizing photobleaching.</figcaption></figure>

## Validation

In this section, we report the numerical assessments performed on MUFASA's ability to replicate the stochastic dynamics of fluorophore emission by comparing it with experimental photon emission traces. As a benchmark, we consider SOFI-Simulation Tools

[16], a widely used simulator based on fluorescence fluctuation analysis. This comparison focuses on two key imaging regimes: Blinking and Fluorescence Fluctuation (FF) microscopy. To quantitatively assess the effectiveness of the simulator in generating data with distributions similar to the observed ones, we used first-order Wasserstein distance [29], which captures differences in the overall shape and structure of the photon count distributions between simulated and experimental traces as detailed in Supplementary Note 9.

## Comparative analysis and fluctuation-based statistical evaluation

SOFI Tools assumes a fixed emission rate during each “on” period, resulting in uniform photon plateaus across molecules. This tool is detailed in Supplementary Note 4. In contrast, MUFASA simulates continuous-time quantum state transitions, allowing emission rates to vary dynamically within “on” regimes.

The figures 7a and 7b display the differences in photon emissions simulated for three molecules. We analyze how the total photon count evolves as the molecules bleach over time. For each simulator, we generated three independent photon time traces over 3,000 frames. The first fluorophore was configured to bleach around frame 1,000, the second around frame 2,000, and the third remained stable throughout the entire acquisition. We then summed the emission traces to yield the total photon count over time, in order to highlight how each model handles photon emission dynamics during the “on” regime. The shaded regions represent the amplitude of emission fluctuations within each plateau.

MUFASA exhibits greater variability compared to SOFI Tools, reflecting its explicit modeling of stochastic transitions between the ground state $S_0$ and the singlet excited state $S_1$. This variability is further illustrated in the violin plots (Figures 7c and 7d), which summarize the photon count distribution.
<figure><img src="image_8.png" /><figcaption>**Figure 7:** MUFASA captures sub-state variability in photon emission dynamics. (a–b) Photon emission traces from three fluorophores under the FF protocol, simulated with MUFASA (a) and SOFI Tools (b). Shaded regions indicate intra-plateau variability. MUFASA reveals greater fluctuations, suggesting more informative sub-state dynamics during the on-regime. (c–d) Violin plots of photon counts per emission plateau. MUFASA displays broader distributions, while SOFI Tools produces narrower, unimodal profiles due to its assumption of constant emission rates.</figcaption></figure>

butions across the three plateaus.

MUFASA’s wider distributions capture the intrinsic randomness of fluorophore dynamics, in contrast to SOFI Tools, which assumes a constant emission rate and yields smoother traces with narrower distributions. This assumption neglects sub-state transitions within the ”on” regime. MUFASA incorporates these transitions, allowing the detection of additional states and producing variable, non-uniform photon emission profiles even during the emitting regime.

## Quantitative assessment of StORM dynamics

To validate MUFASA under the StORM protocol, we compared its simulated photon emission traces to experimental data from [30]. For consistency, the analysis focused on the final 5000 frames of acquisition, where single-molecule events dominate due to widespread bleaching. We note that we validate

MUFASA, under the FF protocol, against experimentally measured data in Supplementary Note [10].

We replicated acquisition conditions based on the metadata associated with the StORM dataset [30], including reported parameters such as excitation wavelength (641 nm), illumination power (~1–2 kW/cm²), and camera specifications (Andor iXon EMCCD, 100 nm effective pixel size). Transition rates not explicitly reported—such as intersystem crossing $k_{1T}$ and triplet decay $k_{T0}$—were empirically tuned to match the temporal dynamics observed in the real data. A detailed summary of the parameter inference and matching procedure is provided in Supplementary Note [8].

Beyond qualitative visual agreement illustrated in the left column of Figure 8, the right column of the same figure shows a qualitative comparison between experimental and simulated traces. MUFASA accurately captures the sparse, burst-like emission events characteristic of StORM imaging. The baseline signal and event amplitude are consistent with the experimental profile, demonstrating MUFASA ability to reproduce the intermittent, high-intensity emission typical of isolated fluorophores.

This alignment is confirmed with a Wasserstein distance ($W_1$ = 32.40) between MUFASA simulations and the experimental data, significantly lower than SOFI Tools ($W_1$ = 170.83).

## Discussion

MUFASA brings together simulation, quantitative validation, and experimental design in a unified framework for the description of fluorescence microscopy dynamics. At its core, it is a continuous-time Markov model that drives both single-molecule and structure-level simulations. This foundation not only enables realistic modeling of fluorophore dynamics (MUFASA-Sim), but also illumination powers MUFASA-Design—our dedicated module for estimating photophysical parameters such as the molar extinction coefficient $\varepsilon$.
<figure><img src="image_9.png" /><figcaption>**Figure 8:** Validation of MUFASA under the StORM protocol: On the left: Each row shows a different time frame: left column corresponds to real experimental acquisitions, and right column displays the corresponding MUFASA-Sim outputs under matched conditions. On the right: Comparison of two photon count time traces between experimental StORM data (left, blue) and MUFASA simulations (right, orange). MUFASA reproduces the burst-like, high-intensity events and stable background observed in real single-molecule datasets. Traces shown correspond to the final 2000 frames of acquisition, where spatial and temporal sparsity dominates.</figcaption></figure>

This mathematical underpinning allows us to infer $ε$ from experimental fluctuation data, a task traditionally hampered by environmental sensitivity and the absence of ground-truth rates. By estimating the excitation rate $k_{01}$ from denoised traces, we sidestep the need for bulk measurements and offer a more direct, physically interpretable estimation pipeline.

MUFASA is also built for flexibility across protocols. It reproduces both sparse emission patterns typical of SMLM methods like StORM and PALM, and the high-density, continuous dynamics seen in FF-based techniques. This is made possible by a modular transition graph structure that adapts to each protocol photokinetics—rather than relying on hand-tuned switching assumptions.

A gap in previous simulation frameworks was the lack of rigorous validation: Virtual SMLM [14] relied primarily on visual matching of StORM im-

ages, without statistical comparison to experimental data, while SOFI Simulation Tools [16] focused on reconstruction quality using SOFI-specific metrics, without assessing the fidelity of raw photon emission dynamics. Here, we introduce the use of the 1-Wasserstein distance to quantitatively compare emission distributions from simulated and experimental datasets. This metric captures not only mean differences but also the geometric shape of distributions—offering a more appropriate evaluation metric than classical ones.

That said, we also acknowledge the limitations. While MUFASA better matches experimental traces than other simulators, it does not achieve perfect modeling. This discrepancy comes from several sources: camera effects in MUFASA are modeled using standard models, without full calibration to specific hardware; the point spread function (PSF) is
<figure><img src="image_10.png" /><figcaption>**Figure 9:** Comparison of simulation fidelity across models using the Wasserstein distance. Wasserstein distances ($W_1$) between experimental photon count distributions and simulation models: MUFASA and SOFI Tools. Lower values indicate closer alignment with the empirical distribution. MUFASA achieves the best match, highlighting its capacity to capture realistic emission variability under the StORM protocol.</figcaption></figure>

over-simplified or borrowed from available datasets, introducing further sources of errors.

More importantly, MUFASA is not designed to replicate individual experiments. It does not simulate molecular diffusion or motion within the sample—fluorophore positions are assumed to be fixed and static. Nor does it capture the full environmental complexity that governs transition rates, such as pH, temperature, or oxygen concentration, which are often not available from experimental datasets. Instead, MUFASA is intended as a pre-experiment planning tool or a data generation engine. Its primary use cases include predicting how experimental parameter choices might impact imaging outcomes. Moreover, it can be used to and generate realistic synthetic data for training deep learning based models such as FluoGAN [31], DeepStORM [32] and DE-CODE [33] for image reconstruction.

Finally, the architecture is modular by design, allowing additional imaging modalities to be incorporated without modifying the core simulation logic. For instance, extension to three-dimensional (3D) structures is straightforward; because photon emissions are modeled independently of spatial arrangement, simulated traces can be assigned to fluorophores distributed in 3D. Similarly, multicolour

imaging can be enabled by assigning distinct excitation and emission properties to different fluorophore types, each governed by its own dynamical parameters. As data-driven methods and large-scale analysis become more prevalent, MUFASA can serve not only as a simulator, but also as a physically grounded prior for inverse modeling, uncertainty quantification, and hypothesis testing. It can be used to generate realistic synthetic datasets for training data-driven models or to constrain machine learning pipelines with physically plausible priors.

By uniting simulation fidelity, statistical validation, and experimental design support, MUFASA lays the groundwork for more interpretable, reproducible, and optimized imaging pipelines.

# Methods

## Photoswitching as continuous-time Markov chain

In fluorescence microscopy, the stochastic absorption and emission of photons cause observable fluorescence fluctuations. Inspired by prior approaches such as as Photo-switching hidden Markov model (PSHMM) [34], Single-Molecule Imaging Simulator (SMIS) [15] and SOFI Simulation Tools [16], we define the molecule's state as a continuous-time stochastic process $X(t)$, with $t \in [0, T]$ representing time over the duration of the experiment $T$ and $X(t)$ assuming values from discrete set of quantum states; more details are given below.

Assuming the Markov property -that future states depend only on the current state- we model, in fact, $X(t)$ as a time-homogeneous *Continuous-Time Markov Chain* (CTMC). In CTMCs, state transitions occur at any time point, with holding times in each state following an exponential distribution. A justification of this assumption can be found in Supplementary Note 1.
In fluorescence microscopy, "On" and "Off" states are commonly used to describe fluorophore behavior. We propose to associate these notions to the notion of 'regimes' to better capture the detailed nature of the transitions. During the "On" regime, a molecule cycles between the ground state $S_0$ and the excited singlet state $T_1$. The emission of one photon corresponds to the transition from $S_1$ to $S_0$. Inspired by the Jablonski diagram [1], we then define the state space of the CTMC model and denote it as $\mathcal{S}$. The state space varies depending on the specific super-resolution protocol. The general state space is thus $\mathcal{S} = \{NA, S_0, S_1, T_1, B\}$, where NA represents the non-activated state, $S_0$ and $S_1$ correspond respectively to the ground and excited singlet states, $T_1$ denotes the triplet state, and B is the irreversible bleached state.

Models in SMIS [15] and PSHMM [34] introduce multiple dark states including triplet or redox states, resulting in $N$ distinct dark states. However, as noted in Remark 1 in [35], in the context of Markov chains, states with identical transition rates can be combined into a single state without losing the Markov property. This allows us to aggregate various dark states into a single regime, termed $T_1$, which is accessible exclusively from $S_1$

Given state space $\mathcal{S}$, we can associate it to a specific transition graph like in figure 10. These diagrams are inspired by prior formulations of fluorophore dynamics using hidden Markov models [34], which we generalize here to accommodate a range of microscopy protocols including StORM, PALM, and FF-based imaging. Note that, depending on the specific experimental conditions or photophysical model under consideration, only a subset of states in $\mathcal{S}$ may be relevant. This modularity enables the model to adapt flexibly to a wide range of protocols and fluorophore behaviors.

The CTMC models we consider are defined in terms of environment-dependent transition rates which depend on specific properties of the considered molecules.

We mainly focus on the absorption/excitation rate $k_{01}$ that governs the transition from $S_0$ to $S_1$. This rate can be linked to specific molecules and the laser properties through:

$$k_{01} = \frac{10^{-6} \ln(10) \varepsilon P \lambda}{N_A h c} \quad (1)$$

where $\varepsilon[M^{-1}cm^{-1}]$ is the molar extinction coefficient of the molecule, $\lambda$ [nm] is the excitation wavelength. Here, $P[W \cdot -2]$ is the illumination power density at $(x, y)$. $N_A$ is the Avogadro number, $h$ is Planck's constant, and $c$ is the speed of light [36]. This equation reveals that excitation becomes more probable for fluorophores with higher $\varepsilon$ or stronger illumination $P$. It thus provides a direct and interpretable relationship between experimental parameters and the molecular transition dynamics that govern fluorescence.

Other transition rates in the CTMC model are expressed through the parameter $\tau$, the mean holding times in given states, with $k = \frac{1}{\tau}$. Specifically, the emission rate $k_{10}$ from $S_1$ to $S_0$ is defined as $k_{10} = \frac{1}{\tau_1}$ where $\tau_1 \approx 10^{-9} s$ is the fluorescence lifetime. The intersystem crossing rate $k_{1T}$ from $S_1$ to $T_1$, is instead given by $k_{1T} = \frac{\phi_{ISC}}{\tau_1}$, where the intersystem crossing quantum yield $\phi_{ISC}$ describes the sensitivity to molecular structures and environment. The triplet state decay is $k_{0T} = \frac{1}{\tau_T}$; $\tau_T$ is significantly bigger than $\tau_1$ due to the spin-forbidden nature of the transition. The photobleaching rate $k_{1B}$ representing the transition from $S_1$ to the irreversible bleached state $B$ is described by $k_{1B} = \frac{\phi_b}{\tau_1}$, where $\phi_b$ represents the inverse of the average number of excitation-emission cycles before bleaching.

Defining precise closed-form expressions for these transition rates is challenging due to complex de-
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10:** Quantum state transition diagrams for fluorophore modeling across microscopy protocols. Each diagram represents a continuous-time Markov model over discrete photophysical states, with arrows indicating allowed quantum transitions. (a) General five-state model incorporating all relevant transitions: activation from a non-active state (NA), cycling between ground $S_0$ and excited $S_1$ states, intersystem crossing to triplet $T_1$, and irreversible bleaching (B). (b) Subset used for StORM and Blinking simulations. (c) PALM-specific model including activation from NA. (d) Model for fluctuation-based (FF) microscopy, assuming continuous excitation without activation or triplet transitions. All states are time-homogeneous and transitions follow exponential dwell-time assumptions.</figcaption></figure>

pendencies on environmental factors (temperature, pH, illumination intensity). Furthermore, the variability introduced by different experimental setups and environmental conditions prevents generalization across different scenarios. We provide general expressions that capture the fundamental behaviors without relying on overly specific parameters. For example, the bleaching rate $k_{1B}$ is computed based on the expected number of excitation cycles a molecule undergoes before photobleaching. The intersystem crossing rate $k_{1T}$ is set as a fixed fraction of the relaxation rate $k_{10}$, reflecting the competition between radiative decay and triplet transition:

$$k_{1T} = \frac{p}{1-p} \cdot k_{10},$$

where $p$ denotes the user-defined probability of triplet transition from the excited state. The return rate from the triplet state, $k_{T0}$, is defined from a user-specified triplet lifetime.

**Emission phenomena as a Poisson process**

Fluorescence microscopy methods modeling photon emission typically assume a Poisson distribution for fluorescence fluctuations [35]. A fundamental justification for this hypothesis was presented by [37], where the authors showed that photon counts during the "on" regime can be modeled as a sum of Bernoulli trials—thus approximating a Poisson distribution in the limit where the number of emission attempts per interval becomes large, the probability of emission per attempt becomes small, and the expected number of emissions remains constant (i.e., $n \to \infty$, $p \to 0$, with $np = \mu$ fixed).

In this work, we provide a complementary justification for the Poissonian nature of fluorescence fluctuations using the continuous-time Markov chain (CTMC) framework described above. Our mathematical derivation is presented in full in Supplementary Note 2; it shows that under the realistic photophysical regime where excitation is much slower than relaxation, the photon emission CTMC process con-
verges to a Poisson distribution without requiring
the short-time approximation.

This connection between CTMC dynamics and
Poisson statistics not only reinforces the modeling
choices made in existing frameworks (see, e.g Fluo-
GAN [31]), but also provides an interpretable phys-
ical mapping: : the reconstructed image intensity x
corresponds to the mean photon count, x = nk₀₁t,
where n is the number of fluorescent molecules and
k₀₁ the excitation rate. This formulation embeds
physical interpretability into the generative model,
linking the learned image directly to the underlying
photophysical parameters.

Optical imaging system

The image formation process in fluorescence mi-
croscopy involves optical blur, spatial sampling, and
the modeling of detector-specific noise. MUFASA
simulates these effects through a modular pipeline,
modeling both the optics and the camera electronics.

In more details, emission at high spatial resolu-
tion is first convolved with a diffraction-limited point
spread function (PSF), modeled for simplicity as a
two-dimensional Gaussian. The resulting blurred
image is then spatially undersampled by a factor
L > 1, reflecting the realistic pixel size of the de-
tector.

MUFASA-Sim further incorporates a detailed
camera model, supporting three commonly
used sensor architectures—CCD, EMCCD, and
CMOS—based on established physical and statisti-
cal models from the literature [38, 39]. A complete
description of the optical and electronic modeling
pipeline is provided in Supplementary Note 3.

Software architecture and performance

All simulations were run on a Linux workstation
equipped with 16 physical CPU cores (22 logical),
30.8 GB of RAM, and Python 3.8. While GPU ac-
celeration was not used, performance remained ef-
ficient thanks to multithreaded CPU execution and

Intel MKL-optimized numerical libraries. Depend-
ing on the protocol and fluorophore density, simula-
tion times ranged from less than 1 to 2 seconds per
frame.

To quantify scalability, we benchmarked the run-
time across protocols and increasing fluorophore
counts. As shown in Supplementary Figure 9, sim-
ulation times grow approximately linearly with the
number of molecules. StORM and PALM protocols
exhibit lower runtime due to rapid photobleaching,
whereas FF and Blinking simulations run longer due
to persistent emission dynamics.

Importantly, MUFASA’s modular simulation ar-
chitecture allows for efficient reuse of previously gen-
erated sequences. Simulated traces and emission
profiles can be saved, recombined, or stacked across
conditions—reducing redundant computations and
enabling fast iteration during model testing or pro-
tocol design.

The graphical interface was developed using the
PySide6 library [40]. The GUI is designed to sup-
port both reproducibility and hands-on exploration
of simulation outputs. While optical effects like PSF
blur are modeled using Gaussian kernels or user-
defined profiles, MUFASA does not aim to replicate
hardware-specific optics or cameras in fine detail.
Instead, it provides a flexible, physically grounded
platform for pre-experiment planning, data genera-
tion, and model testing.

Code availability

The MUFASA simulation framework is openly
available. The complete source code, in-
cluding detailed documentation and a user-
friendly graphical interface, is accessible at
https://github.com/Joywessim/MUFASA_
Fluorescence_Fluctuations_Simulation.
## 814 Acknowledgements

This work is based upon work supported by the ANR JCJC project TASKABILE ANR-22-CE48-0010 and the funding received from the European Research Council (ERC) Starting project MALIN under the European Union's Horizon Europe programme (grant agreement No. 101117133). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. This work has been supported by the French government, through the 3IA Cote d'Azur Investments in the project managed by the National Research Agency (ANR) with the reference number ANR-23-IACL-0001.

Image acquisitions were conducted in the microscopy platform PIM (member of MICA microscopy platform and EMBRC-France), whose French state funds are managed by the ANR-10-INBS-0.

## 834 Competing Interests

The authors declare no competing interests.

## 836 References

[1] Alva, A. et al. Fluorescence fluctuation based super resolution microscopy, basic concepts for an easy start. bioRxiv (2022). URL https://www.biorxiv.org/content/early/2022/05/06/2022.05.06.490863. https://www.biorxiv.org/content/early/2022/05/06/2022.05.06.490863.full.pdf.

[2] Hell, S. W. & Wichmann, J. Breaking the diffraction resolution limit by stimulated emission: stimulated-emission-depletion fluorescence microscopy. Optics Letters 19, 780-782 (1994).

[3] Bailey, B., Farkas, D. L., Taylor, D. L. & Lanni, F. Enhancement of axial resolution in fluorescence microscopy by standing-wave excitation. Nature 366, 44-48 (1993). URL https://doi.org/10.1038/366044a0.

[4] Jacquemet, G., Carisey, A. F., Hamidi, H., Henriques, R. & Leterrier, C. The cell biologist's guide to super-resolution microscopy. Journal of Cell Science 133, jcs240713 (2020).

[5] Rust, M. J., Bates, M. & Zhuang, X. Subdiffraction-limit imaging by stochastic optical reconstruction microscopy (storm). Nature Methods 3, 793-795 (2006).

[6] Betzig, E. et al. Imaging intracellular fluorescent proteins at nanometer resolution. Science 313, 1642-1645 (2006).

[7] Dertinger, T., Colyer, R., Iyer, G., Weiss, S. & Enderlein, J. Fast, background-free, 3d super-resolution optical fluctuation imaging (sofi). Proceedings of the National Academy of Sciences 106, 22287-22292 (2009).

[8] Chen, Y., Müller, J. D., So, P. T. C. & Graton, E. The photon counting histogram in fluorescence fluctuation spectroscopy. Biophysical Journal 77, 553-567 (1999). URL https://doi.org/10.1016/S0006-3495(99)76912-2.

[9] Nieuwenhuizen, R. P. et al. Quantitative localization microscopy: Effects of photophysics and labeling stoichiometry. PloS one 10, e0127989 (2015).

[10] Annibale, P., Vanni, S., Scarselli, M., Rothlisberger, U. & Radenovic, A. Quantitative photo-activated localization microscopy: Unraveling the effects of photoblinking. PloS one 6, e22678 (2011).

[11] Lee, S.-H., Shin, J. Y., Lee, A. & Bustamante, C. Counting single photoactivatable fluorescent molecules by photoactivated localization
microscopy (palm). Proceedings of the National Academy of Sciences 109, 17436-17441 (2012).

[12] Venkataramani, V., Herrmannsdörfer, F., Heilemann, M. & Kuner, T. SuReSim: simulating localization microscopy experiments from ground truth models. Nature Methods 13, 319-321 (2016). URL https://doi.org/10.1038/nmeth.3795.

[13] Lagardère, M., Chamma, I., Bouilhol, E., Nikol ski, M. & Thoumine, O. FluoSim: sim ulator of single molecule dynamics for fluo rescence live-cell and super-resolution imag ing of membrane proteins. Scientific Reports 10, 19954 (2020). URL https://doi.org/10.1038/s41598-020-77094-6.

[14] Griffi è, J. et al. Virtual-smlm, a virtual environment for real-time interactive smlm acquisition. bioRxiv 2020.03.05.967893 (2020). URL https://www.biorxiv.org/content/ early/2020/03/06/2020.03.05.967893. Preprint.

[15] Bourgeois, D. Single molecule imaging simula tions with advanced fluorophore photophysics. Communications Biology 6, 53 (2023).

[16] Girsault, A. et al. SOFI Simulation Tool: A Software Package for Simulating and Testing Super-Resolution Optical Fluctuation Imaging. PLoS ONE 11, e0161602 (2016). URL https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0161602. PMCID: PMC5008722.

[17] Scientific, T. F. The alexa fluor dye series (2024). URL https://www.thermofisher.com/fr/fr/home/references/molecular-probes-the-handbook/technical-notes-and-product-highlights/the-alexa-fluor-dye-series.html. Accessed: 16-Oct-2024.

[18] Sage, D. et al. Quantitative evaluation of software packages for single-molecule localization microscopy. Nature Methods 16, 387-395 (2019).

[19] Liu, Y. et al. Revisiting psf models: unifying framework and high-performance implementation (2025). URL https://arxiv.org/abs/2502.03170. 2502.03170.

[20] Cox, G. & Sheppard, C. J. Shot noise, detection noise, and the limits of the light microscope. Journal of the Optical Society of America A 21, 820-827 (2004).

[21] Banterle, N., Bui, T.-H., Vojnovic, B., Collinson, L. M. & Kaminski, C. F. Camera-based image quality assessment for fluorescence microscopy. Journal of Microscopy 250, 35-49 (2013).

[22] Robbins, M. S. & Hadwen, B. J. The noise performance of electron multiplying charge-coupled devices. IEEE Transactions on Electron Devices 50, 1227-1232 (2003).

[23] Pace, C., Vajdos, F., Fee, L., Grimsley, G. & Gray, T. How to measure and predict the molecular absorption coefficient of a protein. Protein Science 4, 2411-2423 (1995).

[24] Hertzog, M., Wang, M., Mony, J. & Börjesson, K. Strong light-matter interactions: a new direction within chemistry. Chemical Society Reviews 48, 937-961 (2019).

[25] Weight, B. M., Li, X. & Zhang, Y. Theory and modeling of light-matter interactions in chemistry: current and future. arXiv preprint arXiv:2303.10111 (2023).

[26] Kumar, A., Singh, S., Mudahar, G. S. & Thind, K. S. Molar extinction coefficients of some commonly used solvents. Radiation Physics and Chemistry 75, 731-736 (2006).
962

[27] Patterson, G. H., Knobel, S. M., Sharif, W. D., Kain, S. R. & Piston, D. W. Use of the green fluorescent protein and its mutants in quantitative fluorescence microscopy. *Biophysical journal* **73**, 2782–2790 (1997).

963

964

965

966

967

[28] Boudreau, J., Paramasivam, G. & Godin, M. Excitation light dose engineering to reduce photobleaching and phototoxicity. *Scientific Reports* **6**, 30892 (2016).

968

969

970

971

[29] Villani, C. *Optimal transport: old and new*, vol. 338 (Springer, 2008).

972

973

[30] Olivier, N. & Keller, D. STORM Vectashield datasets (Tubulin). https://doi.org/10.5281/zenodo.7620025 (2023). Zenodo dataset.

974

975

976

[31] Cachia, M., Stergiopoulou, V., Calatroni, L., Schaub, S. & Blanc-Féraud, L. Fluorescence image deconvolution microscopy via generative adversarial learning (fluogan). *Inverse Problems* **39**, 054006 (2023). URL https://hal.science/hal-03790156v2.

977

978

979

980

981

982

[32] Nehme, E., Weiss, L. E., Michaeli, T. & Shechtman, Y. Deep-storm: super-resolution single-molecule microscopy by deep learning. *Optica* **5**, 458–464 (2018).

983

984

985

986

[33] Speiser, A. et al. Deep learning enables fast and dense single-molecule localization with high accuracy. *Nature Methods* **18**, 1082–1090 (2021).

987

988

989

[34] Patel, L. et al. A hidden markov model approach to characterizing the photo-switching behavior of fluorophores. *Annals of Applied Statistics* **13**, 1397–1429 (2019).

990

991

992

993

[35] Staudt, T. et al. Statistical molecule counting in super-resolution fluorescence microscopy: Towards quantitative nanoscopy (2020). URL https://arxiv.org/abs/1903.11577. 1903.11577.

994

995

996

997

[36]

Avilov, S. et al. In cellulo evaluation of phototransformation quantum yields in fluorescent proteins used as markers for single-molecule localization microscopy. *PloS One* **9**, e98362 (2014).

998

999

1000

1001

1002

[37]

Aspelmeier, T., Egner, A. & Munk, A. Modern statistical challenges in high-resolution fluorescence microscopy. *Annual Review of Statistics and Its Application* **2**, 163–202 (2015).

1003

1004

1005

1006

[38]

Tian, H. & El Gamal, A. *Noise Analysis in CMOS Image Sensors*. Dissertation, Stanford University, Stanford, CA (2000). Ph.D. dissertation.

1007

1008

1009

1010

[39]

Hirsch, M., Wareham, R., Martin-Fernandez, M., Hobson, M. & Rolfe, D. A stochastic model for electron multiplication charge-coupled devices–from theory to practice. *PLoS One* **8**, e53671 (2013). Epub 2013 Jan 31.

1011

1012

1013

1014

1015

[40]

Company, T. Q. Pyside6 - qt for python. https://doc.qt.io/qtforpython-6/ (2024). Version 6.9.

1016

1017

1018
Supplementary Information for:

MUFASA: A Continuous-Time Stochastic Framework for Realistic
Fluorescence Microscopy Simulation

Wessim Omezzine<fnref n="1" />, Sébastien Schaub<fnref n="2" />, Laure Blanc-Féraud<fnref n="1" />, and Luca Calatroni<fnref n="3" />

This file contains Supplementary Notes, Figures, and Tables referenced in the main manuscript.
We provide also two supplementary videos of two simulations.

Contents

<table><tr><td>1 Supplementary Note 1: Continuous-time Markov Chains</td><td>2</td></tr><tr><td>2 Supplementary Note 2: Approximation of CTMC with |S| = 2</td><td>3</td></tr><tr><td>3 Supplementary Note 3: Optical and Camera System Modeling</td><td>6</td></tr><tr><td>4 Supplementary Note 4: SOFI-Simlulation Tools</td><td>9</td></tr><tr><td>5 Supplementary Note 5: Camera Configurations</td><td>10</td></tr><tr><td>6 Supplementary Note 6: Experimental Protocol Parameters</td><td>11</td></tr><tr><td>7 Supplementary Note 7: Extracting Single Molecule Signals from Experimental Acquisitions</td><td>12</td></tr><tr><td>8 Supplementary Note 8: Experimental Conditions and Parameter Inference for STORM Validation</td><td>14</td></tr><tr><td>9 Supplementary Note 9: Wasserstein Distance W<sub>1</sub></td><td>16</td></tr><tr><td>10 Supplementary Note 10: Quantitative assessments via Wasserstein-distance of FF dynamics</td><td>17</td></tr><tr><td>11 Supplementary Figures</td><td>19</td></tr></table>

<footnote n="1">¹Université Côte d'Azur, CNRS, I3S</footnote>

<footnote n="2">² Sorbonne University, CNRS, LBDV</footnote>

<footnote n="3">³MaLGA Centre, DIBRIS, University of Genoa and MMS, Italian Institute of Technology, Genoa, Italy</footnote>

<footnote>*Correspondence: omezzine@i3s.unice.fr, sebastien.schaub@imev-mer.fr, blancf@i3s.unice.fr, luca.calatroni@unige.it</footnote>
1 Supplementary Note 1: Continuous-time Markov Chains

In fluorescence microscopy, the stochastic nature of photon absorption and emission events leads to observable fluctuations in fluorescence signals. These fluctuations arise because photon-molecule interactions occur randomly over time, even under constant illumination. To model these dynamics, we represent the state of a fluorescent molecule as a stochastic process, denoted as $X(t) \in S$, where $t$ represents continuous time and $S \subset \mathbb{N} \cup \{0\}$ is the state space. Its general form is $S = \{NA, S_0, S_1, T_1, B\}$, where NA represents the non-activated state, $S_0$ and $S_1$ correspond respectively to the ground and excited singlet states, $T_1$ denotes the triplet state, and $B$ is the irreversible bleached state. Given that molecules have a finite set of quantum states, we can assume i.e $|S| < +\infty$.

We now introduce the equivalent in continuous time of the Markov assumption: the probability associated to the future state of the system within a certain time duration depends only on its current state, not on its past history, thus implying a memory-less property [9]. This applies for instance when the system interaction with its environment is minimal and the environment relaxes rapidly, losing its memory faster than the system dynamics, which is the case of molecular state transitions, which can often be effectively modeled as Markovian processes when all relevant states and interactions are accounted for [1].

Building upon the Markov property, we further assume that the process is time-homogeneous, meaning the transition probabilities between states are consistent over time. This leads us to model the fluorescence fluctuations using a Continuous-Time Markov Chain (CTMC) $(X(t))_{t \ge 0}$ with state space $S$, [9]. By definition, the transition probabilities satisfy

$$P(X(s+t) = j | X(s) = i) = P(X(t) = j | X(0) = i) =: P_{ij}(t), \quad \forall s, t \ge 0, \forall i, j \in S.$$

Unlike discrete-time Markov chains, where transitions occur at fixed time steps, CTMCs permit transitions at any continuous time point. The process remains in a given state $i \in S$ for a random duration $T_i$, known as the *holding time*, which—through the memoryless property can be shown to follow an exponential distribution:

$$P(T_i > t) = e^{-k_i t}, \quad t \ge 0,$$

where $k_i > 0$ denotes the total exit rate from state $i$.

The infinitesimal behavior of the CTMC is governed by the *generator matrix* (or rate matrix) $G = (k_{ij})$, whose entries satisfy:

$$k_{ij} \ge 0 \quad \text{for } i \neq j, \quad k_{ii} = \mu_i = -\sum_{j \neq i} \mu_{ij}, \quad \forall i \in S.$$

Thus, each row of $G$ sums to zero (differently from discrete-time Markov chains where transition matrices are stochastic), and the diagonal elements encode the total rate of leaving each state.
## 2 Supplementary Note 2: Approximation of CTMC with $|S| = 2$

The unified CTMC modeling used in MUFASA offers a theoretical insight into the interpretation of such parameters especially in the case of continuous emission before bleaching. We prove here that whenever the cardinality of the state space is equal to two and some assumptions on the transition rates hold, we can in fact approximate the Markov process by a Poisson process with a physical grounded parameter being its mean and variance.

Consider the simple example where $|S| = 2$, with $S = \{S_0, S_1\}$. The generator matrix of a CTMC process in this case reads:

$$G = \begin{bmatrix} -k_{01} & k_{01} \\ k_{10} & -k_{10} \end{bmatrix}$$

and the state transition matrix at time $t$ is given by definition by $P(t) = e^{Gt}$. We now assume

$$k_{01} \ll k_{10}.$$

This is a physically reasonable assumption in fluorescence photophysics, where excitation events are rare and relaxation from the excited state is typically fast [3].

For $t \in [0, +\infty]$, let $N(t)$ be the random process denoting the number of photons emitted in the interval $[0, t]$ and let us define $P(n, T) := \mathbb{P}(N(t) = n), \forall n \in \mathbb{N}^*$ be the probability of emitting exactly $n$ photons until instant $t$. Note that by definition there holds:

$$P(n,t) = \mathbb{P}(N(t) = n, X(t) = S_0) + \mathbb{P}(N(t) = n, X(t) = S_1),$$

which, using similar notations, can be decomposed as

$$P(n,t) = P_0(n,t) + P_1(n,t).$$

We study now the dynamics of the functions $P_i(n, t)$ in order to find its expression. For small time increments $\delta t$, we have:

$$\begin{aligned} P_0(n, t + \delta t) &= \mathbb{P}(N(t) = n, X(t + \delta t) = S_0) \\ &= \mathbb{P}(N(t) = n, X(t + \delta t) = S_0 | N(t) = n, X(t) = S_0) \mathbb{P}(N(t) = n, X(t) = S_0) \\ &\quad + \mathbb{P}(N(t) = n, X(t + \delta t) = S_0 | N(t) = n-1, X(t) = S_1) \mathbb{P}(N(t) = n-1, X(t) = S_1) \\ &= (1 - k_{01}\delta t) P_0(n, t) + k_{10}\delta t P_1(n-1, t) + o(\delta t) \end{aligned}$$

We can thus compute the differential equation governing $P_0(n, t)$, which reads:

$$\frac{d}{dt} P_0(n, t) = -k_{01} P_0(n, t) + k_{10} P_1(n-1, t)$$

Proceeding similarly for $P_1(n, t)$, we get the differential system:

$$\left\{ \begin{array}{l} \frac{d}{dt} P_0(n,t) = -k_{01} P_0(n,t) + k_{10} P_1(n-1,t) \\ \frac{d}{dt} P_1(n,t) = k_{01} P_0(n,t) - k_{10} P_1(n,t) \\ P_0(n,0) = \delta_{n0} \\ P_1(n,0) = 0 \end{array} \right. \qquad (1)$$

where $\delta_{n0}$ denotes the Kronecker symbol. These equations are known as master equations [1]: equations for the probability that a stochastic particle can jump between different states is in one of these states at a time $t$.
To solve the master equations, we introduce the probability generating functions (PGFs):

$$ \forall i = 0, 1 \quad \forall z \text{ s.t. } |z| < 1, \quad G_i(z, t) := \sum_{n=0}^{\infty} z^n P_i(n, t). $$

The PGF compactly encodes the entire desired probability distribution $P_i(n, t)$ into one single analytic expression. In more details, it transforms the infinite sequence of probabilities into a function whose behavior can be studied using differential equations. By multiplying both sides of the master equations (1) by $z^n$ and summing over all $n$, we get:

$$ \left\{ \begin{aligned} \frac{d}{dt} G_0(z,t) &= -k_{01}G_0(z,t) + k_{10}zG_1(z,t) \\ \frac{d}{dt} G_1(z,t) &= k_{01}G_0(z,t) - k_{10}G_1(z,t) \\ G_0(z,0) &= \sum_{n=0}^{\infty} z^n P_0(n,0) = \sum_{n=0}^{\infty} z^n \delta_{n,0} = 1 \\ G_1(z,0) &= 0. \end{aligned} \right. \quad (2) $$

Our goal now is to solve this system for $G(z,t) = G_0(z,t) + G_1(z,t)$, from which we can compute the desired distribution. To do so, we now apply the Laplace transform to both PGFs to turn the system of time-domain differential equations into algebraic equations. The Laplace variable $s \in \mathbb{C}$ is assumed to lie in the domain where the integrals converge, which is typically Re(s) > 0. We obtain:

$$ \tilde{G}_i(z, s) = \mathcal{L}[G_i](z, s) = \int_0^\infty e^{-st} G_i(z, t) dt, $$

so that system (2) becomes:

$$ \left\{ \begin{aligned} s\tilde{G}_0(z, s) - G_0(z, 0) &= -k_{01}\tilde{G}_0(z, s) + k_{10}z\tilde{G}_1(z, s), \\ s\tilde{G}_1(z, s) - G_1(z, 0) &= k_{01}\tilde{G}_0(z, s) - k_{10}\tilde{G}_1(z, s). \end{aligned} \right. $$

Now, since $G_0(z, 0) = 1$ and $G_1(z, 0) = 0$, after some algebraic manipulations we get:

$$ \left\{ \begin{aligned} \tilde{G}_0(z, s) &= \frac{s+k_{10}}{s^2+s(k_{01}+k_{10})+k_{01}k_{10}(1-z)}, \\ \tilde{G}_1(z, s) &= \frac{k_{01}}{s^2+s(k_{01}+k_{10})+k_{01}k_{10}(1-z)}. \end{aligned} \right. $$

We thus deduce that:

$$ \forall z \text{ s.t. } |z| < 1, \quad \tilde{G}(z, s) = \tilde{G}_0(z, s) + \tilde{G}_1(z, s) = \frac{s + k_{10} + k_{01}}{s^2 + s(k_{01} + k_{10}) + k_{01}k_{10}(1-z)}. $$

We can thus find $G(z, t)$ by taking the inverse Laplace transform of $\tilde{G}(z, s)$ so that:

$$ \forall z \text{ s.t. } |z| < 1, \quad G(z,t) = \mathcal{L}^{-1} \left\{ \frac{s + k_{10} + k_{01}}{s^2 + s(k_{01} + k_{10}) + k_{01}k_{10}(1-z)} \right\}. $$

Now, let us express by

$$ s_{\pm}(z) = \frac{-(k_{10} + k_{01}) \pm \sqrt{(k_{10} + k_{01})^2 - 4k_{10}k_{01}(1-z)}}{2}, $$

the two roots of the second degree equation in the denominator above. By now setting $\epsilon := k_{01}/k_{10} \ll 1$ and $\delta := -2\epsilon + 4\epsilon z + \epsilon^2$ we have that:

$$ k_{10}\sqrt{(1+\epsilon)^2 - 4\epsilon(1-z)} = k_{10}\sqrt{1+\delta} = k_{10}(1-\epsilon + 2\epsilon z + O(\epsilon^2)). $$
whence we find:

$$s_{+}(z) = -\epsilon(1-z)k_{10} + O(\epsilon^2 k_{10}) = -k_{01}(1-z) + O(\epsilon k_{01}), \quad s_{-}(z) = -[1+O(\epsilon)]k_{10}.$$

Neglecting the fast-mode contributions. Writing now

$$\tilde{G}(z, s) = \frac{s + k_{10} + k_{01}}{(s - s_{+}(z))(s - s_{-}(z))} = \frac{C}{s - s_{+}(z)} + \frac{D}{s - s_{-}(z)},$$

we have that by direct computations we can compute

$$C = 1 + O(\epsilon), \quad D = O(\epsilon).$$

In the time domain, we thus have

$$G(z, t) = C e^{s_{+}(z)t} + D e^{s_{-}(z)t} = e^{-k_{01}(1-z)t} [1 + O(\epsilon)] + O(\epsilon) e^{-k_{10}t}.$$

Therefore, for any $t$ such satisfying $k_{10}t \gg 1$, we have that $e^{-k_{10}t}$ is exponentially small, hence we may drop the second term in the expression above, thus getting:

$$G(z,t) \approx \exp[-k_{01}(1-z)t] = \sum_{n=0}^{\infty} \frac{(k_{01}t)^n}{n!} z^n,$$

i.e., for $t \gg 1/k_{10}$ the following approximation holds $N(t) \sim \text{Poisson}(k_{01}t)$.

In fluorescence settings where $k_{10} \approx 10^9 \text{s}^{-1}$, the Poisson approximation becomes valid for any observation time more than a nanosecond. Beyond this timescale, the contribution of the fast relaxation mode $e^{-k_{10}t}$ becomes negligible, and the excitation-relaxation cycles effectively behave as independent, memoryless events. As a result, the photon count statistics converge rapidly to a Poisson distribution with rate $k_{01}$, and deviations from this model are only significant at extremely short timescales.
### 3 Supplementary Note 3: Optical and Camera System Modeling

The image formation process in fluorescence microscopy involves optical blur, spatial sampling, and detector-specific noise corruption. MUFASA simulates these effects through a modular pipeline, modeling both the physical optics, with a simplified model, and the camera electronics.

**Optical blur and spatial undersampling.** Photon emission at high spatial resolution is first blurred by the optical point spread function (PSF), modeled here as a 2D Gaussian $G(\sigma)$, where $\sigma$ is a user-defined parameter that defines the spread of the PSF in pixel units. Let $\mathbf{x}^t \in \mathbb{R}^{H \times W}$ denote the high-resolution ground-truth image at time $t$. The blurred image is given by:

$$ \mathbf{X}_{\text{blur}}^t = G(\sigma) * \mathbf{x}^t, $$

where $\mathbf{X}_{\text{blur}}^t \in \mathbb{R}^{H \times W}$ is the 2D blurred image. If needed for downstream processing, this result is then flattened as a vector:

$$ \mathbf{x}_{\text{blur}}^t = \Psi \mathbf{x}^t = \text{vec}(\mathbf{X}_{\text{blur}}^t). $$

To reflect the finite pixel size of the detector, the high-resolution image is spatially undersampled by a factor $L$, corresponding to a pixel area of $L \times L$ in the high-resolution grid. This operation is represented by a linear downsampling operator $U \in \mathbb{R}^{m \times n}$, where $n = H \times W$ is the number of high-resolution pixels and $m = n/L^2$ is the number of detector pixels. Each low-resolution pixel value is computed as the sum of intensities within its corresponding $L \times L$ region of the high-resolution image.

In the context of time-varying molecular emission, the final recorded intensity in each detector pixel integrates the blurred photon contributions from all fluorophores that were active during the acquisition time and spatially located within the corresponding region. This yields the low-resolution intensity:

$$ \tilde{\mathbf{x}}^t = U \Psi \mathbf{x}^t + \mathbf{b}, $$

where $\mathbf{b}$ models background fluorescence coming from out-of-focus regions.

**Photon detection and camera modeling.** ot all emitted photons reach or are registered by the detector due to optical losses, quantum efficiency limitations, and electronic artifacts. This induces the stochasticity of photon detection given by a Poisson distribution. The number of detected photoelectrons $n_{ie}$ is, then, modeled as:

$$ n_{ie} = \mathcal{P}(QE \cdot \tilde{\mathbf{x}}^t + c), $$

where $QE$ is the quantum efficiency and $c$ is the clock-induced charge, a small spurious signal generated during charge transfer in the detector.

MUFASA supports multiple camera types:

* **CCD:** [6] Direct conversion with Gaussian readout noise:

$$ n_{oe} = n_{ie} + \mathcal{N}(\mu, \sigma^2). $$

* **EMCCD:** [6] Stochastic amplification with a Gamma distribution, followed by Gaussian noise:

$$ n_{oe} = \Gamma(n_{ie}, EM_{\text{gain}}) + \mathcal{N}(\mu, \sigma^2). $$
• CMOS: [11] Nonlinear pixel-specific gain function $g_p(\cdot)$, followed by Gaussian readout noise:

$$n_{oe} = g_p(n_{ie}) + \mathcal{N}(\mu, \sigma^2).$$

**Digitization.** The final output is quantized into Analog-to-Digital Units (ADUs):

$$\mathbf{y}^t = \min \left( \left\lfloor \frac{n_{oe}}{e_{ADU}} \right\rfloor + \text{BL}, \text{ADC}_{\max} \right),$$

where:

• $e_{ADU}$ is the electron-to-ADU conversion factor,

• BL is the baseline offset,

• $\text{ADC}_{\max}$ is the digital saturation level (e.g., 65535 for a 16-bit ADC).

This complete modeling chain ensures that simulated data replicates actual real fluorescence microscopy acquisitions.

**Denoising Acquisition Signal** To provide a rough estimation of the underlying photon emission signal $\tilde{\mathbf{x}}^t$ from the observed digitized acquisition $\mathbf{y}^t$, we exploit the known structure of the physical acquisition chain detailed previously. This knowledge allows us to leverage the statistical properties of each noise source—additive Gaussian readout noise, gamma-distributed amplification noise (for EMCCDs), and Poisson-distributed photon counting noise—to denoise the signal in a modular fashion. Each noise component is addressed with methods tailored to its specific distribution.

All camera types include a Gaussian noise component at the readout stage, modeled as:

$$n_{oe} = \hat{n} + \mathcal{N}(0, \sigma^2),$$

where $\hat{n}$ is the intermediate signal prior to quantification. To remove such noise, we use Gaussian denoisers such as wavelet-based methods to estimate $\hat{n}$ and suppress readout noise, yielding a denoised signal $\hat{n}^{(1)}$.

For EMCCD sensors, after Gaussian denoising, we normalize by the known EM gain to approximate the input to the amplifier:

$$\hat{n}_{ie} \approx \frac{\hat{n}^{(1)}}{\text{EM}_{\text{gain}}}.$$

For CMOS cameras, note that a nonlinear, pixel-specific gain function $g_p(\cdot)$ is applied. The calibration data is generally available for CMOS, so we invert this transformation numerically:

$$\hat{n}_{ie} \approx g_p^{-1}(\hat{n}^{(1)}).$$

In all cases, the final signal $\hat{n}_{ie}$ is modeled as a realization of Poisson-distributed photoelectrons:

$$n_{ie} \sim \mathcal{P}(QE \cdot \tilde{\mathbf{x}}^t + c),$$

for which a Poisson denoiser $\mathcal{D}_{\text{Pois}}$, such as PURE-LET [7] can be applied, to estimate a smoother signal:

$$\hat{n}_{\text{denoised}} = \mathcal{D}_{\text{Pois}}(\hat{n}_{\text{ie}}).$$
We retain $\hat{n}_{clean}$ as a denoised representation of the observed photon counts.

This physics-aware, layered denoising pipeline provides a principled way to stabilize and enhance noisy acquisitions, and is used as a preprocessing step throughout MUFASA-Design.
4 Supplementary Note 4: SOFI-Simulation Tools

SOFI Simulation Tools [5] uses a continuous-time Markov model with a state space $|S| = 2$ with two regimes—“on” and “off”—and models transitions between them using exponential holding times. The rates $k_{on}$ and $k_{off}$ are defined by the mean durations in each state:

$$k_{on} = 1/T_{on}, \quad k_{off} = 1/T_{off}$$

During the “on” regime, SOFI assumes a fixed photon emission rate $I_{on}$, so that the total number of photons emitted during a “on” event is:

$$N_{photon}^{SOFI} = I_{on} \cdot t_{on}, \quad \text{with} \quad t_{on} \sim \mathcal{E}(T_{on})$$

This stochastic framework allows simplified simulations of fluorophore blinking, which is valuable for generating synthetic image stacks for testing super-resolution optical fluctuation imaging (SOFI) and balanced SOFI (bSOFI) pipelines. The simulator also includes tunable parameters for choosing fluorophore density, PSF shape, bleaching dynamics, noise statistics (e.g., Poisson photon shot noise, Gaussian camera readout), and supports both SOFI and STORM reconstruction algorithms. Importantly, the bSOFI module applies a linearization step to correct the nonlinear brightness dependence inherent in higher-order cumulants. For further technical details, the full simulation architecture is described in [5].
## 5 Supplementary Note 5: Camera Configurations

The configurations listed in Table 1 are representative examples. Our framework is designed to be modular, allowing any protocol to be simulated with any of the camera models to study the impact of different hardware choices. These configurations correspond to the frames examples given in the main results and the supplementary notes.

<table><caption>**Supplementary Table 1. Camera configurations used for each simulated SRM protocol in MUFASA-Sim.**</caption><thead><tr><th>Parameter</th><th>PALM</th><th>STORM</th><th>FF</th><th>Blinking</th></tr></thead><tbody><tr><td>Camera type</td><td>CCD</td><td>CCD</td><td>EMCCD</td><td>CMOS</td></tr><tr><td>Quantum efficiency (QE)</td><td>0.9</td><td>0.9</td><td>0.9</td><td>0.9</td></tr><tr><td>Read noise σ<sub>R</sub> (e<sup>-</sup>)</td><td>70.58</td><td>70.58</td><td>60.58</td><td>1.58</td></tr><tr><td>Clock-induced charge c (e<sup>-</sup>)</td><td>0.002</td><td>0.002</td><td>0.002</td><td>0.002</td></tr><tr><td>Electron/ADU e<sub>ADU</sub></td><td>10.04</td><td>100.04</td><td>0.74</td><td>20.74</td></tr><tr><td>Baseline (BL) (ADU)</td><td>367</td><td>367</td><td>367</td><td>367</td></tr><tr><td>ADC max</td><td>65535</td><td>65535</td><td>65535</td><td>65535</td></tr><tr><td>EM (Amplification)</td><td>-</td><td>-</td><td>500</td><td>-</td></tr><tr><td>PSF FWHM (nm)</td><td>1200</td><td>1200</td><td>308</td><td>308</td></tr><tr><td>Pixel size (nm)</td><td>100</td><td>100</td><td>100</td><td>100</td></tr><tr><td>Background (photons)</td><td>20</td><td>20</td><td>20</td><td>20</td></tr><tr><td>CMOS gain function g<sub>p</sub>(e)</td><td>-</td><td>-</td><td>-</td><td>0.9e<sub>−</sub>, e &lt; 200<br>180 + 0.4(e − 200), 200 ≤ e &lt; 1000<br>340 + 0.1(e − 1000), e ≥ 1000</td></tr></tbody></table>
156 6 Supplementary Note 6: Experimental Protocol Parameters

<table><caption>Supplementary Table 2. Experimental setup parameters for each protocol.</caption><thead><tr><td>Parameter</td><td>PALM</td><td>STORM</td><td>FF</td><td>Blinking</td></tr></thead><tbody><tr><td>Experiment Duration (s)</td><td>1.00</td><td>1.00</td><td>1.00</td><td>1.00</td></tr><tr><td>Number of Frames</td><td>500</td><td>500</td><td>100</td><td>100</td></tr><tr><td>Frame Length (ms)</td><td>10</td><td>10</td><td>10</td><td>10</td></tr><tr><td>Excitation Power (W/cm²)</td><td>5</td><td>7</td><td>2</td><td>3</td></tr><tr><td>Activation Rate (s⁻¹)</td><td>1</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Excitation Wavelength (nm)</td><td>647</td><td>647</td><td>647</td><td>647</td></tr></tbody></table>

<table><caption>Supplementary Table 3. Molecular parameters used for each protocol.</caption><thead><tr><td>Parameter</td><td>PALM</td><td>STORM</td><td>FF</td><td>Blinking</td></tr></thead><tbody><tr><td>Extinction Coefficient (ε) (M⁻¹cm⁻¹)</td><td>270000</td><td>270000</td><td>270000</td><td>270000</td></tr><tr><td>Excitation Lifetime (τ) (s)</td><td>2 · 10⁻⁹</td><td>2 · 10⁻⁹</td><td>2 · 10⁻⁹</td><td>2 · 10⁻⁹</td></tr><tr><td>Cycles Before Bleaching</td><td>5 · 10⁵</td><td>5 · 10⁵</td><td>5 · 10³</td><td>1 · 10³</td></tr><tr><td>Non-Radiative Decay Rate (α<sub>nr</sub>)</td><td>10</td><td>2</td><td>1</td><td>3</td></tr><tr><td>ISC Rate (α<sub>isc</sub>)</td><td>1 · 10⁻²</td><td>5 · 10⁻⁴</td><td>1 · 10⁻⁵</td><td>5 · 10⁻⁴</td></tr></tbody></table>
# 7 Supplementary Note 7: Extracting Single Molecule Signals from Experimental Acquisitions

To evaluate the accuracy of MUFASA on real data, we considered time-lapse fluorescence microscopy acquisitions of a sample labeled with EGFP (Enhanced Green Fluorescent Protein). The properties of EGFP are summarized in Table 4.

<table><caption>**Supplementary Table 4. Photophysical properties of EGFP used for experimental validation.**</caption><tr><td>Excitation Wavelength (λ<sub>exc</sub>)</td><td>488 nm</td></tr><tr><td>Emission Wavelength (λ<sub>em</sub>)</td><td>507 nm</td></tr><tr><td>Extinction Coefficient (ε)</td><td>55,900 M<sup>-1</sup>cm<sup>-1</sup></td></tr><tr><td>Quantum Yield (QY)</td><td>0.6</td></tr><tr><td>Fluorescence Lifetime</td><td>2.6 ns</td></tr></table>

The experiment consists of recording a 512 × 512 pixel time series of 250 frames using an EMCCD camera. Each frame corresponds to an amplified digital signal. Prior to illumination, a background frame was acquired to estimate the camera intrinsic noise characteristics (Figure 1, left). Subsequent frames captured the fluorescent dynamics of EGFP-tagged molecules under continuous illumination (Figure 1, right).

<figure><img src="image_13.png" /><figcaption>**Supplementary Figure 1. Camera output before and during fluorescence acquisition.**
Left: Background noise frame recorded before illumination. Right: Fluorescence signal after illumination of EGFP-expressing sample at different time frames. (every 2.5s).</figcaption></figure>

We applied a threshold to each frame to isolate pixels with intensities significantly above the background noise, identifying candidate locations likely to contain active fluorophores. For each candidate region, we assumed that the brightest pixel corresponds to the molecule's spatial center, with the surrounding intensity arising from optical blur (PSF spread). We then extracted the full temporal intensity profile from this central pixel across all frames. Signals that did not exhibit the expected fluorescence dynamics were discarded, ensuring that only well-isolated single-molecule traces were retained.
<figure><img src="image_14.png" /><figcaption>**Supplementary Figure 2.** From full-frame fluorescence to single-molecule signal extraction. Top: Raw intensity dynamics from all candidate pixels exceeding the background noise threshold. For each region, the brightest pixel is assumed to correspond to the fluorophore center, and its temporal trace is extracted. Bottom left: Segmented traces of candidate molecules across the entire acquisition sequence (250 frames). Bottom right: Final cleaned dataset of individual single-molecule fluorescence signals, obtained by selecting a time window of 50 frames just before bleaching, during which all chosen emitters remain active. This ensures a consistent set of unbleached, spatially isolated single-molecule trajectories suitable for model validation.</figcaption></figure>

This procedure enabled the extraction of a dataset consisting of 20 single-molecule fluorescence traces across 50 frames each, suitable for validation of simulation pipelines.

174

175
# 8 Supplementary Note 8: Experimental Conditions and Parameter Inference for STORM Validation

To reproduce the experimental conditions for STORM validation, we used data from the Zenodo repository "STORM Vectashield datasets (Tubulin)" [8]. Specifically, we extracted acquisition metadata and sample preparation details from the associated dataset description, which included the following:

* **Excitation wavelength:** 641 nm laser (Coherent, CUBE 640-100C)

* **Excitation intensity on sample:** ~1-2 kW/cm²

* **Camera:** Andor iXonEM+ (EMCCD) operated in conventional CCD mode

* **Pixel size:** Effective 100 nm

* **Objective:** Olympus 100x 1.3NA, with 1.6x magnification lens

All simulations used fluorophores modeled as Alexa Fluor 647 with known spectral and photo-physical properties from datasheets. While excitation rate $k_{01}$ was derived from the given parameters, some transition rates were not directly available, hence they had to be estimated:

* **Intersystem crossing rate $k_{1T}$:** Tuned within a plausible range (based on $\phi_{ISC} \in [0.001, 0.01]$) to match the characteristic dark-state durations observed experimentally.

* **Triplet state decay rate $k_{T0}$:** Adjusted to yield temporal sparsity consistent with recorded emission bursts.

* **Photobleaching rate $k_{1B}$:** Calibrated to reproduce the observed signal decay across 5000 frames.

We stress that the goal was not to fit data exactly, but to match the overall statistics and event dynamics under reported environmental and instrumental conditions. The camera model (including EMCCD noise) used the manufacturer's spec sheet Spec-Sheet_iXon.JPG, as detailed in Supplementary Note 3.

All resulting simulation parameters are available in Supplementary Table 5.
<table><caption>**Supplementary Table 5.** Experimental and photophysical parameters used for MUFASA-STORM simulation.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Description / Source</th></tr></thead><tbody><tr><td colspan="3"><strong>Experimental Setup</strong></td></tr><tr><td>Microscopy protocol</td><td>STORM</td><td>Single-molecule localization microscopy</td></tr><tr><td>Excitation wavelength λ</td><td>641 nm</td><td>Laser from Zenodo metadata [8]</td></tr><tr><td>Laser power density P</td><td>1,000 W/cm²</td><td>Conservative estimate (range: 1–2 kW/cm²)</td></tr><tr><td>Frame duration</td><td>36 ms</td><td>Exposure time per frame</td></tr><tr><td>Number of frames</td><td>2,000</td><td>Full acquisition duration</td></tr><tr><td>Camera</td><td>Andor iXonEM+</td><td>Effective pixel size: 100 nm</td></tr><tr><td colspan="3"><strong>Fluorophore: Alexa Fluor 647</strong></td></tr><tr><td>Molar extinction coefficient ε</td><td>239,000 M<sup>-1</sup>cm<sup>-1</sup></td><td>Thermo Fisher datasheet</td></tr><tr><td>Excited-state lifetime τ<sub>1</sub></td><td>1 ns</td><td>Typical for Alexa 647</td></tr><tr><td>Intersystem crossing coefficient α<sub>ISC</sub></td><td>1.6 × 10<sup>-3</sup></td><td>Calibrated for matching OFF durations</td></tr><tr><td>Triplet decay rate α<sub>nr</sub></td><td>0.1 s<sup>-1</sup></td><td>Tuned to reproduce dark-state returns</td></tr><tr><td>Cycles before bleaching</td><td>10<sup>9</sup></td><td>Effective; negligible bleaching during sim</td></tr></tbody></table>
## 9 Supplementary Note 9: Wasserstein Distance $W_1$

To quantitatively assess the performance of the MUFASA simulations, we consider a validation framework based on 1-Wasserstein distance ($W_1$) [12]. Given the interpretation of our measurements as realisations of stochastic processes, we use this quantity as a measure of discrepancy between the distributions of photon emissions from simulated and experimentally observed datasets.

Specifically, we compute $W_1$ on concatenated photon emission traces from multiple fluorophores. Namely, given $N$ molecules, the emission time series from each molecule (both simulated and experimental) are concatenated into two large arrays:

$$ \mathbf{X}_{\text{sim}} = [x_1^{\text{sim}}, x_2^{\text{sim}}, \dots, x_N^{\text{sim}}] \quad \text{and} \quad \mathbf{X}_{\text{exp}} = [x_1^{\text{exp}}, x_2^{\text{exp}}, \dots, x_N^{\text{exp}}], $$

where, for $i = 1, \dots, N$ each $x_i^{\text{sim}} \in \mathbb{R}^T$ represents the photon count time series for the $i$th fluorophore over $T$ timepoints. The independence of emission events ensures that concatenating photon emission traces is a valid computational procedure capturing the overall statistical behavior across multiple fluorophores.

More precisely, by considering the empirical measures

$$ \mu = \frac{1}{M} \sum_{i=1}^{M} \delta_{x_i^{\text{sim}}}, \quad \nu = \frac{1}{M} \sum_{j=1}^{M} \delta_{x_j^{\text{exp}}}, $$

both having uniform weights, the computation of the 1-Wasserstein distance reduces to consider the optimal-transport linear program:

$$ W_1(\mu, \nu) = \min_{\Gamma \in \mathbb{R}_+^{M \times M}} \sum_{i,j=1}^{M} \Gamma_{ij} C_{ij}, \quad \sum_{j} \Gamma_{ij} = \frac{1}{M}, \quad \sum_{i} \Gamma_{ij} = \frac{1}{M}, $$

where the **cost matrix** has entries

$$ C_{ij} = \|x_i^{\text{sim}} - x_j^{\text{exp}}\|. $$

A lower $W_1$ indicates closer alignment between simulated and experimental emission distributions, thereby validating the accuracy and realism of the simulation framework.

**Practical computation** To compute this distance efficiently, we use the Python Optimal Transport (POT) library [4], which implements a range of optimal transport solvers. Specifically, we employ the exact linear programming solver `ot.emd2`:

```python
import ot
wass_dist = ot.emd2(a, b, C)
```

* a and b are uniform histograms over the rows of $X_{\text{sim}}$ and $X_{\text{exp}}$ respectively.

* C is the cost matrix defined above.
# 10 Supplementary Note 10: Quantitative assessments via Wasserstein-distance of FF dynamics

We evaluated MUFASA's ability to replicate single-molecule emission statistics under the FF protocol by comparing its simulated photon count distributions to experimentally measured data. The experimental dataset was acquired using an EMCCD camera imaging EGFP-labeled samples under continuous illumination, yielding amplified digital signals representative of photon emission time series. Details of molecule preparation, acquisition protocols, and the single-molecule signal extraction pipeline are provided in [10] and in Supplementary Note 7.

To quantitatively assess the simulator's ability to reproduce real-world fluorescence dynamics, we computed the first-order Wasserstein distance, as described in Supplementary Note 9, between the distribution of experimentally observed photon counts per frame and the outputs of three simulation strategies: MUFASA, a standard Poisson process model, and SOFI Tools.

The Poisson process model was included as a baseline, as it represents the simplest counting process: photon arrivals modeled as a homogeneous Poisson process with a constant rate $\mu$, derived from theoretical considerations. Specifically, we used the closed-form excitation rate $k_{01}$ defined in Equation ??, and set $\mu = k_{01} \cdot \Delta t$, with $\Delta t$ being the duration of a frame. This model assumes memoryless photon emission with no sub-state variability, consistent with common assumptions in FF-SRM theory.

Parameter selection for all simulations, including MUFASA and SOFI Tools, was based on the experimental acquisition conditions described earlier. For the FF protocol, no additional parameter tuning was required beyond the known values of illumination power, wavelength, fluorophore properties (e.g., $\epsilon$), and quantum yield. These inputs were sufficient to determine the necessary transition rates for MUFASA via closed-form expressions. For the SOFI Tools baseline, which assumes a two-state Markov model with constant photon emission during the "on" state, we used values extracted from MUFASA simulations to ensure alignment. Specifically, one of its parameters is the emission intensity, which was set equal to the maximum value observed in the MUFASA-generated photon trace. This approach ensures a fair comparison, aligning both simulators with respect to photophysical timescales and emission amplitudes.

The Wasserstein distance provides a meaningful metric for quantifying the difference between two probability distributions, capturing both the distributional shape and the geometry of the data. In our context, we are comparing stochastic photon emission traces. Metrics such as Euclidean distance (e.g., $\ell_2$ norm) are not well suited for this task, as they treat the signals deterministically and penalize pointwise deviations without accounting for underlying statistical variability. Kullback-Leibler (KL) divergence can be undefined or infinite when the two distributions have non-overlapping supports, making it unstable for comparing empirical or noisy data [2]. Photon emission time traces from different simulators may yield distributions where one has a tail or concentration that the other doesn't, especially with noise. KL may diverge or be unstable, but Wasserstein gives a robust, interpretable difference in such cases.

Although MUFASA achieves strong agreement with the experimental data, its Wasserstein distance is slightly elevated. This can be attributed to the explicit modeling of EMCCD cameras, where the stochastic amplification of photoelectrons introduces additional variability. Notably, one can analytically demonstrate that the Wasserstein distance scales linearly with the amplification factor in EMCCD systems.

As shown in Supplementary Fig. 3, MUFASA achieved a Wasserstein distance of $W_1 = 2263$, closely aligning with the experimental data and marginally outperforming the Poisson model ($W_1 = 2269$). By contrast, SOFI Tools produced a substantially higher distance ($W_1 = 3065$),
indicating a poorer fit to the empirical distribution. This deviation reflects the limitations of assuming constant photon emission within on-periods, which fails to capture the intrinsic variability modeled by MUFASA. These results underscore MUFASA's capacity to simulate realistic photo-physical behavior, capturing both the stochastic switching and emission-level fluctuations observed in experimental FF data.

<figure><img src="image_15.png" /><figcaption>**Supplementary Figure 3. Comparison of FF simulations using the 1-Wasserstein distance.** Wasserstein distances ($W_1$) between experimental photon count distributions and three simulation models: MUFASA, Poisson, and SOFI Tools. Lower values indicate closer alignment with the empirical distribution of experimental data. MUFASA achieves the best match, highlighting its capacity to capture realistic emission variability under the FF protocol.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Supplementary Figure 4. Fluorophore grid layout.** Arrangement of 6,498 fluorophores on a 2048 × 2048 grid used for input simulation.</figcaption></figure>
FF Simulation Results

<figure><img src="image_17.png" /><figcaption>**Supplementary Figure 5.** MUFASA-Sim FF protocol output. Selected frames (0, 10, 50, 71, 92) from a 100-frame fluorescence simulation. **Left column:** “Emission” maps showing per-pixel photon emission rate in photons per frame (0–6 photons/frame color-scale). **Middle column:** “Blurred” images after applying the microscope point-spread function, in intensity units. **Right column:** “Observed by camera” output including shot and read noise, in camera digital units (ADU; 0–60,000 ADU scale). The sequence demonstrates continuous and spatially consistent photon emission and detection throughout the simulation protocol.</figcaption></figure>
STORM Simulation Results

<figure><img src="image_18.png" /><figcaption>**Supplementary Figure 6.** MUFASA-Sim STORM protocol output. Selected frames (0, 10, 50, 71, 92) from a 500-frame stochastic blinking simulation. **Left column:** “Emission” maps showing stochastic photon emission events in photons per frame (0–200 photons/frame color-scale). **Middle column:** “Blurred” images after applying the point-spread function, in intensity units (0–200). **Right column:** “Observed by camera” output including shot and read noise, in analog-to-digital units (ADU; 0–400 ADU scale). The output reflects characteristic blinking behavior typical of STORM imaging.</figcaption></figure>
PALM Simulation Results

<figure><img src="image_19.png" /><figcaption>**Supplementary Figure 7. MUFASA-Sim PALM protocol output.** Selected frames (0, 10, 50, 71, 92) from a 500-frame simulation of photoactivated localization microscopy. **Left column:** “Emission” maps showing photon bursts from activated fluorophores in photons per frame (0–1400 photons/frame color-scale). **Middle column:** “Detected” signals after applying the point-spread function, in intensity units (0–2500). **Right column:** “Observed by camera” output including noise and digitization, in analog-to-digital units (ADU; 0–2600 ADU scale). The simulation highlights sparse, stochastic activation consistent with PALM imaging.</figcaption></figure>
Blinking Simulation Results

<figure><img src="image_20.png" /><figcaption>**Supplementary Figure 8. MUFASA-Sim Blinking protocol output.** Selected frames (0, 10, 50, 71, 92) from a 100-frame blinking simulation. **Left column:** “Emission” maps showing frame-wise variation in emission in photons per frame (0–10 photons/frame color-scale). **Middle column:** “Blurred” images after applying the point-spread function, in normalized intensity units (0–14). **Right column:** “Observed by camera” output including camera noise, in analog-to-digital units (ADU; 0–380 ADU scale). The sequence demonstrates stochastic activation and deactivation of emitters across time.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Supplementary Figure 9. MUFASA Simulation runtime as a function of the number of fluorophores across four imaging protocols.** Computational runtimes for simulating PALM, STORM, FF, and Blinking protocols increase approximately linearly with the number of fluorophores. Times were obtained from actual simulations performed with 1,000, 2,000, 5,000, 10,000, and 15,000 molecules for each protocol under identical conditions. PALM and STORM exhibited shorter runtimes due to the early bleaching of the molecules, while FF and Blinking incurred longer runtimes.</figcaption></figure>
[1] Robert Alicki and Karl Lendi. *Quantum Dynamical Semigroups and Applications*. Lecture Notes in Physics. Springer Berlin, Heidelberg, 1 edition, April 2007. Hardcover ISBN: 978-3-540-70860-5, Published: 23 April 2007. Softcover ISBN: 978-3-642-08985-5.

[2] Martin Arjovsky, Soumith Chintala, and Léon Bottou. Wasserstein generative adversarial networks. In Doina Precup and Yee Whye Teh, editors, *Proceedings of the 34th International Conference on Machine Learning*, volume 70 of *Proceedings of Machine Learning Research*, pages 214–223. PMLR, 06–11 Aug 2017.

[3] Jurek W. Dobrucki and Ulrich Kubitscheck. *Fluorescence Microscopy*, chapter 3, pages 85–132. John Wiley & Sons, Ltd, 2017.

[4] Rémi Flamary, Nicolas Courty, Alexandre Gramfort, and et al. Pot: Python optimal transport. *Journal of Machine Learning Research*, 22(78):1–8, 2021.

[5] Arik Girsault, Tomas Lukes, Azat Sharipov, Stefan Geissbuehler, Marcel Leutenegger, Wim Vandenberg, Peter Dedecker, Johan Hofkens, and Theo Lasser. SOFI Simulation Tool: A Software Package for Simulating and Testing Super-Resolution Optical Fluctuation Imaging. *PLoS ONE*, 11(9):e0161602, 2016. PMCID: PMC5008722.

[6] M Hirsch, R.J. Wareham, M.L. Martin-Fernandez, M.P. Hobson, and D.J. Rolfe. A stochastic model for electron multiplication charge-coupled devices–from theory to practice. *PLoS One*, 8(1):e53671, Jan 2013. Epub 2013 Jan 31.

[7] Florent Luisier, Thierry Blu, and Michael Unser. Image denoising in mixed poisson-gaussian noise. *IEEE Transactions on Image Processing*, 20(3):696–708, 2011.

[8] Nicolas Olivier and Debora Keller. STORM Vectashield datasets (Tubulin). https://doi.org/10.5281/zenodo.7620025, 2023. Zenodo dataset.

[9] Sheldon M. Ross. *Introduction to Probability Models*. Academic Press, 11th edition, 2014. Chapter 6, pp. 371–375.

[10] Patrice Royal, Aurélie Andres-Bilbe, Patricia Ávalos Prado, Camille Verkest, Benjamin Wdziekonski, Samuel Schaub, Anne Baron, Florian Lesage, Xavier Gasull, Joshua Levitz, and Guillaume Sandoz. Migraine-associated tresk mutations increase neuronal excitability through alternative translation initiation and inhibition of trek. *Neuron*, 101(2):232–245.e6, January 2019. Epub 2018 Dec 17.

[11] Hui Tian and Abbas El Gamal. *Noise Analysis in CMOS Image Sensors*. Dissertation, Stanford University, Stanford, CA, 2000. Ph.D. dissertation.

[12] Cédric Villani. *Optimal transport: old and new*, volume 338. Springer, 2008.