2.1. SR of microscopy images

While optical SR techniques such as stimulated emission depletion (STED)⁽Reference Hell and Wichmann ^{7 )} and photoactivation localization microscopy (PALM)⁽Reference Betzig, Patterson, Sougrat, Lindwasser, Olenych, Bonifacino, Davidson, Lippincott-Schwartz and Hess ^{8 )} can break the diffraction limit, they are limited by the need for specialized equipment and complex sample preparation. Computational SR methods, on the other hand, while they rely on existing training data, can represent a cheap way to enhance the resolution of images acquired with conventional microscopes. Several approaches were developed to address the problem of computational SR of microscopy images. Deep-STORM (Stochastic Optical Reconstruction Microscopy)⁽Reference Nehme, Weiss, Michaeli and Shechtman ^{9 )} uses an encoder–decoder network to localize emitters in super-resolved images. Content-aware image restoration⁽Reference Weigert, Schmidt and Boothe ^{1 )} uses a U-Net architecture and is trained with low signal-to-noise ratio (input) and high signal-to-noise ratio (target) image pairs. These approaches produced interesting results, but recovered images often lacked high-frequency details due to the mean squared error (MSE) loss. To balance this issue, other approaches such as ANNA-PALM (Artificial Neural Network-based Photo Activated Localization Microscopy)⁽Reference Ouyang, Aristov, Lelek, Hao and Zimmer ^{3 )} are based on GAN. The authors trained a U-Net, and in contrast with Deep-STORM, a combination of pixel-wise reconstruction loss and adversarial losses is used to obtain image reconstructions of better quality. Using a similar network architecture and training loss, Wang et al.⁽Reference Wang, Rivenson, Jin, Wei, Gao, Günaydın, Bentolila, Kural and Ozcan ^{2 )} achieved SR in fluorescence microscopy across different modalities. Overall, using an adversarial loss results in output images that are sharper and of better perceptual quality and these methods demonstrated the potential of deep learning to improve the spatial resolution of fluorescence microscopy images. While the deep learning architectures and applications of these methods differ, they all require hard-to-get paired image data for training and offer a moderate increase in resolution of up to 4x.

It is essential to note that computational SR techniques effectively enhance the resolution of images by optimizing the use of available data in the original images. However, these methods are constrained by the existing information and cannot add details beyond what the original optical systems could capture.

2.2. Style-based generative models

StyleGAN (Style-based Generative Adversarial Network) models⁽Reference Karras, Laine and Aila ^{10 ,} Reference Karras, Laine, Aittala, Hellsten, Lehtinen and Aila ^{11 )} are well known for their ability to generate highly realistic images. The StyleGAN architecture consists of two subnetworks: a mapping network denoted by $ {\mathrm{G}}_m:{\mathrm{\mathbb{R}}}^d\to {\mathrm{\mathbb{R}}}^d $ and a synthesis network consisting of $ L $ layers denoted by $ {\mathrm{G}}_s:{\mathrm{\mathbb{R}}}^{L\times d}\to {\mathrm{\mathbb{R}}}^n $ . Here, $ d $ represents the dimensionality of the latent space. The mapping network takes a sample $ \mathbf{z}\in {\mathrm{\mathbb{R}}}^d $ from a standard normal distribution and maps it to a vector $ \mathbf{w}\in \mathcal{W} $ , where $ \mathcal{W} $ denotes the intermediate latent space. StyleGAN2 introduced path length penalty, which encourages a fixed-size step in $ \mathbf{w} $ to result in a fixed-magnitude change in the image. The regularization term is computed using the Jacobian determinant, and penalizing changes in it promotes the generation of smoother and more realistic images. $ L $ copies of the d-dimensional vector $ \mathbf{w} $ are fed to the $ L $ -layer synthesis network $ {\mathrm{G}}_s $ , with each copy representing the input to the corresponding layer of $ {\mathrm{G}}_s $ . The $ {\mathrm{G}}_s $ network controls the level of detail in the generated image at each layer. Individual modification of these $ L $ layers, by adjusting the latent vector copy of $ \mathbf{w} $ for each layer, extends the latent space into $ {\mathcal{W}}^{+} $ . This extended latent space enhances the model’s capability for accurate image reconstruction, which is vital for SR tasks, offering more nuanced control over the image generation process.⁽Reference Abdal, Qin and Wonka ^{12 ,} Reference Wulff and Torralba ^{13 )}

