## 1.

## Introduction

Optogenetics has become a central tool in neuroscience research.^{1} The ability to remotely stimulate cells with cell-type selectivity and high spatial and temporal resolutions^{2} provides major advantages over earlier experimental perturbation methods, such as electrical stimulation. Optogenetic tools are currently also being developed for applications ranging from vision restoration devices to implantable brain interfaces.^{3}^{–}^{5} Interestingly, computer-generated holographic (CGH) projection systems have emerged as a particularly attractive option for delivering precisely targeted light to optogenetically transduced cells,^{3}^{,}^{6}^{–}^{8} and specifically for the simultaneous stimulation of multiple cells. CGH uses phase-only spatial light modulators (SLMs) for spatially adjusting the wavefront’s phases,^{9} and a lens, which under Fraunhofer diffraction conditions reconstruct a Fourier transform (FT) of the wavefront. The wavefront’s FT is projected onto the optogenetically transduced cells, allowing high rate dynamic generation of distributed patterns for neural stimulation. The use of such holographic optical neural interfaces (HONIs) for photo-stimulation offers substantial benefits over alternative amplitude modulation approaches in terms of power efficiency and maximal irradiance, and ultimately enables safer designs.^{3}^{,}^{6}^{,}^{10} Specifically, for multiphoton stimulation and imaging light efficiency are crucial and HONIs have emerged as the key enabling solution.^{7}^{,}^{8}

As SLMs only modulate the input beam’s phases, CGH requires phase retrieval algorithms for calculating the appropriate phase image to project, based on its target FT reconstruction. For HONIs and many other applications, it is important that the projected pattern will be diffraction efficient, smooth, and uniform. However, conventional iterative CGH algorithms, such as the Gerchberg–Saxton (GS) algorithm and its common derivatives,^{11}^{,}^{12} only constrain the modulus of the projected pattern, allowing its phases to vary randomly. The weighted GS algorithm (GSW), for example, reaches near perfect modulus uniformities and high efficiencies by introducing weights into the projected domain’s constraints but lacks any phase control. These random phases lead to holographic speckles, a phenomenon that arises from the projection of contiguous patterns with varying phases: the spots can interact destructively, which causes high contrast regions to appear in the projected pattern^{13} [illustrated in Fig. 1(a)]. Holographic speckle contrast is worsened by the nonlinearity of two-photon excitation and is recognized as a major challenge throughout the literature on scanless holographic neurophotonic interfaces^{6}^{,}^{7} (that is, excluding systems where single diffraction-limited spots are being projected and mechanically scanned^{8}), thus limiting the performance of the key enabling solution for precision optogenetics. Avoiding holographic speckle has even motivated several key studies exploring generalized phase contrast, an alternative speckle-free diffractive projection method, which is essentially two-dimensional (2-D) and has much lower efficiency.^{14}^{,}^{15}

Several time averaging solutions to the holographic speckle problem were proposed^{16}^{–}^{18} and were shown to provide a reduction of speckle contrast in one- and two-photon holographic projections. However, time averaging requires the projection of multiple holograms to reduce speckle noise, thereby sacrificing the temporal resolution. Alternative approaches to speckle reduction are based on jointly controlling both the amplitude and phase of the projected patterns [Fig. 1(a)], an attribute that has multiple alternative motivations. Physical schemes for realizing this through controlled spatial modulation of the hologram’s phase and amplitude were developed,^{19}^{,}^{20} but their implementations require strict alignment and calibration of the optical components. Algorithmic solutions that achieve phase and modulus control are a possible alternative. Early attempts at speckle reduction focused on applying phase-smoothing constraints in the projected plane^{21}^{,}^{22} but required over-sampling of the pattern, and postiterations with reduced phase freedom, thus increasing the computational load by orders of magnitude. Another potential approach for achieving phase control is to divide the projected plane into constrained and constraint-free regions,^{23} but this results in low efficiencies and uniformities in Fourier-based holographic systems. An important recent contribution^{24}^{,}^{25} addresses this challenge by incorporating phase information into the GS-type iteration process. However, these algorithms do not explicitly maximize the accuracy of the pattern moduli, and require to manually set a predefined rate parameter that affects the diffraction efficiency. Here, we introduce a method for modulus and phase control of projected patterns, which allows projection of speckle-free patterns with high diffraction efficiency. The algorithm performs well given arbitrary one-dimensional (1-D) curved line patterns, with no need for pattern-specific parameter adjustment.

## 2.

## Gerchberg–Saxton with Phase Control

## 2.1.

### Algorithm Description

