© The Author(s), published by EDP Sciences, 2025

1 Introduction

Holographic tomography is a powerful label free, quantitative imaging technique that allows retrieving the 3D Refractive Index (RI) distribution of the probed samples, which represents the most complete information about the sample of interest [1, 2]. By acquiring multiple digital holograms relative to different orientations between the sample and the probing light, it is possible to obtain the RI profile extracting from them the Quantitative Phase Maps (QPMs) and processing them with opportune numerical solvers. To achieve the relative motion between the light source and specimens, two approaches can be used: illumination rotation or sample rotation [2]. In the first, the sample is kept in a stable position while the probing light beam is moved, in the latter the light source is kept in a fixed position while the sample is rotated. The first approach allows to finely control the position of the light and the rotation angles but might reduce the throughput due to the need to physically place the sample on the support, and the impossibility to image cells in suspension conditions. The second approach allows to increase the throughput since there is no need to fix the position of the channel but does not allow to control the rotation angles a priori, which must be retrieved in postprocessing. Sample rotation paradigm has recently been demonstrated for DH in a technique referred to as holo-tomographic flow cytometry [3]; the sample’s rotation is obtained as it flows into a microfluidic channel by creating a gradient in the flow speed. Since both the rotation paradigms discussed above are prone to numerical errors, it is of crucial importance to possess numerical models to be used as ground truth for the development of robust tomographic solvers, denoising, image enhancement, and 3D segmentation algorithms. The problem of lack of universal ground truths has been long known in microscopic imaging, as it is often impossible to have access to the true morphology of the samples of interest [4]. Therefore, often the results have to be cross validated using different imaging techniques or statistics to prove the correctness of the obtained reconstructions [5]. In the case of DH, the verification of experimental results was traditionally made exploiting very simple objects such as uniform micro spheres [6–8] or simple numerical models [9–11], until more sophisticated methods have been proposed to print a test target through photolithography [12–14]. While this method is perfectly suitable for illumination scan rotation approaches, it is more difficult to conjugate it with in-flow tomography based on the rotation of suspended cells, which still strongly relies on numerical phantoms [15]. Great attention must be therefore put into the design of a reliable numerical phantom, able to represent the specimen in the most reliable way and tackle the reconstruction strengths and drawbacks of different solvers, taking into account also the optical properties of the chosen imaging system [16].

In this paper, we propose a method to simulate the numerical phantom of a cell suitable for different cell lines at different levels of complexity. Our objective is to model the 3D Refractive Index distribution of important classes of cells and some internal organelles, with the aim to use this flexible cellular model to tailor algorithms relative to the entire holo-tomographic imaging pipeline. In particular, we apply the method by simulating 1) the RI distribution of a yeast cell with cytoplasmic vacuoles inside, based on the studies reported in [17], 2) of a monocyte with Lipid Droplets (LDs) based on the findings reported in [18], 3) of a monocyte with vacuoles, merging the data from the first two simulations and 4), of a monocyte with lysosomes, based on the studies reported in [19–21]. As will be discussed in the next sections, the choice of this selection of samples is linked to their significance for basic research in biology and diagnostics and is due to the availability of prior information from experimental tests carried out by different groups. The whole method is however applicable to a larger variety of cells and organelles, provided that prior statistical information is available on the morphometry and RI of the sub-compartments to be simulated. We share the simulated dataset of 200 cells in suspension with the aim of providing a tool for the development of phase-contrast tomography algorithms in the absence of a ground-truth.

1.1 Selected case study cells

Yeast cells are important eukaryotic cells used as a model for the study of cytoplasmic vacuoles. At any stage of their life cycle, they are characterized by the presence of a large roundish vacuole along with the potential presence of much smaller vacuoles with low RI values. In mammalian cells, e.g. in monocyte cells of the bloodstream, the presence of various phenotypes of cytoplasmic vacuoles (that can vary largely in number and size) is considered as a biomarker for several disfunctions and pathological conditions, e.g. in lysosomal storage diseases (LSDs) [22], SARS-CoV-2 [23] and other viral infections [24–26], and cancer [27–31], to name a few. Thus, there is a strong interest in phenotyping patterns of cytoplasmic vacuoles and lysosomes inside these cells, some of them challenging the optical resolution of holo-tomographic flow cytometry systems. Since it is possible to experimentally control the volume and the number of cytoplasmic vacuoles, engineering the yeast cells acting on the osmotic conditions of the surrounding buffer, or in the case of mammalian cells by using external drugs to simulate this phenotype, it is important to implement a simulation method capable of describing all these conditions. Accordingly, particular attention was put in the simulation of cytoplasmic vacuoles. Instead, LDs are intracellular compartments mainly located in the cytoplasm [32]. LDs represent small structures with high RI values linked to several regulatory mechanisms in living cells, both in physiological and pathological states. Generally, they present a high numerosity. They were initially recognized only as storage organelles, but it has been discovered that they play multiple roles in maintaining intracellular homeostasis and that they dynamically interact with various organelles, such as mitochondria and lysosomes, thus affecting cellular processes [33]. For this reason, the involvement of LDs has been recognized in several pathologies, such as cancer [34], diabetes [35], or neurodegenerative diseases [36]. In immune cells (e.g., monocytes) LDs are markers of viral infections [37] and inflammations [38], as their number and size increase in response to these stimuli.

