Local aberration control to improve efficiency in multiphoton holographic projections

: Optical aberrations affect the quality of light propagating through a turbid medium, where refractive index is spatially inhomogeneous. In multiphoton optical applications, such as two-photon excitation fluorescence imaging and optogenetics, aberrations non-linearly impair the efficiency of excitation. We demonstrate a sensorless adaptive optics technique to compensate aberrations in holograms projected into turbid media. We use a spatial light modulator to project custom three dimensional holographic patterns and to correct for local (anisoplanatic) distortions. The method is tested on both synthetic and biological samples to counteract aberrations arising respectively from misalignment of the optical system and from samples inhomogeneities. In both cases the anisoplanatic correction improves the intensity of the stimulation pattern at least two-fold.

such distortions. An AO system may consist of three main components: a sensor to measure the aberrations, an adaptive optical element to correct the wavefront distortions, and a controller to feed the information obtained by the sensor to the adaptive element. An alternative to using a wavefront sensor in closed loop is a sensorless approach in which information about the wavefront can be determined indirectly (e.g., from the properties of two-photon excited fluorescence). Both sensor-based and sensor-less AO approaches have been successfully integrated into different microscope designs [8]. In widefield microscopy, a Shack-Hartmann wavefront sensor may be used to measure the distortion on the fluorescent signal using guide stars, or fluorescent landmarks with known shape [9]. Differently, wavefront sensors have been employed in light sheet microscopy to correct the aberrations on the excitation light [10]. The approach most commonly used in laser scanning microscopies, like laser scanning confocal microscopy, is to obtain a single isoplanatic aberration correction for both the excitation light and fluorescence signal over a relatively small FOV [11]. At the same time, sensorless AO has been widely applied to correct excitation and or emission light in different microscopy techniques [12][13][14]. Recently, various research groups proposed anisoplanatic corrections to improve imaging in more heterogeneous samples, such as intact organs [15][16][17][18].
In non-imaging applications, holographic projections are likewise subject to aberrations. Stray photons may have unintended consequences in the case of laser writing or optogenetics. In particular, in the case of optogenetics, chromophores are intentionally designed to perturb the biology of the system rather than simply fluoresce. In fact, photostimulation conditions in optogenetics are often more intense than conditions used in imaging, using longer dwell times (millisecond vs. microsecond), so the effect of off-target photostimulation can be more problematic. Inefficient excitation is also detrimental for photostimulation in living tissues where one wants to use the minimum light needed to be effective.
Here we apply a sensorless AO approach to correct aberrations on holographic projections using two-photon excited fluorescence (TPEF) as a readout. An optical module for photostimulation with CGH is implemented, which we will call 2P-CGH in this paper. A computer generated hologram (CGH) of the desired 3D excitation pattern is calculated and its phase is imposed on a coherent wavefront at the rear pupil plane of an objective through a liquid crystal on silicon spatial light modulator (LCoS-SLM). The SLM serves also as AO element to introduce bias aberration into the stimulation light. The optimal AO corrections are determined through a simple hill-climbing approach where a series of images is acquired under different bias aberrations. To validate the method we demonstrate correction of anisoplanatic aberrations in 2D and 3D using a series of samples with increasing complexity, including tissue. We show that our method can improve the intensity of the stimulation pattern in both synthetic samples and fixed tissue by at least a factor of two.