As schematically shown in Fig. 1(b), GSW-PC follows the general procedure of the GSW algorithm.^{12} Dual FTs between the projected image domain ($V$) and the hologram domain ($H$) are used to satisfy each domain’s constraint in each iteration, while a set of weights ($W$) is used to iteratively correct the pattern moduli in $V$. However, in GSW the projected domain constraints are solely on the modulus of the projected pattern ($|{V}_{\mathrm{init}}|$) while individual pattern phases are randomized to promote convergence. To enforce phase control, GSW-PC adds an additional constraint ($\angle {V}_{\mathrm{init}}$) on the phases of the projected pattern [Eq. (2)]. This step uses a solution suggested by Yuan and Tao,^{24} where the nonpattern (background) pixels of the projected image are not set to zero modulus, so that they preserve their phase information between iterations. As described below, we empirically found that this procedure works very well for 1-D curved line patterns. Following Ref. 24, a control parameter ($\delta $) adjusts how much relative intensity is deflected toward the projected pattern versus the background. Importantly, in GSW-PC, this parameter is allowed to decrease iteratively based on the pattern uniformity levels in each iteration [see Eqs. (7)–(10)]; this maximizes pattern efficiency and provides image-independent parameter optimization.

The algorithm implements the following domain constraints:

With the weights calculated using

## Eq. (3)

$${W}^{0}=J,\phantom{\rule{0ex}{0ex}}{W}^{k}={W}^{k-1}\circ |{\widehat{V}}_{\mathrm{init}}|./|{\widehat{V}}^{-}|,$$## 2.2.

### Three-Dimensional Light Shaping Using the Angular Spectrum Method

Multiplane projection is likely to have an important role in optogenetics for simultaneous stimulation of multiple cells dispersed over a 3-D volume, as well as for optimizing individual cells stimulation by projecting 3-D patterns that wrap around each cell’s membrane and thus activate a large number of photo-sensitive channels [shown in Fig. 2(a)], which adds an additional challenge to the phase retrieval problem. The angular spectrum method (ASM) is an approach used for modeling the propagation of a wave field in free space, based on scalar diffraction theory. The method uses fast FTs to decompose an optical wavefront into multiple plane waves with different spatial frequencies, and to superimpose them in the propagated plane. The general algorithm for the ASM is as follows:

where $H$ is the transfer function of the angular spectrum, and can be expressed as where $z$ is the propagation distance, $\lambda $ is the projected wavelength, $k=\frac{2\pi}{\lambda}$ is the wavenumber, and $u$, $v$ are the spatial frequencies.To implement GSW-PC-based 3-D projection, we followed the general approach described by Wu et al.^{25} [see Fig. 2(b)]. The iteration process starts with projection of the hologram to the projected Fourier plane (i.e., the lens focal plane) using a FFT (step 1). The projected wavefront is then propagated from the Fourier plane to the other axially shifted planes using the ASM (step 2), where the constraints of the individual planes are applied on the wavefield. The wavefield of every plane is then propagated back to the Fourier plane using the ASM (step 3), and subsequently propagated backward to the hologram plane using an IFFT (step 4). The phases of the wavefields are eventually summed, and the process repeats until a quality criterion is met. This process returns a 2-D hologram that produces 3-D structures in its projected domain.

## 2.3.

### Performance Measures

The following quality measures were used to determine the performance of the algorithm (adapted from Ref. 12).

Diffraction efficiency

Uniformity measures (in the case of constant modulus and phase projections)

## Eq. (7)

$$\mathrm{mod}.\text{uniformity}=1-\frac{\mathrm{max}\{|{I}_{p}|\}-\mathrm{min}\{|{I}_{p}|\}}{\mathrm{max}\{|{I}_{p}|\}+\mathrm{min}\{|{I}_{p}|\}},$$## Eq. (8)

$$\text{phase uniformity}=1-\frac{\mathrm{max}\{\angle {V}_{p}+\pi \}-\mathrm{min}\{\angle {V}_{p}+\pi \}}{\mathrm{max}\{\angle {V}_{p}+\pi \}+\mathrm{min}\{\angle {V}_{p}+\pi \}},$$In the case of projection of nonconstant patterns, alternative normalized MSE (NMSE) were used as measures of the phase and amplitude difference between the projected pattern and the target pattern. As defined in Eqs. (9) and (10), the MSEs were normalized by dividing them by the maximal possible MSE value for the given patterns, allowing to set nominal generalized uniformity values (required to update $\delta $) irrespective of the input pattern

## Eq. (9)

$$\mathrm{mod}.\mathrm{NMSE}=1-\frac{\mathrm{mod}.\mathrm{MSE}}{\mathrm{mod}{.\mathrm{MSE}}_{\mathrm{max}}}=\frac{\frac{1}{n}\sum _{p=1}^{n}{(|{\widehat{V}}_{p}^{-}|-|{\widehat{V}}_{{\mathrm{init}}_{p}}|)}^{2}}{\frac{1}{n}\sum _{p=1}^{n}{(\mathrm{max}\{|{\widehat{V}}_{{\mathrm{init}}_{p}}|,1-|{\widehat{V}}_{{\mathrm{init}}_{p}}|\})}^{2}},$$## Eq. (10)