Aside from the specific organelles modeled, the simulation of monocytes is of primary interest in general as, like all the other cells of the blood stream, they naturally live in suspension. Consequently, this class of cells is the one that benefits the most from the use of a holo-tomographic apparatus working in flow cytometry mode, for which the absence of a ground-truth is a non-negligible impairment. The statistics of these simulations were studied and extracted from the literature and from data measured by the holo-tomographic flow cytometry system on suspended cells.

2 Methods

We developed a versatile algorithm for generating numerical phantoms of suspended cells, which we then tailored to simulate the unique characteristics of various cell lines as described in the previous sections. Accordingly, four representative simulations are here discussed: yeast cells with cytoplasmic vacuoles (Simulation 1), monocytes with LDs (Simulation 2), monocytes with vacuoles (Simulation 3) and monocytes with lysosomes and cytoplasmic vacuoles (Simulation 4). For each simulation, different classes of substructures were considered. Three substructures are common to each simulation, in particular: external membrane, cytoplasm and nucleus. In addition, vacuoles are simulated for Simulations 1, 3 and 4, LDs for Simulation 2 and lysosomes for Simulation 4.

2.1 Numerical simulation of a phantom cell

The first simulative step consists in creating a box region of 200 × 200 × 200 pixels, with pixel size 0.125 μm, where the cell will be inserted. To obtain a model as realistic as possible, the dimension and the shapes of the overall cell and of each intracellular structure were varied for each realization and drawn from distributions of real samples [17–21, 39, 40]. Figure 1 sketches the steps to generate the structures common to all cell types. The external shape of the cell is obtained by drawing its volume and sphericity index from a realistic distribution [17, 18, 39, 40] and adapted to the physical dimension of the cells of interest (yeast cell [17] or monocyte [18]), as shown in Fig. 1a. Concerning the yeast cells, they occupy around 1/3 of the modelled volumetric box region. Monocytes are bigger and occupy around 60% of the simulated volume. Once those values are calculated, the cell is initially modeled as an ellipsoid whose semiaxes in a Cartesian coordinate system are labelled as r_x, r_y and r_z. The r_x value is drawn from a Gaussian distribution describing the mean distance of the cell center from the membrane, while r_y and r_z are univocally determined by considering the volume of the cell, the sphericity index and imposing the relation r_x < r_z < r_y. Then, the ellipsoid modelling the cell is numerically altered to obtain an irregular, more realistic external shape, as shown in Figure 1b. This numerical operation is performed generating a second ellipsoid with the same procedure described, rotating it into the 3D space of a randomly chosen angle and performing a Boolean or operation between the two ellipsoids. This allows obtaining tiny protuberances in the external shell which are generally observed in real samples (that are not perfect spheroids). Subsequently, the cytoplasm region is modeled in a similar manner. The values of its three semiaxes are computed multiplying the corresponding values of the overall cell for a factor K which is randomly determined for each semiaxes and for each realization of the simulation, falling into the interval [90%, 98%], as shown in Figure 1c. Setting the cytoplasmatic semiaxes to a smaller value compared to the ones of the overall external cell allows obtaining a leftover border modelling the external membrane. Concerning the nucleus, its center is randomly placed into the cytoplasmic region and varied over the different realizations of the simulated dataset, sketched in Figure 1d. Its dimension is determined by first choosing its sphericity index, values of minor and major principal axes, and nucleus to cell ratio in accordance with the known values. The latter parameter is particularly important for yeast cells as they are characterized by important changes in nucleus dimensions depending on the cell life stage. Consequently, the volume of the nucleus is calculated in accordance with the volume of the cell through the nucleus to cell ratio. The third principal axis is then univocally determined. Finally, after defining the geometry of the internal structures, their RI values are assigned according to the reference values labeled in the relative literature. The procedure described above is used to model the three structures common to all the simulations. The geometrical modelling of cytoplasmic vacuoles, LDs and lysosomes will be discussed in detail in the next sections.

Fig. 1 Conceptual steps common to all the proposed simulations. a), Schematization of the modelling of the external shape. V = Volume, SI = Sphericity Index, rx, ry and rz: semiaxes of the simulated ellipsoid. b), creation of the irregular external shape. c) Simulation of the cytoplasm. rx,cyto, ry,cyto and rz,cyto: semiaxes of the ellipsoid simulating the cytoplasm. d), Simulation of the nucleus. rx,nuc, ry,nuc and rz,nuc: semiaxes of the ellipsoid simulating the nucelus, NCR = Nucleus to Cell Ratio.

2.2 Simulation of cytoplasmic vacuoles

Each vacuole is modeled as an ellipsoid whose center is randomly allocated into the cell, respecting the experimentally observed statistics of vacuole placement. Indeed, depending on the cell stage and on their sizes, vacuoles tend to occupy different positions in the cytoplasm of the cell. Then, the sphericity index of the vacuole is drawn from the statistics in [17]. The most important factor used for the parametrization of the vacuole is its equivalent diameter. By selecting this value and combining it with the sphericity index, the three semiaxis of the ellipsoid vacuole r_xv, r_yv and r_zv are calculated.

