Probing ordered lipid assemblies with polarized third-harmonic-generation microscopy

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés. Probing ordered lipid assemblies with polarized third-harmonic-generation microscopy Maxwell Zimmerley, Pierre Mahou, Delphine Débarre, Marie-Claire Schanne-Klein, Emmanuel Beaurepaire


I. INTRODUCTION
Ordered lipid assemblies are physiologically important structures participating in functions such as skin barrier and neuronal signaling.Although the supramolecular organization of lipid arrays is critical for proper function, it cannot be studied in situ and at the subcellular scale, however, using techniques such as x-ray scattering (XRS) and nuclear magnetic resonance (NMR).Recently, laserbased, coherent nonlinear optical microscopy (NLOM) has emerged as a possible approach for investigating molecular order in intact biological tissues.In addition, there has been an increased interest in characterizing the polarization-based optical response of ordered structures by using second-harmonic [1][2][3][4][5][6][7] and third-harmonic [8][9][10] generation (SHG and THG, respectively) and four-wave mixing [11] (FWM), which are all sensitive to the electronic properties of the material, as well as vibration-sensitive coherent-anti-Stokes-Raman-scattering (CARS) microscopy [12][13][14][15] and stimulated-Ramanscattering (SRS) [16] microscopy.
As in conventional polarized-light microscopy, polarization control in NLOM (P-NLOM) can be used to reveal material anisotropy and molecular ordering by determining the polarization dependence of the nonlinear optical (NLO) signal.This approach couples the biocompatibility, threedimensional (3D) signal localization, and specificity inherent in NLOM with the ability to probe molecular order, holding promises for advanced imaging studies of molecular arrays in tissues.In addition, while the theory describing multiscale, multiphoton light-matter interactions is complex, experimentally the typical NLOM setup need only be slightly modified by the simple incorporation of wave plates and polarizers.
In this study, we establish the sensitivity of P-THG to the orientation and ordering of lipids and thoroughly characterize this novel NLOM contrast for biomedical applications.This verification requires calculating the third-order nonlinear optical properties at each scale, followed by simulating the nonlinear light-matter interactions near the microscope focus.To that end, we first develop and validate a model for calculating the nonlinear optical response of individual lipid molecules and ordered lipid structures.We compare the calculated results to P-THG imaging experiments of multilamellar lipid vesicles (MLVs), establishing that polarization sensitivity correlates with molecular ordering, and validating the theoretical model.Finally, we present a biomedical application of P-THG toward monitoring lipid order in the outermost layer of human skin [called the stratum corneum (SC)] at the submicrometer scale of the intercorneocyte multilamellar matrix.

II. MODEL
On the molecular scale, the polarization-dependent, nonlinear optical response of third-order processes such as THG and CARS is described by the hyperpolarizability . is an 81-element tensor that maps three incident electric fields to the induced molecular dipole moment er ð3Þ according to [17] where e is the electron charge, i; j; k; l are any of the three Cartesian indices in the molecular frame, E loc ð! n Þ is the local electric field oscillating at frequency !n , and For the THG process, ! 1 ¼ ! 2 ¼ ! 3 ; for the CARS process, ! 2 < 0, and ! 1 þ ! 2 matches a vibrational eigenfrequency of the molecule.
In the approximately 1-m 3 excitation volume typical for NLOM, the detected signal reflects the averaged nonlinear optical response of the entire molecular ensemble located at the focus.At this scale, the induced nonlinear polarization P ð3Þ is where the nonlinear susceptibility tensor ð3Þ is dependent on the molecular geometries and at the focus, and I; J; K; L are any of the Cartesian indices in the laboratory frame.The symmetry properties of ð3Þ can be determined experimentally from the magnitude and directionality of the P-NLOM response, revealing information on the local molecular organization.For example, when tuned to match the CH 2 stretch frequency, P-CARS images of lipid bilayers show a maximum nonlinear response for excitation fields oriented parallel to the bilayer surface [18], that is, parallel to the CH 2 moiety.The CH 2 vibrational modes determine the symmetry properties of the molecular vibrational hyperpolarizability vib which, when averaged over the orientations found in a bilayer, produces an asymmetric vibrational ð3Þ .
Conversely, we note that the carbon-carbon (CC) bonds forming the lipid backbone show an average directionality parallel to the bilayer normal.According to Kajzar and Messier [19], the average electronic hyperpolarizability of CC bonds h CC i is approximately 6 times larger than that of CH bonds, hinting that this may be the source of a previously unnoticed asymmetric electronic ð3Þ that can be probed with P-THG.We later present experimental data on model lipid systems and human skin samples demonstrating that this is indeed the case.Before we do that, however, we derive the P-THG response from multilamellar lipid structures using a multiscale approach, summarized visually in Fig. 1.First, we calculate the electronic hyperpolarizability of a single lipid (1,2-Dimyristoyl-sn-glycero-3phosphorylcholine [DMPC]) molecule DMPC from an energy-optimized structure calculated by Krishnamurty et al. [20].Next, we use DMPC to determine the electronic nonlinear susceptibility tensors ð3ÞML of liquid-and gelphase multilamellar (ML) DMPC arrays.Then, with the determined values of ð3ÞML , we perform numerical calculations to predict the expected P-THG signal at lipid interfaces.