Computer generated hologram calculation
The optical layout to project computer generated holograms is based on Fourier holography, as shown in Fig. 1(a). The target field U o (z i ) is defined by user-selected points in a virtual representation of the real space FOV, describing the 3D pattern of points at depths z i . The hologram field U h is related to the field at the focal planes through a discrete Fourier transform. The phase of the hologram field is then addressed to the SLM and the hologram is reconstructed optically by a spatial Fourier transformation by the objective lens. Holograms are calculated with the Compressive Sensing Weighted Gerchberg-Saxton (CS-WGS) algorithm [19]. This algorithm can generate holograms of point clouds arbitrarily distributed in the 3D FOV of the optical system. As reported in [19], CS-WGS provides comparable results in terms of spot brightness and uniformity as Gerchberg-Saxton and Weighted Gerchberg-Saxton approaches, while reducing drastically the computational time. For a CGH with N spots, the field at the SLM is calculated as the superposition of the fields generating each individual spot: where u o,n is the amplitude of the field associated with the n th spot at the focal plane, ϕ n is a constant term and ϕ POS,n is a phase describing the position. The position ϕ POS,n is determined by combination of tip, tilt and defocus localizing each of the N points in 3D space: where x n , y n , z n are the coordinates of each point, x ′ , y ′ are the SLM coordinates, f is the focal length of the optical system and λ is the wavelength of the light source. Any iterative calculation of the hologram field aims at maximizing the quality of the hologram by optimizing the value of the phase term ϕ n . The approach here described allows to include in the hologram field (Eq. (1)) a phase term, ϕ AO,n , representing the full-pupil phase applied to the sub-hologram for each point n to compensate local optical aberration [15]. Hence, the CGH phase encoding for a 3D distribution of points with their corresponding aberrations correction is computed as : Here the aberrations are encoded in the CGH using a modal representation. Specifically the aberration phase is expanded as a linear combination of Zernike polynomials as follows: where Z j is the Zernike polynomial function, j is the index of the expansion and w j is the coefficient or weight. Zernike modes are a set of basis functions defined on a circle, in this case corresponding to the objective back aperture. We follow the ANSI order where the index j carries information about the radial (k) and angular (m) indices of the Zernike polynomials. Zernike modes with the same radial index k belong to the same aberration order [20]. These polynomials are widely used [21] and easily applicable. The low order modes are directly linked to identifiable aberrations such as defocus and spherical aberration allowing for clear interpretation of the results. The coefficients w j of the Zernike polynomials are calculated in units of the wavelength and are hereby reported in µm. The 3D CGH can include an isoplanatic phase term ϕ AO,n which is the same across all the points in the field of view. This case is depicted in Fig. 1(b) where a uniform vertical coma aberration (Z 8 ) with coefficient w = 1 µm is superimposed to the CGH on the SLM. The left image shows the phase of the coma aberration addressed to the SLM. The right shows a simulation of a 16-point grid pattern extending over a FOV of 80 µm at focal plane z = 7 µm in the sample, here all the points are affected by the same aberration.
Alternatively, as shown in Fig. 1(c) anisoplanatic aberrations are encoded in the CGH. Here distinct aberration phases are targeted to different location of the FOV. The left image shows different aberration phases, from oblique astigmatism (Z 3 ) to horizontal secondary coma (Z 18 ) with the same coefficient w = 1 µm, for each point of the same grid shown in Fig. 1(b). The right image displays the image of the grid simulated at the focal plane where each point is degraded by different aberrations. The following equation describes the calculation to obtain the simulated images: where e i zπ λ 2 f 2 ·(x ′2 +y ′2 ) is a propagator factor to have the focal plane at a specific z plane, F is the Fourier transform applied to propagate the field from the pupil plane to the focal plane, γ is used to apply a gamma correction to enhance the features of the spots. The simulated images in Fig. 1(b) and 1(c) were generated with γ = 0.5.

Imaging system
The optical setup shown in Fig. 1(d) uses short pulsed near-infrared (NIR) light emitted by a Ti:Sapphire laser mode-locked (Coherent, Mira 900-F) tuned to a central wavelength of 800 nm with 76 MHz repetition rate. Power at the sample is controlled by adjusting the angle between a Glan-Laser Calcite Polarizer (GLP, Thorlabs, GL5B) and a half-wave plate (HWP). The beam is magnified by a telescope (lenses L1, f=25.4 mm and L2, f=150 mm) to fill the chip of the LCoS-SLM (Meadowlark Optics, P1920-600-1300-HDMI, pixel size: 9.2 µm) placed in a plane where the wavefront is nominally flat. The SLM is illuminated obliquely with an angle of 4.6 • to be able to separate the modulated beam from the incident one. A mechanical shutter (S, Uniblitz, LS2S2Z1) installed at the focal point of lens L1 temporally gates the NIR illumination with a response of 2 ms. A second telescope (lenses L3, f=250 mm and L4, f=500 mm) magnifies the beam to fill the back aperture of the water immersion objective lens (O, Olympus XLUMPLFLN, 20X/NA 1.0). In the focal plane of lens L3, an inverse pinhole (IP, custom tungsten deposition, diameter 1.3 mm, on glass window of thickness 0.5 mm) blocks the zero-order diffraction spot. The focal lengths of lenses L3 and L4 are chosen to match the dimension of the SLM image at the back aperture of the microscope objective and the SLM chip. In this case the lateral resolution is only limited by the objective NA. TPEF is collected through the objective, transmitted through a low-pass dichroic mirror (Semrock, FF01-720/SP-25) and finally the image is formed by a tube lens (L5, f=300 mm). The image is relayed to the sCMOS camera (Andor, Zyla 4.2) by a 1:1 telescope (L6 and L7, f=150 mm) and an electrically-tunable lens (ETL, Optotune EL-16-40-TC-VIS-20D) lens which scans the signal at different depths. Brightfield (BF) images are acquired by transmission of a white light LED (Thorlabs, LEDW25E) placed at the opposite side of the objective.
As reported previously [22] the precision to address a CGH to a specific target depends on the number of pixels in the SLM and gray scale values available. We use a LCoS-SLM with 1920 x 1152 pixels so that the number of pixels in the shorter axis permits a spatial accuracy better than 1 µm and also a high number of degrees of freedom when the SLM is used as adaptive element to correct for high-order aberrations. Regions of interest (ROIs) for generating the CGH can be located anywhere within the FOV which is a combination of lateral and axial FOV. The lateral FOV, defined the maximum tilt that can be applied on the first order of diffraction, is 391 µm. The theoretical axial FOV is 166 µm, where the SLM is the limiting factor for generating the CGH. Nonetheless, a trade off between the available lateral and axial FOV exists since increasing the axial displacement of a feature decreases the accessible amount of lateral displacement. Limitations in the applicable corrective patterns arise from aliasing as discussed in Supplement 1 (Section S1).

