Enabling second harmonic generation as a contrast mechanism for optical projection tomography (OPT) and scanning laser optical tomography (SLOT)

Volumetric imaging of connective tissue provides insights into the structure of biological tissue. Second harmonic generation (SHG) microscopy has become a standard method to image collagen rich tissue like skin or cornea. Due to the non-centrosymmetric architecture, no additional label is needed and tissue can be visualized noninvasively. Thus, SHG microscopy enables the investigation of collagen associated diseases, providing high resolution images and a field of view of several hundreds of μm. However, the in toto visualization of larger samples is limited to the working distance of the objective and the integration time of the microscope setup, which can sum up to several hours and days. A faster imaging technique for samples in the mesoscopic range is scanning laser optical tomography (SLOT), which provides linear fluorescence, scattering and absorption as intrinsic contrast mechanisms. Due to the advantages of SHG and the reduced measurement time of SLOT, the integration of SHG in SLOT would be a great extension. This way SHG measurements could be performed faster on large samples, providing isotropic resolution and simultaneous acquisition of all other contrast mechanisms available, such as fluorescence and absorption. SLOT is based on the principle of computed tomography, which requires the rotation of the sample. The SHG signal, however, depends strongly on the sample orientation and the polarization of the laser, which results in SHG intensity fluctuation during sample rotation and prevents successful 3D reconstruction. In this paper we investigate the angular dependence of the SHG signal by simulation and experiment and found a way to eliminate reconstruction artifacts caused by this angular dependence in SHG-SLOT data. This way, it is now possible to visualize samples in the mesoscopic range using SHG-SLOT, with isotropic resolution and in correlation to other contrast mechanisms as absorption, fluorescence and scattering.


Introduction
For imaging biological samples, like skin or cornea, second harmonic generation (SHG) has become a standard method [1]. Non-centrosymmetric components, as collagen and myosin can be visualized without the need of an external stain, due to the noninvasive nature of this process. Furthermore, the penetration depth in scattering media is increased by using wavelength in the near infrared, compared to linear fluorescence imaging, where usually light in the visible range is used. This way 3D imaging of collagen in turbid media becomes possible [1,2]. SHG is a parametric process that creates no energy deposition inside the sample which reduces photo damage and eliminates effects like photo bleaching [3]. This way SHG imaging has been used successfully for diagnostics by imaging cornea [4], atherosclerotic arterial tissue [5], skin [6,7] or cancer [8][9][10][11]. Even more structural information can be gained by using polarization resolved SHG imaging. Due to the polarization dependence of the SHG signal, the 3D orientation and structural changes of collagen can be identified [12][13][14][15][16][17][18]. Thereby, most imaging setups are designed for high resolution imaging, with a small field of view using high NA objectives, with relatively small working distances. This limits the in toto visualization to small samples or to small subunits of larger samples. To visualize samples in the mesoscopic range, objectives with larger working distances can be used, but still suffers from long acquisition times.
A method, developed for 3D imaging of samples in the mesoscopic range, has been shown for the first time in 1992 by Brown et al. [19]. This method is based on the principle of computed tomography using visible light instead of x-rays and arose from the work of Kawata et al. [20][21][22], who applied the method in 1985 to microscopic samples. In 2002 Sharpe et al., who were unaware of the work of Brown and Kawata, published and patented the optical computed tomography and called it optical projection tomography (OPT) [23]. Projection images are acquired during sample rotation and 3D information can be reconstructed by the inverse radon transform. Scanning laser optical tomography (SLOT) is a further development of OPT, with improved collection efficiency [24]. SLOT allows the usage of scattering, linear fluorescence and absorption as intrinsic contrast mechanisms. With this method the successful imaging and segmentation of the murine lung and cochlea was realized in toto after optical clearing [25][26][27]. Furthermore, surfaces of non-transparent samples have been imaged, such as biofilm on metal surfaces and implants [28].
A valuable extension of the SLOT setup would be the integration of SHG as a new contrast mechanism. This way all the benefits of SHG mentioned above would complement the SLOT technique and SHG images could be directly correlated to linear fluorescence, scattering and absorption. Furthermore, larger samples in the mesoscopic range, can be visualized in toto with a shorter acquisition time and isotropic resolution. However, the dependence of the SHG signal on sample orientation and laser polarization, which is the basis for polarization resolved SHG, becomes a major challenge when 3D reconstruction is performed. SLOT uses angular projections of a full 360 • rotation to reconstruct volumetric images. The underlying computed tomography reconstruction process (the inverse Radon transform) requires, that the intensity of each microscopic volume element (voxel) is independent on the rotation of the sample. This condition is fundamentally violated by SHG. This is because the SHG signal of a scatterer strongly depends on the relative orientation of the polarization of the illuminating light and the scattering objects axis of symmetry. During sample rotation this orientation will constantly change and typically result in a fluctuating intensity of each scatterer. Reconstruction artifacts will result, if no measures are taken to compensate the rotational dependence of the SHG signal.
In this paper we simulate the SHG intensity, that is generated by collagen, for different laser polarizations and collagen orientations during sample rotation. Furthermore we simulate the according resulting artifacts in the reconstructed image using an area of circle as a phantom. The extent of these artifacts are measured and displayed for all polarization angles and sample tilts. Subsequently, we performed SHG-SLOT measurements of a rat tail fascicle and a rat lymph node and compare these with the simulations. Finally, we provide a solution for the usage of SHG as a contrast mechanism in radon based imaging techniques as SLOT and OPT.