A. Calculating the molecular hyperpolarizability
We calculate the elements of DMPC using the bondadditivity model.We therefore assume that the nonlinear response of the molecule is a direct sum of the individual bond hyperpolarizabilities and that each bond hyperpolarizability is independent of the local electronic structure.To simplify the calculation, we relegate each bond type to having either a net isotropic (i) or polarization-sensitive (p) effect, and split up the hyperpolarizability [see Fig. 1(a)]: Examination of the molecular structure of DMPC allows us to identify the PO, CN, and CO bonds as those contributing to i .In the head group of the DMPC molecule, the four PO and four CN bonds are oriented in a tetrahedral geometry around the P and N atoms, respectively.For the CO bonds, the flexibility of the head group, the low cumulative directionality of the single CO bonds, and the CCto-double-CO-bond ratio of 29:2 all indicate that the anisotropic contributions of the CO bonds are negligible.Finally, while the number of CH bonds is approximately double the number of CC bonds, the orientation and low magnitude of h CH i compared to h CC i indicate that the CH bonds only increase the isotropy of the DMPC molecular hyperpolarizability, allowing us to relegate all non-CC bonds to i .We may thus calculate the quantity i as the direct sum of the various non-CC-bond hyperpolarizabilities h AB i: where N AB is the number of bonds of type AB in the molecule.For clarity, we present the average bond hyperpolarizabilities h ðABÞ i in Table I, and refer the interested reader to Appendix A for the derivation of the exact values used.We briefly note here, however, that we lack experimental values for PO corresponding to single and double bonds.We thus use the value for SiO single bonds, calculated from glass, and assume that the CC-to-PO bond ratio will minimize error (also potentially mitigating the error from assuming CH isotropy).Using these values for an isotropic DMPC liquid, we calculate ð3ÞDMPC ¼ 3:68 Â 10 À22 m 2 V À2 , i.e., a correct order of magnitude compared to the values 2.58 and 2.71 experimentally measured for triglycerides and vegetable oil, respectively [21].While the agreement on absolute average values is only semiquantitative, we show later that this model accurately predicts the P-THG response observed in imaging experiments of multilamellar lipids.
To evaluate p , we must first calculate the tensor elements of CC .For this, we use the bond-charge model [22][23][24][25] in which the bonding electrons are viewed as point charges centered between two positively charged nuclei.In this model, the symmetry of the CC bond results in nine nonzero tensor elements of CC , only two of which are unique: k and ? .For the CC bond, we calculate (in Appendix A) that k ¼ 0:69 and ?¼ 0:14 Â 10 À50 m 5 V À2 .W en o w calculate p as a sum of the projections of all tensor elements of each CC bond onto the molecular frame: where the first sum is over all CC bonds and the second sum is over the elements of CC .Figure 1(a) depicts the polarization dependence of the nonlinear dipole jer ð3Þ ð3!Þj calculated with Eq. (1) using either p (blue), i (green), or their sum, DMPC (red, with values given in Appendix A).It can be seen from the blue and red meshes that the symmetries of p result in lobed patterns in the induced nonlinear dipole which follow the orientation of the DMPC molecule.
We now focus on determining the averaged macroscopic, nonlinear optical susceptibility ð3Þ of DMPC arrays from DMPC .

B. Calculating the nonlinear susceptibility
As discussed previously, the nonlinear susceptibility of multilamellar lipid arrays ð3ÞML is related to the averaged molecular hyperpolarizability hi, which depends on the orientation of the lipid molecules in the assembly.We define the molecular z axis as the line from the midpoint of the two terminal aliphatic carbons to the midpoint of the two carbonyl carbons [Fig.1(a)] and define the laboratory X axis to be parallel to the surface normal of the multilamellar lipid structure [Fig.1(b)].Each lipid molecule in the array may freely rotate about both its own z (c ) and the laboratory X () axes.Neglecting birefringence, and assuming that the molecule adopts a Gaussian distribution of angles relative to the X axis [Fig.1(b)], we can calculate the elements of ð3ÞML as [17] ð3ÞML IJKL ðhi;Þ¼ where V is the average volume of DMPC [26]  The symmetry of stacked lipid bilayers (specifically, C 1V ) leaves 21 nonzero tensor elements of ð3ÞML , only three of which are independent: ð3Þ XXXX , ð3Þ XXYY , and ð3Þ YYYY .Importantly, the ð3Þ elements parallel and perpendicular to the lipid surface are linked to the orientation distribution parameters (hi, ) through Eq. ( 6); later, we determine how changes in these tensor elements affect the P-THG response from lipid-vacuum [27]and lipid-water interfaces [Fig.1(c)].We calculate ð3Þ elements for liquid-phase (hi¼0 ) and gel-phase (hi¼30 ) DMPC arrays spanning a range of disorder parameters ¼ 5:7 to ¼ 74 [Fig.2(a) and 2(b)].Of course, some of the states for which the calculation is done are not thermodynamically stable, but what we are particularly interested in here is the dependence of the optical modulation on the disorder parameter.In both systems, the tensor elements show a nonlinear relation with , which is decreasing for the element parallel to the bilayer normal XXXX and the cross term XXYY , and increasing for the perpendicular YYYY .Note that, with higher , the systems approach isotropy; in this case, the ð3Þ tensor elements have the relation ð3Þ XXYY .Comparison between the two systems shows that disorder in the gel phase induces an optically isotropic system for lower disorder parameters than the liquid phase; however, both systems approach the same average value of h ð3ÞDMPC i¼3:68 Â 10 À22 m 2 V À2 .
Using the calculated ð3Þ , we can visualize the nonlinear polarization jP ð3Þ j for ordered lipid systems [Fig.1(b)].The elongation along the X axis shows that the ordering information is retained even after averaging over the molecular configurations.While the displayed trace shows the behavior of a liquid-phase lipid system, a similar, slightly more spherical shape would be expected for the gel phase.

