1. Introduction

Nonlinear optical (NLO) phenomena arise from multiple interactions between optical fields and matter. For example, in second-harmonic generation (SHG), an incident field at the fundamental frequency ω interacts twice with matter to give rise to a field at the doubled frequency 2ω [1], and in third-harmonic generation (THG), three interactions give rise to a field at the tripled frequency 3ω [2,3]. These coherent processes conserve energy, and momentum considerations lead to phase-matching effects [1,2]. Such phenomena can provide new contrast mechanisms for three-dimensional (3D) label-free microscopy [4,5].

A majority of biological molecules are chiral, i.e., they lack mirror symmetry and occur in two enantiomers that are mirror images of each other. Unfortunately, chirality affects the linear optical properties of the molecules only weakly and the resulting optical activity effects, such as circular dichroism typically have a relative strength of only ~10⁻³ [6]. The relative weakness is due to the fact that such effects are forbidden in the electric-dipole approximation of the light-matter interaction [6]. NLO phenomena also provide new opportunities to study chiral molecules [7]. Nonlinear optical activity (NOA) effects can be allowed even in the electric-dipole approximation, and can have relative strengths approaching and even exceeding unity [7,8]. These properties have motivated the development of several NLO techniques to study chiral systems [9–13].

Optical processes can be described using susceptibility χ tensors, where the symmetry of the material structure dictates the structure of the tensor, i.e., the non-zero components, their interdependences, and their transformation properties [6]. For example, as a second-order process, SHG is forbidden in centrosymmetric media in the electric-dipole approximation, but can always occur at interfaces or surfaces [1]. The χ tensors can also be complex-valued, where the imaginary parts are related to absorptive and/or retarded responses [14,15]. For example, circular dichroism and second-harmonic generation circular difference (SHG-CD) responses from thin samples require that the related χ tensor is complex-valued [14,16], although the imaginary parts are often neglected in SHG microscopy studies. Considering only the real parts seems reasonable when the excitation and emission frequencies are far from any material resonances, since the χ tensor should be approximately real-valued as expected for an energy conserving process. But such approximation does not always seem to be valid, as recent experimental observations of SHG-CD responses from collagenous tissues have shown [8,17].

To understand the origins of the measured NLO signals in biological tissues, it is therefore important to develop new modeling tools. A majority of formulated approaches are based on a simplified version of the discrete dipole approximation [18–20], where interparticle coupling effects are neglected, hence the term uncoupled dipole model (UDM) [21–24]. Recently, this approach has also been generalized to model NLO microscopy [25,26]. Also, a full NLO discrete dipole approximation formulation has been developed [27]. But as biological tissues scatter light more weakly than for example metallic nanoparticles, the computationally more efficient UDM provides a good starting point for studying how NLO responses could arise from complex 3D biological materials.

Here, we use UDM to numerically investigate the origin of NOA signals from macroscopically organized 3D biological materials. We perform numerical SHG and THG microscopy on structures resembling collagen fibers and multilamellar vesicles (MLV), which were chosen as model systems due to their biomedical relevance [28–31]. The numerical SHG results elucidate the recently measured strong SHG-CD responses from type-I collagen fibers [8], where the 3D macroscopic organization of the fibers together with the complex-valued response tensors give rise to the strong NOA signals. Complementarily, the numerical THG results from MLVs predict the potential of a novel contrast mechanism in imaging using third-harmonic generation circular difference (THG-CD) responses, which have been recently measured from lipid droplets in mammalian cells [32].

2. Theory

2.1. Uncoupled oscillator model and Green’s function approach

The nonlinear process of Nth order can be described using the respective susceptibility χ^(Ν)(ω = ω₁ + ω₂ + ... + ω_Ν, r) tensor consisting of 3^N+1 components. The structure of the tensor is dictated by the symmetry of the medium and the particular NLO process considered. Note that we allow the tensor to vary spatially to account for variations in molecular orientation within the sample. Here, we concentrate on SHG and THG processes using a single excitation beam at the fundamental frequency ω, for which the susceptibilities are χ⁽²⁾(2ω,r) and χ⁽³⁾(3ω,r), respectively. More detailed descriptions of the susceptibilities follow in later sections.