2.3. GAN-based high-resolution image reconstruction

The problem of obtaining a SR image of dimension $ n $ from a low-resolution (LR) image of dimension $ m $ is ill-posed as, for a non-invertible forward operator $ \mathrm{D} $ with $ m<n $ , there are infinitely many HR images that match a given LR image. Thus, the reconstruction procedure must be further constrained by prior information to better define the objective and lead to a stable solution. One such prior consisted of considering the reconstructed HR image to be part of a given domain. First, the distribution of HR images is learned in an unsupervised fashion, thanks to a GAN, and then the latent space of this trained GAN is searched to find a latent vector producing an HR image that, once downscaled, is the closest to the LR image input. GAN prior-based image reconstruction was first introduced by Bora et al.⁽Reference Bora, Jalal, Price and Dimakis ^{14 )} and further improved using StyleGAN⁽Reference Karras, Laine and Aila ^{10 ,} Reference Karras, Laine, Aittala, Hellsten, Lehtinen and Aila ^{11 )} by Menon et al. in PULSE (Photo Upsampling via Latent Space Exploration)⁽Reference Menon, Damian, Hu, Ravi and Rudin ^{15 )} by constraining the search to remain on the image manifold. To this end, PULSE and two other studies⁽Reference Wulff and Torralba ^{13 ,} Reference Zhu, Abdal, Qin, Femiani and Wonka ^{16 )} use an invertible transformation of the intermediate latent space $ {\mathcal{W}}^{+} $ _, which includes a leaky rectified linear unit (rectified linear unit (ReLU))⁽Reference Goodfellow, Bengio and Courville ^{17 )} followed by an affine whitening transformation, so that transformed latent vectors approximately followed the standard Gaussian distribution $ \mathcal{N}\left(\mathbf{0},{\mathbf{I}}_d\right) $ . Sampled vectors are then constrained to lie around a hypersphere with radius $ \sqrt{d} $ hypothesizing that most of the mass of a high-dimensional Gaussian distribution is located at or near $ \sqrt{d}{\mathcal{S}}^{d-1} $ , where $ {\mathcal{S}}^{d-1} $ is the $ d $ -dimensional unit hypersphere. Constraining samples to lie in a dense area of the StyleGAN style distribution resulted in increased realism of the generated images.

Although the above approach showed major improvements over previous work, it also presents important caveats that in practice led to image artifacts. Here, we show that transforming the intermediate latent space in this way does not lead to an accurate standard Gaussian distribution, and so prevents proper regularization based on this hypothesis. Moreover, as the search is strictly limited to the spherical surface $ {\mathcal{S}}^{d-1} $ , this limitation may prevent the search from reaching the closest reconstruction of the HR image.

In this work, we regularize the search in the latent space for a latent code located in “healthy” regions of the latent space. In this way, the system is constrained to produce images that belong to the original image domain StyleGAN was trained on. To do so, we take advantage of normalizing flow to Gaussianize the latent style sample distribution, which leads to a much closer standard Gaussian distribution. We then use this revertible transformation to regularize the search in $ {\mathcal{W}}^{+} $ such that it remains in a high-density area of the style vector distribution. We then show experimentally that the latter produces reconstructed images that are not only realistic but also more faithful to the original HR image.

4.1. Experimental settings