Sensorless AO correction of the holographic stimulation
AO using direct measurement of the wavefront, e.g., with a Shack-Hartmann sensor, is generally limited to a single isoplanatic patch. For correction over multiple isoplanatic patches in the CGH, we instead chose a modal-based sensorless approach. Each aberration mode is imparted as a trial correction by the AO element and an image of the sample is acquired. The result is a collection of images for each of the optical aberration with a different bias coefficient. The best coefficient is chosen measuring a metric that quantifies image quality.
Sensorless methods require some a priori information about the system, such as a mathematical model describing how the corrective phase influences the chosen metric. For TPEF, the intensity of the fluorescence is known to depend quadradically on the peak intensity of the excitation light. Since we use widefield microscopy to capture fluorescence excited by the reconstructed CGH, we calculate the TPEF intensity of each spot in the pattern as the sum of gray scale pixel values within a bounding region, In the ideal case of diffraction-limited spots in the CGH, the resulting TPEF spots should also be minimized in dimension. Bearing these two criteria in mind, we propose two different metrics M 1 and M 2 , that combine the TPEF intensity with a feature that quantifies the spot shape. The metric M 1 combines the intensity I F and the diameter D, determined as the maximum dimension of a single-spot CGH as shown in the following equation: where the diameter D is taken to be the larger of the vertical or horizontal axes of the single-point image at the focal plane. On the other hand metric M 2 combines the intensity the intensity I F and the area A of a single-spot CGH as shown in the following equation: Both area and diameter are measured after a threshold of 1 e · I max is applied. By dividing the intensity I F with either the diameter or the area of the spot, we ensure that these combined metrics reach the maximum when the spot is more confined. Although the metric M 1 is sensitive to the variation of the spot shape on the lateral direction, the metric M 2 takes into account also any variation of the spot on oblique directions. For such reason metric M 2 is more robust to correct for aberrations that severely impair the shape of the points.
The metric is integrated in an optimization algorithm scheme. The optimization scheme affects both speed and efficiency of the correction. Those are equally important factors to correct aberrations in turbid samples. Here we use an hill-climb algorithm that sequentially loops over the Zernike polynomial base and explores a pre-defined range of bias coefficients for each Zernike polynomial. The corrections found are sequentially applied to each Zernike term before moving to the next one. The optimal coefficient is found by fitting the metric data points with a combination of a Lorentzian function with a line (Eq. (10)). Piston, tip, tilt and defocus are excluded in the optimization process because these are defined in the positioning of the points (Eq. (2)) and may lead to loss of fidelity in the CGH as they cause displacement of the ROIs [23].

Isoplanatic correction
Isoplanatic corrections represent an average correction for the whole FOV. This correction may be optimal for points near the middle of the FOV, but sub-optimal for points near the edges. The pseudo-code in Algorithm 1 and Fig. 2 highlight the main steps to implement such corrections.

Anisoplanatic correction
The pseudo-code for Algorithm 2 gives the key steps of the algorithm for anisoplanatic AO. Being an extension of the previous case to multiple point-like features in the FOV, the anisoplanatic algorithm provides local corrections. This feature is an important factor to compensate spatially Algorithm 1. Isoplanatic correction pseudo-code Zernike polynomial. The corrections found are sequentially applied to each Zernike term before 202 moving to the next one. The optimal coefficient is found by fitting the metric data points with 203 a combination of a Lorentzian function with a line (Eq. 10). Piston, tip, tilt and defocus are 204 excluded in the optimization process because these are defined in the positioning of the points 205 (Eq. 2) and may lead to loss of fidelity in the CGH as they cause displacement of the ROIs [23]. Isoplanatic corrections represent an average correction for the whole FOV. This correction may 208 be optimal for points near the middle of the FOV, but sub-optimal for points near the edges. The 209 pseudo-code in Algorithm 1 and Figure 2 highlight the main steps to implement such corrections. The pseudo-code for Algorithm 2 gives the key steps of the algorithm for anisoplanatic AO. Being inhomogeneous aberrations induced by turbid media such as biological tissues. A drawback compared to the isoplanatic correction process is that the CGH is recalculated after optimization of each aberration mode to encode the corrections just found. This extra step increases the total optimization time by a factor 2.6. To compare the TPEF images before and after AO correction, the intensity of ROI sub-images are 220 normalized with respect to the corrected ROI sub-image ( ) with the following normalization 221 formula: where is the intensity of pixel in the image to normalize and and are pixels with, where is the amplitude, 0 the center of the curve, the width of the curve and the slope 228 of the line. The uncertainty on each fit parameter is given by the square root of the diagonal 229 elements of the covariance matrix [24]. The same function is also used to fit the measured metric 230 values.