$$\text{phase}\text{\hspace{0.17em}}\mathrm{NMSE}=1-\frac{\text{phase}\text{\hspace{0.17em}}\mathrm{MSE}}{\text{phase}\text{\hspace{0.17em}}{\mathrm{MSE}}_{\mathrm{max}}}=\frac{\frac{1}{n}\sum _{p=1}^{n}{(\angle {\widehat{V}}_{p}^{-}-\angle {\widehat{V}}_{{\mathrm{init}}_{p}})}^{2}}{\frac{1}{n}\sum _{p=1}^{n}{\pi}^{2}}.$$## 3.

## Numerical Results

The following subsections show numerical simulation results examining GSW-PC performance in different core scenarios (algorithms were implemented in MATLAB).

## 3.1.

### Uniform-Intensity Patterns

Figure 3 shows the process of optimizing the different performance measures during the iterations of the GSW-PC algorithm for different input images. During the preadaptation period, the $\delta $ value is high (initially set to 0.1), allowing the algorithm to rapidly converge to almost perfect modulus and phase uniformities but with a suboptimal efficiency. After these high uniformities are reached, adaptation commences and $\delta $ begins to decrease (by a factor of 0.9 each step), causing the efficiency to rise; $\delta $ continues to decrease until the uniformity drops, wherein the iteration process is stopped.

Next, we compared the algorithm to related alternatives including: (1) a nonadaptive (constant $\delta $) version of GSW-PC, (2) GSW, and (3) Yuan and Tao’s “phase and amplitude” algorithm^{24} (referred to as YT). The comparative results (Fig. 4) are somewhat expected: GSW, which has no control over the projected pattern’s phases, leads to null phase uniformity, whereas YT achieves good phase uniformity and efficiency, but low modulus uniformity. By adjusting the pattern weights, the nonadaptive GSW-PC dramatically improves the modulus uniformity compared with YT. Finally, the complete GSW-PC (adaptive) algorithm reaches high phase and modulus uniformities and optimal efficiency adapted to each given input 1-D pattern.

## 3.2.

### Nonuniform Patterns

Projection of nonuniform patterns is an important feature in many applications, e.g., in optogenetics it can be used for compensating for nonuniformities caused by variable expression levels of the photo-sensitive channels as well as by the inherent nonuniformity of the projection system.^{3} Figure 5 shows results achieved with the GSW-PC method in projecting patterns that are nonuniform in both intensity and phase: amplitude and phase nonuniformities were introduced through ${V}_{\mathrm{init}}$. The results clearly indicate that even in this more challenging scenario, the method is generally able to combine very high pattern controllability with high efficiencies (approaching 50%).

## 3.3.

### Three-Dimensional Holography

Figure 6 shows results of testing the multilayer CGH capability of the 3-D GSW-PC algorithm by projecting 3-D pattern of rings with varying diameters. This pattern can be designed to match the membrane’s geometry for efficient multiphoton optogenetic stimulation of single cells. The adjacent rings were calculated at a distance of $z=2\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ from each other, with a projected wavelength of $\lambda =920\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{nm}$. The performance measures results for each plane show that the solution reaches near perfect phase and modulus uniformities, with efficiencies that are relatively high for multiple layer projection. In comparison, generating the same 3-D pattern using the algorithm presented in Ref. 25, led to a highly nonuniform pattern ($<70\%$ modulus uniformity values for generally similar efficiencies, results not shown).

## 4.

## Experimental Results

A two-photon holographic projection system [Fig. 7(a)] was used to test the performance of the GSW-PC algorithm. An expanded beam from an ultrafast titanium-sapphire laser (MaiTai WB, Spectra-Physics, set to $\lambda =920\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{nm}$) had its polarization rotated using a $\lambda /2$ waveplate, before hitting a phase SLM (XY Phase, Boulder Nonlinear Systems) and a Fourier lens. The zero-order spot was blocked in the projection plane, and this plane was imaged through a custom microscope (objective: $40\times $, 0.8 NA, Apo NIR; Nikon) onto a sample containing a Fluorescein solution. The fluorescence image resulting from two-photon excitation was captured using a CCD camera (GC1380H, Prosilica).

We used the system to project holograms calculated using GSW and GSW-PC for several different patterns (circles, curved lines, and spirals). As can be clearly seen in Fig. 7(b), images obtained using GSW are strongly speckle contaminated, whereas the GSW-PC patterns do not exhibit this phenomenon. Quantitatively, the GSW patterns had 50% to 100% higher speckle contrast^{26} than the corresponding GSW-PC projections [Fig. 7(c)], using the following measure for contrast:

## Eq. (11)

$$\text{contrast}=\frac{1}{n}\sum _{i=1}^{n}\frac{\mathrm{std}({p}_{i})}{\text{mean}({p}_{i})},$$## 5.

## Discussion

We have presented a CGH method for controlling both the phase and amplitude of the projected pattern. The GSW-PC algorithm allows projection of arbitrary 1-D curved line patterns with high efficiency and uniformities, and thus provides speckle-free projections using a single hologram in either 2-D or 3-D (e.g., using the ASM adaptation). The algorithm was tested with a 2P projection system, and demonstrated a marked reduction of speckle contrast compared with the GSW method, at the cost of only a moderate increase in computational load (calculation of the quality measures in each iteration, and some extra iterations for algorithm convergence). Numerical simulations of GSW-PC in projection of 3-D target patterns with upward of 20 planes have also been tested, and have not shown any systematic drop in the algorithms performance with increasing number of planes (data not shown).

We anticipate that GSW-PC will provide significant benefits for a variety of CGH applications, due to different aspects of its ability to allow uniquely precise and flexible control over properties of the projected patterns. The likely benefits for patterned neurophotonic holographic interfaces, already extensively used for stimulation and imaging, are potentially substantial: rendering excitation patterns to be speckle-free does not only allow for more highly controlled and smooth, less variable, stimulation intensity profiles, but also avoids hot spots where intensities could be high enough to cause bleaching, membrane poration,^{27} or other forms of tissue-damage. Collectively, the solution thus addresses issues that both limit stimulation efficacy by constraining the stimulation intensities used (to avoid damage) and lead to less reproducible results during extended precision optogenetic stimulation experiments. Rendering the patterns to be speckle-free does not only allow for more highly controlled and smooth stimulation intensity profiles but also avoids “hot spots” where intensities could be high enough to cause bleaching, membrane poration,^{27} or other forms of tissue-damage, which collectively may lead to less reproducible results during extended precision optogenetic stimulation experiments. Avoiding laser speckle is also highly desirable and has been intensely researched in holographic laser machining, processing and polymerization^{28}^{–}^{30} and in holographic 3-D display systems.^{31}^{,}^{32} Furthermore, GSW-PC’s ability to directly and flexibly control the pattern’s phase could be used in optical tweezers setups for manipulating particles along arbitrarily shaped phase ramped 1-D patterns^{20}^{,}^{23} (optionally together with flexible amplitude variations).

One noted limitation of the proposed approach is its ability to generate only 1-D projection patterns (i.e., curved lines, rings, spirals, etc.) with high efficiencies. To directly examine this issue, we evaluated the different performance measures as a function of the projected pattern width [see Fig. 8(a)], finding that efficiency drops very strongly with pattern width. A similar drop in efficiency was observed for the projection of straight lines (not shown). We believe that the lack of randomness in the projected pattern’s phase is the cause for this behavior due to inherent properties of the complex FT, which decomposes a complex image with randomized phase into a relatively uniformly distributed wide range of frequencies. This allows algorithms such as the GSW to converge to high efficiencies when imposing a uniform modulus in the hologram domain. On the other hand, a complex image with regions of constant phase is decomposed by the FT into a nonuniform distribution of frequencies, with high-energy residing in its low frequencies. This nonuniform distribution prevents the GSW-PC algorithm from converging to high efficiencies when projecting 2-D constant phase patterns. Nevertheless, Fig. 8(b) shows that regardless of the efficiency, 2-D wide patterns projected with the GSW-PC algorithm continue to display speckle-free results, compared with patterns projected with the GSW algorithm, due to high uniformities level conservation [see Fig. 8(a)].

## Disclosures

The authors declare no relevant financial or competing interests, nor other potential conflicts of interest.

## Acknowledgments

The research was financially supported by the European Union’s Horizon 2020 research and innovation program under ERC grant agreement 648927, and by the U.S. National Institutes of Health (Award No. 1U01NS090498-01).

## References

## Biography

**Tal Aharoni** received his BSc degree in biomedical engineering from the Technion-Israel Institute of Technology and his MSc degree in autonomous systems and robotics in 2014 and 2017, respectively. Currently, he leads the algorithms team at Medic Vision imaging Solution Ltd.

**Shy Shoham** received his BSc degree in physics from Tel Aviv University and his PhD in bioengineering from the University of Utah, and postdoctoral training at Princeton University. He is a professor of ophthalmology/neuroscience at NYU Langone Health and of biomedical engineering at the Technion. His lab develops photonic, acoustic and computational tools for spatiotemporal interfacing with neural circuits. He serves on the editorial boards of the Journal of Neural Engineering and of Translational Vision Science and Technology.