Theory
For the following calculations and measurements, a reference coordinate system needs to be defined with the following three angles. The polarization of the incident electromagnetic field ì E is tilted by polarization angle α to the x-axis in the x-y-plane. The coordinate system of the collagen strand (cx, cy, cz) is tilted by the tilt angle θ around the y-axis and by the rotation angle ϕ around the x-axis (see Fig. 1). Light propagates in z-direction and the sample rotation occurs around the x-axis. E (red arrow) and collagen orientation (green rod). Light propagates in z-direction. The laser polarization is placed in the x-y-plane and is tilted by the polarization α to the x-axis. The orientation of the collagen fiber is specified by the tilt angle θ and the rotation angle ϕ. The sample rotates around the x-axis during SLOT measurements. To obtain a better overview, all relevant axes and angles are listed in the right side of the figure.

The second order susceptibility of collagen
A second harmonic signal can only be generated, if the second order susceptibility tensor of the sample χ (2) i jk is non-zero, which is only the case for noncentrosymmetric media. Due to intrinsic permutation symmetry the last two indices of χ (2) i jk can be freely interchanged [29] which enables the usage of the contracted description of the nonlinear susceptibility. This is given by d is with i = {1;2;3} and s = {1;2;...6} , where 1 = (xx), 2 = (yy), 3 = (zz), 4 = (yz), 5 = (xz) and 6 = (xy) [14]. This reduces the number of independent coefficients from 27 to 18. Furthermore, cylindrical symmetry can be assumed for collagen with cx being the axis of symmetry (see Fig.  1). Under this condition, the only remaining non vanishing components are d 11 , d 12 = d 13 and d 26 = d 35 [17]. For interactions far from resonance, the susceptibility is independent of the frequency which leads to free permutation of all indices i, j and k (known as Kleinmann symmetry [14]). This allows the following simplification: (1) The ratio ρ = d 11 /d 12 depends on the architecture of the collagen and varies in literature between 1.3 and 1.4 [30-32]. The following simulations where all performed for ρ = 1.3. Varying this value in the mentioned range results in the same simulation results that vary only slightly in amplitude. Furthermore the first component (d 11 ) was set to 1.

SHG signal strength in respect to the laser polarization and collagen orientation
The induced SHG polarization ì P SHG depends on the susceptibility tensor d and the incident electro-magnetic field ì E c in the coordinate system of the collagen strand (cx, cy, cz): where 0 is the vacuum permittivity [29]. To calculate the electro-magnetic field in the coordinate system of the collagen strand, the incident electromagnetic field ì E needs to be transformed from the lab frame work (x, y, z) into the coordinate system of the collagen fiber (cx, cy, cz). To do so, the rotation matrix M = R y (θ) · R x (ϕ) composed of a rotation by the tilt angle θ around the y-axis and a rotation by the rotation angle ϕ around the x-axis needs to be multiplied with ì E: Combining equations (2) and (3), the resulting induced SHG polarization can be calculated for any combination of α, θ and ϕ. The total radiated SHG intensity is proportional to | ì P SHG | 2 : where c is the speed of light, k the SHG-wavenumber, V the sample Volume, and 0 the vacuum permittivity [33]. All these variables are independent on α, θ and ϕ and were combined to (ck 2 V 2 )/(12π 0 ) = a = 1 for simplification.