A critical aspect of our study was the necessity to work with high downscaling factors, such as 16x and 32x, which are central to our application scenario. These factors are not commonly used in the literature, and to the best of our knowledge, there are no experimental datasets available with such high downscaling factors. Most existing datasets typically focus on moderate downscaling factors of 2x or 4x, which is not aligned with the needs of our study, where we aim to address more extreme cases of resolution enhancement. Therefore, we had to generate our own LR images.

As for implementation details, we began by training a StyleGAN2-ADA (StyleGAN2 with Adaptive Discriminator Augmentation)⁽Reference Karras, Aittala, Hellsten, Laine, Lehtinen and Aila ^{23 )} model on a subset of the BBBC021 dataset,⁽Reference Caicedo, McQuin, Goodman, Singh and Carpenter ^{24 ,} Reference Ljosa, Caie and Ter Horst ^{25 )} which comprises cells treated with drugs acquired following the cell painting assay.⁽Reference Caie, Walls, Ingleston-Orme, Daya, Houslay, Eagle, Roberts and Carragher ^{26 )} That is, cells were fluorescently stained with markers for F-actin, B-tubulin, and deoxyribonucleic acid (DNA), as described in.⁽Reference Caie, Walls, Ingleston-Orme, Daya, Houslay, Eagle, Roberts and Carragher ^{26 )} The dataset comprises wide-field epi-fluorescence images, captured using the automated ImageXpress imaging platform. From this dataset, we extracted images of size $ 128\times 128 $ centered around each cell nucleus and used 400 images per compound treatment (approximately 28,000 in total) to train StyleGAN2-ADA.⁽Reference Karras, Aittala, Hellsten, Laine, Lehtinen and Aila ^{23 )} We then proceeded to generate 100000 random samples $ \mathbf{z}\sim \mathcal{N}\left(\mathbf{0},{\mathbf{I}}_d\right) $ and used their associated style vector $ \mathbf{w}={\mathrm{G}}_m\left(\mathbf{z}\right) $ to train the normalizing flow model. We opted for MAF, as it tends to perform better than RealNVP (Real-valued Non-Volume Preserving) for density estimation tasks. Our normalizing flow model comprised five flow blocks with all hidden dimensions set to 1024. For super-resolving the images using our RLS algorithm, we employed an Adam optimizer over 200 iterations with a learning rate of 0.5 and initialized the search using the mean of 10,000 randomly generated latent vectors.

To fine-tune the regularization parameters $ {\lambda}_w $ and $ {\lambda}_c $ , we empirically evaluated a spectrum of values, eventually setting them to 5e-5 and 0.01, respectively. This process revealed that the algorithm’s performance remained relatively stable across a broad range of these parameters, indicating a lack of sensitivity as long as the values fell within a specific boundary. Precisely, we observed optimal performance when $ {\lambda}_w $ was between 1e-6 and 5e-4, and $ {\lambda}_c $ ranged from 0.005 to 0.05. These findings suggest that while the exact values of $ {\lambda}_w $ and $ {\lambda}_c $ are flexible, maintaining them within these determined ranges ensures the algorithm functions effectively.

For evaluation, we used 100 random samples from various compound treatments, ensuring samples were not included in the training of the generative network. We simulated degraded images from these HR images using bicubic downsampling, which is a common choice in the literature. We then compared the results of our algorithm with those of Pix2Pix and two state-of-the-art unsupervised image reconstruction methods based on StyleGAN inversion: PULSE⁽Reference Menon, Damian, Hu, Ravi and Rudin ^{15 )} and BRGM.⁽Reference Marinescu, Moyer and Golland ^{18 )} For the baseline methods, we used the same parameter settings reported by the original papers.

4.2. Results

4.2.1. RLS recovers high-quality cell images

We reconstructed images using RLS and compared them to images generated with baseline methods. The results of this comparison are presented in Figure 2. At these high upscaling factors, due to the lack of proper regularization, the competing methods fail to produce reasonable cellular details and tend to produce images that can be accurately downscaled to the LR image at the cost of generating distorted cellular structures. In contrast, our approach produces highly realistic images for most examples and can reproduce cell phenotypes induced by compound treatment even if the fine-grain details are not similar. For instance, when reconstructing the images of cells treated with nocodazole, a known microtubule destabilizer, RLS captures the typically fragmented microtubule phenotype, resulting in a more accurate texture. Reconstructing images of cells at higher resolution obviously provides access to finer-grain measurable features.