Image processing
To compare the TPEF images before and after AO correction, the intensity of ROI sub-images are normalized with respect to the corrected ROI sub-image (I n ) with the following normalization formula: where I i is the intensity of pixel i in the image to normalize and I n min and I n max are pixels with, respectively, minimum and maximum value in the corrected image. The intensity improvement is quantified at each ROI in the normalized image through a line profile across the pixel of maximum intensity. The line profile data points are then fitted with a combination of a Lorentzian function with a line as follows: where a is the amplitude, x 0 the center of the curve, σ the width of the curve and b the slope of the line. The uncertainty on each fit parameter is given by the square root of the diagonal elements of the covariance matrix [24]. The same function is also used to fit the measured metric values.

Sample preparation
Anisoplanatic correction of CGH was performed on both synthetic and biological samples. The synthetic sample was a fluorinated ethylene propylene (FEP, Bola, S2022-04) tube with an internal diameter of 0.8 mm and an outer diameter of 1.6 mm. The FEP tube is attached to a needle (B. Braun, 100 Sterican, 21G x 1 1 2 ′′ ) and to a syringe (Braun, Omnifix F Solo 1 ml) and filled with a diluted solution of fluorescein (Invitrogen, F 1300) dissolved in demineralized water. The colour of the solution was vibrant dark yellow with an optical density of 0.056 in 0.8 mm path length measured through a spectrometer (BioDrop, µLITE). The same type of FEP tube is used to hold biological samples.
The biological sample was made from fixed tissue of zebrafish embryos. Zebrafish embryos were obtained by inbreeding adult fish from the casper strain [25]. Larvae were euthanized at 5 days past fertilization (dpf), fixed in a 1:1 solution of embryo water (60 µg/ml) with 0.0002% methylene blue and 4% para-formaldehyde in phosphate saline buffer (PBS) for 2 hours, rinsed 3x in PBS, and stored in PBS at 4 • C until the experiment. For the experiment they were mounted into a straightened FEP tube in a solution of 1.5% low-melting temperature agarose (Sigma-Aldrich, A4018) in embryo water and fluorescein.

Comparison of sensorless AO metrics for widefield TPEF
We first tested each of the metrics in the isoplanatic correction algorithm (Fig. 2). From each TPEF image, a 40x40 pixel sub-image was retrieved around a single ROI selected for isoplanatic correction ( Fig. 2(a), left). The metric I F was calculated from a 10x10 pixel patch centered around the ROI (Fig. 2(a), middle), while the area to compute the metric M 2 is evaluated on the binary image ( Fig. 2(a), right). This image was obtained by applying a threshold to the original sub-image. A range of 15 bias coefficients in the interval [−3, 3] µm was explored for each Zernike polynomial. Those trial aberrations were sequentially applied to the point chosen for the correction and an image collected as shown in the progression of 40x40 sub-images in Fig. 2(b), top. Once a progression of images as function of the bias aberration was acquired, we measured different metrics using Eq. (6) -Eq. (8). As seen in Fig. 2(c), each metric gives a different optimal coefficient to correct spherical aberration (Z 12 ).