As described in [17], the vacuole region is characterized by a higher value of RI in its membrane region, which we refer to as perivacuolar region, and then by a value of RI gradually decreasing up to the center of the vacuole. To model this, we first randomly generate the value of RI of the perivacuolar region RI_perivac and of the central vacuolar region RI_vac considering the reported statistics. Then, the internal area is modelled through concentric ellipsoidal shells. Once the number of shells, N, is selected, depending on the equivalent radius of the vacuole, their geometry is described as reported below.

• The most external shell (n = 1) coincides with the generated ellipsoid representing the entire vacuole. In formulas, r_x1 = r_xv, r_y1 = r_yv and r_z1 = r_zv.

• The semiaxes of the most internal shell (n = N) are equal in size to the 10% of the ones of the most external shell. In formulas, r_xN = 10%r_xv, r_yN = 10%r_yv and r_zN = 10%r_zv.

• An equally spaced vector of shells’ semiaxes is generated for each dimension, ranging from the calculated outermost and the innermost ones, and then assigned to each shell. In formulas s_x = linspace(r_xN, r_x1, N), s_y = linspace(r_xN, r_x1, N), s_z = linspace(r_xN, r_x1, N), where s_x is the vector containing the x semiaxes of each shell, s_y the vector containing y semiaxes and s_z the vector containing the z semiaxes. The MATLAB^® function linspace [41] returns a vector of N elements equally spaced from the first argument (r_xN in this case) to the latter (r_x1 in this case).

• The value of RI_n assigned to each shell n comes from a Gaussian distribution of mean RI_mean,n and standard deviation std_n.

An equally spaced array of mean RI value is created for each shell, ranging from the RI assigned to the peri-vacuolar region to the one assigned to the internal vacuole. In formulas, A = {RI_mean,N,…, RI_mean,1} = linspace(RI_perivac, RI_vac, N)., where A is an array storing all the mean RI values for each shell.

• A normal distribution with mean equal to 1 and standard deviation equal to 0.03 is used to generate a standard deviation value for each shell. In formulas, B = {std_N…std₁}, where B is an array storing all the standard deviation values for each shell.

• The final RI value for each shell n RI_n is drawn from a Gaussian distribution with mean RI_mean,n and standard deviation std_n.

• Finally, a smoothing operation is applied to the entire vacuolar region to smoothen the sharp variations between the shells followed by a volumetric erosion to reduce the width of the peri-vacuolar region. This erosion is performed through the MATLAB^® function strel [42] with disk-shaped structuring element of radius r. r was manually tuned and set as 22% of the equivalent radius of the vacuole.

Figure 2 sketches the geometrical structure of the proposed vacuole before the smoothing operation.

Fig. 2 schematical representation of the geometric structure used for modelling the cytoplasmic vacuole. N shells are modeled, each of them characterized by x and y dimensions rx and ry and a Refractive Index value RIn.

2.3 Simulation of lipid droplets and lysosomes

LDs and lysosomes are modeled in similar way. Therefore, for the sake of simplicity, in this section we describe the procedure for LDs only unless differently specified. As for the case of the vacuoles, LDs are modeled setting their equivalent radius as the principal parametric parameter. Like all the organelles discussed, LDs are modeled as ellipsoids. Combining the value of equivalent radius with the sphericity index and the volume, as discussed in [18], the three semiaxes are determined. The number of LDs simulated is chosen considering the ratio between the LDs aggregate volume and the volume of the overall cell. Knowing the average volume of each LD, the chosen number of LDs inserted in the cell is 4. Differently, lysosomes are generally more numerous in the cells, so a number between 10 and 13 is simulated here. As a last geometric consideration, it was noted in ref. [18] that the LDs tend to aggregate in a small area which is displaced from the centroid of the cell. Therefore, this specification is incorporated in the placement of the centroid of each LD. Concerning the assignation of RI for LDs, its mean value and standard deviation were reported in ref. [18]. The statistics for lysosomes are reported in refs [19–21]. However, it was proven that the RI profile is not constant in the LD volume but increases going towards the center. Accordingly, we exploited a simulation mechanism similar to the one used for the vacuoles (sketched in Fig. 2), where three shells were exploited to model the RI profile:

• The most external shell coincides with the generated ellipsoid representing the entire LD.

• The semiaxes of the most internal shell are equal in size to the 30% of the ones of the most external shell.

• The semiaxes of the median shell are calculated as the mean value between the internal and the external.

• The RI value of the most internal shell is drawn from a Gaussian distribution with mean RI equal to the one of the overall LD and standard deviation = 0.001.

• The RI value of the most external shell is equal to the 70% of the one of the inner shell.

• The RI value of the middle shell is the mean value between the internal and the external ones.

• Finally, a smoothing operation is applied to the entire LD region to smoothen the sharp variations between the shells.

Differently from the vacuole simulation, a smoother profile of RI is simulated, as the presence of a sharper LD membrane was not detected in ref. [18], thus explaining the presence of only three shells and the absence of the erosion needed to simulate the peri-vacuolar membrane. Finally, the mean value of the LD’s RI profile is checked to ensure it falls into the expected ranges.

3 Results