Figure 2. Qualitative comparison of SR reconstructions on the BBBC021 dataset: visualizing the performance of RLS against baseline methods in reconstructing cellular structures and phenotypes under negative control (DMSO) and various treatment conditions at a 16x upscaling factor.

To quantitatively assess this gain in performance, we evaluate the SR images using Frechet inception distance (FID)⁽Reference Heusel, Ramsauer, Unterthiner, Nessler and Hochreiter ^{27 )} and kernel inception distance (KID)⁽Reference Bińkowski, Sutherland, Arbel and Gretton ^{28 )} to measure the discrepancy between the real HR images and the reconstructed one. As expected, the scores of FID and KID listed in Table 1 show that with both upscaling factors, our method significantly improves realism.

Table 1. Quantitative evaluation of RLS and baseline methods for SR on the BBBC021 dataset at 32x and 16x upscaling factors

Note: The best and the second best are highlighted in bold and italic, respectively.

4.2.3. RLS preserves interpretable features

RLS achieves a balance between realism and fidelity when super-resolving microscopy images. However, as any deep learning method for SR, it cannot generate details it would not have seen during training. Therefore, it cannot be used in a context where novel events can be expected. However, it could possibly be used to perform measurements to quantify expected phenotypic changes. This is the case for many assays used in basic research in biology. It is also the case in diagnostics such as parasitemia for instance where the tool must assess if a parasite is present or not. Instead of relying on a blind classification or regression of LR images, these tasks on LR images in such a context could be decomposed into two steps. A first step based on our deep learning approach would consist of reconstructing a HR image, while a second step would use a dedicated handcrafted analysis to quantify a phenotypic feature, making the analysis explicitly interpretable.

To evaluate to what degree the information conveyed by the SR images can be used for such quantitative assays and maybe later for interpretable diagnostics, we applied it to two assays that allowed straightforward quantification of “interpretable features” on HR images but not on LR images. Here, interpretable features refer to the measurable and explainable phenotypic changes crucial for biological assessment, such as the nucleocytoplasmic ratio in response to TNF- $ \alpha $ treatment and alterations in the Golgi apparatus morphology in reaction to nocodazole. Quantitative analysis was conducted on the SR images to assess the reconstruction of these features, and results were compared against those obtained from HR images as well as baseline methods, including BRGM, PULSE, and an unregularized latent search (“w/o Regu.”).

The first assay aimed to track the location of the nuclear factor kappa B (NF- $ \kappa $ B) protein within the cell. Upon treatment with TNF- $ \alpha $ , a pro-inflammatory cytokine, the protein moves to the nucleus, causing a shift in fluorescence signal from the cytoplasm to the nuclear area, resulting in bright green nuclei. We observed that the nucleocytoplasmic fluorescence ratio measured on 1000 HR-treated and HR-untreated images could be replicated when computed from super-resolution (SR) images (Figure 5a). We also provided some visual examples of super-resolving the LR images, that is, the first step in Figure 3, which shows that our method can reconstruct images of wild-type and treated cell images obtained with this common assay.

Figure 3. Visual examples of super-resolving translocation assay LR images at a 32x upscaling factor under negative control (DMSO) and tumor necrosis factor (TNF)- $ \alpha $ treatment conditions: left: LR image, middle: SR reconstruction, and right: GT.

The second assay we conducted aimed to monitor changes in the morphology of the Golgi apparatus. When cells are exposed to nocodazole, microtubules disassemble, causing the Golgi, originally located near the center of the cell, to break up into smaller stacks. Figure 4 shows that standard assays such as nocodazole-induced Golgi scattering (green) were reproduced with SR images. Moreover, as illustrated in Figure 5b, when computed only from SR images, a straightforward average spot size difference measured on 1000 HR-treated and 1000 HR-untreated images could be retrieved.

