Multi-scale fractal Fourier Ptychographic microscopy to assess the dose-dependent impact of copper pollution on living diatoms

Accumulation of bioavailable heavy metals in aquatic environment poses a serious threat to marine communities and human health due to possible trophic transfers through the food chain of toxic, non-degradable, exogenous pollutants. Copper (Cu) is one of the most spread heavy metals in water, and can severely affect primary producers at high doses. Here we show a novel imaging test to assay the dose-dependent effects of Cu on live microalgae identifying stress conditions when they are still capable of sustaining a positive growth. The method relies on Fourier Ptychographic Microscopy (FPM), capable to image large field of view in label-free phase-contrast mode attaining submicron lateral resolution. We uniquely combine FPM with a new multi-scale analysis method based on fractal geometry. The system is able to provide ensemble measurements of thousands of diatoms in the liquid sample simultaneously, while ensuring at same time single-cell imaging and analysis for each diatom. Through new image descriptors, we demonstrate that fractal analysis is suitable for handling the complexity and informative power of such multiscale FPM modality. We successfully tested this new approach by measuring how different concentrations of Cu impact on Skeletonema pseudocostatum diatom populations isolated from the Sarno River mouth.


Materials and methods
Sampling, preparation and exposure of diatom probes to Cu Skeletonema pseudocostatum was isolated on June 2020 from water samples collected from the Sarno River mouth (40.7291N, 14.4698E) by hand pipetting, correctly identified by PCR amplification and sequencing of 18S and 28S genes (Fig. 1A), and treated as described in Supporting Information.Skeletonema pseudocostatum was then exposed to scalar Cu concentrations to obtain the following final concentrations: 0 (e.g.control conditions), 5 µM, 10 µM, 15 µM, 25 µM, 35 µM and 50 µM (see Supporting Information).Concentration was daily assessed at different time intervals (0, 24, 48 and 72 h), hereafter referred to as T 0 , T 1 , T 2 , T 3 , respectively.Experimental samples were all acquired as described in Bianco et al. 30 at the same cell density (around 190,000 cells/ml), obtained by appropriate dilution in f/2 or centrifugation (4 °C, 2700 g, 10 min, Thermo Scientific SL16R) and resuspension of the cells in the desired volume.For each acquisition, 900 µl of lugol-fixed sample were placed in a Petri dish (glass bottom WillCo-dish® GWST-5040, 40 mm aperture diameter, 7 mm height) to form a very thin film of liquid to obtain images of in-focus diatoms.3 FoVs were acquired for each sample.In particular, we focused on the exposure times T 2 and T 3 , since the exposure times T 0 and T 1 are expected to be too short to observe significant differences among diatoms subjected to different Cu doses.As regards the exposure time T 2 , we recorded PCMs for each Cu concentration under test.The same procedure was repeated in the T 3 case, with the sole exception of the highest concentration tested (50 μM), that caused a consistent reduction of cell density and was considerably higher than half maximal effective concentration.