C. P-THG as a probe of molecular ordering
In a microscope, THG signal can be induced (i) within birefringent media [8,28] or (ii) at an interface between two materials (isotropic or ordered) with either different values for ð3Þ or refractive indices.the signal due to ð3Þ mismatches between the two materials is proportional to their difference squared [29]: where indices 1 and 2 refer to the two materials, and, for simplicity, this expression neglects refractive-index variations and dispersion.If one of these samples is ordered, the polarization of light can be varied to probe different ð3Þ tensor elements, e.g., an X-polarized field produces an X-polarized THG signal I X /ðÁ ð3Þ XXXX Þ 2 .Thus, a linearlypolarized electric field in the XY plane can be used to probe the three tensor elements of ð3ÞML according to FIG. 2. Calculated nonlinear optical properties of multilamellar lipid structures and interfaces.The ð3ÞML tensor elements of (a) liquid-phase ðhi¼0 Þ and (b) gel-phase ðhi¼30 Þ lipids as a function of disorder .Higher ordering results in greater ð3ÞML anisotropy for the liquid phase; in the gel phase, optical isotropy is achieved at lower disorder.The elements are also shown to approach the same value of h ð3ÞML i (dotted line) as the system approaches isotropy.(c),(d) The anisotropy coefficients X , Y characterize the magnitude of the anisotropy-induced P-THG signal oscillation.Generally, for a pair of coefficients, X > 0 and Y < 0, and their magnitude is larger for liquid-phase interfaces (gray curves) relative to the gel-phase interfaces (purple curves).For an interface between two isotropic media, X ¼ Y ¼ 0; thus, the gel-phase curves [corresponding to the tensor elements in (b)] cross at % 45 , and the liquid-phase curves [corresponding to (a)] cross at % 75 .Comparing the (c) vacuum-lipid and (d) water-lipid interfaces shows that ð3Þ matching of the second, isotropic material can potentially be used to enhance the P-THG contrast.
where is the direction of the field relative to the X axis, and ? ¼ðÁ ð3Þ YYYY Þ sinðÞ 2 þ 3ðÁ ð3Þ XXYY Þ cosðÞ 2 and k and ?are defined relative to the laboratory X axis.
Because of the relationship between tensor elements discussed above, the differences between the elements of two isotropic materials are related as Á XXXX ¼ Á YYYY ¼ 3Á XXYY .Therefore, deviations from this ratio can be used to characterize anisotropy: where X and Y are termed the interface ''anisotropy coefficients.''Importantly, X and Y depend on the ð3Þ of both materials.The disorder dependence of X and Y is plotted in Fig. 2(c) and 2(d) for lipid-vacuum and lipid-water interfaces.As can be seen, X > 0 and Y < 0 for most disorder parameters, but they approach zero as the system approaches isotropy.Comparing the two interface types shows that the anisotropy coefficients are larger for better ð3Þ -matched materials; we show later that this effect can greatly increase the sensitivity of P-THG to molecular ordering.The crossing of X and Y for larger disorder parameters reflects the change in the relative magnitude of the tensor components, indicating that disorder in the CC-bond orientation might reverse the P-THG contrast.
Rewriting Eq. ( 8) using the anisotropy coefficients, the THG signal can be expanded as where we see the nonlinear, oscillatory dependence of the P-THG intensity on the material anisotropy coefficients.This equation provides a good starting point to describe the extraction of ordering coefficients from the P-THG response of lipid interfaces; however, a more rigorous treatment is needed for an analysis of P-THG microscope images, as we explain in the next section.