Unlike in the full discrete dipole approximation approach, the scattering effects at the fundamental and NLO frequencies are neglected in UDM [27], which means that the total and incident fields are equal E(ω,r) = E₀(ω,r). This condition is also known as Rayleigh-Gans approximation (RGA) [33,34], and can be justified by assuming that the refractive indices of all media are close to each other. RGA provides a simple framework to understand the underlying physics. For simplicity, we take all refractive indices equal to unity, and neglect the effects of optical nonlinearities to the refractive indices. We use RGA for simplifying the scattering problems occurring both at the fundamental frequency as well as at the harmonic frequencies Nω [25,27]. Other more rigorous approaches to consider the scattering and birefringence effects have been presented elsewhere [26,34,35].

When RGA and electric-dipole approximations are valid [1], the incident field E(ω,r) gives rise to NLO polarization distributions as

Due to the scaling of the NLO response with high powers of the incident field, the responses are appreciable only close to the focus. Therefore, the source volume V can be restricted to the vicinity of the focal plane of the focusing objective. In addition, the vectorial aspects of light become more important when light is strongly focused [36]. Recently, we have demonstrated that SHG emissions from isotropic surfaces with tightly focused beams can considerably deviate (~5-fold difference) from the predictions of scalar approaches, which was interpreted in terms of effective multipoles due to the vector excitation [37]. Therefore, the excitation E(ω,r) is considered fully vectorial and is calculated using the angular spectrum representation [16,21,36].

In the calculations, the following parameters were used: fundamental wavelength of 1 μm, focal length f = 2.625 mm, beam waist w₀ = 4 mm (at the back aperture of the lens), and apodization factor f₀ = 1.9. The corresponding NA was 0.8, and similar NA was assumed for the collection lens. Here, only the harmonic emission in the transmitted direction was considered. Note that in the NLO microscopy simulations, E(ω,r) for the given excitation polarizations needs to be calculated only once and therefore does not limit the overall calculation speed. Now we have all we need to numerically calculate how the NLO signal originates from the sample volume V, described by some inhomogeneous χ^(Ν)(Νω,r) distribution. We wrote the codes to simulate scanning NLO microscopy of complex 3D samples using MATLAB, and studied how NOA effects could arise from the samples of our two example cases resembling collagen fibers and MLVs.

2.2. SHG and THG susceptibilities for collagen fibers and MLVs

The tropocollagen molecules consisting of three polypeptide helices are usually considered to belong to C_∞v, C_∞, or C₃ symmetry groups [8,26]. The tropocollagen molecules are approximately 300 nm long and a few nanometers in diameter, which in collagenous tissues are often further organized into fibrils with ~200 nm of diameter or fibers with a few micrometer diameter [26]. Therefore, for our discretization accuracy, the SHG-active collagen structures correspond to a portion of a collagen fibril, which we can assume to have cylindrical symmetry and belong to either C_∞v or C_∞ symmetry group, which is higher symmetry than the C₃ used sometimes for the individual tropocollagen molecules.

For third-order processes, the symmetry properties of the sample do not have such a special role as for second-order processes, and for example THG is allowed also from centrosymmetric media [1]. The non-zero susceptibility components for THG from isotropic media are ${\chi }_{ijkl}^{(3)}={\chi }_{xxyy}^{(3)}({\delta }_{ij}{\delta }_{kl}+{\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk})$, where is the Kronecker delta.

3. Numerical results for polarized NLO microscopy

3.1. SHG-CD microscopy of collagen fibers