In this section we report the results for the performed simulations and the parameters used for the creation of the synthetic datasets. For all the simulations, Matlab^® 2024b was used, running on a conventional Desktop Computer equipped with 32GB RAM. Concerning the chosen RI values, in Simulation 1, since an accurate statistic was available in the literature for each of the internal structures and for each stage of division of yeast cells [17], the values of the modeled RI where varied for each realization, drawn from opportune distributions. On the other hand, in Simulations 2, 3 and 4 the main interest was set into LDs, vacuoles and lysosomes, therefore the RI values of the external membrane, cytoplasm and nucleus where not varied in the reported examples. Table 1 summarizes the type of cells and substructures included, and the RI values simulated. The simulation time for a single cell of this geometry was evaluated in Simulation 4, which was the most time-consuming due to the monocyte’s larger size and the number of internal structures simulated compared to the yeast cell. To simulate 10 substructures (5 vacuoles and 5 lysosomes), the average time required was 60 seconds. Simulating 15 substructures (5 vacuoles and 10 lysosomes) increased the average time to 2 min.

In Figure 3, the simulation for two different yeast cells is reported, (Simulation 1). As can be noticed by the scalebar, yeast cells are small with respect to the majority of mammalian cells. In Figures 3a–3d, a yeast cell with one single vacuole of large dimensions (i.e. with a volume occupying approximately 10% of the cell’s volume) is depicted. The cell has a quasi-spherical external shell. In Figures 3a and 3c the central slice (cuts of the simulated cell in the plane yz and xy respectively) of the simulated 3D RI distribution are reported. The corresponding isolevel representation is shown in Figure 3d. In Figure 3b it is possible to visualize the cytoplasmic vacuole with the peri-vacuolar membrane, defined as in [17], whose structure can be observed in detail in the RI profile along the red dotted line shown in Figure 3a. The panels in Figures 3e–3h report the simulation of a yeast cell characterized by a moderately large roundish vacuole (i.e. with a volume occupying approximately 7.5% of the cell volume) and a multitude of smaller vacuoles (i.e. with a volume occupying approximately 0.3% of the cell volume each) in the peripheral region (Fig. 3h). The external shell of the cell is non symmetric. The difference in size and the internal structures can be noticed from peripheral yz and xy slices of Figures 3g and 3e respectively, and the corresponding RI profile along the red dashed line, shown in Figure 3f.

Fig. 3 simulation results for (a–d) a yeast cell with a large vacuole inside (Visualization 1) and (e–h) a yeast cell with a large vacuole and a multitude of smaller vacuoles (Visualization 2). (a, g) Visualization of a slice in the yz plane from the 3D RI profiles; (b, f) Cut along the red dotted line in subfigures (a, e), highlighting the RI profile into the vacuole; (c, e) Visualization of a slice in the xy plane from the 3D RI profiles; (d, h) Isolevel visualizations. EM = External membrane, C = Cytoplasm, N = Nucleus, V = Vacuole. For each cell, the isolevels’ thresholds for the visualization are set with respect to the maximum value of the simulated RI profile.

Figures 4a and 4b reports the isolevel visualization of two different yeast cells in the late-G1 stage, according to the nomenclature reported in [17] that describes the different stages of the cell budding process. This simulation is aimed at showing the possibility of including the case of a cell that does not match exactly the common approximation of spheroidal external shell. Hence, eventual asymmetries in the external shell are simulated. Accordingly, the simulation is performed setting the simulation in agreement with the statistics reported for the set of experimental results obtained in the case of late-G1 stage [17]. An ensemble of 50 cells of late-G1 yeast cells is simulated and histograms of relevant statistics are computed. In particular, Figure 4c shows the histograms relative to global parameters such as the overall volume of the cell and its equivalent diameter; the average value of RI and the value of the three principal axis of the cell. In addition, the nucleus’ volume, equivalent diameter and average RI histograms are reported, considering that, as mentioned in Section 2, the nucleus’ geometry is parametrized starting from the available value of sphericity index for a late-G1 cell. Finally, histograms of the vacuole’s sphericity index, average RI and average RI of the perivacuolar membrane are reported, considering in this case that the vacuole geometry is parametrized over the equivalent diameter, which is fixed for the late-G1 class of cells, thus providing a constant equivalent volume. As it can be inferred from the set of histograms, the simulated statistics, also summarized in Table 2, are in good accordance with the tables reported in [17], which confirms the reliability of the simulation.

Fig. 4 simulation results for late-G1 yeast cells. (a, b) isolevels of two different simulated cs; (c) Histograms of some relevant statistics evaluated over a population of 50 cells. For each cell, the isolevels’ thresholds for the visualization are set with respect to the maximum value of the simulated RI profile. The statistics are in good accordance with the ones in the literature.

Figure 5 reports a similar analysis for the case of monocytes with LDs (Simulation 2). In this case the simulation is performed setting the RI values for the external membrane, cytoplasm and nucleus as described in Table 1 and the ones for the LDs as in [18], Table 1. As can be observed in the reported isolevel representation (Figs. 5a and 5b), LDs tend to aggregate into a reduced portion of the cell’s volume, a behavior observed in [18] and that we included in our modelling. From this visualization, their volume can be quickly compared to the volume of the overall cell (each LD occupies about 0.006% of the whole cell volume). The slice-by-slice visualization of Figure 5b allows visualizing the spatially variating RI distribution of each LD. Also in this case, statistics of relevant parameters relative to LDs are evaluated for a population of 50 simulated cells (sphericity index, average RI and RI standard deviation, reported in Fig. 5c) and show a good accordance with the values reported in [18].