FPM principle
All the FPM acquisitions consist of illuminating the object with plane waves distinct in angular inclinations while simultaneously capturing the corresponding intensity images.This procedure can be clearly understood considering that fine structures of the sample transmit through the low-pass filtering imaging system when the sample is probed with oblique illumination, Fig. S1.Sequential-illumination FPM measurements of a thin object o(r) with a transversal vector r = (x,y) are often accomplished with an LED array placed far enough from the sample; thus, the probing can be considered a locally-coherent plane-wave illumination of a central wavelength λ.Consequently, each of the N used LEDs contains the spatial carrier frequency ν n = (ν xn , ν yn ) , n = 1,…,N.The light behind the sample is straightforwardly considered a product of the sample and the oblique illumination optical fields as o(r)exp(i2πν n • r) , which generates a shift of the object's Fourier spectrum O(v − v n ).Here O(ν) = FT{o(r)} represents the original-sample Fourier spectrum with FT{} being the 2D Fourier transform operator.The optical field subsequently passes the imaging system serving as a low-pass filter of a pupil function P(ν) .The truncating function P(ν) is commonly considered a circle of radius NA MO /λ, where NA MO is numeri- cal aperture of the used microscope objective (MO); thus, the field arising at the pupil of the MO is expressed as the product of O(ν − v n ) and P(ν) .The low-resolution intensity image created in the detector plane is then calculated as I n (r) = |IFT{O(ν − v n )P(ν)}| 2 , where IFT{} represents the 2D inverse Fourier transform operator.The primary motivation of the complex-amplitude retrieval process is to jointly estimate the object function O(ν) and the pupil function P(ν) , the latter containing the optical system aberrations that may distort the image quality 42 .The whole sequence of N low-resolution intensity snapshots I n (r) , n = 1,…,N, serves as the input of the retrieval algorithm, where the up-sampled version of the central low-resolution bright-field image acts as an initial guess of the object, and the pupil function is initialized as a circle of radius NA MO /λ and the constant phase.All the captured images are addressed sequentially from n = 1 to n = N and considered in this order, with both the pupil function and sample spectrum being updated in each of the j-th loop (j = 1,…,J) 42 .
In the j-th iteration of the retrieval algorithm, functions P (j) (ν) and O (j) (ν) are created using their estimates obtained in the previous (j−1)-th loop.The output field of the object illuminated by the n-th LED in the pupil plane of the used microscope objective is calculated as ϕ (j−1) (ν) = O (j−1) (ν − v n )P (j−1) (ν) , and the correspond- ing low-resolution complex-amplitude image in the detector plane is simulated by its inverse Fourier transform as � (j−1) (r) = IFT(ϕ (j−1) (ν)) .Subsequently, based on the Gerchberg-Saxton-Fienup retrieval approach, the amplitude of such complex distribution is replaced by the square root of the measured low-resolution intensity image; thus, �′ (j−1) (r) = (I n (r)) 1/2 � (j−1) (r)/|� (j−1) (r)| and the updated pupil plane distribution are calculated according to ϕ′ (j−1) (ν) = FT(�′ (j−1) (r)) .Following the gradient-descent principle, both estimates of O (j−1) (ν) and P (j−1) (ν) are further updated based on previous calculations as and respectively.The functions G (j−1) and H (j−1) vary during each iteration due to their dependency on the estimated pupil function, object's Fourier spectrum, and additional parameters as where * represents a complex conjugation, max(�) is maximum of , and δ 1 , δ 2 are regularization constants ensuring numerical stability 43 .Simultaneous retrieval of the sample spectrum and the pupil function enhance the accuracy of the complexamplitude recovering process.These updates are repeated until all n = 1,…,N captured low-resolution intensity images are used, which can be considered a single complete iteration that recovers the high-resolution spectrum.This single complete iteration is repeated J times to improve convergence, and in turn to get the steady-state pupil function and the sample spectrum estimates P (J) (ν) and O (J) (ν) .Finally, the obtained high-resolution spectrum O (J) (ν) is inverse Fourier transformed to obtain the high-resolution complex-amplitude image of the investigated specimen o (J) (r) .More details about the used algorithm are included in the reference 43 .