We modeled collagen fibrils or fibers as simple cylinders as seen in Fig. 1(a).The single fiber had a radius of ~200 nm and was tilted from the focal plane since the 3D orientation of the fiber affects the SHG-CD signals [8]. The side view of the sample is shown in Fig. 1(b) where the tilt angle was arbitrarily chosen to be 11.5 degrees. Scanning SHG microscopy simulations using left-handed circular polarization (LHCP) and right-handed circular polarization (RHCP) inputs for three fibers with different second-order susceptibilities were performed [see Fig. 1(c)-1(h)]. The results for the first fiber with susceptibility belonging to C_∞v are shown in Fig. 1(c) and 1(f). The non-zero${\chi }_{ijk}^{(2)}$components for the first fiber were taken as z’z’z’ = 1, z’x’x’ = z’y’y’ = 0.63 and x’x’z’ = y’y’z’ = x’z’x’ = y’z’y’ = 0.48, where z’ is along the long axis of the fibril [8]. The second fiber was assumed to belong to C_∞, and had an otherwise identical susceptibility than the first fiber, except it had also an additional non-zero chiral components x’y’z’ = x’z’y’ = -y’x’z’ = -y’z’x’ = −0.14. The results for the second fiber are shown in Fig. 1(d) and 1(g). The third fiber was otherwise identical to the second fiber, except the chiral components were complex-valued x’y’z’ = −0.14-0.19i. These values were used as they were previously found to provide qualitative agreement between SHG-CD simulations and experiments [8]. The results for the third fiber are shown in Fig. 1(e) and 1(h). Note that the SHG signals are maximized (minimized) at the center (edges) of the images, which is due to the oblique orientation of the fiber with respect to the focal plane. As the excitation field is scanned in the focal plane (where z = constant), it is no longer totally focused on the fiber at the edges of the images.

 

Fig. 1 Numerical SHG microscopy of a collagen fiber. Schematics showing the relative size and position of the modeled collagen fibril as viewed from the top (xy plane) (a) and side (xz plane) (b). The molecular coordinate frame (x’,y’,z’) is shown with dotted arrows. SHG images using (c-e) LHCP and (f-h) RHCP input polarizations were used to calculate the corresponding SHG-CD images for fibers with: (i) achiral and real-valued, (j) chiral and real-valued and (k) chiral and complex-valued${\chi }_{ijk}^{(2)}$susceptibilities. The normalized maximum SHG values are shown at the upper right of the SHG images and the maximum and minimum SHG-CD values are shown at the upper right of the SHG-CD images. Video showing SHG and SHG-CD images while varying the fiber orientation in the focus for the complex-valued case is given in (Media 1). Scale bars are 800 nm.

Download Full Size | PDF

We then calculated the SHG-CD signals pixel-by-pixel as $SHG-CD=2(I_{LHCP}-I_{RHCP})/(I_{LHCP}+I_{RHCP})$and plotted the resulting image [see Fig. 1(i), 1(j) and 1(k)]. In all the results, pixels with normalized SHG values smaller than 0.1 were excluded from the analysis to remove possible image artifacts due to the side lobes of the excitation field, which would not be present in most of the experiments due to experimental noise.

We then studied how the 3D orientation and position of the fibers affect their SHG and SHG-CD responses. A series of simulated SHG and SHG-CD scans while changing the relative orientation of the fiber at the focal volume with the complex-valued susceptibility is shown in Media 1. The results for fibers with real- and complex-valued and chiral susceptibilities while changing the relative z-position of the fiber with respect to the focal plane (z = 0 μm) are shown in Fig. 2.

 

Fig. 2 Numerical SHG-CD microscopy of a collagen fiber as a function of the relative z-position of the fiber with respect to the focal plane (dotted line). (a) Side views show the z-position of the fiber, where z = 0 μm corresponds to fiber being at the focal plane. SHG-CD images using (b) real-valued and (c) complex-valued chiral susceptibilities. The maximum SHG-CD values are shown at the upper right of the images. Complex-valued susceptibility leads to larger SHG-CD values, which are dominantly positive. Scale bars are 800 nm.

Download Full Size | PDF

3.2. THG-CD microscopy of MLVs