Fig. 5 results for the simulations of a monocyte with LDs (Visualization 3). (a) Isolevel. From this visualization, the size and the positioning of the LDs can be seen, compared to the overall cell; (b) yz cut of the simulated 3D RI distribution. The non uniform distribution of RI into the LD can be noticed; (c) Histograms of the sphericity index, average RI and standard deviation values are evaluated for the LDs in a population of 50 monocytes. μ = mean value. σ = standard deviation. EM = External membrane, C = Cytoplasm, N = Nucleus, LD = Lipid Droplets. For each cell, the isolevels’ thresholds for the visualization are set with respect to the maximum value of the simulated RI profile.

The results for the case of a monocyte with simulated cytoplasmic vacuoles (Simulation 3) are presented in Figure 6. In this case, small vacuoles are simulated, in a numerosity ranging from 7 to 10 for each cell. The equivalent radius of the vacuole was drawn from a Gaussian distribution with mean 0.8 μm and standard deviation 0.8 μm, since those structures are small compared to the overall volume of the cell. Concerning their mean RI value, it was drawn from a Gaussian distribution with mean 1.360 and standard deviation 0.008. Figure 6a presents an isolevel and a peripheral yz slice from the cell, where the asymmetry of the cell’s external shell can be noted, together with the placement of the vacuoles and of their non-uniform RI distribution.

Fig. 6 Results for the simulations of monocytes with cytoplasmic vacuoles and monocytes with cytoplasmic vacuoles and lysosomes. Top panel, (a), an isolevel visualization and a yz slice from a simulated cell, highlighting the positioning of the vacuoles and their spatially-varying RI distribution (Visualization 4). Bottom panel, (b), an isolevel visualization and an xz slice from a simulated cell, where the high RI and numerosity of lysosomes can be appreciated (Visualization 5). EM = External membrane, C = Cytoplasm, N = Nucleus, V = Vacuoles, L = Lysosomes. For each cell, the isolevels’ thresholds for the visualization are set with respect to the maximum value of the simulated RI.

The results for the last case (Simulation 4) are reported in Figure 6b. In this case monocytes with 5 vacuoles and multiple lysosomes were modelled (10–13). The RI values of the vacuoles are set up as in Simulation 3 and their equivalent radius is 1 μm. The RI values of the lysosomes are drawn from a Gaussian distribution with mean 1.46 and standard deviation 0.008 as reported in [19–21], while their equivalent diameter is drawn from a Gaussian distribution of mean 0.6 and standard deviation 0.3, as reported in [19–21]. Figure 6b depicts an isolevel and an xz slice of a simulated monocyte with 13 lysosomes, where the difference in size between lysosomes and vacuoles can be noted.

3.1 Simulation of digital holograms

As mentioned in the previous sections, having a tool allowing to simulate a ground truth distribution of refractive index is useful to tune all the steps of an optical imaging pipeline. In particular, it is possible to numerically reproduce the digital holograms by simulating the interference between an object beam propagating through the simulated RI distribution [43] and a reference beam. The results are presented in Figure 7, where the hologram was simulated for an example taken from Simulation 1, a yeast cell with a round vacuole inside. An off-axis digital hologram was simulated. The wavelength is 495 nm, the pixel size of the camera is 4.5 μm, the numerical aperture of the objective is 0.95 and the lateral magnification 76.4. The size of the simulated hologram is 5124 × 5124 pixels. The inclination of the fringes is of 45° in the hologram plane. Through this simulation it is possible to visualize also the Fourier spectrum containing the three holographic diffraction orders.

Fig. 7 Left: Numerical digital hologram obtained starting from a simulated yeast cell with a vacuole (Simulation 1). The interference fringes are highlighted in the yellow box. Right: corresponding simulation of the amplitude Fourier spectrum, in which the holographic diffraction orders are clearly visible.

3.2 From the simulated cell to its simulated reconstructed tomogram

In this section, example use-cases for the presented cell models are reported. Concerning the nomenclature, we herein refer to the 3D RI cellular model proposed as “simulated cell”, which is the high-fidelity representation of the 3D RI distribution of the sample of interest. Instead, the tomogram is the simulated reconstruction of the 3D RI distribution, after the cell has been probed through an arbitrary optical system and reconstructed through an arbitrary inversion algorithm.

The tomogram hence embeds the effects of the optical setup, noise, and the specific algorithms exploited for each reconstruction step, including the tomographic inversion algorithms. All of these aspects determine a difference between the “simulated cell” and its relative “tomogram”, which is a degraded version of the cell after the passage through a non-ideal holo-tomographic measurement-reconstruction system. It is worth pointing out again that the aim of our work was to propose a high-fidelity ground truth cellular model which could be used to tune different algorithms relative to the holo-tomographic acquisition-reconstruction pipeline. Depending on the specific holographic setup or tomographic paradigm employed, the demonstrative and simplified use-cases discussed in this section can be modified to suit the one of interest, for example embedding the measured PSF of the setup or exploiting different tomographic solvers to reconstruct the tomogram.