Numerical simulations
For a complete SLOT measurement, the sample needs to be rotated by 360 • orthogonal to the optical axis. This results in a variation of rotation angle ϕ from 0 • -360 • . The SHG intensity was calculated numerically as a function of ϕ using equations (2), (3) and (4) for different contributions of the polarization angle α and the tilt angle θ. Numerical simulations were performed in GNU Octave [34] (see Fig. 2). Depending on the polarization angle α and the tilt angle θ, the SHG intensity varies or stays constant during sample rotation. Only for α = 0 • or θ = 0 • the SHG intensity remains constant for a full revolution of the rotation angle (ϕ), which is essential for a successful 3D reconstruction. Since the tilt angle θ is random for a biological sample, it can not be controlled and the laser polarization needs to be aligned in parallel to the axis of rotation (α = 0 • ). Only this way the SHG intensity can be kept constant during sample rotation and the inverse radon transformation can be applied successfully. Also equation (3) shows directly, the independence of ì E c on the rotation angle ϕ for α = 0, since all ϕ dependent components vanish to zero after matrix multiplication. However, depending on the tilt angle θ of the collagen fiber the amplitude of the SHG-signal differs. This means the tilt angle θ is an additional variable that influences the SHG intensity, next to phase matching and the number of scatterer. As a consequence no direct statement about the number of scatterer inside the sample is accessible by measuring the SHG intensity.
The effects of the intensity modulation during sample rotation are shown using the example of a collagen phantom. An area of circle was used as a phantom, that can be seen as a homogeneous cross section through a collagen fascicle (see Fig. 3 (A) left).
The x-axis of the reference coordinate system is directed into the image plane of the phantom. The orientation of the corresponding collagen filament and its axis of symmetry are specified by Fig. 1. To simulate the SHG-SLOT measurements, the Radon transform (RT) was applied on the phantom and the resulting sinogram was subsequently multiplied by the SHG intensity modulation (SHG IM) for the cases α θ = 0 • 45 • and α θ = 90 • 45 • (see Fig. 3(A), center). This SHG sinogram is subsequently reconstructed using the filtered back projection algorithm, which is the standard algorithm for tomographic reconstruction [35], and shown as the reconstructed phantom in Fig. 3(A) right. For α = 0 • no artifacts are visible in the reconstructed image. In case of α = 90 • , the intensity outside the circle is not flat but shows variations. To point out these differences intensity profiles along the indicated lines are shown in Fig. 3(B) (horizontal profile) and C (vertical profile). The profiles for α = 0 • show an equivalent shape with an increased amplitude as for the non-modulation case (compare green dotted and black solid line in Fig. 3). The profile for α = 90 • shows a broadening of the reconstructed circle in horizontal direction (see orange dashed line in Fig. 3 (A)) and a drop in intensity below zero in horizontal direction (see orange line in Fig. 3(B)). To visualize the extent of the reconstruction artifacts for many different combinations of the polarization angle α and the tilt angle θ, two artifact-values were defined. One is the geometrical distortion of the structure, which is measured by the averaged pixelwise amount differences between the reconstructed SHG-image (px SHG ) and the reconstructed image without amplitude modulation (px). This is reduced to the area outside the the circle: ∆ = |px SHG − px| . It is important to note, that this artifact-value strongly depends on the phantom size inside the image. As for all simulations the same input phantom is used, a relative comparison of this value is valid. The absolute value, however, is not significant. The second artifact value is the change in the maximum amplitude of the reconstructed circle with SHG modulation (A SHG ) and the image without amplitude modulation (A). This is expressed by the ratio of both: η = A SHG /A. 0.5 • steps and visualized in a false color image (see Fig. 4). In Fig. 4(A) the extent of the artifacts outside the circle (∆) is shown for all combinations of α and θ. As already seen in Fig. 2, there are no reconstruction artifacts for α = 0 • and θ = 0 • . All other combinations suffer from the intensity modulation and result in reconstruction artifacts outside the sample. The amplitude ratio η is only for some combinations of α and θ equal to one and can not be controlled by external parameters (see Fig. 4(B)).