We modeled a MLV as a sphere consisting of a core filled with isotropic molecules and a thin outer layer with denser concentration of isotropic molecules as seen in Fig. 3. As previously shown, THG can occur in nonlinear medium either due to inhomogeneities in the refractive index or in the nonlinear susceptibility [2]. In our approach, we neglect the effects of changing refractive indices, but note that they might play a role in experiments. The 21 non-zero susceptibility components were taken as i’j’k’l’ = x’x’y’y’(δ_i’j’δ_k'l’ + δ_i’k’δ_j’l’ + δ_i’l’δ_j’k’), where x’x’y’y’ = 0.3. This assumption is based on the relatively well-known structure of MLVs exhibiting several lipid crystalline bilayers that encloses a core filled with aqueous solutions, surfactants, biomolecules, polymers, drugs or nanoparticles [31]. Here, the relative concentration c of the outer layer was taken to be 10-fold as compared to the inner isotropic part of the MLV, but even the outer layer was modeled as isotropic. This model is clearly highly simplified in the sense that we neglect the dependence of the refractive index on the density of molecules and the fact that the outer layer is usually oriented. However, the point here is that even the simplified model will lead to a rich variety of interesting unexpected effects without possible signal artifacts due to the ordered outer layer. The radii of the outer and inner parts of the MLV were taken as 0.5 μm and 0.47 μm, respectively [see Fig. 3(a)]. We then performed numerical THG microscopy of the MLV (see Fig. 3), while varying the relative z-position of the MLV (see Media 2). We used RHCP, LHCP and three linear polarizations for the excitation. The three linear polarizations were oriented 0°, 45° and 90° from the x-axis.

 

Fig. 3 Numerical THG microscopy of a MLV with a homogenous core. Schematics showing the relative size and position of the modeled MLV as seen from the top (xy plane) (a) and side (xz plane) (e). Scanning was performed in the xy-plane. THG images taken with (b) LHCP and (c) RHCP inputs were used to calculate the (d) THG-CD image. (f-h) Microscopy images taken using linear input polarizations (0°, 45° and 90° with respect to the y-axis) was also modeled predicting a ~10 000-fold THG signal enhancement compared to THG using CP inputs. The normalized maximum THG values are shown at the upper right of the THG images and the maximum THG-CD value is shown at the upper right of the THG-CD image. Video showing THG images while varying the MLV z-position in the focus is given in (Media 2). Scale bars are 500 nm.

Download Full Size | PDF

Then we introduced a smaller but ordered spherical volume, with a radius of 100 nm, to the center of the otherwise isotropic MLV core [see Fig. 4(a)].Such smaller spherical volume can be linked with the inclusion of closely-packed molecules inside a MLV. The ordered spherical volume was associated with THG-active structures that were assumed to be anisotropic and possess C_∞v symmetry. The 21 non-zero susceptibility components, of which only three are independent, were taken as z’z’z’z’ = 1, y’y’z’z’ = x’x’z’z’ = 0.1 and x’x’y’y’ = y’y’x’x’ = 0.05. The orientation of the ordered structures was arbitrarily taken to point at 68° from the z-axis and 24° from the x-axis towards y-axis. The concentration of the outer MLV layer was chosen to be 10-fold as compared to the inner isotropic part of the MLV. Numerical THG microscopy of this MLV containing a partly ordered core was then performed, while changing the relative concentration c2 of the ordered part with respect to the isotropic inner MLV part (see Fig. 4 and Media 3).

 

Fig. 4 Numerical THG microscopy a MLV with a partially ordered core. Schematics showing the relative size and position of the MLV as seen from the top (xy plane) (a) and side (xz plane) (e). Scanning was performed in xy-plane. THG images taken with (b) LHCP and (c) RHCP inputs were used to calculate the (d) THG-CD image, which shows a clear THG-CD response. (f-h) Microscopy images taken with linear input polarizations (0°, 45° and 90° with respect to the y-axis) was also modeled predicting decreased THG signal enhancement compared to THG using CP inputs (now only ~1800-fold). The normalized maximum THG values are shown at the upper right of the THG images and the maximum THG-CD value is shown at the upper right of the THG-CD image. Video showing images while varying the relative concentration c2 of the ordered MLV center is given in (Media 3). Scale bars are 500 nm.

Download Full Size | PDF

4. Discussions

First, we consider the computational aspects of the approach. We discretized the integral of Eq. (3) over volume V into 101 × 101 × 101 cubical sub-units V_i and denote the discretized NLO polarization inside V_i as $P^{(N)}\left(N\omega ,r_i\right)=P_i^{(N)}$. Changes in the local concentrations of the samples were implemented by multiplying the local susceptibilities, which is justified in the UDM approach where local-field effects are neglected. The edge lengths of the cubical sub-units were ~40 nm for SHG and ~20 nm for THG simulations. These discretizations were smaller than the traditional rule of thumb λ/10 for discrete dipole approximation. The integration of Eq. (4) over area A was numerically calculated by discretizing A into M × M sub-units A_j. The convolutional integration of Eq. (3) was sped up using fast convolution algorithms based on either fast Fourier transform or singular value decomposition depending on which one was faster for given input parameters [19,20,38]. For singular value decomposition approach we divided the 3D G(Nω,R_j-r_i) into 101 (M × M) stacks. These 2D stacks were then separable matrices and singular value decomposition could be used to divide the 2D convolution into two 1D convolutions [38].

The calculation times depended on how sparse the P_i array was, of the value of M and of the simulated NLO process. For SHG calculations (with sparser P_i) the singular value decomposition approach was faster for all the considered M values of 26, 52 and 101. The fast Fourier transform approach was faster for THG calculations for M = 101. Already with M = 26 the related discretization error was estimated to be smaller than 5%, and was seen to cause the slight discrepancy in Figs. 3(f-h), where physically the three different linear input polarizations should lead to identical THG images. But as the error was already quite small, discretization with M = 26 was deemed to provide sufficiently accurate results for the present study.

The images shown in Figs. 1 and 2 illustrate how a 3D macroscopic object can give rise to non-zero SHG-CD responses, which is perhaps surprising since the previous theoretical NOA studies have concluded that complex-valued response functions are needed for NOA effects [14–16]. But unlike this work, the previous works have considered thin films where field retardation effects present in thick media could be neglected. For fibers parallel to the focal plane, the maximum SHG-CD values were very small (~0.001). But when the fibers were slightly tilted, large non-zero SHG-CD responses (~0.25) were predicted also for fibers consisting of either achiral or chiral fibril structures with real-valued susceptibilities [see Fig. 1(i) and 1(j)]. This implies that the retardation effects can play a role in the SHG-CD responses, which can be qualitatively understood by considering the light-matter interactions occurring at the focal volume, where the overall symmetry of the interaction is due to both the excitation and the sample. Therefore, when an achiral but anisotropic sample is moved away from the center of the excitation, the overall symmetry of the macroscopic environment decreases and NOA effects, such as SHG-CD, become allowed [Fig. 1(i)]. Interestingly, for the real-valued susceptibilities we found zero SHG-CD responses when the fiber was centered at the beam focus, and non-zero SHG-CD when the sample was off-centered. In addition, the SHG-CD responses were seen to be anti-symmetric with respect to the long axis (x-axis) of the fiber for the achiral susceptibility as is seen in Fig. 1(i). The SHG-CD responses changed sign upon reversing the tilt of the fiber. Therefore, our results imply that non-zero SHG-CD responses could occur also purely due to the macroscopic geometry of the 3D sample. The overall symmetry of the interaction is further broken when the sample is chiral also at the level of the fibrils, in which case the SHG-CD responses are no longer perfectly anti-symmetric [Fig. 1(j)]. This is seen as a small positive offset of 0.035 in the SHG-CD values.

The strongest SHG-CD responses (~0.65) were predicted for the fiber composed of chiral fibrils with complex-valued susceptibilities [see Fig. 1(k)]. Interestingly, for the slightly tilted fiber and complex-valued susceptibilities, a strong positive offset of ~0.3 was predicted due to which the SHG-CD responses of the whole fiber are positive. This agrees qualitatively very well with the recent SHG-CD microscopy measurements of type-I collagen fibers [8]. Therefore, despite the possibility for non-zero SHG-CD responses due to field retardation effects even with real-valued susceptibilities, we believe that the imaginary parts of the susceptibility tensors could play a surprisingly large role in the SHG-CD responses of collagens. We also believe that it is very important to develop NLO characterization techniques to provide such complex-valued information of the material responses with sufficiently high spatial resolution, and in the future perform spectral SHG-CD experiments to provide further evidence of the complex-valued response functions.

Similar to the numerical SHG-CD results of collagen fibers, we found that non-zero THG-CD signals can also arise from 3D isotropic MLVs excited by oppositely-handed circular polarizations (Fig. 3). We note that similar THG-CD results have been recently measured from lipid droplets in mammalian cells [32]. We found that the THG-CD effects near the outer boundary of the MLV are relatively weak attaining a maximum (~0.2) when the center of the MLV is located very near the focal plane. Interestingly, we found that THG signals using circularly polarized inputs become stronger for MLVs that contains ordered domains [Fig. 4(b) and 4(c)], when compared to the THG signals using linear input polarizations [Fig. 4(f)-(h)]. Unlike for the linear input polarizations, the THG peak signals with circular input polarizations also become more localized to the vicinity of the ordered core. Although small anisotropy in the THG signals using linear polarizations is visible, the effects are small. In addition, THG signals with linear input polarizations showed visible polarization dependences even with rotationally symmetric isotropic samples [see Fig. 3(f)-(h)], which is attributed to the non-zero longitudinal field component present at the focus. This suggests that dependence of THG on linear input polarizations is not necessarily related to ordering [32]. Therefore, circular input polarizations seem to provide better contrast towards ordered structures. This is strengthened by the fact that THG with linear polarization is allowed from the bulk of an isotropic medium near interfaces [2], leading to a strong background that overwhelms the signal from the ordered core. THG with circular polarizations, on the other hand, is forbidden in the bulk of isotropic materials, providing thus superior sensitivity to any deviations from isotropy. Of course, these rules of thumb are modified in the present case by the vectorial character of the focal fields but the sensitivity of circular polarizations to ordering is still maintained.

In order to provide more quantitative measures towards ordering inside MLVs, we calculated also the THG-CD responses for a MLV containing ordered domains in the core [Fig. 4(d)]. The THG-CD responses are in general stronger compared to the isotropic MLV with an isotropic core, and more importantly the maximum THG-CD values are correlated with the relative concentration of the ordered domain inclusion. This is seen in Media 3, which shows THG microscopy results while varying the relative concentration ratio c2 of the ordered core and the surrounding isotropic constituents, reaching maximal THG-CD signals of ~1.1 when the concentration ratio c2 is 5-fold.

Similarly promising measure of the ordering would be to calculate the averaged ratio between the THG signals using linear and circular input polarizations $\delta =(I_0+I_{90})/(I_{LHCP}+I_{RHCP})$. As the THG signals are normalized to the maximum THG signal of the circular input polarizations, the maximum δ values are approximately the maximum THG values of linear input polarizations shown in the upper right corners of Figs. 3(f) and 4(f). Similar to THG-CD responses, also the δ values correlate with the concentration ratio of the ordered domain c2. For purely isotropic MLVs, the δ ≈10 400 and decreases down to δ ≈340 when c2 is 5-fold (see Media 3).

In our case, the ordered domain was assumed to possess quite high symmetry of C_∞v and we show that visible changes in the THG-CD responses could occur due to symmetry breaking in the focal volume caused by the domain. Similar results could be expected also for ordered domains with lower overall symmetries. Therefore, we propose that both the THG-CD responses and the δ values could be used to provide more quantitative imaging contrast mechanisms to detect ordered domains in 3D biological systems than what is possible using only linearly polarized input polarizations. We believe that these proposed contrast mechanisms could provide new easy-to-use and clinically relevant tools for lipid vesicle studies [31,32].

Finally, we checked the importance of using full vector approach for our results by repeating some of our simulations while forcing the longitudinal field component of the excitation field to be identically zero. We found out a general trend of roughly ten-fold (five-fold) decrease in the simulated SHG-CD (THG-CD) responses. This result agrees well with our previous computational work on NOA effects from surfaces [16]. In addition, we have earlier measured approximately five-fold differences between transmitted and reflected SHG emission from isotropic thin film sample, where traditional scalar approach does not predict any difference in emission strength [37]. Therefore, we believe our approach of using full vector excitation instead of paraxial approximation is both well justified and relevant.

5. Conclusion

We have used the uncoupled-dipole method to numerically study how strong nonlinear optical activity (NOA) effects can arise from three-dimensional morphology of biological materials. We have modeled second-harmonic generation (SHG) and third-harmonic generation (THG) microscopy of three-dimensional samples resembling collagen fibers and multilamellar vesicles. We studied the possible molecular and morphological origins of NOA effects, and propose that NOA effects could be utilized to provide morphological contrast for various nonlinear microscopy modalities. Our SHG microscopy results from samples resembling collagen fibers help elucidating the recently measured strong SHG circular dichroism signals [8]. In particular, our results show that the complex-valued material response tensor seems to play a vital role in reported SHG responses. We also show that non-negligible THG circular dichroism effects could occur in biological materials, such as multilamellar vesicles. More importantly, we propose that such THG circular dichroism effects could provide information of local ordering in otherwise homogeneous materials, and could therefore be an interesting new nonlinear microscopy modality.