Results
Samples were analysed using a custom FPM inverted microscope, sketched in Fig. 1B and described in the Supporting Information (Fig. S1).Here we used FPM to image in the same large FoV hundreds of diatoms, used as bio-probes after being exposed to controlled Cu doses.Indeed, the low 4×g permits accessing a wide 3.3 mm 2 FoV, while the FPM synthetic aperture principle allows super-resolved quantitative phase-contrast imaging, with 0.5 μm lateral resolution.This unique feature of FPM allows imaging in the same FoV large numbers of diatoms on a millimetric scale while preserving the micron-scale characteristics of single diatom elements (e.g., see the insets in Fig. 1B).An example of FPM processing is reported in Fig. 2 for one of the control acquisitions.We report (a) the full FoV of one of the 177 captured images that corresponds to the Low Resolution (LR) brightfield intensity obtained by switching on the central LED, (b) the full FoV reconstructed High Resolution (HR) amplitude and (c) the full FoV reconstructed HR phase-contrast map.The HR phase-contrast map is the one we use for the further fractal analysis.In particular, in Fig. 2d we show some zoom-in details corresponding to the areas marked by yellow boxes in Fig. 2c, where the effect of super-resolution processing and phase retrieval is apparent from the comparison between the LR bright-field and the HR phase-contrast.Figure 3 shows one characteristic example of diatoms imaged by FPM for each exposure time and Cu dose.The phase-contrast maps at the single diatom level are shown.In some cases, especially for the highest Cu doses, diatoms experience morphological variations and/or cytoplasm extrusion after membrane lysis.In these cases, the variation of the phase contrast is made noticeable by observing the phase-contrast map.However, there are cases (especially for the sub-lethal doses) where the morphological variations are not so evident for the single diatom and membrane lysis did not occur.
In order to fully exploit the multi-scale property of FPM, it is convenient to look for the analysis methodology in the framework of multi-scale math categories.So far, fractal geometry has been demonstrated the most powerful instrument able to provide an insight of nature deeper than its Euclidean counterpart.Indeed, its main credit lies in the ability of breaking the mold of the classical geometry that frames the topological dimension of an object in the realm of integer numbers, thus pushing towards a fascinating non-integer vision of reality 44 .To accomplish this goal, we inspected the samples by a multi-scale fractal analysis that sinks its roots on the box-counting concept.Among the possible fractal parameters that could be considered, here we focused on lacunarity [44][45][46] .Lacunarity can be thought as the distribution of hole sizes inside a single object 41 or within a group of objects imaged in the same FoV, e.g. a tissue or confluent cell layers 46 .Lacunarity is expected to provide a reliable solution to the problem we are tackling, since it is a quantitative descriptor very sensitive to the tiny dissimilarities between the PCMs of the bio-sentinels exposed to different Cu doses.We introduced two lacunarity-based parameters (Fig. 1C), here defined as the global lacunarity (GL) and the local lacunarity (LL), described in detail in the Supporting Information section.The global lacunarity GL is an ensemble descriptor of the hole sizes distribution computed at different scales and therefore it considers the global context of the PCM, i.e. both the diatom probes and the background medium.Instead, the local lacunarity LL measures the hole sizes distribution at a fixed scale, properly selected here in order to consider only the local interplay between the bioprobes and the medium at the single-diatom level.The following calibration experiments and the proposed test are based on the Multi-Scale Lacunarity (MSL) parameter, defined here as the ratio between the global and the local lacunarities and used to measure the Cu concentration the imaged diatoms were exposed to.

Calibration
We employed the described pipeline of FPM and multi-scale fractal analysis to record and characterize liquid samples made of diatoms exposed to different Cu doses for both the exposure times T 2 and T 3 .The two exposure times have been treated separately, as shown in Fig. 4A,B.For each Cu dose, we measured the MSL from the PCMs and exploited it to calibrate the system.As desired, and also expected on the basis of the previous considerations, the corresponding dots reported in Fig. 4A,B follow monotonic trends with the Cu dose.Therefore, we fitted them by exponential curves (solid lines), used here as reference calibration curves.The comparison between the T 2 and T 3 calibration curves in Fig. 4C is remarkable.Indeed, the two curves start from the same MSL value at 0 μM (control sample), despite they were obtained independently of each other from different experimental sets and acquired in different days in order to avoid any experimental bias.
Moreover, they have comparable MSL values at low doses, and strongly diverge at high doses.Measuring the same starting point underlines the robustness and repeatability of this strategy and serves as a negative control for the sensor, since the diatoms in the control images were not exposed to any Cu dose, therefore no difference (4)  between the T 2 and T 3 cases should have emerged.Instead, the divergence at high doses is a first validation of this system on the basis of biological considerations.In fact, because of the longer exposure time, the T 3 curve grows faster with the Cu dose than the T 2 curve, since diatoms are more impacted by the toxic effects due to the longer Cu internalization, and therefore they start exhibiting discrepancies with respect to the control at lower Cu doses.Furthermore, the steepest trend of the T 3 curve leads to a greater sensitivity of the proposed sensor in detecting the correct Cu dose intervals with respect to the T 2 curve, as larger MSL differences can be measured in correspondence of the same Cu dose interval.Hence the curves shown in Fig. 4 put in evidence the capability of the system and the overall measurement and analysis protocol of sensing this difference of impact on the bio-sentinels.