Figure 4. SR of images from a Golgi assay at a 16x upscaling factor under negative control (DMSO) and nocodazole treatment conditions. The left column shows the LR images, the middle column shows the SR reconstructions, and the right column shows the GT images.

Figure 5. Making interpretable measurements from LR images: first increasing the resolution by a SR method and then measuring a handcrafted interpretable feature. Each pair in the boxplots displays the distribution of handcrafted interpretable measurements, where the solid box represents the negative control (DMSO) and the dotted box signifies the positive controls (TNF- $ \alpha $ for translocation and nocodazole for Golgi), across various SR methods including RLS, BRGM, PULSE, and “w/o Regu.” alongside with the HR images for benchmarking. (a) Translocation ratio measurement: The y-axis quantifies the translocation ratio, an interpretable metric indicating TNF-induced NF- $ \kappa $ B translocation (green). The translocation ratio can be differentiated between two conditions not only from real HR images but also from SR images. (b) Mean spot area measurement: The y-axis quantifies the mean spot area, an interpretable metric indicating nocodazole-induced Golgi spreading (green), distinguishable between two conditions not only from real HR images but also from SR images reconstructed by our method.

Furthermore, we utilized cells treated with DMSO (Dimethyl Sulfoxide) as negative controls and cells treated with nocodazole (or TNF in the case of the translocation assay) as positive controls. We developed a classifier trained on HR images to distinguish between the two phenotypes and evaluated its performance on both HR and super-resolved images, ensuring that the datasets for GAN training and classifier training were distinct to maintain the integrity of our results. Additionally, we trained a separate classifier on LR images for phenotype differentiation, and its testing was confined to LR images. The results as shown in Table 2 indicate that the SR operation does not degrade the discriminative signal contained in LR images, and the classification accuracy is about the same on SR and LR images. Furthermore, we show evidence that this signal difference can be retrieved when quantifying the same interpretable feature variation present in HR images using handcrafted image analysis (see Figure 5).

Table 2. Comparison of classification accuracy for identifying phenotypic changes between negative control (DMSO) and positive control (TNF- $ \alpha $ and nocodazole and nocodazole conditions) in the translocation and Golgi assays, respectively, using super-resolved images and HR images as a benchmark

It is worth noting that in Figure 5 the similarity in performance between the “w/o Regu.” model and RLS during the translocation assay reveals an interesting insight. It suggests that, for this specific assay, the additional regularization does not significantly change the outcome. However, it is important to note that this assay was designed to evaluate strong phenotypic changes where the expected changes are substantial and easy to detect. In this case, the nonregulated version can also accurately capture the phenotypic changes, since the differences between the images are large. This is evidenced by the comparative classification accuracies in Table 2 for the translocation assay, where RLS achieves similar performance to the “w/o Regu.” version and the baseline methods. The advantages of regularization in RLS are more apparent in assays with subtler phenotypic variations, where detecting nuanced differences is more challenging. For example, in the Golgi morphology assay, where detecting changes requires more detailed analysis, RLS’s regularization reconstructs realistic images that more faithfully represent biological structures required for detailed quantitative evaluations. Furthermore, considering the overall performance, RLS demonstrates consistent and enhanced performance compared to both the baselines and the non-regularized version, as indicated by the superior classification results in Table 2. This suggests that RLS is a more robust method for super-resolving microscopy images, especially in assays where the phenotypic changes are subtle. Overall, the experiments confirm that the quality of super-resolved images is adequate for further analysis. Our method enables a two-stage process: employing deep learning for the challenging task of SR followed by a handcrafted, interpretable method for the subsequent quantitative measurements.

4.3. Robustness