D. P-THG microscopy modeling
In practice, THG microscopy usually involves focusing a beam into the sample of interest with a high-numericalaperture (NA) objective lens.The simple field description assumed in the previous section is therefore insufficient.Importantly, the tight focusing of a linearly polarized laser produces a transverse asymmetry of the excitation field.This asymmetry in the focal electric field results in an additional P-THG modulation, expected even for interfaces between two isotropic media [9,30].In simulations and imaging experiments of interfaces between isotropic lipids and water, we observe that this extra modulation is well fitted with a cosðÞ 2 term with a magnitude related to the NA.
To ascertain that the measured response retains ordering information, we have performed numerical calculations [9] of P-THG induced by a linear beam focused at 0.8 NA onto lipid-water and lipid-vacuum interfaces, varying .(For full details, see Appendix B.) In addition, from the P-THG angular profile, we define the modulation M ¼ ðI max À I min Þ=I max , which, as we will show, allows us to characterize the molecular ordering even in the presence of focus asymmetry.
We have calculated ð3ÞML ðhi;Þ and the modulation for both liquid-phase (hi¼0 )a n dg e l -p h a s e( hi¼30 ) DMPC arrays as a function of .Figures 3(a) and 3(b) depict the general trend predicted by the calculations for liquid-phase and gel-phase lipid arrays with lipid-water and lipid-vacuum interfaces: The optical modulation shows a nonlinear, monotonically decreasing relationship with the disorder parameter for all cases simulated.Comparing the lipid-vacuum interfaces (gray curves) and lipid-water interfaces (blue curves) shows the effect of ð3Þ matching on polarization contrast; as Á ð3Þ XXYY decreases in Eq. ( 9), the anisotropy terms X , Y and the relative modulation increase.Thus, while ð3Þ matching usually reduces contrast in THG imaging [21], it increases contrast in P-THG microscopy of anisotropic structures.We also note that, as the system approaches isotropy, the modulation is nonzero due to the focus asymmetry as described above; however, we can interpret modulation higher than approximately 0:15 as reflecting the molecular arrangement of the lipids within the bilayer.
Also of interest is the full form of the THG intensity profiles as a function of .The values obtained from the simulations are displayed as points in Figs.3(c) and 3(d), color-coded to match Fig. 2(a).We find that a fitting function of the form of Eq. ( 10) (where A and B are the fitting parameters) provides an accurate approximation of the angular dependence of the P-THG signal.[See the continuous curves in Figs.3(c) and 3(d).]Because of the relative magnitude of X , Y , 2 X , 2 Y , we find that a fit, [see curves in Figs.2(c) and 2(d)], provides an accurate approximation of the numerical simulation [points in Figs.2(c) and 2(d).We thus use this equation to fit experimental P-THG profiles from MLVs in the next section.
Because the NA-related modulation f NA is present for any interface, we may assume that it is proportional to the material signal (that is, f NA =ðÁ ð3Þ Þ 2 ¼ const for all materials).In addition, we make the assumption that the PROBING ORDERED LIPID ASSEMBLIES ... PHYS.REV.X 3, 011002 (2013) 011002-5 material and NA-related anisotropies are independent of one another.(See Appendix C for further discussion.)First, this allows us to rewrite Eq. ( 9) from the last section as Second, we may separate the material-anisotropy and focus-asymmetry effects in the constant Bðhi;;NAÞ of our fitting equation as follows: Bðhi;;NAÞ¼ Bðhi;Þþf NA .Using the definition of the modulation with both the fit equation (10) and Eq. ( 11) above, we can derive the relationship between X , Y and A, B as and the factor f NA can be determined from the modulation at the interface of two isotropic materials ( X , Y ¼ 0): Importantly, Eq. ( 10) can be used to fit P-THG microscopy profiles from lipid interfaces, and Eqs. ( 11)-( 14) can be used to quantify molecular ordering.

A. Microscopy
Multiphoton images of MLVs and skin are made with a custom-built, scanning, upright microscope (Fig. 4) utilizing a 20x 0.8 NA water-immersion objective (Olympus).Microscope resolution for THG imaging is estimated as the full width at half maximum of the THG profile from glasswater interfaces; we have determined the transverse (XY) and axial (Z) resolutions to be 0:5 m and 4:3 m, respectively.For MLVs, 8 images taken at a pixel rate of 400 kHz (pixel dwell time 2:5 s) are averaged for analysis.For the skin-biopsy samples, 4 images with 200 kHz pixel rate (pixel dwell time of 5 s) are averaged, taken approximately 30 m from the skin surface.THG imaging is performed using 1180 nm excitation pulses provided by a KTP-based OPO (APE, Germany) externally compressed to about 100 fs using dispersive prisms.CARS images are generated using longer pulses (150 fs pump, 400 fs Stokes) from a second Ti:S/OPO laser chain (Chameleon-OPO, Coherent, Inc.) with frequency difference tuned to the CH 2 stretch frequency (2840 cm À1 ).For P-THG and P-CARS experiments, a polarizer and achromatic near-infrared half wave plate are placed just before the objective for polarization control.Polarization profiles of MLVs and skin samples are generated from 31 images covering a range of 0 -180 .

B. Ordering in multilamellar lipid vesicles
THG and CH 2 -resonant CARS images of MLVs are shown in Fig. 5(a).As can be seen immediately, the two techniques provide complementary information: Not only are the polarization dependencies perpendicular to one another, but the THG contrast is highest at the water-MLV interfaces and weaker in the bulk, whereas the CARS signal is weaker on the surface and strongest in the bulk.This polarization complementarity may explain discrepancies between P-CARS imaging experiments of MLVs [12,13], which utilized the ratio of resonant-tononresonant CARS to quantify molecular order.As we have demonstrated here, DMPC is oriented parallel to the hydrocarbon chains (and thus perpendicular to vib-CH 2 ).We point out that nonresonant CARS and THG are both sensitive to the electronic nonlinear optical response; normalization to either of these signals would add modulation artefacts to any analysis.Conversely, we can predict the polarization dependence of THG from lipid interfaces.In addition, as we show later in skin samples, THG is a background-free and high-contrast imaging modality for interface detection, easily combined with SHG and/or 2PEF microscopy [31] (which is the gold standard in tissue fluorescence microscopy).
The MLVs generate THG by the two contrast mechanisms discussed earlier: (i) Away from the edge, phase matching is possible due to MLV birefringence, and signal is observed locally when the excitation field is oblique to each symmetry axis of the liquid crystal; (ii) At the MLV-buffer interface, the differences in ð3Þ and refractive indices both contribute to signal generation.Polarizationresolved-detection THG images [Fig.5(a) and movie S1 in Ref. [32]) further show that the interface signal retains the same polarization as the excitation, while birefringencerelated emission from the bulk is perpendicularly polarized, in agreement with Ref. [8].Here, we focus on contrast mechanism (ii), which is most often met in biological imaging, and analyze the P-THG response at the MLV-water interface.
P-THG intensity profiles from the water-MLV interface are plotted in Fig. 5(b), along with the intensity profile of a plant oil droplet, in which lipids are randomly oriented.The difference in modulation between signals from the oil droplet (0:13 AE 0:02) and the MLV (0:73 AE 0:075) confirms that the polarization response in THG is sensitive to lipid ordering.The three MLV profiles compare different sampling sizes in the same region of the upper MLV interface.The experimental fits show relatively little change with an increase in sampling.Indeed, for the experimental conditions used in this study, we predict only a small gain in precision when averaging over more pixels.Our analysis (see Appendix D for a full discussion) indicates that we extract A, B, and 0 with good precision even down to the limit of the excitation volume, with a detection of about 5000 photons per P-THG pixel (distributed over the 180 profile).Moreover, the agreement between the predicted modulation and the experimentally observed modulation validates our theoretical framework for constructing the ð3Þ of the lipid arrays.The nonuniform modulation measured along the surface of the MLV probably reflects ordering heterogeneity due to imperfect lipid packing.This interpretation is corroborated by our detecting heterogeneity within the MLV in the THG images.
By fitting the profiles with Eq. ( 10), we generate corresponding average interface angle 0 and modulation (M) images on a per-pixel basis [Figs.5(c) and 5(d)].Figure 5(c) is a hue-saturation-value rendering of the MLV in Fig. 5(a) where the hue is the extracted 0 , saturation is 1, and the value or brightness depends on the relative error of the fit, where we have limited the processing to a sample-dependent intensity threshold.As can be seen, the fit accurately determines the angle of the bilayer normal which points radially outward.This sensitivity can be used to detect lipid orientation with sub-m spatial resolution in complex samples, such as the lipid aggregate shown in Fig. 5(e) and movie S2 in Ref. [32].

C. Quantitative comparison with XRS
To compare P-THG to other methods for quantifying order, we correlate the MLV modulation amplitude, depicted in Fig. 5(d), to the disorder parameter .The modulation of the MLV signal in dilute buffer corresponds to % 20 AE 8 from the ordered lipid-water interface profile [Fig.3(b)].To relate this to x-ray scattering measurements, using our data, we calculate hP 2 i¼0:77 [33,34] assuming hi¼0 (see Appendix E), and compare the result to hP 2 i%0:4 determined by Pan et al. [35] for DMPC (which correlates with % 38 ).This discrepancy of hP 2 i is likely due to the overestimation of h ð3ÞDMPC i¼ 3:68 Â 10 À22 m 2 V À2 calculated using the bond-additivity model.As noted earlier, the approximately 2:6 Â 10 À22 value for h ð3ÞLipid i measured by our group in a previous [21] more closely matches that of water (1:83 Â 10 À22 ) which, as we have shown, results in a higher modulation.If we uniformly scale the tensor elements of ð3ÞML ðhi¼0; ¼ 38 Þ by a factor of 0.63, we can calculate the experimentally observed modulation of 0.73 (Appendix E).This calculation results in h ð3ÞDMPC i¼ 2:3 Â 10 À22 m 2 V À2 , in very good agreement with the expected value.

D. THG modulation in human skin stratum corneum
Ordered lipid assemblies play a key role in the hygroscopic nature of the skin SC [36], and their disruption is associated with different skin diseases [37].The SC, depicted in Fig. 6(a), can be viewed as a ''brick-and-mortar'' assembly [38] where dried, dead corneocytes form large, flat disks surrounded by stacked lipid bilayers (approximately 10 bilayers per intercorneocyte stack, with a total approximate thickness of 50 nm [39]) consisting of cholesterol, free fatty acids (FFAs), and ceramides [40].Figures 6(b)-6(e) show the overlay of a THG and a twophoton excited fluorescence (2PEF) image of the skin stratum corneum along a fold, with THG signals anticorrelated relative to the endogenous fluorescence observed from the corneocytes.These maps explicitly show the physical separation of the signals and localize the THG signal to the ordered-lipid-disordered-corneocyte interface.
Figures 6(f) and 6(g) show P-THG images of the skin, taken with two orthogonal polarizations of the excitation field, demonstrating the structure-dependent variation of the signal.From a series of P-THG images, we construct the directional hue-saturation-value map of the lipid orientation of the intercellular skin lipids, shown in Fig. 6(h).In Fig. 6(i), it is interesting to note that the layers closer to the surface have the lowest modulation, approximately 0.6 compared to 0.75 for the layers further from the surface.This observation may reflect variations in lipid ordering for different layers of the SC.Importantly, as is evident in movie S3 in Ref. [32] (recorded at a depth of 35 m), P-THG still detects molecular order with high specificity even within complex tissue environments.Indeed, strong modulation is observed only from the stacked lipids in the folds, whereas disordered structures produce only small, geometry-related modulation.
While we observe excellent qualitative agreement between the stratum corneum and the MLVs, quantitative correlations between modulation and lipid disorder require further characterization.The density, length, and angular distribution of CC bonds within SC bilayers is different than in the MLVs.In addition, the geometry of the lipids in the skin is better approximated by a thin slab within the focal volume, requiring a more accurate theoretical treatment.Still, because the THG polarization dependence in the SC is strong and specific, and because THG microscopy is background free, this method should be a sensitive probe to disruptions of the lipid layers integrity (e.g., from solvents or potentially from lipid-related diseases).
We finally note that scattering may alter polarizationbased optical measurements in environments that exhibit large-scale anisotropy such as tendon [6].However, such artifacts should generally not be present in standard tissues where scattering is not expected to cause major distortion in the focal-field distribution [41].Ratiometric P-THG should therefore be a valid method for probing lipid order in various tissue contexts.

IV. CONCLUSION
We have identified a novel mechanism of nonlinear optical contrast sensitive to molecular ordering in lipid assemblies.We have presented a systematic, multiscale method for calculating the electronic nonlinear susceptibility of lipid arrays.Furthermore, we have derived a Note how the folds that are oriented perpendicular to the excitation polarization appear brighter than the ones with parallel orientations.As with the MLV, the P-THG signal can be fitted to extract the stacked lipids orientation (h) and modulation (i), to gain insight into the lipid ordering within the interfaces.Even within the tissue, a strong modulation is observed exclusively from ordered phases and not from other structures.See also the movie S3 in Ref. [32].
simple expression relating ð3Þ tensor elements to the P-THG response, physically relevant for interfaces between isotropic materials and any multilamellar lipid structures (or other arrays with C 1V symmetry).
The general agreement between numerical calculations and P-THG imaging experiments of MLVs validates the use of P-THG as a method to probe lipid ordering and facilitates the extraction of molecular distribution parameters with high precision even at the diffraction limit.P-THG imaging of human skin biopsy samples has led to the discovery of a remarkably strong probe of biological lipid ordering, and it is easily coupled to other multiphoton imaging modalities such as SHG and 2PEF.Importantly, P-THG is an experimentally simple technique requiring a single excitation beam, and skin is optically accessible [42] and amenable to in-vivo microscopy, potentially leading to biomedical applications.

APPENDIX A: DETERMINING THE PARAMETERS USED IN THE MULTISCALE MODEL
It is difficult to find comprehensive multiscale data and calculations in the literature, so we present a general guideline for determining the nonlinear optical susceptibility of dynamic systems.The calculation of ð3ÞML is broken into parts, each pertaining to a progressively larger scale: First, calculation and determination of the bond hyperpolarizabilities AB on the A ˚ngstrom scale; next, calculation of the molecular hyperpolarizability DMPC from the molecular structure on the nm scale; and, finally, averaging over molecular orientations to calculate the nonlinear susceptibility of multilamellar lipid structures ð3ÞML on the 0:1-1-m scale.From experiment, we already have ð3Þvac ¼ 0 and

The CH, CN, CO, and PO hyperpolarizabilities
We start from the work of Kajzar and Messier [19] for the experimental values of AB that we use in this study.In their work, they measured the ð3Þ of various transparent liquids relative to the fused silica ð3Þglass , from which they calculated h mol i for each liquid.By systematically measuring similar structures, they were able to estimate AB for several bond types of interest to our study of DMPC; specifically, they presented values for CC , CH for single bonds and CO for both single and double bonds.
At this point, we lack values for CN and PO single and double bonds.We obtain an estimate for CN from Kajzar and Messier's experimentally determined molecular hyperpolarizability of N,N-dimethylformamide DMFA , composed of seven CH bonds, one CO double bond, and three CN bonds.Because we could not find a tabulated reference value for PO , we use the value for SiO , calculated as half the molecular hyperpolarizability of SiO 2 [19].While not ideal, we believe that this approximation is sufficient due to the ratio of PO to non-PO bonds in the DMPC molecule and the similarity between the electronic structures of Si and P.
While these numbers are a valuable starting point, more recent, similar THG experiments from our group have used a more accurate and significantly different value for ð3Þglass [43].Because all of the AB by Kajzar and Messier were calculated relative to ð3Þglass , we need to update the values for this study.Fortunately, since both groups measured the same ratio ð3ÞH 2 O = ð3Þglass , we may simply multiply each value of AB by the ratio of the measured ð3ÞH 2 O to calculate the updated values presented in the main text.

CC hyperpolarizability
We make use of the bond-charge model to determine the ratio k = ? of the elements of the CC-bond second hyperpolarizability.Originally, Van Vechten [22] derived simplified relations between the linear susceptibility ð1Þ and the band gap of a material E g .The energy E g has ionic C and covalent E h components, each of which is dependent on the position of the electrons relative to the nuclei and was measured for many semiconductors.Based on this work, Levine [23] later noted that an induced polarization changed the electron positions (Ár) and thus the bond energies.By expanding the expression for the electricfield-induced polarization to first order in Ár, Levine was able to calculate the first hyperpolarizability from atomic properties and crystal structure.
Following a similar path, Chemla et al. [25]d e r i v e d expressions for the second-order expansion of the energy and for calculating the two independent tensor elements of AB .Taking the z 0 axis as the bond axis, the nine nonzero elements of AB are related to these two values according to k ¼ z 0 z 0 z 0 z 0 ; ?¼ i 0 i 0 i 0 i 0 ¼ 3 i 0 i 0 j 0 j 0 ¼ 3 i 0 j 0 i 0 j 0 ¼ 3 i 0 j 0 j 0 i 0 ; i 0 ;j 0 ¼ x 0 ;y 0 ; i 0 Þ j 0 ; which shows that fields parallel (perpendicular) to the bond axis induce only parallel (perpendicular) dipoles er ð3Þ ð3!Þ, respectively.While the bond-charge model itself provides an order-of-magnitude estimate of bond hyperpolarizabilities, the ratio of i to p determines the observed polarization dependence of the THG signal.Therefore, we limit the use of the bond-charge model to MAXWELL ZIMMERLEY et al.
PHYS.REV.X 3, 011002 (2013) 011002-10 calculating the ratio k = ?and later couple this ratio to the experimentally measured h CC i.
For brevity, we present the calculation of k = ?specific to the CC bond.For bonds in a diamond-type crystal (i.e., no ionic components, C ¼ 0), where we explicitly use the notation set out by Van Vechten [22] for the homopolar energy E h and the notation by Chemla et al. [25] for the factors g 1 and h 2 .This model was developed using the Gaussian system of units; here, we perform the calculation in this system and convert to le Syste `me international d'unite ´s (SI) at the end of the calculation.For a CC bond: where a 0 is the Bohr radius; e ¼ 4:803 Â 10 À10 esu is the electron charge; r C and Z C are the covalent radius and number of valence electrons on carbon, respectively; the exponential constant s ¼ 2:48 has been determined experimentally [22]; b is related to the average coordination number of the two atoms (and is generally between 1 and 2); k s (k f ) is the Thomas-Fermi screening constant (wave number), respectively; N is the number of bonding electrons per diatomic volume for a compound with molar mass m and density ; and N A is Avogadro's number.
Using the values presented in Table II, we calculate the ratio k = ?¼ 5:08.
The relationship between the average, isotropic hyperpolarizability hi of a collection of randomly oriented molecules and the molecular hyperpolarizability tensor elements is We may apply this relation to a theoretical liquid made up of noninteracting single CC bonds to obtain where we use ?¼ xxxx ¼ yyyy ¼ 3 xxyy , and iizz ¼ 0. With k = ?¼ 5:08 and with h CC i¼0:21 Â 10 À50 m 5 V À2 [19], we obtain k ¼ 0:69, and ?¼ 0:14 Â 10 À50 m 5 V À2 .

DMPC hyperpolarizability
The determination of the tensor elements of the molecular hyperpolarizability is performed according to Eq. ( 5)in the main text.The tensor elements of DMPC calculated in this study are summarized in Table III.Since we assume only Kleinman symmetry, there are 81 nonzero elements, 15 of which are independent.The independent tensor elements are related by the number of times a particular Cartesian axis appears in the subscript.Thus, ð2xÞyz refers to all elements such as xxyz , xzyx , etc., which are all equal.The only unique value of note is zzzz , since an arbitrary rotation about the molecular z axis (defined in the main text) will change the elements with x and y indices, while reproducing the same macroscopic ð3ÞML .

Lipid nonlinear susceptibility
The tensor elements of ð3ÞML are determined by averaging DMPC over all assumed molecular orientations.This calculation is performed according to Eq. ( 6) in the main text, where the factors L ! , L 3! , and N are [17] and the factors T nm are obtained according to where m is in the reference frame of the molecule and n is in the laboratory frame.We ignore the effects of both birefringence and dispersion in this study and use n X ¼ n Y ¼ n Z and n !;3! ¼ 1:43.

APPENDIX B: DETAILS OF P-THG MICROSCOPY SIMULATION
Once the ð3Þ tensor elements for all materials have been calculated, the P-THG microscopy simulation involves (1) calculating the 3D electric field near the focus, (2) simulating the induced nonlinear polarization P ð3Þ within the focal region, and (3) propagating it to the far field and coherently summing all the elementary contributions in the detection plane, taking into account their relative phase.

Field at the focus
We use the fairly straightforward approach laid out by Leutenegger et al. [44] for fast calculation of the focal-field distribution.Briefly, we simulate a linearly polarized Gaussian beam of light at 1180 nm incident on a 0.8-NA water-immersion objective with a filling factor 10=6:5.The focal volume is sampled at a size spanning AE1:17 m in the (lateral) X and Y directions and AE3:94 m in the (axial, parallel to beam propagation) Z direction, discretized into a 201 Â 201 Â 200 grid.In our simulation, we assume a refractive index of 1.33 for the entire focal volume.Thus, we ignore the refractive index and birefringence effects related to THG signal generation and investigate solely the Á ð3Þ effects upon signal generation, as stated in the main text.

Sample geometry
For the geometry of the sample in focus, we create a grid the same size as the sampled focal volume and assign the localized ð3Þ for each 3D grid point.The geometry we choose emulates an idealized lipid-water or lipid-vacuum interface between two semi-infinite slabs.The interface is located at the X ¼ 0 plane where all grid points with X<0 described by ð3ÞML and X !0 by ð3ÞH 2 O or ð3Þvac .T o simulate different field-interface angles (), we can recalculate either the sample geometry and ð3ÞML or the focal-field distribution for each angle.Because of the simple geometry assumed, we find that recalculating the sample geometry is faster than recalculating the focal volume, so we choose the former.The nonlinear polarization P ð3Þ is then calculated at each grid point according to Eq. ( 2) in the main text; here, because we assume Kleinman-symmetry conditions are fulfilled, we use the d-matrix formalism as in Ref. [25].

Far-field THG signal
Plane waves oscillating at the TH wavelength (393.3 nm) are propagated to a far-field plane 10 cm from the focus and 25 Â 25 cm (41 Â 41 grid points) in size using Green's functions [9,30] with a refractive index of 1.36.The intensity over the whole far field is integrated to obtain the THG signal.The THG signal is calculated with steps of size Á ¼ 11:25 (=16) to determine the analytical form of the angular dependence shown in Figs.2(c) and 2(d) of the main text.

APPENDIX C: FOCUSING RELATED P-THG MODULATION
In the main text, we assume that the polarization dependence of the focus transverse asymmetry (described by the factor f NA ) can be separated from the material effects.Here, we justify that this assumption is reasonable for the DMPC arrays simulated here.To do so, we compare the modulation determined from the simulations FIG. 7. Calculation of the signal modulation using the fast, approximate method.(a) Modulation as a function of the disorder parameter, similar to the predicted trend from the full numerical simulations.(b) Relative error between the two predictions, showing that the fast calculation is consistently within 10% (although generally even lower) of the numerical simulations.presented in Fig. 3(b) (which we assume are the accurate or ''true'' modulations) with the modulation calculated according to Eq. ( 11) and using f NA =ð3Áð3Þ XXYY Þ 2 ¼ 0:15.Figure 7(a) shows the predicted microscope modulation without the use of a full simulation as a function of disorder for liquid and gel phase, lipid-vacuum and lipid-water interface types.The relative difference 2ðM Ap À M Sim Þ=ðM Ap þ M Sim Þ is plotted in Fig. 7(b).As is readily apparent, the anisotropy of the lipids prevents this approximation from perfectly predicting the P-THG modulation; however, the fast, approximate calculation is still within 10% of the full simulation for all values, and within 5% for the majority.
In addition, we determine an additive correction factor F NA to the approximated modulation of Eq. ( 10) that may be used to calculate the expected modulation in the microscope: (C1) Importantly, this calculation allows us to predict the expected modulation in the microscope without performing the full P-THG simulation.

APPENDIX D: FITTING ERROR AS A FUNCTION OF SIGNAL STRENGTH
Here, we discuss the precision with which we can determine optical parameters from the P-THG measurements, using an approach to error analysis similar to the one found in Ref. [45].The fitting function defined in the main text [Eq.(10)] was used because of the simple relationship between the nonlinear susceptibility elements and the fit parameters A, B, and 0 .Here, we use a different fitting function f for the P-THG signal, which can be related to the fit parameters discussed in the main text but which allows for a simpler derivation of measurement precision: where A ¼hIi is the signal average amplitude, À1 < M<1 is the modulation depth of the signal as a function of polarization, and 0 is identical to 0 defined in Eq. (10).The function is measured at N different points f i g i¼1...N , and the measurement yields the values fI i g i¼1...N .First, we note that the residual error for a least-square fit if there is zero noise on the signal and an error A, M, 0 on the parameters can be written as assuming that 0 ( , which is reasonable if we investigate the measurement error, and considering only the first order in infinitely small quantities A, M, 0 .I f we also assume that the N points f i g i¼1...N are evenly spaced between 0 and , the residual error reduces to Second, we now evaluate the residual error 0 purely due to shot noise in the signal.Here, we assume that the fitting parameters are perfectly known (A ¼ B ¼ 0 ¼ 0), so that the difference between the measured values fI i g i¼1...N and the fit values fIð i Þg i¼1...N is due only to the shot noise.Because this noise is on average equal to the square root of the intensity, the statistical mean of the total error can be expressed as where hIi is the average number of photons over the N measurements.
Last, the fitting parameter precision can finally be estimated by assuming that the residual error due to the fitting parameters is comparable to the residual error due to the shot noise, that is, ¼ 0 .For instance, for the statistical error in the fitting parameter A, we solve the equation [45] The quantities M m , 0m are the values that minimize in Eq. (D3); specifically, @M m =@A, @ 0m =@A are determined from the partial derivative of with respect to M and 0 , respectively.Here, we assume that there is no change in A or M with respect to changes in 0 (and vice versa).Following this scheme, we obtain the following values for the statistical error in the fitting parameters: As expected, the errors for M and 0 decrease as the square root of the signal (a typical behavior in the presence of shot noise) and with the modulation depth M. (The greater the modulation, the more information in the measured values.)For the curve presented in Fig. 3(b) of the main text, the numbers are as follows: typical maximum signal of approximately 1 photon per 10 pulses (i.e., twice below photon-counting saturation); pixel dwell time of 2:5 s; 8 images per orientation; and 31 orientations, resulting in a typical total signal of about 5000 photons per P-THG pixel.The modulation depth M as defined above is around 0.6, which yields

FIG. 1 .
FIG. 1. Visual summary of calculations.(a) In the bond-additivity model, each bond contributes mol .The non-CC bonds contribute isotropically while the CC-bond directions give rise to a directed ð3Þ .The nonlinear dipole moment jer ð3Þ ð3!Þj is calculated as a function of excitation angle [Eq.(1)] for p (blue mesh), i (green mesh), and DMPC (red mesh).The black lines are the same length and indicate the molecular z axis.(b) The lipids within the bilayer rotate freely about c , and adopt a Gaussian distribution of angles relative to the bilayer normal (laboratory X axis).Averaging of DMPC about these axes results in ð3ÞML .The grey mesh shows the lipid nonlinear polarization jP ð3Þ ð3!Þj, and the black line indicates the bilayer normal.(c) Left: Simulated focal volume half-filled with ordered lipid (dark gray) and either vacuum or water (clear).Right: Focal distribution of induced Pð3!Þ in each material at the focus (NA ¼ 0:8, ¼ 1180 nm, n !¼ n 3! ¼ 1:33, scale bar 500 nm).

FIG. 3 .
FIG. 3. Calculated P-THG response as a function of lipid order parameters.(a) The components of ð3ÞML are determined by the physical distribution of lipid angles relative to the X axis.A low angle and low correlate with the highest P-THG modulation.The arrow points from the highest-modulation parameters to lowest.(b) Calculated modulation amplitude as a function of lipid parameters and interface type.Modulation is generally higher at lipid-water (blue curves) than at vacuum-lipid (grey curves) interfaces, showing the increased sensitivity to anisotropy due to the partial ð3Þ matching.The modulation is also consistently larger for the liquid-phase lipids (closed symbols) relative to the gel-phase (open symbols) lipids for a particular interface type.Experimental values for the modulation in the MLV and lipid droplet (LD) are indicated with dashed lines.(c),(d) The calculated angular profiles for four multilamellar lipids, color-coded to match (a).The P-THG signal is lowest when the excitation field is parallel to the bilayers.The fitting curves are of the form I ¼ A þ B cosðÞ 2 .

FIG. 4 .
FIG.4.Setup used for the various imaging experiments.OBJ ¼ Objective, PMT ¼ photomultiplier tube.The excitation path is represented by the red beam; blue beams indicate both the forward and the epi (backward) emission pathways.THG and CARS images are usually detected in the forward direction, while the 2PEF images are epidetected using a dichroic mirror.The general form of this setup is typical for any nonlinear optical microscope, with the inclusion of a polarizer and a half wave plate to guarantee accurate control of the polarization state at the focus.

FIG. 5 .
FIG. 5. P-NLOM imaging of ordered and amorphous lipid structures.(a) THG and CARS images of MLVs generated with two orthogonal polarizations (scale bar 15 m).The maximum CARS signal is obtained with excitation fields polarized parallel to the CH 2 moieties, i.e., to the MLV surface.Two types of THG signals are detected.The interface THG signal (analyzed in this study) is maximized for an excitation field polarized parallel to the CC bonds (i.e., to the surface normal) and retains the excitation polarization.The bulk THG signal is due to MLV birefringence and is polarized perpendicular to the excitation.(See also movie S1 in Ref.[32].)(b) Comparison between the polarization dependence of THG from the surface of a lipid droplet (isotropic lipids) and a MLV (with approximately 10 3 bilayers), clearly showing the sensitivity of P-THG to molecular ordering.The MLV data analyses indicate that pixel averaging is not necessary to retrieve accurate estimates of the interface angle and modulation amplitude.(c)-(e) Extraction of order parameters.On a per-pixel basis, the fits can be used to construct maps of (c) the lipids orientation and (d) modulation amplitude based on the interface THG signal (scale bar 15 m).The same fitting procedure can accurately extract physical parameters in more complex systems, such as the large lipid aggregate in (e) (scale bar 15 m; see also movie S2 in Ref.[32]).

FIG. 6 .
FIG. 6. Nonlinear optical imaging of samples from human skin biopsies.(a) Within the skin stratum corneum (labeled ''SC''), stacked bilayers exist between the dead corneocytes (labeled ''c'').Vertically oriented folds extending tens of microns deep into the epidermis (labeled ''e'') provide an optimal geometry for P-THG imaging (dashed blue line).(b) THG images of the skin highlight cell nuclei (labeled ''n'') as well as the folds (''f'') (scale bar 50 m).(c)-(e) Subregion of the rectangular area indicated in (b).In human skin biopsies, endogenous 2PEF from the corneocytes alternating with THG signals (see white arrowheads) localizes the THG response to the intercellular lipid interfaces.(f),(g) THG images of a human skin biopsy approximately 25 m from the skin surface generated using two different polarizations, indicated by the straight doubleheaded arrows (scale bar 50 m).Note how the folds that are oriented perpendicular to the excitation polarization appear brighter than the ones with parallel orientations.As with the MLV, the P-THG signal can be fitted to extract the stacked lipids orientation (h) and modulation (i), to gain insight into the lipid ordering within the interfaces.Even within the tissue, a strong modulation is observed exclusively from ordered phases and not from other structures.See also the movie S3 in Ref.[32].

TABLE I .
Values used to compute DMPC , ð3ÞML .