In the reported examples, we consider the effect of an isotropic system PSF operating as a 3D low-pass filter on the cell. This is reasonable in all the cases of in-flow holo-tomography where isotropic resolution is achievable. Then, we consider the effects of speckle noise as well as the use of the tomographic inversion algorithms. Two different classes of inversion algorithms are considered here, namely the filtered backprojection [44] and the filtered backpropagation with Rytov approximation [45].

After the ground-truth 3D RI distributions are simulated, the system PSF, speckle noise, and the influence of the employed inversion algorithm have to be considered. Here, a simplified general method is considered, which allows to simulate and test each step of the holo-tomographic acquisition-reconstruction pipeline [3, 46, 47]. This is summarized in the sketch of Figure 8a. The starting point is the cell simulation we propose in this work. Convoluting this model with the PSF of the optical system in use, the simulated tomogram model is obtained. Reprojecting this model along the available tomographic angles allows simulating the QPMs, which can be eventually corrupted by noise [48]. Finally, composing these QPMs with a suitable tomographic solver [48] allows computing the reconstructed, simulated tomogram. Figure 8b represents each step of the considered pipeline applied to our in-flow tomographic flow cytometry optical setup for a sample simulated cell.

1. A RI distribution of interest, C, of size n × n × n is simulated according to the presented numerical models. Figure 8b1 sketches the central slice of C.

2. The PSF of the optical system is modeled through a 3D Gaussian kernel G of size 5 × 5 × 5 pixels and standard deviation 2. These values where experimentally chosen for this specific example to obtain a loss in the resolution visually similar to the one observed in real reconstructed tomograms. Figure 8b2 depicts the central slice of the Gaussian kernel.

3. The simulated tomogram T of size n × n × n is obtained as the convolution between the modeled PSF and the RI distribution C. In formulas, T = C ⨂ G, where ⨂ is the convolution operator. Figure 8b3 sketches the central slice of T.

4. The QPMs are obtained by reprojecting T along the available projection angles. In formulas Q = ψ{T, theta}, where ψ{…} is the operator modelling the projection, theta is the set of projection angles with size 1 × t and Q is the stack of resulting QPMs of size n × n × t. In our case, the operator ψ{…} is obtained through the MATLAB^® radon function [49]. Figure 8b4 depicts one of the obtained QPMs.

5. The noisy QPMs, QPM_noise, are obtained by adding a suitable noise distribution N over the uncorrupted QPMs. In our case, the noise is modeled as a distribution of multiplicative speckle noise with mean 1 and variance 5e–6. This was made resorting to the imnoise MATLAB^® function with “speckle” specification [50].

6. The simulated reconstructed tomogram $T_{\mathrm{rec}}$ is obtained by inverting the operator ψ{…}, composing all the available projections and the associate projections angles with a suitable physical criterion. In our case, this step is performed resorting to the MATLAB^® iradon routine in the first example [44] and using a custom filtered backpropagation routine using the Rytov approximation in the second example [45]. In formulas, $T_{\mathrm{rec}}= {\psi }^{-1} \left\{{\mathrm{QPM}}_{\mathrm{noise}}, \mathrm{theta}\right\}$. Summing up, the simulated reconstructed tomogram is obtainable from the simulated cell as:

Fig. 8 General sketch of a holo-tomographic acquisition-reconstruction pipeline. (a) Most important steps are schematically depicted; (b) Images corresponding to each step particularized to our optical setup and reconstruction pipeline are discussed (Visualization 6). For the sake of visualization clarity, a 13 × 13 × 13 Gaussian kernel is shown. In this case an equally spaced angular sequence was used, with the angles ranging from 0 to 360° with a step of 5°, recalling an in-flow tomographic reconstruction pipeline.

$$T_{\mathrm{rec}}= {\psi }^{-1} \left\{\psi \left\{C⨂G, \mathrm{theta}\right\}+N, \mathrm{theta}\right\}$$(1)

Figure 8b6 depicts the central slice of the simulated reconstructed tomogram T_rec.

Figure 9 depicts the application of this pipeline to the cells simulated and discussed in this work. The top row shows an yz slice from the simulated cells in Figures 3, 5 and 6a, the central row shows the resulting simulated tomogram reconstruction obtained through the filtered backprojection algorithm, and the bottom row shows the reconstructions obtained through the filtered backpropagation with Rytov approximation. Each step of the pipeline is important and can be flexibly used to tailor suitable algorithms of interest according to the specific needs of the user and the features of the employed optical system. For example, we can investigate the effects of limited projection angles and the impact of errors in projection angle recovery, particularly in in-flow setups. In addition, denoising procedures can be experimented to study their effects on QPMs or even on the tomogram.

Fig. 9 application of the introduced pipeline to the simulated cells displayed in this work. Top row, slices from the simulated cells; Middle row, corresponding slices of the reconstructed simulated tomogram through the filtered backprojection routine. Bottom row, corresponding slices of the reconstructed simulated tomogram through the filtered backpropagation routine. Columns 1a and 1b refer to Simulation 1 with 1 and 7 cytoplasmic vacuoles respectively, columns 2, 3 and 4 refer to Simulation 2, 3 and 4. An equally spaced angular sequence was used, with the angles ranging from 0 to 360° with a step of 5°.