Testing and sensitivity enhancement
In order to test the ability of the calibrated system in measuring a certain Cu dose on the basis of the MSL, in   μM, and [42.5,∞) μM.Instead, the greater slope of the T 3 case displayed in Fig. 5B allows revealing four Cu dose intervals, i.e. [0,17.5]μM, [17.5,22.5]μM, [22.5,30] μM, and [30,∞) μM.On the basis of the reported analysis, we conclude that the exposure time T 3 is the optimal solution to determine the unknown Cu doses due to the greater sensitivity of the system.This is particularly important for discriminating the Cu doses in the lowest concentrations range, where the current biological procedures are not able to catch very small changes among the different Cu effects.Moreover, the exposure time T 3 also allows an enhancement of the system sensitivity in finding the lowest Cu doses within the [0,17.5]μM interval.Toward this scope, we characterized each PCM through a more conventional entropy-based texture parameter, i.e. the range GLCM entropy ∆s described in the Supporting Information and shown in Fig. S4O.Entropy is well-known as a metric of the disorder of a system.Hence, entropy is expected to grow with the Cu dose because, as discussed in the previous section, the sample goes from a clean medium with embedded healthy elements to a dirtier medium with elements undergoing membrane lysis that provokes cytoplasm leaks, i.e. a situation with a higher disorder level.This hypothesis is verified in Fig. 5C, where the entropy-based texture parameter corresponding to the calibration doses follow again a monotonic trend, but with more marked differences among them in the [0,15] μM interval.Therefore, we used a cubic fitting to create the low-dose calibration curve (solid line in Fig. 5C).Finally, as shown by the test PCMs (black dots in Fig. 5D) overlapped to the entropy-based calibration curve, three Cu sub-intervals can be found in correspondence to the lowest doses range, i.e. [0,2.5]μM, [2.5,7.5]μM, [7.5,12.5]S4O, even if the T 3 curve has a monotonic trend, the sole entropy-based texture parameter cannot replace the MSL in the first step of the detection system since the basic congruence between the T 2 and T 3 data discussed in Fig. 4C is missing when the Cu dose is higher than 15 μM.In the Supporting Information and Fig. S4, other conventional features have been tested, but none of them shows the congruence property exhibited by the fractal analysis (see Fig. S4F,H,K,L).Some features do not show a monotonic trend even in the T 3 case (see Fig. S4A-E,G,I,J,M,N).This congruence property is instead pivotal as it gives robustness to the proposed Cu detection system, since the output of the fractal analysis is able to replicate what is expected from the biological point of view, i.e. it is able to provide a parameter that follows a monotonic trend in both the T 2 and T 3 cases but with different growth rates due to the different exposure times, while starting from the same value corresponding to the negative control.

Fractal-FPM timing analysis
The time needed for detecting the Cu dose can be considered as the sum of three terms.
The first term is related to the FPM recording time.In fact, 177 low-resolution images are acquired in about 4 min.However, it has been demonstrated that the number of employed LEDs can be reduced without losing resolution by employing coded illumination schemes.The second term is the computational time requested by the phase retrieval algorithm for obtaining a highresolution FPM map from the recorded 177 low-resolution intensities.We have demonstrated in a previous work that, using deep learning, the phase retrieval can be greatly speeded up, thus taking about 1 min for reconstructing the entire 7000 × 9500 (i.e., 3.3 mm 2 ) FoV 33 .The third term is the computational time requested by the fractal analysis for computing the global lacunarity, which is about 5-6 min for each 7000 × 9500 (i.e., 3.3 mm 2 ) FoV.However, due to the structure of the fractal analysis algorithm, parallel computing along with the use of GPUs could significantly speed up this process.