Experimental validation of anisoplanatic AO
To validate the anisoplanatic AO method, we calculated a test pattern with specific aberrations simulated computationally. We then experimentally applied either the isoplanatic and anisoplanatic algorithm to optimize the reconstructed hologram based on its TPEF image. The test pattern consisted of a 16-point grid covering an area of 150 x 150 µm, where spherical aberration (Z 12 ) was introduced at point 1 and point 7 with weights w = -2 µm and w = 2 µm, respectively. This CGH was projected into a chamber filled with fluorescein solution, where sample-dependent aberrations were expected to be minimal. Figure 3 shows the results of both isoplanatic and anisoplanatic algorithms using the intensity-diameter metric, M 1 . For both algorithms, the correction procedure started with Z 12 and then cycled through Zernike coefficients up to the third aberration order (from Z 3 to Z 9 ). This order of modes was chosen a priori because first correcting for dominant modes improves the accuracy of the final correction [26]. In the initial TPEF image of the hologram (Fig. 3(a), left), points 1 and 7 are visibly distorted due to the added spherical aberration. As expected, when isoplanatic correction was performed for point 7, the intensity of point 7 improved by a factor of 10.6 whereas all the other points were degraded by this correction (Fig. 3(a), middle). In contrast, when anisoplanatic correction was performed on all points, the intensity of all spots improved by an average factor of 2.7 ( Fig. 3(a), right).
Analysis of the retrieved Zernike coefficients, shown in Fig. 3(b), showed that the spherical aberrations (Z 12 ) retrieved for points 1 and 7 under anisoplanatic algorithm (orange bars) were w = 2.34 µm and w = −1.87 µm, respectively, while contributions for other orders were non-zero. For comparison of the retrieved and simulated aberrations, we also obtained a baseline anisoplanatic correction of the system using a 4x4 test pattern with no added spherical aberration (blue bars). A slight positive spherical aberration was found in the baseline at the locations of points 1 and 7, respectively, w = 0.26 µm and w = 0.15 µm. This can explain why the retrieved coefficients for spherical aberration is slightly higher than w = 2 µm for point 1, and slightly higher than w = −2 for point 7. We estimated the error on the Zernike coefficient as the standard deviation in the lower orders (Z 3 -Z 9 ) recovered from the three independent corrections for point 7 (baseline, isoplanatic, and anisoplanatic) to be 0.10 µm. Therefore, these data show that the anisoplanatic method accurately retrieves the applied Z 12 coefficient. The anisoplanatic corrective phase patterns applied to the SLM are depicted in Fig. 3(c). Closer inspection of the xy and yz projections of the TPEF image of point 1 (Fig. 3(d)) illustrates the improvement in intensity and shape after the anisoplanatic correction in both the lateral and axial direction. Both images are normalized with respect to the corrected image and the aberrated images are shown with a reduced dynamic range for the pixel values for visualization purposes. The intensity improvement and the spot size are quantified on the intensity profiles (Fig. 3(e)) taken from a line profile across the pixel of maximum intensity. The dashed lines on the images in Fig. 3(e) reveal that spherical aberration displaced the point laterally along the y direction as well as axially. The axial displacement is compensated in the corrected image where intensity is also recovered. The small lateral displacement is most likely due misalignment between the rear pupil of the objective and the SLM plane.
The shape and size of the TPEF spots improved for all the points in the grid. Table 1 shows the reduction in mean lateral and axial dimensions determined as the full width half maximum (FWHM) of a Lorentzian fit to the measured line profile. The FWHM of points focused in the dye solution is an indication of the resolution of the holographic projections because the TPEF image of point-like features are the convolution between the excitation and emission point spread functions (PSFs). Since the sample is a homogeneous solution, TPEF can be generated throughout the geometrical focus. We can take the dimensions of the squared intensity distribution I 2 ex as the upper bound for the spot diameter expected from an ideal Gaussian excitation focus with waist w o = λ πNA and Rayleigh range z R = nw o NA , where n represents index of refraction. The lower bound for the FWHM of the spots detected is limited by the dimension of the Airy disk of fluorescence generated by an infinitesimally small TPE focus, given by R xy ≈ 2 3 1.22λ NA , and depth of field R z = 2nλ NA 2 , where λ em = 520 nm. Noting that the size of the optimized TPEF spots is still larger than either the theoretical beam waist or PSF, we can conclude that the AO optimization of the excitation beam alone does not produce a diffraction-limited image. This can be expected because the system does not correct for aberrations in the pathway of the fluorescence.