Sample preparation
Tendon was isolated from rat tail by surgical instruments. The extracted tendon was subsequently rinsed with PBS (phospahate buffered saline) and incubated in 70% ethanol for 2 h and in 100% ethanol for 4 h to remove all water from the sample. Subsequently the tendon was stored in benzyl benzoate (n = 1.568) over night to achieve optical clearing. Afterwards the tendon was dissected in single fiber bundles, which are thinned to single fascicles to the end of the bundle [36]. The rat lymph node was cleared using CRISTAL (Curing Resin-Infiltrated Sample for Transparent Analysis with Light), whereby clearing is performed with monomeres, which are subsequently polymerized to obtain an optically cleared sample in a rigid polymer block [37]. The lymph node was incubated in 70% ethanol for 24 h and in 100% ethanol for 3 d. Afterwards the lymph node was stored in xylene over night and subsequently transferred into the monomer (NOA 68, Norland Products, USA) for 3 d (protected from light). To polymerize the monomer, the sample was exposed to UV light (375 nm) for 2 h resulting in a polymer with a final refractive index of Figure 3. The effects of the angular intensity modulation of the SHG signal on the reconstruction are shown. An area of circle is used as a phantom and shg reconstruction is compared to the non-modulated case. (A) The procedure of the simulation is shown. An area of circle is used as an input phantom, which can be seen as a cross section of a tendon. This phantom is radon transformed (RT), as it is the case during a SLOT measurement. In case of SHG-SLOT, the sinogram is affected by the intensity modulation (IM) that is due to the angel dependence of SHG. This modulation also depends on the input polarization α . Reconstructing this SHG sinogram results in a reconstruction with artifacts (for good visualization the contrast of the images is adjusted individually). The profiles of these images are shown as green dotted lines and orange dashed lines in B and C, respectively. The black solid line shows the profile for a reconstruction without SHG affected intensity modulation. (B) Profiles of the reconstructed images in horizontal direction are shown. For α = 0 • (green dotted line) the shape of the profile is equal to the non-modulated case (black solid line). Only the amplitude is increased. Also for α = 90 • (orange dashed line) the amplitude is slightly increased. Additionally the shape of the profile is affected by a broadening of the structure. (C) Profile of the reconstructed images in vertical direction. For α = 0 • (green dotted line) the shape of the profile is equal to the non-modulated case (balck solid line). Only the amplitude is increased. Also for α = 90 • (orange dashed line) the amplitude is slightly increased. Additionally the shape of the profile is affected by a drop to negative values right next to the structure. n = 1.54. All steps were performed at room temperature and during gentle agitation.