Conclusions and discussion
High doses of the Cu bioavailable fraction can play a negative role in the morphology, growth, survival rate, and reproduction of primary producers.Sensors have been proposed to detect and quantify the concentration of HMs in water.SEM has been already employed in experiments performed on Skeletonema spp.exposed to metals, but it was mainly employed to study and assess the distribution of metal nanoparticles in the culture media rather than to observed eventual impairments in the algal morphology 47 .Moreover, metal-treated cultures can be characterised by cell aggregations and metal sorption on the cell wall which make difficult a cell-by cell observation at SEM or TEM microscopes 48,49 .This problem could be partially overcome by diluting samples before the observation, but this strategy further limits the number of cells which can be observed, rendering their number non-statistically significant.TEM micrographs were acquired in a previous work on the model species P. tricornutum exposed to copper, and images include only 2-4 diatoms 22 .Another bottleneck related to the use of Electron Microscopy is the time needed for sample preparation, that requires sample fixation with toxic agents and several dehydration steps with alcohols.In light of this, we have opted for an experimental design that includes the estimation of cell concentration as broad indicator of a dose-dependent effect of toxicity of copper, and the FPM as image descriptor system.Our work followed a functional approach, i.e. we proposed to use an all-optical microscopy paradigm to characterize the effects HMs may produce on aquatic communities, which intrinsically depend on the concentration of their bioavailable fraction.In particular, we introduced a test that relies on three main concepts: i) The use of microalgae sensitive to the presence of HMs.ii) FPM as a multi-scale QPI method.iii) Multi-scale fractal analysis of FPM-QPI maps at the aim of assessing stress condition in the algal population when still capable of sustaining a positive growth.The use of FPM provides the unique opportunity to image the population and the liquid under test as a whole over a wide FoV, while returning super-resolved phase-contrast information at the single-cell level and in label-free mode.Potential alternatives to the employed FPM apparatus within the QPI framework are in-line digital holography systems or different FPM configurations employing efficient source coding schemes.In-line digital holography by lensless configurations can ensure a quite wide FoV, i.e. as wide as the active area of the employed sensor 50 .However, resolution is inherently limited by the pixel pitch and it is strictly dependent on the adopted sensor.To keep the required resolution across the multiple scales and thus to permit the fractal analysis, we preferred to use FPM schemes.In our case, samples are fixed in lugol, so that time resolution was not a requirement.Whenever time resolution is important, or simply one wants to reduce the acquisition times for wide data collection campaigns, efficient source coding schemes implementing illumination multiplexing should be adopted 42,[51][52][53] .Within the FPM framework, it has been demonstrated that the number of images can be reduced to 72 using annular bright-field illumination 52 and up to 42 if one considers symmetrical dark-field illumination schemes 53 .This system improvement will be object of future investigation from our group.
We found that fractal analysis well manages the complexity of multi-scale FPM maps and can determine dose-dependent stress conditions in the microalgae populating the liquid under test.As a testbed for the proposed method, we use Skeletonema Pseudocostatum isolated from the Italian Sarno River mouth to probe the dose-dependent effects of the Cu bioavailable fraction.We conducted a study to investigate the optimal sample treatment protocol, (see Fig. 1A) and test modalities to optimize the sensitivity of the system to tiny dose variations that can provoke a non-negligible impact on diatoms.In these experiments, the initial cell density was relatively high, and algal cultures reached the late stationary phase within 72 h.This time frame is usually employed in several ecotoxicological assays aimed at assessing the effect of copper on diatoms and, in general, on model microalgal species [54][55][56] , allowing an immediate comparison between standardized and validated methods and the new protocol described in the present paper.We compared, in terms of system sensitivity, experiments conducted after exposing diatoms for 48 h and 72 h to different Cu doses, (i.e.T 2 and T 3 ).Although both conditions allowed us detecting the presence of Cu, we found that T 3 is an optimal exposure to maximize the sensitivity.Instead, preliminary experiments performed at T 0 and T 1 showed poor system sensitivity, especially in the low and medium doses range, which means those exposure times are too short to obtain a non-negligible effect on this diatom species.
The proposed method is applicable in practice by sampling natural water from the location point under test, applying water filtering and/or digestion to remove the largest objects/particles, then to put S. pseudocostatum diatoms (from lab cultivation) and to let them in the liquid under test for the desired amount of time corresponding to T 3 for maximum sensitivity.In this way one could be certain about the time diatoms are exposed to Cu and can use them as bioprobes.In the conditions of our experiments, 72 h were sufficient to reach the late exponential phase in algal cultures.Additional experiments with less concentrated inocula of S. pseudocostatum could allow to evaluate the efficiency of this method to assess the effect of a longer-term exposure (7-10 days) of this species to Cu.Moreover, the efficiency on long timescales (e.g. 1 month) could be also assessed by employing other microalgal species -such as Chlorophyceae and Eustigamtophyceae-exhibiting slower growth rates (and thus longer growth curves) as test-organisms.
From the mere image analysis standpoint, results presented here are a first proof that multi-scale microscopes can reach their highest informative potential when coupled to analysis methods that consider the image features on multiple scales.Although this is a concept generalizable to a wider set of microscopy modalities and mathematical descriptors, we believe that fractal geometry is the most suitable candidate and can provide the optimal framework within this scope.It is worth to mention that FPM microscopes can be realized with non-demanding hardware components 31 and in principle any benchtop optical microscope could be turned into a FPM setup by using a commercial LED matrix 28,30,31,42 .We believe this is an advantageous feature to broadly promote a new type of dose-dependent bioimpact-based ecotoxicological assays in biology labs.