Correction on synthetic sample: FEP tube filled with fluorescein
Next we characterized anisoplanatic AO in an aberrating sample, a FEP tube filled with a fluorescein solution (Fig. 4(a)). The FEP tube has an index of refraction of n = 1.344 at 589 nm, which is near to that of water (n = 1.3324) so introduces relatively weak aberrations compared to sample containers made of other common materials (e.g. borosilicate glass) [27]. We expected the cylindrical geometry of the capillary to affect vertical astigmatism [28]. A CGH was calculated for a 4x4 grid with a spacing of 25 µm between points. The TPEF image of the CGH grid projected into the sample is showed in Fig. 4(b), left. On this image we performed the anisoplanatic correction with the metric M 2 , correcting up to the eighth aberration order, corresponding to Zernike coefficients Z 3 -Z 44 . The correction process took 30 minutes, corresponding to calculation of 43 unique CGH and acquisition of 630 images testing 15 weights per order for 42 orders. The resulting corrected image is reported in  Fig. 4(b), right. The corrective phase maps shown in Fig. 4(c) provide a full characterization of the sample-induced aberrations across the central 100 µm of the FOV. As anticipated, vertical astigmatism (Z 5 ) was one of the dominant aberrations given the cylindrical geometry of the sample under examination. (See Supplement 1, Section S2.) Moreover, we also found a prominent contribution for spherical aberration (Z 12 ). As visible from the corrective phase maps (Fig. 4(c)), specific ROIs also had contributions from higher-order aberrations, as in the case of 3 µm coefficient for Z 36 at point 11. The lines profiles across the pixel of maximum intensity for the raw (orange) and corrected (blue) spots represented in Fig. 4(d) reveal an average intensity improvement factor of 2.07 ± 0.33.
We then asked whether the AO corrections obtained for one hologram would be sufficient to optimize a novel point cloud pattern extending over the same FOV. We calculated a second CGH of 12 points randomly distributed over the same 100 x 100 µm area, shown in Fig. 5(b), left, together with the location of reference points from the grid CGH (red dashed lines). For each point in the random pattern, the corrective phase was calculated as a distance-weighted average of the local corrective phases found for all points in the reference grid. The weight coefficient was inversely proportional to the distance between the point to correct and each point in the reference grid CGH, as described in Fig. 5(a) and by the following formula: (11) where w i is the coefficient of a given Zernike mode for a point i in the random pattern, w k are the coefficients of a given Zernike mode for all the k points in the reference pattern, and d ki is the distance between point i and point k. Comparison of the TPEF images of the random points before and after AO correction (Fig. 5) indicates that the distance-weighted anisoplanatic correction improves the spot intensity. The resulting corrective phase maps in Fig. 5(c) are more homogeneous between the random ROIs compared to the ROIs in the reference grid. The weighted average of points from across the whole FOV can be seen as diluting each local correction. Nonetheless, by applying aberrations based on the distance-weighted average correction, the random points increased their average TPEF signal by a factor of 2.07 ± 0.17 on average, as depicted in (Fig. 5(d)). Since the average intensity improvement is comparable to the one obtained on the reference grid (Fig. 4), we can conclude that the reference set of anisoplanatic corrections sampled at 25 µm intervals are sufficient to correct for different patterns covering the original FOV. In a sample with greater complexity, the needed sampling distance can be expected to change. However, for points that are closer together we observed a coupling of the corrections even if the ROIs belong to different anisoplanatic patches. (See Supplement 1, Section S3.) The same approach can be extended to a 3D pattern. To test this experimentally, we first retrieved anisoplanatic corrections for a reference 3D CGH made of 10 points randomly distributed over a volume of 100 x 100 x 50 µm 3 where the axial direction was sampled at a spacing of ∆z = 7 µm (Fig. 5(f)). Here we optimized up to the fourth aberration order, corresponding to Zernike modes in the interval Z 3 -Z 14 , testing 9 bias coefficients in the interval [−2, 2] µm for each mode. This smaller number of optimization parameters still provided good intensity improvement for this sample, while also reducing the number of measurements overall. In this case the correction process took 26 minutes. In three dimensions the correction process was longer compared to the two dimensional case since it was necessary to acquire an image at each depth by changing the ETL voltage. Hence, the optimization time increases linearly with the number of the z depths explored.
As in the 2D case, all points in the 3D reconstruction improved in TPEF intensity at the focal plane. On average the improvement factor was 1.76 ± 0.25, as reported in Fig. 5(g). Figure 5(h) shows that the ROIs after correction are displaced with respect to their positions in the original volume. Displacements were calculated as the distance between the centroids of spots retrieved from the starting image and centroids of spots retrieved from image following AO. The average lateral displacement was 0.6 µm while the average axial displacement was 0.9 µm. Lateral displacement can be introduced by Zernike polynomials (Eq. (4)) with odd k index and m = ± 1, such as coma aberration, while the axial shift can be addressed to radial Zernike modes with even k and m = 0.

Correction on fixed biological sample
Finally, to validate our method on highly aberrated samples we applied the anisoplanatic correction to a CGH projected through approximately 140 µm of tissue in the tail of a 5 dpf zebrafish embryo embedded in a mixture of 1.5% low-melting temperature agarose and fluorescein. As shown in Fig. 6(a), top, each point targets a different anatomical location in the tail. The lateral points are projected through muscles areas whereas the middle point passes through the notochord. The CGH extends over a lateral FOV of 115 µm. The TPEF signal generated behind the tissue, is degraded as shown in Fig. 6(a), second row. The middle point is brighter compared to the lateral ones revealing the inhomogeneity of the tissue.
To compensate for the aberrations we initially applied to each point the local corrections found for the same pattern projected into the FEP tube filled with fluorescein. Those corrections accounting for FEP-induced aberration already improve our pattern projected into the tissue as shown in Fig. 6(a), third row. Further improvement was achieved when we corrected for tissue-induced aberration by using as starting point the pattern corrected for the FEP-induced aberrations. This final result is shown in Fig. 6(a), bottom. To retrieve those corrections we optimized up to the eighth aberration order, corresponding to Z 3 -Z 44 , using the intensity-diameter combined metric M 1 (Eq. (7)). The correction process took 24 minutes. Analysis of the xy profile of TPEF for each ROI, shown in Fig. 6(b), reveals that all three points improve in terms of intensity and shape. However, the correction has a stronger impact on the middle point than the lateral points, as confirmed by the intensity profiles in Fig. 6(c). The average intensity improvement was 3.09 ± 0.26. Figure 6(d) shows the contribution to each local corrective phase in the case of FEP-induced (top) and tissue-induced (bottom) aberrations. In this experiment we observed that the improvement of the metric M 1 as function of Zernike order is different than in the case of FEP-induced aberrations.
Here, as expected, we have an increase of about 2-fold when we use more than 10 Zernike modes and after this increase the correction reaches a plateau with little additional improvement from higher orders (see Fig. 6(e), top). On the other hand, the correction of tissue-induced aberrations starting from the FEP parameters, led to greater contributions from higher-order modes. From the bar plot of the Zernike coefficients in Fig. 6(f), we observed that vertical astigmatism (Z 5 ) has a prominent contribution, as well as higher order aberrations (Z 27 , Z 35 and Z 36 ).
Inspection of the full sequence of bias corrections tested during optimization, shown in the Supplemental document (See Supplement 1, Section S4, Fig. S4), casts doubt whether the hill-climbing algorithm finds a global maximum value across all Zernike modes. For comparison, the same experiment was also carried out by correcting for aberrations directly on the tissue without performing the extra step to compensate for the FEP tube. This approach provided a lower average improvement in the TPEF intensity, limited to a factor of 1.59 ± 0.04, as reported in Supplemental document (See Supplement 1, Section S4, Fig. S5).