Finally, tomographic solvers can be studied and fine-tuned, to compare the ideal tomogram to the reconstructed one. In the specific examples shown in Figures 8 and 9, differently from the cell, the tomogram shows artifacts in the external region and residuals due to noise, which can be detected in comparison with the proposed ground truth model and could guide the refinement of the denoising and tomographic algorithms.

3.3 Experimental use of the simulated models

Here we present an example of how our simulated model can be embedded into an experimental context. Our simulation was used to tune a Convolutional Neural Network (CNN) to replenish the reconstruction obtained via sparse view filtered backprojection [16, 51]. The experimental data of interest were THP1 monocytes, acquired using the setup described in [16] and processed accordingly. The simulated model was used to create a population of simulated monocytes, which were then convoluted with a PSF modelling our optical system to obtain the corresponding simulated ground-truth tomograms, as depicted in Figure 10. The SSIM between the experimental tomograms and the simulated ones is 0.9977, suggesting a high-fidelity representation by the model. The simulated QPMs were then computed by projecting the simulated tomograms along chosen angular directions and used to compute the sparse-view tomograms, which contain artifacts. The sparse view tomogram, along with the ground-truth data, was used to train the CNN to remove existing artifacts. Again, the elevate SSIM value between the real QPMs and the real ones (0.9951) suggests the good quality of the model. Training the CNN over a simulated dataset allowed us to extend its cardinality without limitations, to overcome the risk of overfitting. In addition, the simulated dataset allowed us to monitor the training parameters, evaluate our models and assess the performance of the reconstruction technique. Finally, we were able to inference the CNN trained using simulated-only data over experimental ones, as shown in Figure 8c. The CNN’s ability to generalize to unseen experimental data confirmed that our simulated model was sufficiently representative of the real-world data.

Fig. 10 application of the proposed simulated model to an experimental context. a), comparison between the simulated QPM and tomogram and the experimental ones. b), training of a CNN using a simulated dataset. c), reconstruction of an experimental tomogram via filtered backprojection (FBP) and via CNN.

4 Conclusions

In this paper, a simulation method has been proposed for the high-fidelity modelling of the 3D RI distribution of suspended cells. As a benchmark, we simulated yeast cells as an example model for eukaryotic cell lines and human monocytes as an example model for mammalian cell lines. We simulated the 3D RI distribution accounting for the presence of internal organelles/sub-compartments. Thus, we created a dataset of cells with various levels of complexity. As a proof of concept, we simulated the presence of various phenotypes of cytoplasmic vacuoles inside yeast cells and monocytes, and the presence of LDs and lysosomes inside monocytes. All these compartments are of great interest since they are biomarkers associated with various dysfunctions, inflammations, stress conditions and diseases in cells of the human bloodstream. Besides the reported examples, the proposed approach is general enough to simulate realistic RI distributions of various organelle complexes provided that prior statistical information is available from the literature and/or measured tomographic data. Here, the simulations are based on experimental statistics present in the literature and in some cases on data gathered by using a holo-tomographic flow cytometry apparatus. Once the RI distributions of the cells are simulated, a simplified example is shown of how our proposed cell model can be applied to tune the steps of a general holo-tomographic reconstruction pipeline. In this way, the loss of resolution due to the passage through the system PSF, the tomogram quality worsening in terms of SNR due to the presence of speckle noise, and the numerical artefacts due to the non-ideal inversion algorithm (e.g., the Inverse Radon Transform or the filtered backpropagation with Rytov approximation) can be considered. Results show the versatility of the simulation method for the realization of numerical phantoms which are of primary importance for the development of holo-tomographic flow cytometry reconstruction, denoising, 3D segmentation and analysis algorithms in the absence of a ground-truth. A consideration is due on the status of the literature and research efforts in the field of additive manufacturing [52, 53]. Indeed, extremely realistic phantoms have been developed that simulate the 2D and 3D RI distribution of cells [54, 55]. So far, these have been developed for simulating cells and more complex objects in adhesion, but a promising perspective is envisaged towards the realization of 3D phantoms of suspended cells. In the future, the presented method can be easily extended to represent different classes of cells or organelles whose statistics are available in the literature and could complement the abovementioned research efforts. For example, the recently introduced methods to add specificity about cell sub-compartments of suspended cells [17, 20] can be applied to study the RI distribution of different organelles such as the nucleus, and then embed them into our cellular model, or modify the statistics representing the general shape of the cell to describe pathological conditions. Results show the versatility of the simulation method for the realization of numerical phantoms which are of primary importance for the development of holo-tomographic flow cytometry reconstruction, denoising, 3D segmentation and analysis algorithms in the absence of a ground-truth.

Funding

This work was supported by a project PRIN 2022 – Label-free cytoplasmic vacUoles pheNotyping plAykit (LUNA) Prot. 960, 30th June 2023 – funded by the Italian Ministry of University & Research in the framework of the European Union program NextGenerationEU (2022BKYM22 – Project CUP: B53D23002490006).

Conflicts of interest

Authors declare no conflict of interest associated to this work.

Data availability statement

The whole simulated dataset is shared publicly and can be downloaded at the following link: https://doi.org/10.6084/m9.figshare.27226464.

Author contribution statement