Figure 1 .
Figure 1.Analysis protocol to use diatoms as bio-probes of Cu doses in water.(A).Sampling and preparation pipeline to create the sensor calibration curves.(B).Sketch of the multi-scale FPM imaging device.The output is the high-resolution phase-contrast map that we use for downstream fractal analysis.(C).Multi-scale analysis of FPM maps of marine diatoms based on fractal geometry.In the inset we report example plots of global and local lacunarity parameters.As a result of the analysis, each water sample is associable to one interval of Cu doses out of seven.Figure 1A: Created with BioRender.com.
Figure 1.Analysis protocol to use diatoms as bio-probes of Cu doses in water.(A).Sampling and preparation pipeline to create the sensor calibration curves.(B).Sketch of the multi-scale FPM imaging device.The output is the high-resolution phase-contrast map that we use for downstream fractal analysis.(C).Multi-scale analysis of FPM maps of marine diatoms based on fractal geometry.In the inset we report example plots of global and local lacunarity parameters.As a result of the analysis, each water sample is associable to one interval of Cu doses out of seven.Figure 1A: Created with BioRender.com.

Figure 2 .
Figure 2. FPM processing results.(a) Full FoV LR bright-field intensity, central LED.(b) Full FoV reconstructed HR amplitude.(c) Full FoV reconstructed HR phase-contrast map.(d) Zoom-in details corresponding to the areas marked by yellow boxes in (c).

Figure 3 .
Figure 3. details of FPM phase-contrast maps obtained from experiments at different Cu doses and exposure times.
Fig.5A,B we report as black dots the MSL values obtained from PCMs not considered during the calibration step, overlapped to the T 2 and T 3 calibration curves, respectively.As for the exposure time T 2 (Fig.5A), due to the lower slope of the calibration curve, only three Cu dose intervals can be detected, i.e. [0,30] μM,[30,42.5]

Figure 4 .
Figure 4. Calibration of the system for detecting the Cu dose.A,B Multi-scale lacunarity of the calibration data (dots) about the FPM-PCMs with diatoms recorded respectively after T 2 and T 3 exposure times at different Cu doses.The equations of the best exponential fittings (solid lines) are reported.C Comparison between the T 2 and T 3 calibration curves.

Figure 5 .
Figure 5. Testing the system for detecting the Cu dose.A,B MSL calibration curves (solid lines) with overlapped the test data (dots) about the FPM-PCMs with diatoms recorded respectively after T 2 and T 3 exposure times at different Cu doses.The detectable Cu ranges are highlighted.C Entropy-based texture parameter of the calibration data (dots) about the FPM-PCMs with diatoms recorded after T 3 exposure time at different Cu low-doses.The equation of the best cubic fitting (solid lines) is reported.D Entropy-based calibration curve (solid line) with overlapped the test data (dots) about the FPM-PCMs with diatoms recorded after T 3 exposure time at different Cu low-doses.The detectable Cu ranges with higher sensitivity at the low doses are highlighted. https://doi.org/10.1038/s41598-024-52184-3 μM, and[12.5,17.5]μM.In summary, as reported in Table1, the proposed test for detecting the impact of different Cu concentrations is made of two steps.The former exploits the fractal geometry to compute the MSL, which value yields four possible Cu dose intervals, i.e. [0,17.5]μM,[17.5,22.5]μM,[22.5,30]μM,and[30,∞)μM.The latter is performed only in the case the MSL measured at the end of the first step is below the 17.5 μM value.In such a case, the system sensitivity is enhanced by computing a PCM entropy-based texture parameter, thus providing four further Cu low-dose sub-intervals, i.e. [0,2.5]μM,[2.5,7.5]μM,[7.5,12.5]μM,and[12.5,17.5]μM.It is worth noting that, as reported in Fig.

Table 1 .
Thresholds of the MLS Fractal Parameter and Entropy-based Texture Parameter for detecting different Cu Ranges inside a water sample.