Discussion
Anisoplanatic adaptive optics has already been shown to improve image quality in laser scanning confocal microscopy where the anisoplanatic approach is used to correct the one-photon excitation light of scanning lattice microscope using confocal detection [15]. There the anisoplanatic corrections were encoded, together with the spots forming the scanning lattice, as a CGH on a SLM. Differently, in this paper we applied anisoplanatic AO to arbitrary CGH projections using two-photon excited fluorescence as an indirect measure of peak pulse intensity. We demonstrated that correction of local aberrations in point-cloud CGHs projected into either synthetic and biological samples improved the TPEF signals by at least 2-fold. We also characterized the limitations of the approach using widefield camera-based imaging. While the dimensions of features in the corrected point cloud CGHs, quantified as FWHM of the TPEF spots, improved compared to the case before correction, diffraction limited spots were not obtained due to the residual aberrations in the fluorescence emission path of our optical system. In fact, the TPEF images depend on the convolution between the stimulation and the emission PSFs. Since the method described here addresses only aberrations in the stimulation path, images obtained following AO correction still appear aberrated, (e.g., the quality of AO-corrected spots projected into tissue in Fig. 6(a) is lower than spots projected into less complex sample of Fig. 4(b)). Remaining distortions could be compensated by installing an adaptive element along the detection path [15].
In this paper anisoplanatic adaptive optics is performed in a sensorless scheme, hence it does not require any additional component to measure the wavefront. This scheme supports a very simple optical setup where the adaptive optics element is the same SLM used to project and correct CGHs. Aberration corrections are simply included as a modification of a target CGH. This aspect highlights an advantage of sensorless approaches over direct wavefront sensing. In the latter, the wavefront distortion is directly measured (e.g. with a Shack-Hartmann sensor) on the light emitted by a guide star or on light back-scattered from the sample. Measured distortions are corrected through an adaptive element in a plane conjugated with both the objective rear pupil and the wavefront sensor. As demonstrated by Wang, et al. [29], direct wavefront sensing enables to recover diffraction limited resolution in living zebrafish embryos over large volumes (240 µm 3 ) in both laser scanning confocal and two-photon excited fluorescence microscopy. Nonetheless, wavefront sensors add an extra layer of complexity to the existing microscope and can introduce non-common-path errors. Scattering could make a wavefront sensor solution unfeasible when the scattering components prevail over the ballistic components in the wavefront measurement. Most importantly, Shack-Hartmann sensors measure the aberration of only one isoplanatic patch at a time [10] making it difficult to implement for anisoplanatic field corrections.
Notably, a recent work by Ancora and colleagues [30] presents a new method to measure locally the wavefront in a widefield microscope pupil and to apply anisoplanatic correction of aberration through a multi-region deconvolution approach. The aberration detection is based on a spinning subpupil aberration measurement (SPAM) where a subaperture is scanned across the microscope's pupil via a motorized device and it is able to measure local PSFs. This module can be integrated in the detection path of any fluorescent microscope as described in [31] where it is integrated into a light sheet microscope. Here, the SPAM module operates in closed loop with a deformable lens to correct the aberrations in the center of the FOV while the residual distortions at the edges are compensated with the multi-region deconvolution approach.
Sensorless approaches have been realized in many flavors, from pupil segmentation zonal methods to modal-based AO. In pupil segmentation strategies, the objective pupil is divided into sub-regions to measure and correct aberrations. Those sub-regions are sequentially illuminated to measure the local wavefront slope from the image shift. Ji, et al. [8] implemented this technique in a two-photon microscope and corrected for complex aberrations in mouse cortical slices at 400 µm depth to provide near diffraction-limited resolution imaging. A generalized framework for sensorless approaches by Booth, et al. [32] concluded that in general pupil segmentation zonal methods, whereas they provide better performances for high order aberration, they are less robust to noise than modal methods.
As shown in a recent work by May, et al. [18], sensorless sample-conjugated AO enables to correct scattering on multiple anisoplanatic patches distributed over large FOVs for application on two-photon imaging on both fixed and living microglia in mouse hippocampal tissue. Here local corrections are provided within short time scales on the order of 10 s. Such correction speed is achieved thanks to the algorithm used for the correction search, which is a fast converging phase retrival method called Dynamic Adaptive scattering compensation Holography (DASH) [33]. Moreover the use of the photomultiplier tube (PMT) detector enables to collect fluorescent signals at a faster rate than an sCMOS camera.
In the work presented here, the speed of the optimization is constrained by the number of images acquired to assess the optimal correction as well as by the complexity of the aberration. When the sample introduces severe distortions an higher number of Zernike modes has to be included increasing the optimization time. The speed factor is particularly limiting when dealing with samples where the fluorescent signal comes from a limited pool of proteins. Here the results can be impaired by sample drift, photobleaching and other phenomena resulting in loss of intensity during the optimization process, as reported in Supplemental document (See Supplement 1, Section S5, Fig. S6). Moreover a faster optimization translates in a decreased exposure of the sample. This is a critical aspect when handling living samples where photodamage, both thermal and photochemical, could impair the normal physiology of the sample under investigation [34]. To speed up the correction process, we could retrieve an anisoplanatic correction on a reference CGH that samples the FOV and then apply this reference correction to other patterns, as demonstrated in Section 3.3. Applied to the context of holographic photostimulation in living tissue, the reference corrections would best be retrieved on the TPEF signal from guide-star fluorophores (e.g., beads injected in the tissue). The caveat of this solution is that the guide-stars need to be adequately distributed over the same FOV of the target ROIs to provide sufficient sampling of the local aberrations.
Another improvement to increase the speed of the correction process on complex samples, was to perform a step-wise correction of first the sample holder (e.g., FEP tube) and use these values as a starting point for tissue correction. We found this approach improved the overall performance of the optimization process compared to a single-step correction of the full sample.
The AO algorithm itself may be sped up in several ways. For example, the calculation of CGHs can be accelerated by use of a GPU, as recently proposed by Pozzi, et al. [35]. During the optimization process an updated CGH is calculated after each Zernike mode correction. This step takes several seconds and this dead time could be reduced in the order of milliseconds by introducing the GPU calculation. Moreover, a non-iterative algorithm based on convolutional neural network with unsupervised learning is also available to compute CGHs at video rate reducing at least 10 times the computational time with respect to standard iterative algorithms as Gercheberg-Saxton [36]. A further possibility to speed up the whole correction algorithm is to replace the sequential search of Zernike modes, which is quite inefficient, with a smarter search algorithm. As shown in Supplemental document (See Supplement 1, Section S4, Fig.  S4), most Zernike coefficients tested in the hill-climbing search algorithm yielded low metric values, whereas it may be more efficient to vary more than one coefficients value in combination, or in non-sequential order. For instance, in the work of Verstraete, et al. [37], the range of the coefficients explored for each Zernike mode is based on the previous measurement with a random perturbation to ensure an adequate exploration of the bias coefficients space.
Another limiting factor of our method is cross-talk between different modes. This drawback arises from lack of orthogonality in the polynomial base used to model the aberrations. In fact, orthogonality is a necessary property to express the metric as a quadratic function of the chosen base [23]. In our case, Zernike polynomials are not an orthogonal basis with respect to either of our combined metrics M 1 and M 2 , and the order in which the modes are optimized and applied matters. Other aspects to consider are the accuracy of the correction, which is affected by different parameters such as noise and cross-talk between the modes. The major factors contributing to noise in such experiments are Poisson photon noise and changes of the sample over the optimization time, mainly due to sample drift and photobleaching. Noise in the image can induce the algorithm to find a wrong maximum. This effect is mitigated when the maximization of the metric is realized through a fitting process [38], but the number of data points per mode also affects the accuracy of the fit and, ultimately, the performance of the AO. As discussed in [26], the number of measurements for each point needs to be greater than the number of fit parameters.
Finally, we note some displacement of the spots after correction. For instance, with the distance-weighted average anisoplanatic correction to a 3D CGH, we observed that the corrected points were displaced laterally and axially from their original position. The shift is likely due to the polynomial base used, which does not guarantee invariance of the geometrical center of the PSF during the correction. This is a common problem in non-linear excitation schemes where the the PSF is non-linearly proportional to the light intensity. A solution to this problem could be to develop a polynomial base that conserves the first moment of the square intensity or implement an experimentally-calibrated shift-less basis set [39].

Conclusions
Sensorless adaptive optics techniques have been extensively employed to improve signal and optical resolution in the context of imaging though turbid samples. In this paper we apply sensorless AO techniques to correct local aberrations in the context of projection of multi-point holograms through turbid media. We demonstrated an anisoplanatic AO approach that can compensate for both system and sample-induced aberrations in synthetic and biological samples. Our strategy does not require additional elements to an existing optical setup for holographic projections, but rather uses optimization of two-photon excited fluorescence as a guide star.
We can envision applications in optical trapping, parallel laser writing through two-photon polymerization and parallel photostimulation of optogenetics actuators in living cells or intact animals. If integrated with a light sheet microscope as a setup for optogenetics, the stimulation module presented can simultaneously excite multifocal points targeted to regions of interest within the sample, while the light sheet microscope provides flexibility to readout fluorescence signals, e.g., image calcium activity from small samples, such as larval zebrafish.