The SHG-SLOT imaging setup
For measuring an SHG signal efficiently, the SLOT setup needed to be modified in a way, that enables the detection of the generated signal in forward direction. The initial SLOT setup, which was limited to the detection of linear contrast mechanisms, has been described in detail elsewhere [24]. Briefly, the light of a laser diode was spatially cleaned by a single mode fiber and subsequently adjusted in diameter by a zoom lens (ZL) (see Fig. 5). The beam diameter determines the numerical aperture of the imaging setup and thus the optical resolution. For a successful reconstruction of the generated data, the Rayleigh range of the focused beam needs to be adjusted to the sample thickness. This simultaneously limits the resolution of the imaging setup. The beam is subsequently scanned by a x-/y-scanner (SM) over the sample while absorption and fluorescence/scattering can be detected by a photodiode (PD) and a photo multiplier tube (PMT), respectively. For this study several major additions to the imaging setup were necessary to enable SHG as a contrast mechanism. A fs-laser (Chameleon Ultra II, Coherent, USA) was integrated to the setup, which can be adjusted in power by a λ/2-plate (WP) and a polarizing beam splitter (PBS) (see Fig. 5). For this study a laser power of 250 mW at 800 nm on the sample was used, which results in good image quality while limiting the data acquisition to less than 6 h. To adjust the light polarization angle α on the sample, subsequently another λ/2-plate is inserted in the beam path. The beam diameter of the fs-laser can also be adjusted by a zoom lens (BE052-B,Thorlabs Inc, USA) and is subsequently overlayed by a dichroic miror with the path of the cw-laser. The detection path has been extended by a second PMT (R10699, Hamamatsu Photonics K.K., Hamamatsu City, Japan), that also detects in orthogonal direction. For SHG detection the third PMT in forward direction is essential. Behind the cuvette, a collection lens (#66-008, Edmund Optics, USA), two filters (F75-680 and F39-390, AHF, Germany) for additional removal of the excitation light and a focusing lens (LB1027-A, Thorlabs, USA) are placed to detect the forward scattered SHG signal with PMT3. Additionally, another dichroic mirror (F38-699, AHF, Germany) is placed in the signal path in forward direction, to enable the detection of the absorption signal by the photo diode.The following measurements were all performed in forward direction by PMT3.
The field of view was reduced to the thinned part of the fascicle or the lymph node respectively, which were connected to the stepping motor, which enables the rotation of the sample by the angle ϕ. The sample is held in a cuvette (c) filled with benzyl benzoate for the rat tail tendon and with silicon oil (n = 1.54) for the lymph node. The fluid bath for the sample in the rectangular cuvette prevents the deformation and deflection of the scanning laser beam due to refractive index changes. This preserves the parallel displacement of the laser beam inside the sample. Projection images were acquired for every rotational step of 0.45 • with a pixel resolution of 2.9 µm/px. The beam diameter was set to 5 mm for the rat tail tendon and to 1 mm for the lymph node to ensure that the length of the excitation point spread function (ω Z ) covers the full thickness of the sample (ω Z (rat tail tendon) = 133 µm and ω Z (lymph node) = 3.33 mm, calculations according to [38]). This results in an optical resolution of 4.3 µm for the rat tail tendon and 21.7 µm for the lymph node. Further image processing, as contrast adjustment and image centering, was performed in Fiji [39,40]. Reconstruction was performed by the filtered back projection algorithm using the IMOD tomographic reconstruction package [41].

Experimental validation
To verify the simulations above, measurements on a rat tail fascicle, were performed. The tendon was prepared and measured as described in section 4. The average SHG intensity was measured for every rotational step with an input polarization parallel (α = 0 • ; green) and orthogonal (α = 90 • ; orange) to the axis of rotation (see Fig. 6 (A)). In agreement with the simulations in section 3, there is a strong intensity modulation over the rotation angle ϕ for α = 90 • , whereas the intensity is almost constant for α = 0 • . The remaining intensity variation for α = 0 • , which results in an oval shape of the plotted data, is due to the different interaction lengths of laser light and sample. During rotation of a sample with elliptical cross section, the projection of the SHG signal occurs once along the long axis of the sample and once along the short axis. Due to the coherent nature of SHG, the intensity can add up coherently, which leads to a quadratic dependence of the signal strength on the sample thickness and number of scatterer N if phase matching is fulfilled (W SHG ∝ N 2 ) [29]. This seems to be valid here, too, even though the coherence length of SHG in benzyl alcohol amounts to only 2.4 µm, which is shorter than the thickness of the sample. Due to that, quasi phase matching is assumed to be valid here for this imaging geometry, as it was shown in periodically poled cristals and also in collagen albeit for a tightly focused beam [29,42,43]. By calculating the square root of the intensity value of each pixel in the acquired projection images, the ellipticity of the angle dependence in Fig. 6(A) (green) can be eliminated and a linear dependence on the number of scatterer is maintained ( √ W SHG ∝ N) as required by the inverse Radon transform (see Fig. 6(B)). Additionally the simulation is shown as a red dashed line for Figure 5. A schematic drawing of the SLOT setup is shown. The light of a laser diode (cw-Laser) is coupled into a single mode fiber (SMF) and subsequently collimated and adjusted in diameter by a zoom lens (ZL). The beam of a fs-pulsed laser (fs-Laser), also adjustable in diameter (ZL), is overlayed with the cw-beam path by a dichroic mirror (DM). By scanning (SM) the beam over the sample (S), which is placed inside a cuvette (C), the fluorescence/scattering (yellow) or the SHG signal (green) can be generated. The sample is connected to a stepping motor, which can rotate the sample by the rotation angle ϕ. Fluorescence/scattering or SHG can be detected (separated by a filter (F)) from both sides of the sample by PMT1 and PMT2 (PMT: photomultiplier tube) or in forward direction by PMT3, which was used in this study, exclusively. The laser light that transmits the sample is reflected by a dichroic mirror (DM) behind the sample (S) and detected by a photo diode (PD). both polarization angles and θ = 35 • . The tilt angle θ of the collagen strand was previously measured in the projection data. Experiments and simulations are in good agreement. The reconstructions of one single plane of the acquired data are shown in Fig. 6(C). The images show similar features as the simulations in Fig. 3. For the case of α = 0 • , there are no artifacts outside the sample. For α = 90 • , the sample is surrounded by a brighter and a darker area. To obtain a closer insight, the profiles in horizontal (see Fig. 6(D)) and vertical (see Fig. 6(E)) direction are plotted. The profiles for α = 90 • (orange dashed line) show the broadening and the drop in intensity as previously simulated (compare to Fig. 3). The profiles for α = 0 • (green dotted line) do not show this behavior, which is in agreement with the simulations. Also the reduced amplitude for the reconstruction of the data for α = 90 • (orange dashed line) is in accordance with the simulations.
Furthermore a rat lymph node was imaged using the SHG-SLOT setup with α = 0 • . The reconstructed 3D data is shown in Fig. 7. To compensate for the quadratic intensity dependence on the sample thickness, again the square root of the intensity value of each pixel in the acquired projection images was calculated, as discussed above for the rat tail fascicle. If this phase matching assumption is transferable to the lymph node needs to be investigated in more detail. A different correction procedure might be more accurate, however this is beyond the scope of this paper. Due to the large field of view of the imaging setup, the whole lymph node (length: 5.3 mm) can be imaged at once. This way it is possible to visualize the full distribution of the collagen inside the lymph node. Due to the dense occurrence of collagen in the capsule (Ca), the 3D shape of the lymph node is depictable. Furthermore the trabeculae into the center of the lymph node (*) are visible. The shell of vessels (V) (artery and vein) provide a strong SHG signal, appearing as a hollow structures in the cross sections (see Fig. 7(A) and (B)) and as two thick bundles in the maximum intensity projection (see Fig. 7(C)). Furthermore the isotropic resolution of SLOT provides good capability for 3D rendering of the reconstructed data (see Fig.  7(D)). The acquisition time amounts less than 6 h. Using a conventional multiphoton microscope, the acquisition time would be increased to multiple days, under preservation of full sampling and comparable resolution in axial direction. It must be noted, however, that a conventional multiphoton microscope has a higher resolution in lateral direction, compared to the SHG-SLOT data.

Conclusion
In this paper have proven for the first time in simulation and in experiment the capability of radon based imaging techniques, as OPT and SLOT, to visualize samples in 3D using SHG as a contrast mechanism. Due to the angular dependence of the SHG signal on sample and light polarization, the SHG intensity fluctuates during sample rotation. By choosing a light polarization parallel to the axis of rotation this fluctuation can be eliminated. Only the amplitude of the reconstructed structure is affected by the tilt angle of the sample itself, which prohibits the quantitative evaluation of the number of scatterer inside the sample. All geometrical distortions that are due to the angle dependence of the SHG signal on the reconstructed image are eliminated. This way, it is possible to reconstruct SHG data in 3D, that has been acquired using SLOT (or similar techniques as OPT). This enables now the fast visualization of samples in the mesoscopic range with SHG as a contrast mechanism, as shown here for a rat lymph node. More investigation needs to be addressed to the phase matching behavior in this imaging setup. To exclude, that unfulfilled phase matching generates additional artifacts inside the image, it is the next step to investigate this in more detail. This addresses specifically the correction procedure of the projection data that was realized here by calculating the square root of each pixel value in the projection data. This might not be sufficient for all kind of samples.
Furthermore, the SHG data can be directly correlated to other contrast mechanisms as absorption, scattering and fluorescence that are available in the very same imaging setup. This way structural information about the sample that is accessible via different contrast mechanisms (as specific antibody labeling via fluorescence) can be overlayed and correlated to each other.