Investigation: P.F., L.M.; Supervision: V.B.; Methodology: F.B., D. P.; Software: F.B., D.P., P.M.; Validation: L.M.; Visualization: F.B., D.P.; Formal analysis: F.B., P.M., V.B.; Funding acquisition: V.B.; Writing – original draft: F.B.; Writing – review & editing: V.B., L.M., P.M., P.F.

References

All Tables

All Figures

Fig. 1 Conceptual steps common to all the proposed simulations. a), Schematization of the modelling of the external shape. V = Volume, SI = Sphericity Index, rx, ry and rz: semiaxes of the simulated ellipsoid. b), creation of the irregular external shape. c) Simulation of the cytoplasm. rx,cyto, ry,cyto and rz,cyto: semiaxes of the ellipsoid simulating the cytoplasm. d), Simulation of the nucleus. rx,nuc, ry,nuc and rz,nuc: semiaxes of the ellipsoid simulating the nucelus, NCR = Nucleus to Cell Ratio.

In the text

	Fig. 2 schematical representation of the geometric structure used for modelling the cytoplasmic vacuole. N shells are modeled, each of them characterized by x and y dimensions rx and ry and a Refractive Index value RIn.
In the text

Fig. 3 simulation results for (a–d) a yeast cell with a large vacuole inside (Visualization 1) and (e–h) a yeast cell with a large vacuole and a multitude of smaller vacuoles (Visualization 2). (a, g) Visualization of a slice in the yz plane from the 3D RI profiles; (b, f) Cut along the red dotted line in subfigures (a, e), highlighting the RI profile into the vacuole; (c, e) Visualization of a slice in the xy plane from the 3D RI profiles; (d, h) Isolevel visualizations. EM = External membrane, C = Cytoplasm, N = Nucleus, V = Vacuole. For each cell, the isolevels’ thresholds for the visualization are set with respect to the maximum value of the simulated RI profile.

In the text

	Fig. 4 simulation results for late-G1 yeast cells. (a, b) isolevels of two different simulated cs; (c) Histograms of some relevant statistics evaluated over a population of 50 cells. For each cell, the isolevels’ thresholds for the visualization are set with respect to the maximum value of the simulated RI profile. The statistics are in good accordance with the ones in the literature.
In the text

Fig. 5 results for the simulations of a monocyte with LDs (Visualization 3). (a) Isolevel. From this visualization, the size and the positioning of the LDs can be seen, compared to the overall cell; (b) yz cut of the simulated 3D RI distribution. The non uniform distribution of RI into the LD can be noticed; (c) Histograms of the sphericity index, average RI and standard deviation values are evaluated for the LDs in a population of 50 monocytes. μ = mean value. σ = standard deviation. EM = External membrane, C = Cytoplasm, N = Nucleus, LD = Lipid Droplets. For each cell, the isolevels’ thresholds for the visualization are set with respect to the maximum value of the simulated RI profile.

In the text

Fig. 6 Results for the simulations of monocytes with cytoplasmic vacuoles and monocytes with cytoplasmic vacuoles and lysosomes. Top panel, (a), an isolevel visualization and a yz slice from a simulated cell, highlighting the positioning of the vacuoles and their spatially-varying RI distribution (Visualization 4). Bottom panel, (b), an isolevel visualization and an xz slice from a simulated cell, where the high RI and numerosity of lysosomes can be appreciated (Visualization 5). EM = External membrane, C = Cytoplasm, N = Nucleus, V = Vacuoles, L = Lysosomes. For each cell, the isolevels’ thresholds for the visualization are set with respect to the maximum value of the simulated RI.

In the text

	Fig. 7 Left: Numerical digital hologram obtained starting from a simulated yeast cell with a vacuole (Simulation 1). The interference fringes are highlighted in the yellow box. Right: corresponding simulation of the amplitude Fourier spectrum, in which the holographic diffraction orders are clearly visible.
In the text

Fig. 8 General sketch of a holo-tomographic acquisition-reconstruction pipeline. (a) Most important steps are schematically depicted; (b) Images corresponding to each step particularized to our optical setup and reconstruction pipeline are discussed (Visualization 6). For the sake of visualization clarity, a 13 × 13 × 13 Gaussian kernel is shown. In this case an equally spaced angular sequence was used, with the angles ranging from 0 to 360° with a step of 5°, recalling an in-flow tomographic reconstruction pipeline.

In the text

Fig. 9 application of the introduced pipeline to the simulated cells displayed in this work. Top row, slices from the simulated cells; Middle row, corresponding slices of the reconstructed simulated tomogram through the filtered backprojection routine. Bottom row, corresponding slices of the reconstructed simulated tomogram through the filtered backpropagation routine. Columns 1a and 1b refer to Simulation 1 with 1 and 7 cytoplasmic vacuoles respectively, columns 2, 3 and 4 refer to Simulation 2, 3 and 4. An equally spaced angular sequence was used, with the angles ranging from 0 to 360° with a step of 5°.

In the text

	Fig. 10 application of the proposed simulated model to an experimental context. a), comparison between the simulated QPM and tomogram and the experimental ones. b), training of a CNN using a simulated dataset. c), reconstruction of an experimental tomogram via filtered backprojection (FBP) and via CNN.
In the text