In opposition to supervised methods that are sensitive to the input image domain, this approach is not restricted to a particular degradation operator that is used during training. To evaluate this aspect, we applied additional degradation operators such as Gaussian noise, salt and pepper, and Gaussian blur to a bicubic downscaled image DS before reconstruction. As depicted in Figure 6, the reconstruction closely matches the DS image. This result validates our choice of using the bicubic downscaling operator during training instead of more complicated specific degradation.

Figure 6. Evaluation of RLS performance under various degradation conditions: Images downscaled by bicubic method are further altered with Gaussian noise, salt and pepper noise, and Gaussian blur to assess the stability of the proposed method across a range of image perturbations (at a 16x upscaling factor).

4.4. Ablation study

Figure 7 demonstrates ablation experiments highlighting the impact of the components of our image prior. First, “w/o Regu.” searches the latent space without any regularization for the image that, once downscaled, matches the LR image. The second variant is denoted “w/o $ {\mathcal{P}}_w $ _,” that is, the image prior does not include the prior term $ {\mathcal{P}}_w $ . Similarly, “w/o $ {\mathcal{P}}_{\mathrm{cross}} $ ” refers to the suppression of $ {\mathcal{P}}_{cross} $ .

Figure 7. Ablation study showcasing the impact of regularization components on RLS performance, with qualitative results in the left column and quantitative results in the right column. The variants include “w/o Regu.” (searching the latent space without any regularization), “w/o $ {p}_w $ ” (the image prior does not include the prior term $ {p}_w $ ), “w/o $ {p}_{\mathrm{cross}} $ ” (the image prior does not include the prior term $ {p}_{\mathrm{cross}} $ ), and “RLS” (the full RLS model) (at a 16x upscaling factor).

To evaluate the three variants, we use the same set of parameters described in Section 4.1. One can see that searching the latent space without any regularization produces images that do not necessarily belong to the original image manifold and therefore do not appear realistic. It also tends to generate images that are accurately downscaled to the LR image but at the cost of generating distorted images when $ {\mathcal{P}}_{\mathrm{cross}} $ or $ {\mathcal{P}}_w $ is discarded. This implies that both priors play an important role in generating realistic details.

4.5. Uncertainty

An important challenge of the SR task is that it is an ill-posed problem. Although we can improve this aspect by using an image prior constraint, several closely related highly resolved images could still be consistent with a single LR image. To generate $ n $ realistic SR images, we sample $ n $ latent codes denoted as $ {\mathbf{w}}_1,{\mathbf{w}}_2,\dots, {\mathbf{w}}_n $ . We assume that their distribution follows a Gaussian distribution $ \mathcal{N}\left(\mu, \sigma \right) $ , where the parameter $ \sigma $ follows an inverse gamma distribution. Using the Bayes rule, we estimate the distribution’s parameters by

$$ \underset{\mu, \sigma }{\max}\left[\sum \limits_{i=1}^n\log p\left({\mathbf{w}}_i|y\right)-\frac{n}{2}\log {\sigma}^2-\frac{n}{2}\log \left(2\pi \right)-\frac{1}{2{\sigma}^2}\sum \limits_{i=1}^n{\left({\mathbf{w}}_i-\mu \right)}^2+\log p\left(\sigma \right)\right] $$

Here, the first term is the log-likelihood of the posterior distribution of $ p\left(\mathbf{w}|\mathbf{y}\right) $ , which is defined in Eq. (5); the second term is the regularization term, which penalizes large values of $ \sigma $ _; and $ \log p\left(\sigma \right) $ is the log prior distribution of $ \sigma $ . Figure 8 displays multiple solutions for a given LR image we can obtain with our approach. Sampling close variations of HR images from a single LR input can be used to enhance the robustness of the estimation of an image-based quantitative feature by reducing the effects of noise or artifacts that may exist in the input LR image.

Figure 8. Visualizing the ill-posed nature of the SR task and the uncertainty associated with SR reconstruction. Five distinct SR images (SR₁ to SR₅) are generated for each LR image by sampling five different latent codes from the latent space, alongside the GT (at a 16x upscaling factor).