# Multiscale characterization of myelin distribution with polarized THG microscopy

Josephine Morizet, Nicolas Olivier, Marie-Stephane Aigrot, Pierre Mahou, Elodie Martin, Anne Desmazieres, Bruno Stankoff, Chiara Stringari, Emmanuel Beaurepaire

► To cite this version:

Josephine Morizet, Nicolas Olivier, Marie-Stephane Aigrot, Pierre Mahou, Elodie Martin, et al. Multiscale characterization of myelin distribution with polarized THG microscopy. Optica, 2025, 12 (7), pp.1122. 10.1364/optica.562091. hal-05166303

HAL Id: hal-05166303

https://hal.science/hal-05166303v1

Submitted on 16 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Multiscale characterization of myelin distribution with polarized THG microscopy

JOSEPHINE MORIZET,<fnref n="1" />,<fnref n="2" /> NICOLAS OLIVIER,<fnref n="1" /> ID MARIE-STEPHANE AIGROT,<fnref n="3" /> PIERRE MAHOU,<fnref n="1" /> ID ELODIE MARTIN,<fnref n="3" /> ANNE DESMAZIERES,<fnref n="3" /> BRUNO STANKOFF,<fnref n="3" /> CHIARA STRINGARI,<fnref n="1" />,<fnref n="4" /> ID AND EMMANUEL BEAUREPAIRE<fnref n="1" />,*, ID

Received 12 March 2025; revised 12 June 2025; accepted 16 June 2025; published 15 July 2025

Myelin is essential for axonal conduction and metabolic support. To better understand its role in health and disease, it is necessary to establish accurate methods for *in situ* mapping of myelin at scales ranging from submicrometer to centimeters. Third-harmonic generation (THG) microscopy has recently been proposed as an efficient label-free method to visualize myelin in thick and living tissue. However, the contrast mechanism of THG from myelinated axons is complex and poorly described, which has limited the development of THG as a quantitative probe of myelin distribution. Here, we present a systematic characterization and modeling of polarization-resolved THG (pTHG) signals from individual axons as a function of their diameter and myelin thickness, and we show that pTHG can be used to derive myelin scores in several biological systems. First, we confirm the sensitivity and specificity of the THG contrast for myelinated axons in mouse brain tissue and its ability to detect isolated micrometer-sized axons oriented both in-plane and out-of-plane. We then present a detailed characterization of the pTHG contrast of small and large axons in live zebrafish larvae at different developmental stages, and we demonstrate that pTHG detects early axon development *in vivo*. We show that classical models of coherent multiphoton microscopy fail to reproduce pTHG profiles of axons because they neglect myelin-induced optical aberrations, and we establish a numerical strategy based on finite-difference time-domain calculations that can accurately relate pTHG signal profiles to axonal diameter and myelin thickness. Finally, we illustrate the relevance of pTHG microscopy for characterizing myelin distribution at different scales in fixed mouse and human brain tissue.

Published by Optica Publishing Group under the terms of the [Creative Commons Attribution 4.0 License](https://creativecommons.org/licenses/by/4.0/). Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

https://doi.org/10.1364/OPTICA.562091

## 1. INTRODUCTION

The vertebrate nervous system comprises a complex network of neurons whose axons are surrounded by a multilayered lipid membrane, the myelin sheath. The structure of myelin enables essential functions such as rapid neuronal communication and axon metabolic support [1–3]; therefore, myelin plasticity dynamically shapes neuronal function and brain connectivity [4–10]. Specifically, neuronal conduction velocity is related to myelin thickness, which is often characterized by the ratio of the radius of the axon without myelin to the radius of the myelinated axon (g-ratio) [11,12].

Given the vital importance of myelin in the nervous system, its pathological degradation is the cause of serious neurological manifestations. For example, the autoimmune attack that occurs in the pathology of multiple sclerosis causes myelin loss and axonal injury in the central nervous system (CNS) [1,13]. This structural

alteration leads to major functional changes and is associated with severe impairment. The ability to detect changes in both the large-scale distribution and local structure of myelin is critical for further understanding of myelin physiopathology [14–16]. Magnetic resonance imaging has long been the standard method for *in vivo* whole-brain imaging [17,18], while light microscopy has been largely used for small-scale exploration of myelinated tissue with micrometric resolution [19]. Confocal fluorescence microscopy was established as an effective technique for sensitive imaging of stained myelinated axons. Two-photon excited fluorescence (2PEF) was then used to overcome the limited penetration depth of confocal microscopy in tissues. However, because fluorescent labeling of tissues can be challenging and may interfere with cellular function, the development of label-free methods with high sensitivity for myelin is also beneficial for non-invasive studies. In recent years, several label-free optical modalities with micrometer-scale resolution have been applied to myelinated tissue imaging,

<footnote n="1">¹Laboratory for Optics and Biosciences, CNRS, INSERM, École Polytechnique, Institut Polytechnique de Paris, 91128 Palaiseau, France</footnote>

<footnote n="2">²Institut de Biologie de l'ENS, École Normale Supérieure, CNRS, INSERM, Université PSL, 75005 Paris, France</footnote>

<footnote n="3">³Paris Brain Institute, Sorbonne University, Pitié Salpêtrière Hospital, INSERM, CNRS, Paris, France</footnote>

<footnote n="4">⁴chiara(stringari@polytechnique.edu</footnote>

<footnote>*emmanuel.beaurepaire@polytechnique.edu</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** THG contrast from myelinated axons in fixed mouse nervous tissue. (a)–(e) Nonlinear microscopy setup and principles. (a) Microscope layout. See Section 5. (b) Multimodal THG-2PEF acquisition sequence. (c) Basic principles of THG contrast. The THG intensity profiles (blue lines) depend on the size and geometry of the object (gray) in relation to the focal volume of the excitation beam (red). Small objects yield the maximum signal when their size is comparable to the excitation focus. THG signals peak when the beam is focused on vertical (middle) and horizontal (right) optical discontinuities. (d) Axon and myelin structure. (e) Sample mounting and geometry. gm, gray matter;wm, white matter. (f) Multiscale THG imaging of myelin in a spinal cord section. df, vf, lf, dorsal, ventral, lateral funiculus; dh, dorsal horn; vh, ventral horn. Scale bars = 100, 50, and 10 µm. (g) THG-2PEF imaging in the cortex of plp-GFP mouse brain tissue, illustrating the visibility of individual axons both in-plane (layer II–III) and out-of-plane (layer I). ctx, cortex; cc, corpus callosum. Scale bars = 100 and 20 µm. (h) THG-2PEF imaging in the corpus callosum and cortex of plp-GFP mouse brain tissue. THG and 2PEF signals can both be used to distinguish domains with different myelin densities, albeit with a different contrast. Scale bar = 100 µm.</figcaption></figure>

ranging from linear techniques such as polarization microscopy [20], optical coherence tomography (OCT) [21,22], and spectral confocal reflectance (SCoRe) [23] to nonlinear techniques including coherent Raman scattering methods [24] and third-harmonic generation (THG) microscopy [25].

In particular, THG is a nonlinear process in which light of wavelength $\lambda$ is converted to $\lambda/3$ by a medium under pulsed laser illumination [Figs. 1(a) and 1(b)]. Under the tight focusing conditions used in microscopy, THG signals are observed near discontinuities in the refractive index $n$ and/or the third-order susceptibility $\chi^{(3)}$ [Fig. 1(c)] [26–28]. Therefore, THG images highlight interfaces and submicrometer heterogeneities, such as lipid droplets, due to the optical difference between lipid and water [29], while being compatible with other multiphoton contrast modalities. THG emission occurs mainly in the direction of the excitation beam, but signal epi-detection is also possible in thick tissue by taking advantage of THG light backscattering [29,30]. For these reasons, THG microscopy has become a valuable tool

for 3D morphological imaging of tissues [31–38]. Few studies have investigated the suitability of THG microscopy for myelin imaging. Farrar *et al.* [25] first demonstrated THG sensitivity to CNS myelin in *ex vivo* and *in vivo* mouse and zebrafish models using 1043 nm excitation. Witte *et al.* [39] confirmed THG contrast on myelinated axons as well as additional structures such as blood vessels in living mice with 1200 nm excitation. Lim *et al.* showed that THG microscopy allows the detection of myelin domains in large-diameter sciatic nerves of the peripheral nervous system [40], as well as the characterization of myelinated axon distribution in cortical tissue *in vivo* [41]. Although these studies and a few others [42–45] have demonstrated the sensitivity of THG imaging to myelin, no model has been developed to explain the coherent signal produced by a myelinated axon. In addition, the dependence of axon THG signals on the polarization of the excitation beam has not been investigated. It is generally accepted that the THG contrast observed in axons is due to discontinuities in the optical properties between the myelin sheath and its internal
and external environment. Early numerical analyses simply considered a horizontal lipid slab immersed in water [25]. However, recent studies of THG from vertical interfaces in the presence of index discontinuity [46,47] suggest that axon edges may induce focal field distortions that strongly affect THG signals. Myelin birefringence [48–51] due to the packing of lipid polar molecules [52,53] [Fig. 1(d)] could also have an effect on the THG signal. Thus, a more complete model that includes both geometrical and optical properties is needed to understand the THG contrast observed in axons. Finally, the potential of THG to probe the state of myelination of an axon should be explored.

In this work, we introduce the use of polarization-resolved THG (pTHG) [54,55] for mapping myelin distribution at submicrometer-to-tissue scales, and we develop a complete model that successfully accounts for the complex pTHG contrast patterns observed on individual axons as a function of axon radius and myelin thickness. In the first part, we present multiscale THG images of myelinated axon networks in fixed mouse tissue. We illustrate the visibility of submicron cortical myelinated axons with different orientations by comparing 2PEF and THG contrasts in transgenic fluorescent mouse lines. We then perform a systematic analysis of pTHG contrast in reticulospinal fibers of developing zebrafish larvae, and we demonstrate the suitability of pTHG for probing early stages of axon myelination *in vivo*. We then establish a numerical modeling strategy capable of accurately reproducing the pTHG signals observed in myelinated axons as a function of axon radius and myelin thickness. Finally, we illustrate the potential of THG imaging to study myelin distribution under normal and pathological conditions in mouse and human tissues.

## 2. RESULTS

### A. Myelinated Axons Are a Major Source of THG Contrast in Fixed Mouse Brain Tissue

For all the work presented in this paper, we performed THG imaging with an excitation wavelength of ≈1100 nm to avoid possible three-photon resonance from blood in THG signals [56,57]. Following previous studies [25,41], we first investigated the specificity and sensitivity of THG contrast for myelin in fixed mouse nervous tissue (Fig. 1). Figure 1(f) shows multiscale THG images of a fixed transverse section of mouse spinal cord tissue. Large-area THG imaging shows tissue morphology with intensity differences between white matter (ventral, lateral, and dorsal funiculus) and gray matter (ventral and dorsal horns). Higher-resolution THG imaging resolves individual myelinated axons near the ventral horn and funiculus. Axons from white matter tracts orthogonal to the imaging plane appear as contrasting rings surrounding a dark region. Individual in-plane axons appear as parallel lines.

To characterize the THG contrast of submicrometer myelinated axons in different orientations, we acquired xyz images in the cortical region of a mouse brain slice [Fig. 1(g)]. We used a transgenic mouse expressing green fluorescent protein (GFP) associated with the proteolipid protein (plp) of axonal myelin to allow identification of myelinated axons [58]. We acquired high-resolution THG-2PEF images to a depth of 100 µm in the cortical region near the midline. In-plane and out-of-plane myelinated axons were equally visible in the tissue, highlighting the sensitivity of THG contrast to submicron diameter tubular structures, regardless of their orientation. Specifically, cortical layer I showed

a majority of transverse myelinated axons visible on an XZ projection, while layers II and III showed a majority of myelinated axons lying in the XY plane. These observations are consistent with the work of Redlich and Lim [41], who described THG contrast across cortical layers *in vivo*. To compare the sensitivity of lower-resolution THG and 2PEF images to myelin density, we imaged areas of the cortex (gray matter) and the corpus callosum (white matter) on a coronal mouse brain slice [Fig. 1(h)]. The two regions were easily identifiable using 2PEF contrast, and the overall similarity of THG with 2PEF across the field of view confirms that myelin is a dominant source of THG contrast here. Cortical areas showed lower signal intensity and heterogeneous contrast with both modalities, whereas the corpus callosum showed brighter signals with more homogeneous contrast. Figure 1(h) shows histograms of THG and fluorescence signals recorded in the cortex and corpus callosum. These measurements confirm that THG contrast can be used to distinguish areas with different myelin densities. However, THG intensity is not proportional to the corresponding fluorescence signals, and therefore to myelin concentration.

Overall, this first set of multiscale THG imaging experiments confirms several key observations. (i) THG highlights small and large axons, as well as myelin-rich regions. (ii) THG highlights both in-plane and out-of-plane axons. This is true even for small cortical axons and is a distinguishing feature of THG imaging over optical imaging modalities based on reflectance. (iii) The THG contrast is not trivial and does not directly report on myelin density (as assessed by plp-GFP fluorescence levels). This observation highlights the need to accurately characterize and model the axon THG signal in order to develop quantitative THG-based measurements.

### B. pTHG Imaging of Axons in the Spinal Cord of Live Zebrafish Larvae

We then investigated the possibility of using polarization-resolved THG (pTHG) contrast to detect changes in the morphology of individual axons during myelination. Experimentally, we performed a systematic characterization of the pTHG contrast in the spinal cord of 2–5 days post-fertilization (dpf) zebrafish larvae using a fast polarization-switching scheme (see Section 5) [55]. Numerically, we developed a finite-difference time-domain (FDTD) approach to predict the pTHG contrast as a function of the dimensions and optical properties of an individual myelinated axon, taking into account the local field distortions encountered by the excitation beam as it is scanned across the axon. In the following sections, we first discuss the pTHG contrasts observed *in vivo* in the zebrafish larval spinal cord (Fig. 2) during development, and then present the numerical approach and compare it with the experimental observations.

Although the THG contrast highlights myelinated axons, it is not chemically specific to myelin or lipids but reveals interfaces and heterogeneity in general, as well as strongly absorbing structures [57] and biominerals [54,59]. This is illustrated in Figs. 2(a) and 2(b) with a large-scale reconstructed THG image showing a transverse plane in a live 5 dpf zebrafish larva oriented laterally. Several structures are visible, such as the boundaries of the fins, the notochord vacuoles, and the spinal cord. THG images from the spinal cord [Fig. 2(c), Visualization 1] show large Mauthner axons (ma) with up to 10- to 15-fold enhanced contrast compared to nearby cell structures and a tubular appearance reminiscent of previously
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** pTHG imaging of spinal cord axons in live zebrafish embryos. Larvae from the mbp-mCherry-CAAX line were imaged at 5 dpf. (a) Schematic representation of a zebrafish larva in coronal orientation and an enlarged diagram of the spinal cord. (b) Reprojection in the same orientation of a 3D large-scale THG-2PEF image. sc, spinal cord; nc, notochord; df, dorsal fin; vf, ventral fin. (c) Zoomed-in image of the spinal cord. ma, Mauthner axons; d, dorsal axons; v, ventral axons. (d), (e) Implementation of the polarization-resolved THG acquisition. (d) Experimental principle: an electro-optical modulator (EOM) rotates the polarization direction of the excitation beam between acquiring two successive lines of a pTHG image. Static half-wave plates (HWP) and quarter-wave plates (QWP), installed around the EOM and before the objective (Obj), minimize polarization ellipticity. (e) pTHG acquisition sequence: each line is acquired several times with different polarization states $P_n$ (4 shown here and 18 for most experiments) in an X-P-Y sequence, to ensure a temporal precision in the ≈20 – 40 ms range for the polarimetric measurement, depending on the size of the image in the x direction. XYP images are then reconstructed. (f) Example of pTHG contrast in the spinal cord of a 5 dpf embryo. Left image: region of interest. Other images: THG with an incident polarization parallel (THG₀) and orthogonal (THG₉₀) to the axon axis, and 2PEF fluorescence from labeled myelin. (g) Definition of the pTHG modulation amplitude and angle. (h) Corresponding maps for the images shown in (f). See also Visualization 1, Visualization 2 and Visualization 3.</figcaption></figure>

described myelinated axons in the mouse spinal cord [Fig. 1(f)]. Smaller dorsal (d) and ventral (v) axons are also visible, although they produce a signal only 3 to 4 times higher than nearby cells. To confirm that these signals originate from axons, we used the fluorescent transgenic line mbp-mCherry-CAAX [6], which expresses the fluorescent protein mCherry under a myelin basic protein (mbp) promoter. We then acquired a single-plane pTHG series

on a Mauthner and ventral axon bundle oriented parallel to the imaging plane. Rotating the excitation polarization over an angular range of 180° resulted in significant changes in THG intensity [Fig. 2(f), Visualization 2 and Visualization 3]. We characterized the THG intensity variations using the modulation parameter M, defined as $M = 4a_1/(a_1 + 2a_0)$, where $a_0$ and $a_1$ are the two first orders of the pTHG signal Fourier transform. This modulation
parameter is equal to $(THG_{max} - THG_{min})/THG_{max}$ for a $A + B \cos^2(\phi - \phi_1)$ dependence on the polarization orientation $\phi$ as contrasts obtained at vertical interfaces. As shown in Fig. 2(h), a pTHG modulation amplitude of $M = 60-80\%$ is consistently observed for myelinated axons, while nearby structures such as cell boundaries generally exhibit $M = 20-50\%$. This strong dependence on polarization orientation can be used to isolate axons by applying appropriate filtering (see Section 5 and Figs. S1 and S2 in Supplement 1). These experiments also lead to two interesting observations: (i) On large 5 dpf Mauthner axons with thick myelin sheaths, where both lateral axon boundaries are visible, the THG profile across an axon boundary exhibits a single peak when the incident laser polarization is parallel to the axon (THG₀) and a double peak when the polarization is orthogonal to the axon (THG₉₀). (ii) THG intensity from these axons is higher when the incident polarization is parallel to the axon, whereas THG intensity from 5 dpf dorsal axons is higher when the incident polarization is perpendicular to the axon (in other words, the modulation is phase-shifted by 90° in small axons compared to large axons). Lastly, as a first approximation, regions of stronger pTHG modulation correspond to the location of 2PEF signals from labeled myelin [Figs. 2(f)-2(h)].

### C. pTHG Profiles of Axons at Different Stages of Development

We further analyzed the pTHG response of dorsal and Mauthner axons at different developmental stages between 2 and 5 dpf. Variations in axonal radii and myelin thickness based on the g-ratio have been documented in the zebrafish larva during this period (see Table 1), as well as the onset of myelination, which occurs before 3 dpf [60-62]. It has been reported [60] that myelin sheath thickness in the developing Mauthner axon is regulated according to increasing axon diameter (resulting in a constant g-ratio), whereas in dorsal and ventral axons, myelin thickness increases while axon diameter remains constant, resulting in a variable g-ratio (see Table 1 and [60]). Figures 3(a) and 3(b) show representative pTHG intensity, modulation amplitude and angle maps and corresponding reference 2PEF from labeled myelin (mbp-mCherry-CAAX) displayed on small domains around a 4.5 dpf dorsal axon and 2-to-5 dpf Mauthner axons. We observed that the Mauthner axon is visible in THG images as early as 2 dpf, whereas smaller axons are visible only at the onset of myelination (2.5 or 3 dpf). 2PEF and pTHG profiles were obtained from 2PEF and pTHG images averaged over a few tens of microns along the axons after applying coregistration and correlation operation to straighten the axons (see Section 5 and Fig. S1 in Supplement 1). When plotting the modulation profiles, we used the convention that $M$ is positive when $THG_{90} > THG_0$ to account for the modulation phase shift. Figure 3(c) and Visualization 4 also show

<table><caption>**Table 1.** Axon Diameter and Myelin Thickness for Different Neurons in Zebrafish Spinal Cord Measured by Electron Microscopy [60]</caption><thead><tr><th>Axon Type</th><th>Axon Diameter (µm)</th><th>Myelin Thickness (nm)</th><th>g-Ratio</th></tr></thead><tbody><tr><td>Dorsal, ventral</td><td>~ 1</td><td>50-250</td><td>0.8-0.95</td></tr><tr><td>Mauthner 3 dpf</td><td>3-4</td><td>150-200</td><td>0.9</td></tr><tr><td>Mauthner 5 dpf</td><td>4-5</td><td>200-300</td><td>0.9</td></tr></tbody></table>

pTHG data from a mouse Purkinje axon in cultured mouse cerebellum tissue for comparison. Figure 3(d) shows the statistical analysis of the amplitude of the spatial pTHG modulation variation defined as $\Delta M = M_{out} - M_{in}$ [see Fig. 3(c)], as well as THG and 2PEF peak intensities in Mauthner and dorsal axons between 2 and 5 dpf (n = 5 at 48, 52, 56, 60 hpf and 5 dpf; n = 10 at 3 and 4 dpf). These measurements provide much original information and require a number of comments. First, the pTHG profile across an individual axon can vary significantly with the development of the size of the axon and its degree of myelination [Figs. 3(a)-3(c), column “pTHG”]. In this regard, the Mauthner axon has a remarkable development starting at 48 hpf. At this stage, when the axon is not myelinated and its diameter is less than 3 µm [60], the pTHG profiles at 0° and 90° show a unique central maximum. During the early stages of myelination (60 hpf-3 dpf) [60], lateral peaks appear on the THG₉₀ profile. As the Mauthner axon development continues between 3 and 5 dpf and axonal diameter increases along with myelin thickness [60], the pTHG profiles become more complex, resulting in a tubular appearance with double peaks at the edges when the incident polarization is orthogonal to the axon. Note that at 5 dpf, the Mauthner axon diameter is ~4.5 µm, and the myelin thickness is ~220 nm [60]. The measurements summarized in Fig. 3(d) show that the myelination of the Mauthner axon as well as that of small dorsal axons can be detected at very early stages by an increase in the pTHG modulation gradient $\Delta M$. In dorsal axons, $\Delta M$ increases from 0.52 ± 0.05 to 0.72 ± 0.05 between 3 and 5 dpf. In Mauthner axons, $\Delta M$ increases from 0.48 ± 0.07 to 0.76 ± 0.07 between 2 and 3 dpf. In comparison, the low level of expression of mbp-mCherry fluorescent reporters before 3.5 dpf makes this detection difficult from the 2PEF signal.

As a consequence of the complex dependence of pTHG profiles on axon morphology, the angular modulation maps also show a non-trivial evolution. Depending on the characteristics of the axon, some regions of the image give more THG for incident polarization parallel to the axon (“positive” modulation) and other regions for polarization orthogonal to the axon (“negative” or 90°-shifted modulation), all of which evolve as the axon develops.

### D. Advanced Numerical Modeling of pTHG Contrast from Myelinated Axons

Since no model has been proposed in the literature that can account for the previous observations, we developed a suitable numerical framework to predict the pTHG profiles of axons as a function of their optical and physical properties. A critical implementation decision was to work with a numerical FDTD approach that can simulate the distortions experienced by the excitation beam when focused on an axon. Recent work from our group has shown that THG profiles across index-mismatched vertical interfaces (e.g., lipid vesicles in water [46] or stretched micropipette tips [47]) exhibit a strong polarization dependence related to the distortions experienced by the excitation beam near its focus. These profiles can be accurately predicted using a FDTD strategy, in contrast to more classical coherent nonlinear models [28,63,64]. Myelin forms a hollow cylinder with a thickness of tens to hundreds of nanometers and has a high refractive index compared to surrounding structures. It can be expected that the THG profile observed on a myelinated axon is influenced by beam distortions near the focus. Therefore, we extended our FDTD approach to model a myelinated axon with micrometric dimensions. We considered
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** pTHG profiles from axons at different development stages. (a) pTHG contrast from dorsal axons in 4.5 dpf zebrafish embryos. Scale bar: 1 µm. (b) Evolution of pTHG contrast on Mauthner axons in 2–5 dpf embryos. Scale bar: 3 µm. (c) pTHG contrast from Purkinje axons in an organotypic mouse cerebellum slice. See also Visualization 4. Scale bar: 3 µm. (a)–(c) THG profiles were measured across individual axons for different directions of incident polarization. The profiles are shown in colors ranging from green (THG₀, polarization parallel to the axon) to blue (THG₉₀, polarization perpendicular to the axon). pTHG modulation profiles are shown as solid black lines. Dashed blue lines show the corresponding THG₉₀ profiles for comparison. In the rightmost column, the 2PEF profiles are shown as solid red lines, and the corresponding THG₉₀ profiles are shown as dashed blue lines for comparison. Horizontal axis labels correspond to *x* in µm; vertical axis values are normalized. (d) Measurement of pTHG modulation gradient Δ*M*, peak THG intensity, and peak 2PEF intensity in zebrafish Mauthner and dorsal axons of 2–5 dpf mbp:mCherry-CAAX zebrafish larvae.</figcaption></figure>

the simplified geometry of a cylinder characterized by two geometrical parameters (axon radius and myelin sheath thickness) and three media with different optical properties (axon, myelin, and external) [Fig. 4(a)]. One difficulty in constructing such a numerical model is that it requires realistic values for both the linear and nonlinear susceptibilities of the three media, for which there are few reported measurements in the literature. For the linear indices,

we chose values close to those used in [65], slightly lowered to account for the longer excitation wavelength used in THG imaging: $n_{axon} = 1.37$; $n_{myelin} = 1.43$; $n_{ext} = 1.39$. For the nonlinear susceptibilities, we chose values based on experimental measurements in liquids [66] [Fig. 4(a)]. The value for $\chi^{(3)}_{myelin}$ was taken close to the value measured for lipids, i.e., higher than $\chi^{(3)}_{axon}$ and $\chi^{(3)}_{ext}$, which were chosen closer to the value measured for aqueous
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** FDTD modeling of axon pTHG. (a) Model geometry and parameters. (b) Calculated THG profiles in the equatorial plane of a myelinated axon for a polarization parallel (THG₀) and perpendicular (THG₉₀) to the axon, for different axon radii and myelin thicknesses. Horizontal axis labels correspond to *x* in µm. Vertical yellow bands indicate the position of the myelin sheath. Vertical axis labels are in arbitrary units. Additional profiles are shown in Figs. S5–S7 in Supplement 1.</figcaption></figure>

solutions [66]. One implementation choice was to set $\chi_{\text{axon}}^{(3)} < \chi_{\text{ext}}^{(3)}$, consistent with the observation that the interior of axons appears less dense than the extracellular medium in both published electron microscopy data [6,60,67] and in CARS microscopy images, which depend directly on $\chi^{(3)}$ (Fig. S3 in Supplement 1 and [68]). As discussed later, this choice turned out to be important for the shape of the simulated profiles. An additional assumption was to neglect the anisotropy and the wavelength dependence of the susceptibilities in a first approximation in order to limit the number of model parameters. We therefore chose $(\chi_{\text{axon}}^{(3)}, \chi_{\text{myelin}}^{(3)}; \chi_{\text{ext}}^{(3)}) = (1.7; 2.5; 1.9) \times 10^{-22} \, \text{m}^2 \, \text{V}^{-2}$. To model imaging conditions similar to those of the experiments, we made the following additional assumptions: an excitation wavelength of 1.2 µm, an absence of a resonant term in the nonlinear susceptibility, and a numerical aperture of 1.0. Using these parameters, we calculated theoretical THG intensity profiles in the equatorial plane of in-plane axons for incident polarizations parallel and orthogonal to the cylinder axis (Fig. S4 in Supplement 1), for axon radii ranging from 1.0 to 2.5 µm, and for myelin thicknesses ranging from 30 to 800 nm [Fig. 4(b) and Figs. S5–S7 in Supplement 1]. Remarkably, these calculations predicted pTHG intensity and modulation profiles that exhibited a wide variety of shapes. In this geometry, where characteristic sample sizes are comparable to the dimensions of the focal volume, the far-field pTHG profiles result from interference

between three contributions and are affected by polarization-dependent field distortions caused by the nearby vertical interfaces. In particular, the FDTD simulation was able to recapitulate our experimental observations for corresponding axon dimensions (Fig. 5). Figure 5(a) displays THG profiles for an excitation parallel (“THG₀,” in green) and perpendicular (“THG₉₀,” in blue) to the cylinder for four different geometries. As myelin thickness and axon radius increase, the simulations predict a progressive transition from a single-peak profile with weak polarization dependence to a multi-peak profile with strong polarization dependence. Characteristic shapes and modulation inversions observed experimentally (Fig. 3) are correctly predicted. We have verified that the shape of the simulated profiles is preserved when the focusing plane is shifted by $\pm 0.5$ µm for the [r = 1250, t = 100] geometry (Fig. S8 in Supplement 1). We also found that the shape of complex experimental profiles was less well reproduced when reversing the $\chi^{(3)}$ values between inside and outside the axon (Fig. S9 in Supplement 1).

The above results confirm the suitability of the FDTD approach for modeling the pTHG profiles observed on a myelinated axon and demonstrate that the complex shapes of pTHG profiles observed on axons can be explained by interferences originating from three media stacked on a subwavelength scale and by field distortions induced by vertical interfaces. It is then particularly interesting to use the model to explore the sensitivity of
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Numerical and experimental axon pTHG profiles. (a) Experimental and calculated THG profiles in the equatorial plane of a myelinated axon for a polarization parallel (THG₀) and perpendicular (THG₉₀) to the axon, for different axon radii and myelin thicknesses. Horizontal axis labels correspond to *x* in µm. Vertical axis values are normalized. (b) Calculated maps of THG₀ intensity and pTHG modulation as a function of axon diameter and myelin thickness. The stars correspond to the profiles shown in (a). (c) Calculated pTHG intensity profiles as a function of increasing myelin thickness for a 1.75 µm radius axon. Vertical yellow bands indicate the position of the myelin sheath.</figcaption></figure>

THG profiles to myelin sheath thickness for a given axon diameter. We considered a model neuron with an axon radius of 1750 nm and different myelin thicknesses. Figure 5(c) shows the calculated THG profiles for two orthogonal polarizations. As the amount of myelin increases, THG profiles generated with excitation polarization parallel to the axon (green curves) show a single peak profile of increasing intensity, while THG profiles generated with excitation polarization perpendicular to the axon (blue curves) gradually evolve from a single peak to a double peak profile and show a smaller increase in intensity.

Figure 5(b) and Fig. S10 in Supplement 1 summarize the dependence of pTHG profiles on axon radius and myelin thickness for all the simulations shown in Figs. S5–S7 in Supplement 1. The map on the top shows that the maximum intensity of the THG₀ profile increases in a nearly quadratic manner with myelin thickness for the range of axon sizes considered. The map on the bottom shows that the pTHG modulation gradient Δ*M* has a more complex dependence on myelin thickness and axon radius and is maximal for myelin thicknesses in the range of 100–200 nm. In summary, these analyses confirm that (i) the FDTD numerical strategy is valid for explaining the pTHG signals observed experimentally on single axons, and that (ii) pTHG imaging is able to detect early stages of axon maturation and myelin thickness changes in the submicrometer range.

E. Relative Effects of Axon Geometry and Myelin
Anisotropy

The above results showed that the pTHG images of axons result
from the sample geometry at the subwavelength scale. We then
numerically and experimentally evaluated the possible contribu-
tion of optical anisotropy of the myelin sheath to the observed
pTHG profiles. Figures S11 and S12 in Supplement 1 show the
changes in the predicted profiles for five different geometries when
linear (birefringence) and nonlinear anisotropy are included in
the simulation. Myelin birefringence was assumed to be positive
(n<sub>e</sub> = 1.445, n<sub>o</sub> = 1.435, and Δn = 0.01) with a higher index
for a polarization parallel to the lipid tails, similar to experimental
data reported in the literature [69]. For the nonlinear anisotropy
Δχ<sup>(3)</sup>/χ<sup>(3)</sup>, we considered several values between 0 and 12%,
since no direct measurement is currently available [70]. Figures S11
and S12 in Supplement 1 show that anisotropy can alter pTHG
profiles and modulation. This effect is generally larger for the
THG<sub>90</sub> profile, i.e., when the incident polarization is parallel to the
lipid tails. However, changes remain limited (<20%) for myelin
thicknesses less than 300 nm. The effects on the pTHG profiles
are non-trivial due to the complex phase matching involved in the
axon geometry, but overall, the profile shapes are preserved. These
simulations suggest that the experimentally observed pTHG con-
trast depends primarily on sample geometry, i.e., myelin thickness
and axon radius, and marginally on myelin anisotropy.
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Effect of index mismatch and material anisotropy on pTHG in simple experimental geometries. (a) pTHG modulation when the excitation beam is focused on a vertical glass-water or glass-TDE interface after 7 µm propagation. The amplitude and sign of the measured modulation are strongly dependent on the index discontinuity and can be modeled using a FDTD approach. (b) Measured pTHG modulation at the vertical edge of amorphous oil droplets (LDs, n = 5) and organized MLVs (n = 5) embedded in agarose. pTHG is sensitive to lipid ordering, but index discontinuity accounts for most of the pTHG modulation.</figcaption></figure>

In a complementary manner, we examined experimental pTHG profiles from simple model samples (Fig. 6). We considered (i) a simple vertical interface in the presence of an adjustable index discontinuity to evaluate the influence of index mismatch (as done in [47]), and (ii) oil droplets and multilamellar lipid vesicles (MLVs) embedded in agarose to evaluate the influence of molecular ordering and anisotropy. In the vertical interface experiment, the beam was focused on a quartz-liquid interface after 7 µm of propagation. The liquid medium used was water and TDE at different concentrations. As we previously reported [47], pTHG in this geometry is governed by the beam distortion along the interface. pTHG modulation can take a wide range of values from +0.45 to -0.3 depending on the index mismatch, which can be reproduced by FDTD calculations. In the droplets experiment, we imaged the vertical edge of amorphous oil droplets (LDs, n = 5) and organized MLVs (n = 5) embedded in agarose. This mounting protocol ensured that the sphericity of the droplets was preserved, in contrast to previous experiments where LDs were spread on a glass coverslip [70]. We measured a pTHG modulation of 57 ± 6% on oil droplets and 72 ± 3% on MLVs. This experiment confirms that pTHG is sensitive to lipid ordering, but that index discontinuity, rather than material anisotropy, accounts for most of the pTHG modulation.

## F. Measuring Axon Diameters Using pTHG

Having established the origin of the pTHG contrast in myelinated axons, we now explore the possibility of measuring relevant morphological parameters based on THG signals. When the axon diameter exceeds 2–3 µm, THG profiles exhibit well-resolved peaks that have been previously interpreted as a signature of the myelin sheath or its boundaries and used to estimate the axon

diameter [25,40,41]. We can revisit this type of measurement using our novel numerical framework.

An interesting prediction of the FDTD calculations is that, contrary to what was assumed in previous work [40,41], axon THG peaks do not generally colocalize with myelin boundaries. In particular, for large axons (1.5 µm radius or more), THG₀ maxima are localized within the inner region and tend to reduce the apparent axon diameter [Fig. 7(a)]. In order to establish guidelines for accurate axon size measurements, we therefore measured this offset in simulated THG₀ THG₉₀, and 2PEF profiles [Fig. 7(b)]. Simulations were run for axon radii of 1, 1.5, and 2 µm, and myelin thicknesses ranging from 0 to 600 nm. The black dots in Fig. 7(b) indicate the theoretical positions of the center of the myelin sheath for the different cases, for reference. As expected, 2PEF signal maxima provide a relatively accurate estimate of the myelin position, although they exhibit a small bias (2–7%) toward the axon interior [see orange stars in Fig. 7(b)]. THG₀ maxima (green squares) systematically underestimate axon dimensions by amounts that can reach 30% and cannot be used directly for size measurements. Finally, the outer peaks of the THG₉₀ profiles (blue triangles) generally colocalize well with myelin and provide a good probe of axon diameter. To explore the possibility of characterizing axon growth during development, we estimated the Mauthner axon diameter as a function of rostro-caudal distance in zebrafish larvae at different stages (five larvae/stage, six areas imaged/larva) based on the outer peaks of the THG₉₀ profiles [Fig. 7(c)]. These measurements [Fig. 7(d)] illustrate that this approach recovers the diameter of the Mauthner axon as a function of time and space, with smaller diameters toward the tail and larger diameters as the fish develops.

## G. Multiscale pTHG Imaging of Fixed Mouse and Human Tissue

Finally, as a complement to the previous analysis of the contrast mechanism of polarized THG in single axons, we present preliminary explorations of the utility of pTHG microscopy for more general characterization of myelin organization in biomedically relevant contexts.

### 1. pTHG Contrast in a Model of Hypomyelination

We investigated the potential of pTHG contrast to study myelin disorganization in a dysmyelination model tissue. *Shiverer* mice are myelin basic protein (mbp)-deficient mice with altered or absent myelin sheaths in central nervous system axons, resulting in tremor behavior [71,72] (Fig. 8). We compared the polarimetric THG response obtained from white matter tracts in fixed normal (n = 24) and *shiverer* (n = 32) cerebellar folia. Three regions of interest were analyzed in each folium. Under our sampling and resolution conditions, not all fibers could be resolved individually. We found that THG intensity in myelinated areas was 0.4 times lower in *shiverer* mice compared to normal mice in agreement with previous observations [25]. In addition, we measured that pTHG modulation in the central axon bundle was 33 ± 3% in wild type and 26 ± 3% in *shiverer* tissue (Fig. 8). These observations suggest that pTHG provides relevant contrast even when individual axons are not resolved. Polarization-sensitive imaging adds a dimension to the differentiation between *shiverer* and wild type and provides increased sensitivity to myelin disruption or disorganization compared to simple intensity-based imaging.
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Monitoring Mauthner axon growth with THG. (a) FDTD simulation of the THG signal from an axon with a 2 µm radius and 200 nm myelin thickness, for two orthogonal excitation polarizations (THG₀ in green and THG₉₀ in blue). (b) Suitability of 2PEF and THG peaks for estimating axon radius (FDTD simulations). Simulations were run for axon radii of 1, 1.5, and 2 µm and for myelin thicknesses ranging from 0 to 600 nm. The black dots indicate the positions of the center of the myelin sheath for reference. The 2PEF signal maxima (stars) and the outer peaks of the THG₉₀ profiles (blue triangles) colocalize well with the myelin and effectively probe the axon diameter. THG₀ maxima (green squares) cannot be used directly for size measurements. (c), (d) Experimental measurement of the apparent size of Mauthner axons as a function of position in the fish (distance from the hindbrain) and as a function of age using the second peak of the THG₉₀ signal.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 8.** pTHG contrast in a model of hypomyelination. (a) Myelin distribution in cerebellum folia. (b) THG images of normal (wild-type) and hypomyelinated (*shiverer*) mouse cerebellar tissue. Scale bar: 100 µm. (c) pTHG images recorded in the regions indicated by the red stars in (b). Scale bar: 15 µm. (d) Measured pTHG intensity and modulation in wild-type and *shiverer* folia.</figcaption></figure>

## 2. Multiscale THG Imaging of Fixed Human Tissue

We also investigated the feasibility and potential relevance of multiscale pTHG microscopy for the characterization of *post-mortem* human brain tissue samples. We acquired THG images of a 20 µm × 2.0 cm × 1.6 cm slice of fixed human cerebellar tissue using a tiled acquisition scheme. Figure 9 shows a stitched 60 × 40-tiles mosaic recorded with micrometer-scale sampling. Large-scale THG imaging reveals the tree structure of the white

matter fibers, which extend from the main trunk down to the individual folia. Granular and molecular layers show lower THG signals. In zoomed views, THG contrast reveals the complete structure of individual folia, including the central axon bundle, the heterogeneous granular layer with bright structures corresponding to myelinated fibers, and the molecular layer producing low THG signals. We note the similarity of the THG contrast obtained on human and mouse cerebellar folia. Finally, pTHG recordings
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Multiscale THG imaging of fixed human cerebellum tissue. Tiled high-resolution THG imaging of *post-mortem* human brain tissue. Large-scale imaging shows global myelin density, while THG contrast in zoomed views shows myelin distribution at micrometer resolution. The lower right panel shows pTHG contrast obtained from a single in-plane Purkinje axon.</figcaption></figure>

on smaller areas containing individually resolved Purkinje axons yield profiles very similar to those obtained on mouse Purkinje axons and mature zebrafish Mauthner axons (Fig. 2). This suggests that the modeling shown in Fig. 4 also applies to human axons. Overall, this proof-of-principle experiment shows that polarized THG imaging is a potentially very informative contrast modality for characterizing myelin distribution in mammalian tissues at micrometer to centimeter scales.

### 3. DISCUSSION

The study presented in this manuscript introduces the use of polarization-resolved THG microscopy for imaging myelinated axons and nervous tissue. Our primary goal was to describe and explain the pTHG contrast observed on axons of different sizes and myelin thicknesses to allow quantitative interpretation of these images. Our observations confirm that the THG contrast observed in mouse brain and nerve tissue reveals myelinated regions with good specificity when the excitation wavelength is around 1100–1150 nm. It should be noted, however, that this specificity depends

on the excitation wavelength. We have recently shown that when the excitation wavelength is shifted into the 1200–1350 nm range, the THG contrast changes and highlights red blood cells due to a strong three-photon resonance in hemoglobin [57]. Away from this resonance, the sensitivity of THG microscopy for myelinated axon detection is remarkable, as it can highlight individual cortical axons of micrometric size whatever their orientation relative to the imaging plane. If their diameter is greater than a few microns, axons appear as hollow cylinders in THG. This ability to detect axons of any orientation distinguishes THG from techniques based on the detection of reflected light, such as OCT or confocal reflectance microscopy, which essentially detect axons parallel to the imaging plane. In addition, we have shown that polarization-resolved imaging enhances the specificity of THG for myelinated axons and enables the detection of very early stages of axonal maturation in live zebrafish embryos (Fig. 2).

The pTHG contrast observed on an individual axon is particularly complex. It results from phase-matching conditions in a tightly focused geometry and in the presence of excitation field distortions along the axon's inner and outer vertical walls.
However, these pTHG profiles can be reproduced and explained with a FDTD numerical approach modeling propagation in a heterogeneous medium in the vicinity of the focal point (Fig. 4). The pTHG profiles can then be related to relevant morphological parameters, namely axon diameter and myelin sheath thickness. We have thus illustrated how THG signals should be used to reliably measure the radius of large axons (Fig. 7). A possible follow-up to this work would be to develop an automatic estimate of myelin thickness from experimental pTHG profiles based on the predictions of our numerical model.

It is interesting to note that the distortion of the excitation field due to index discontinuities in the focal region appears to be a major source of THG and pTHG contrast. Material anisotropy also influences this contrast but is not necessarily the dominant factor (Fig. 6, Figs. S9 and S10 in Supplement 1). Variations in polarization-resolved signals cannot therefore be interpreted solely in terms of sample anisotropy, particularly in the case of vertical interfaces. In the case of axons, birefringence and anisotropy have a more important influence on the pTHG contrast for myelin thicknesses beyond 100–150 nm. Thus, a follow-up to the modeling part could be to refine our numerical approach by performing a more systematic investigation of the effect of anisotropy in this size range. From a general perspective, contrast in THG microscopy has sometimes been viewed as a band-pass filter, highlighting spatial variations in optical properties such as interfaces and heterogeneities [30,33,73]. Recent findings ([46,47], this work) suggest that THG can also be considered as a detector of focal volume distortions.

Beyond the scale of the individual axon, we have also illustrated that THG contrast provides information related to myelin density at different scales and enables regions to be identified according to their richness in myelinated tracts (Figs. 1, 8, and 9). We have shown that polarization information adds a dimension enabling us to distinguish different myelin distributions even when axons are not individually resolved (Fig. 8). One perspective of this work is, therefore, to model the pTHG response of a bundle of myelinated axons and off-plane oriented axons. Finally, we have shown that large-scale, high-resolution imaging can provide insight into myelin organization from submicrometer to centimeter scale, even in post-mortem human tissue (Fig. 9).

pTHG microscopy, therefore, offers an interesting approach for characterizing myelin distribution in brain and nerve tissue at high resolution and without labeling. A distinctive advantage of multiphoton techniques (THG, CARS, and SRS) is their ability to detect out-of-plane fibers, unlike linear reflectance approaches. Although THG contrast is, in principle, less specific than vibrational signals (CARS and SRS), we have shown that its sensitivity to early stages of axonal myelination is excellent (Fig. 2), potentially superior to what has been reported with CARS microscopy [68]. An additional benefit of THG imaging is that it is directly compatible with other femtosecond pulse-based multiphoton microscopy modalities, such as 2PEF, SHG, and 3PEF microscopy. In particular, the recent development of OPA sources optimized for 3PEF microscopy [74,75] enables THG images to be recorded in epi-detection at millimeter depths in brain tissue [38,76,77], simultaneously with fluorescence imaging.

Finally, it should be noted that implementing polarimetric detection in a THG microscope involves several constraints. First, detecting polarized THG information requires at least two images, making it slower than simple intensity imaging with the

same excitation power. Second, introducing a polarization control system increases the microscope's complexity. However, we note that, in recent years, compact alternatives to the electro-optical modulator-based setup used in this work have been explored for pTHG imaging. These alternatives include the use of motorized waveplates [70] or the multiplexing of orthogonally polarized pulses [78].

## 4. CONCLUSION

To summarize, in this work we explored the relevance of polarization-resolved THG microscopy for probing the distribution and morphology of myelinated axons in fixed and living tissue. We first confirmed the sensitivity of the THG contrast to myelin density in *ex vivo* mouse brain tissue. In particular, THG highlights both in-plane and transverse micrometer-sized myelinated axons in cortical regions. We demonstrated the ability of polarimetric THG to monitor early stages of axonal myelination in live zebrafish embryos and characterized the complex pTHG profiles experimentally observed on Mauthner axons as a function of their developmental stage. We established a new numerical model based on FDTD calculations to relate the pTHG profile to axon morphological parameters. Finally, we illustrated how pTHG data can be used to obtain morphometric information at different scales in zebrafish, mouse, and human tissues. This work establishes pTHG microscopy as a new imaging method for mapping myelinated axons and paves the way for the interpretation of THG signals in terms of morphological parameters. As a perspective, it would be interesting to automate the determination of axonal diameters and myelin content at different scales in the context of axonal regeneration research. In addition, by taking advantage of the recent development of efficient laser sources and microscopes optimized for three-photon microscopy [38,44,74,75,77,79,80], this contrast can, in principle, be adapted to *in vivo* deep-tissue microscopy, and is directly compatible with fluorescence imaging.

## 5. METHODS

### A. Multiphoton Microscope

A custom-built microscope optimized for 2PEF and polarized THG contrast collection was used for the experiments [55]. Excitation was supplied by an OPO femtosecond laser source with a repetition rate of 80 MHz, a pulse duration of 100 fs, and a wavelength of $\lambda_{exc} = 1100$ nm (Insight DS+, SpectraPhysics, USA). Focusing was achieved with a $25 \times 1.05$ numerical aperture objective XLPLN25XWMP2, Olympus, Japan). Signals were detected using photon-counting PMTs (Sens-Tech, UK) and lab-designed MHz-rate counting electronics. 2PEF imaging was performed in the epi configuration, and signals were selected with a 625/25 (resp. 525/50) bandpass filter (Semrock) for mCherry (resp. GFP) labels. THG contrast was collected in transmission with a 377/50 filter (Semrock) on the detection path. Lateral and axial resolution of the THG measurements were evaluated to be 0.4 and 1.6 µm, respectively, measured as the FWHM of $(xy)$ and $(xz)$ THG intensity profiles retrieved from submicrometer polystyrene beads in agarose.
## B. Polarization-Resolved THG Imaging

Polarization-resolved THG imaging is based on an acquisition scheme where each line is scanned 18 times with different linear polarization states (X-P-Y mode). This approach provides a temporal resolution in the 50–100 kHz range for polarimetric measurements [55]. Polarization state switching was performed using an electro-optic modulator driven by a high-voltage amplifier (models 350-150 and 302 RM, Conoptics, CT, USA) and set between a half waveplate and a quarter waveplate. An additional quarter waveplate was placed at the entrance pupil of the objective in order to compensate for polarization distortions caused by the microscope optics and the scanning system, resulting in an ellipticity less than 0.09 for all incident linear polarizations. The X-P-Y acquisition was implemented by synchronizing galvanometer scanning and EOM voltage control using LabView-written software and an I/O board (National Instruments, USA).

## C. CARS Imaging

CARS imaging was performed on the same microscope as P-THG imaging, using the two outputs of the laser at 802 nm (pump) and 1045 nm (Stokes), temporally synchronized by a motorized delay line (ODL220/M, ThorLabs, USA). Pulse durations were ≈150 fs for the pump (spectral bandwidth ≈7 nm) and ≈200 fs for the Stokes beam (spectral bandwidth ≈5 nm). The CARS signal was emitted at ≈650 nm and detected in transmission with an appropriate emission filter (647/57, Semrock).

## D. Zebrafish Embryos Mounting and Imaging

The Tg(mbp:mCherry-CAAX) transgenic zebrafish line [6] was used to image the onset of myelination in this study. Adult male and female zebrafish with proven reproductive capacity were placed in adjacent compartments overnight and allowed to mate in the morning. The resulting eggs were placed in an incubator at 28°C until they reached the stages of development chosen for imaging (48, 52, 56, and 60 hpf; 3.0, 3.5, 4.0, 4.5, and 5.0 dpf). The chorion enveloping the youngest embryos was eventually mechanically removed to avoid chorion-induced aberration effects during multi-photon imaging. Immediately prior to imaging, five embryos were anesthetized by immersion in a 0.016% (0.16 mg/mL) solution of tricaine, prepared from a stock solution (E10521, MS-222, Sigma-Aldrich) diluted in culture water. Anesthetized embryos were then embedded in 2% low melting point agarose (UltraPure Low Melting Point Agarose, Cat. No. 16520-100, Invitrogen). Zebrafish embryos were transferred to two 35 mm glass bottom dishes, where they were oriented laterally and parallel to each other before the agarose solidified. The Petri dish was finally filled with the above-mentioned tricaine solution to maintain anesthesia during imaging. For imaging, a glass bottom dish with mounted embryos was placed under the immersion objective of the microscope at ambient temperature. The dish was rotated such that the Mauthner axons of each embryo were approximately parallel to the vertical edge of the imaged field of view (FOV) to facilitate subsequent data analysis. The axial position of the embryo was then adjusted to visualize the Mauthner axons in a segment of the spinal cord near the yolk onset, which was used as a reference region for imaging all specimens at different developmental stages. The z-position of the sample along the optical path was chosen to image the equatorial plane of the upward-facing Mauthner axon (or the core of the upper dorsal axon bundle). In both planes, 2PEF and

pTHG contrasts were acquired sequentially by two separate PMTs with a laser power of 8–20 mW (resp. 80 mW) under the objective and a 625/25 nm (resp. 377/50 nm bandpass filter for 2PEF (resp. pTHG) contrast. Typical pixel dwell times were 10 µs/polarization state and typical signal levels from Mauthner axons were 13 ± 3 detected photons at 48 hpf and 40 ± 10 detected photons at 5 dpf.

## E. Mouse Models

The care and use of mice were in accordance with institutional policies and guidelines (UPMC, INSERM, French and European Community Council Directive 86/609/EEC). The following mouse strains were used: C57bl6/J (Janvier Labs), Shiverer (gift from Dr Anne Baron Van Evercooren, Paris Brain Institute). Shiverer (Shi/Shi) is an autosomal recessive mouse mutation (in the myelin basic protein gene) that produces a shivering phenotype associated with CNS hypomyelination in affected mice [71,72]. Mice were housed under standard conditions of 12 h light/12 h dark cycles with ad libitum access to dry food and water in the animal facility of the Paris Brain Institute.

## F. Preparation and Imaging of Fixed Mouse Cortex, Cerebellum, and Spinal Cord Slices

Adult mice (P60) were perfused intracardially with PFA 4%, and the brain and spinal cord were dissected and postfixed in PFA 4% overnight at 4°C. Tissues were maintained in PBS supplemented with 0.005% sodium azide (NaN₃), embedded in agarose 4%, and sectioned with a vibratome (100 µm thick sections). Sections were stored in PBS-NaN₃ until imaging. Imaging conditions were as follows: pixel rate = 100 kHz; power after objective = 70–90 mW (resp. 20–25 mW) for THG (resp. 2PEF) imaging.

## G. Preparation of Organotypic Mouse Cerebellum Slices

*Ex vivo* sagittal mouse cerebellum slices were prepared and then maintained in an incubator until the imaging days at 8 days *in vitro*, as previously described in reference [81]. Briefly, mouse cerebella at postnatal days 8 and 10 (P8 to P10, respectively) were dissected in ice-cold Gey's balanced salt solution supplemented with 4.5 mg/mL d-glucose and penicillin-streptomycin (100 IU/mL, Thermo Fisher Scientific). The cerebella were then cut into 250 µm parasagittal slices using a McIlwain tissue chopper, and the slices were placed on Millicell membranes (3–4 slices per membrane, 2 membranes per animal; 0.4 µm Millicell, Merck Millipore). These slices were cultured in a medium consisting of 50% BME (Thermo Fisher Scientific), 25% Earle's balanced salt solution (Sigma), and 25% heat-inactivated horse serum (Thermo Fisher Scientific), supplemented with GlutaMax (2 mM, Thermo Fisher Scientific), penicillin-streptomycin (100 IU/mL, Thermo Fisher Scientific), and d-glucose (4.5 mg/mL, Sigma). Cultures were maintained at 37°C under 5% CO₂, with the medium changed every 2–3 days. Imaging was performed at 8 DIV to obtain a fully myelinated folia in the cerebellar tissue.

## H. Human Cerebellum Sections

Human tissues were obtained after informed consent (healthy donors, Brainbank Neuro-CEB, Pitié-Salpêtrière Hospital, Paris, under the ethical approval of the National Research Ethics
Committee, as snap-frozen blocks. Sections were cut using a Leica CM 1950 cryostat (20 µm thick sections), collected on Superfrost slides, and stored at -80°C until used. Human cerebellar sections were first rehydrated with PBS prior to THG imaging.

# I. Human Tissue Tiled Large-Area Acquisitions

Large-area human tissue images were obtained by acquiring 65 × 45 THG stacks of 230 × 230 pixels over 10 z-planes (Power after objective = 50 mW, pixel rate = 100 kHz, pixel size = 2.3 µm, accumulation = 1, Δz = 2 µm). After edge removal (4 pixels up and left, 15 pixels right and down), planes were selected from the z-stack images to retrieve images from the same depth in the tissue, thus compensating for the non-horizontality of the sample (less than 10 µm over 1.5 cm). Due to the large distribution of THG intensities throughout the imaged volume, two field intensity profiles were generated by averaging THG profiles from two subsets of individual THG images (subset 1: mean signal between 2 and 8 detected photons/pixel and subset 2: mean above 8 detected photons/pixel), and THG images were then corrected with one or the other filter depending on their mean. Stitching was then performed using the “Grid/collection stitching” plugin in FIJI. A mean filter (2 pixels) was applied to the mosaic. A binary filter obtained by applying a threshold (>2 counts/pixel) to a Gaussian-filtered mosaic image (sigma = 100 pixels) was applied to remove THG signals from debris around the tissue. Higher-resolution tissue images were collected from 6 × 6 THG stacks of 1140 × 1140 pixels over 3 z-planes (Power after objective = 50 mW, pixel rate = 200 kHz, pixel size = 0.46 µm, accumulation = 1, Δz = 3 µm). Similar post-processing was applied, including edge removal, field correction, followed by stitching. The images were smoothed using a median filter (2 pixels).

# J. pTHG Fourier Analysis

Briefly, $I(\phi)$ the THG intensity at each pixel as a function of polarization angle $\phi$ with $\phi \in [0 : \pi/2n : (\pi - \pi/2n)]$ (experimentally $n = 9$) is Fourier transformed using 1D FFT as a vector $\mathcal{F}(I) = [c_0, c_1, \dots, c_n, c_{-n}, \dots, c_{-1}]$ with $c_k = a_k e^{i\psi_k}$, $\psi_k \in [0 : 2\pi]$, and $a_k \in \mathbb{R}^+$. We then keep only the first two terms for the analysis: $a_0$ corresponding to the total signal and $a_1$ corresponding to a polarization dependence in $\cos(2\phi)$, from which we can calculate the pTHG modulation $M$ and the angle $\phi_1$ producing the maximum THG signal (see Fig. 2) as follows:

$$M = 4a_1/(a_1 + 2a_0),$$

$$\phi_1 = -\psi_1/2.$$

# K. Analysis of pTHG Data from Zebrafish Axons

An analysis pipeline was developed in MATLAB to extract metrics describing axon morphology from 2PEF and pTHG intensity stacks. First, pTHG stacks were corrected for excitation power variations between the different polarization states due to the dependence of optics transmission on the polarization orientation. To do this, each polarization-resolved image of a pTHG stack was corrected by the averaged intensity collected with the same polarization at the interface between the objective immersion water and a glass coverslip and normalized to the maximum value obtained for all polarizations.

Second, a region of interest (ROI) with a dimension of 30 × 50 pixels (6.9 × 11.5 µm) around the Mauthner axons and 10 × 30 pixels (2.3 × 6.9 µm) around the dorsal axon bundles was manually delineated. Potential deviation of the axons from the axis parallel to the vertical edge of the FOV (averaging axis) was corrected by realigning each line to the top line using a correlation operation. This procedure is further described in Fig. S1 of Supplement 1.

Pixel filtering of the maps was performed as in [55,59] by using a threshold on both the intensity and the $R^2$ parameter between the pTHG intensity and the first two Fourier orders corresponding to $a_0 + 2a_1 \cos(2(\phi - \phi_1))$. Typically, the value of the $R^2$ threshold was chosen as 0.6 (Fig. S2 in Supplement 1).

Scores were derived from 1D profiles of 2PEF and pTHG intensity, allowing thus for more robust quantification. These profiles were obtained by averaging 2PEF and pTHG contrasts on each line of the ROI. Semi-automatic peak detection was used to locate maxima in the 2PEF profile, and maxima (positive) and minima (negative) in the pTHG modulation profile, identified at the outer edge ($M_{out}$) and on the inner edge ($M_{in}$), respectively. Maxima and minima were evaluated on both the left and the right sides of the profiles, resulting in two positions for each extremum:

$$(x_{2\text{PEF}}^{\text{left}}, x_{2\text{PEF}}^{\text{right}}, x_{M_{\text{in}}}^{\text{left}}, x_{M_{\text{in}}}^{\text{right}}, x_{M_{\text{out}}}^{\text{left}}, x_{M_{\text{out}}}^{\text{right}})$$

2PEF and THG$_0$ intensity and pTHG modulation evaluated at these locations yielded the following scores: (i) maximum 2PEF intensity, calculated as the mean between 2PEF intensity at $x_{2\text{PEF}}^{\text{left}}$ and $x_{2\text{PEF}}^{\text{right}}$; (ii) maximum THG$_0$ intensity, calculated as a mean between THG$_0$ intensity at $x_{M_{\text{out}}}^{\text{left}}$ and $x_{M_{\text{out}}}^{\text{right}}$; and (iii) the maximum outer modulation $M_{\text{out}}$ (resp. $M_{\text{in}}$), calculated as a mean between the modulation $M_{\text{out}}$ at $x_{M_{\text{out}}}^{\text{left}}$ and $x_{M_{\text{out}}}^{\text{right}}$ (resp. at $x_{M_{\text{in}}}^{\text{left}}$ and $x_{M_{\text{in}}}^{\text{right}}$). We introduce an additional score, the “modulation gradient” $\Delta M$, as the difference between $M_{\text{out}}$ and $M_{\text{in}}$, i.e., the amplitude of the modulation variation along the radial axis.

Statistics show respectively (i) 2PEF intensity, (ii) THG$_0$ intensity, and (iii) $\Delta M$ at 7 time points: 48, 52, 56, and 60 hpf ($N=5$), at 3.0 and 3.5 dpf combined, 4.0 and 4.5 dpf combined ($N=10$), and 5.0 dpf ($N=5$). Comparison of pTHG intensities between all embryos at different developmental stages, despite the variability in thickness and transparency, was made possible by normalizing the ROI-restricted pTHG stack by the mean intensity from a 50 × 50 px (11.5 × 11.5 µm) domain chosen at the same longitudinal position in the spinal cord on the same image. Comparison of 2PEF intensities was simply corrected by dividing the intensities by the square of the laser power used for imaging.

In a similar manner, Mauthner axon diameter was measured as the distance between $x_{M_{\text{out}}}$ and $x_{M_{\text{in}}}$.

# L. Analysis of pTHG from Mouse Cerebellar Folia

THG images of cerebellar folia were acquired in mouse tissue from wild-type and *shiverer* lines at a depth of 15 µm below the tissue surface (power after the objective = 50 mW; pixel rate = 100 kHz; pixel size = 0.3 µm; accumulations = 4). Images were corrected from the field intensity profile obtained by THG images from a z-stack acquired from a homogeneous area of cerebellar tissue. For comparison between wild-type and *shiverer* tissues, pTHG analysis was performed in the dense bundle at the base of the folia. Three ROIs (pixel size = 0.23 µm) were acquired along fiber bundles of 3–5 folia from 6 to 7 fixed cerebellar slices (wild-type: 26 folia;
shiverer: 31 folia). Pixel selection was based on $R^2$ calculation between the pTHG signal and the function derived from Fourier coefficients: $a_0 + 2a_1 \cos(2(a-p1)) + 2a_2 \cos(4(a-p2)) + 2a_3 \cos(6(a-p3))$ with $R^2 = 0.5$. The illustrations in Fig. 8(a) were inspired by a BioRender (Toronto, Canada) illustration.

## M. FDTD Simulations

We performed FDTD calculations using a commercial implementation (Lumerical Device Suite, Ansys Inc, Canonsburg, Pennsylvania, USA) by elaborating on the approach described in [46,47]. As in these previous works, we neglected dispersion by holding the linear and nonlinear indices constant, and considered only polarization along *x* and *y* due to the use of diagonal tensors. We considered an incoming focused (NA = 1.0) Gaussian beam with a central wavelength of 1.2 µm and a bandwidth of 50 nm, and performed calculations over a focal region spanning 18 × 18 × 13 µm³ discretized over 20–40 nm steps with a time resolution of ≈0.5 fs. Within this region, we defined a central nonlinear region of size 6 × 6 × 10 µm³ composed of different materials with non-zero third-order nonlinear susceptibility. The THG signal of a given geometry was measured by creating a 2D array of point detectors at the end of the simulation domain and calculating the integrated intensity around 1/3 of the incident wavelength (385–415 nm). Figure S4 in Supplement 1 summarizes the simulation geometry and parameters. The FDTD calculations were performed on a Dell Precision 7920 (2× Intel 6140 CPUs, 384 GB RAM, DDR4 2666 MHz) running Lumerical version 2019b, 2020a, or 2021R2.4 with a typical computation time of ≈50 min per condition.

## N. Statistics

Error bars are standard deviations. Numbers of samples for each experiment are specified in the main text and figure legends. In Fig. 8(d), tests for normal distributions were done using the MATLAB function `kstest`, and comparisons between the two populations were done using `ttest2`.

**Funding.** Agence Nationale de la Recherche (ANR-EQPX-0029, ANR-10-INBS-04, ANR-15-CE11-0012); Fondation Bettencourt Schueller (Brain Harmonics).

**Acknowledgment.** We thank Júlia Ferrer Ortas, Bahar Asadipour, Matthew Gleeson, and Marie-Claire Schanne-Klein for discussions on THG and polarized multiphoton microscopy. We thank Delphine Débarre for discussions on THG and aberrations. We thank Laurent Bourdieu for discussions on myelin imaging. We thank David Lyons for providing us with the Tg(mbp:mCherry-CAAX) zebrafish line.

J. Morizet, N. Olivier, C. Stringari, and E. Beaurepaire designed the study. J. Morizet acquired and analyzed the imaging data. N. Olivier performed the numerical simulations. P. Mahou designed and optimized the multiphoton microscope. M. S. Aigrot, E. Martin, A. Desmazieres, and B. Stankoff designed the mouse and human tissue studies and provided the samples. J. Morizet, N. Olivier, C. Stringari, and E. Beaurepaire wrote the manuscript with input from all other authors.

**Disclosures.** The authors declare no conflicts of interest.

**Data availability.** Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request. The scripts used to run the simulation are available in Ref. [82].

**Supplemental document.** See Supplement 1 for supporting content.

## REFERENCES

1. C. Stadelmann, S. Timmler, A. Barrantes-Freer, et al., “Myelin in the central nervous system: structure, function, and pathology,” *Physiol. Rev.* **99**, 1381–1431. (2019).

2. J. E. Rinhol and L. H. Bergersen, “Neuroscience: the wrap that feeds neurons,” *Nature* **487**, 435–436 (2012).

3. K. Nave, “Myelination and the trophic support of long axons,” *Nat. Rev. Neurosci.* **11**, 275–283 (2010).

4. G. Tomassy, D. Berger, H. Chen, et al., “Distinct profiles of myelin distribution along single axons of pyramidal neurons in the neocortex,” *Science* **344**, 319–324 (2014).

5. E. Gibson, D. Purger, C. Mount, et al., “Neuronal activity promotes oligodendrogenesis and adaptive myelination in the mammalian brain,” *Science* **344**, 1252–304 (2014).

6. S. Mensch, M. Baraban, R. Almeida, et al., “Synaptic vesicle release regulates myelin sheath number of individual oligodendrocytes in vivo,” *Nat. Neurosci.* **18**, 628–630 (2015).

7. C. Sampaio-Baptista and H. Johansen-Berg, “Plasticity in the adult brain,” *Neuron* **96**, 1239–1251. (2017).

8. A. Saab and K. Nave, “Myelin dynamics: protecting and shaping neuronal functions,” *Curr. Opin. Neurobiol.* **47**, 104–112 (2017).

9. M. Swire and C. Ffrench-Constant, “Seeing is believing: myelin dynamics in the adult CNS,” *Neuron* **98**, 684–686 (2018).

10. L. Osso and E. Hughes, “Dynamics of mature myelin,” *Nat. Neurosci.* **27**, 1449–1461 (2024).

11. W. Rushton, “A theory of the effects of fibre size in medullated nerve,” *J. Physiol.* **115**, 101–122 (1951).

12. T. Chomiak and B. Hu, “What is the optimal value of the g-ratio for myelinated fibers in the rat CNS? A theoretical approach,” *PLoS One* **4**, e7754 (2009).

13. K. J. Smith, W. F. Blakemore, and W. I. McDonald, “Central remyelination restores secure conduction,” *Nature* **280**, 395–396 (1979).

14. N. Moll, A. Rietsch, S. Thomas, et al., “Multiple sclerosis normal-appearing white matter: pathology-imaging correlations,” *Ann. Neurol.* **70**, 764–773. (2011).

15. A. Petiet, M. Aigrot, and B. Stankoff, “Gray and white matter demyelination and remyelination detected with multimodal quantitative MRI analysis at 11.7T in a chronic mouse model of multiple sclerosis,” *Front. Neurosci.* **27**, 491 (2016).

16. W. Rushton, “Towards in vivo g-ratio mapping using MRI: unifying myelin and diffusion imaging,” *J. Neurosci. Methods* **348**, 108990 (2021).

17. M. Filippi, M. Rocca, N. De Stefano, et al., “Magnetic resonance techniques in multiple sclerosis: the present and the future,” *Arch. Neurol.* **68**, 1514–1520. (2011).

18. C. Louapre, “Conventional and advanced MRI in multiple sclerosis,” *Rev. Neurol.* **174**, 391–397. (2018).

19. E. Romanelli, C. Sorbara, I. Nikić, et al., “Cellular, subcellular and functional in vivo labeling of the spinal cord using vital dyes,” *Nat. Protoc.* **8**, 481–490. (2013).

20. N. Blanke, A. J. Gray, R. E. Robinson, et al., “Practical considerations for birefringence microscopy of myelin structure: microscope design and tissue processing for effective imaging,” *Imaging Neurosci.* **2**, 1–22 (2024).

21. J. B. Arous, J. Binding, J.-F. Léger, et al., “Single myelin fiber imaging in living rodents without labeling by deep optical coherence microscopy,” *J. Biomed. Opt.* **16**, 119802 (2011).

22. N. Blanke, S. Chang, A. Novoseltseva, et al., “Multiscale label-free imaging of myelin in human brain tissue with polarization-sensitive optical coherence tomography and birefringence microscopy,” *Biomed. Opt. Express* **14**, 5946–5964 (2023).

23. A. J. Schain, R. A. Hill, and J. Grutzendler, “Label-free in vivo imaging of myelinated axons in health and disease with spectral confocal reflectance microscopy,” *Nat. Med.* **20**, 443–449 (2014).

24. H. Wang, Y. Fu, P. Zickmund, et al., “Coherent anti-stokes Raman scattering imaging of axonal myelin in live spinal tissues,” *Biophys. J.* **89**, 581–591 (2005).

25. M. J. Farrar, F. W. Wise, J. R. Fetcho, et al., “In vivo imaging of myelin in the vertebrate central nervous system using third harmonic generation microscopy,” *Biophys. J.* **100**, 1362–1371 (2011).

26. R. W. Boyd, A. L. Gaeta, and E. Giese, “Nonlinear optics,” in *Springer Handbook of Atomic, Molecular, and Optical Physics* (Springer, 2008), pp. 1097–1110.
27. Y. Barad, H. Eisenberg, M. Horowitz, et al., "Nonlinear scanning laser microscopy by third harmonic generation," Appl. Phys. Lett. 70, 922-924 (1997).

28. J.-X. Cheng and X. S. Xie, "Green's function formulation for third-harmonic generation microscopy," J. Opt. Soc. Am. B 19, 1604-1610 (2002).

29. D. Débarre, W. Supatto, A.-M. Pena, et al., "Imaging lipid bodies in cells and tissues using third-harmonic generation microscopy," Nat. Methods 3, 47-53 (2006).

30. D. Débarre, N. Olivier, and E. Beaurepaire, "Signal epidetection in third-harmonic generation microscopy of turbid media," Opt. Express 15, 8913-8924 (2007).

31. C.-K. Sun, S.-W. Chu, S.-Y. Chen, et al., "Higher harmonic generation microscopy for developmental biology," J. Struct. Biol. 147, 19-30 (2004).

32. R. Genthial, E. Beaurepaire, M.-C. Schanne-Klein, et al., "Label-free imaging of bone multiscale porosity and interfaces using third-harmonic generation microscopy," Sci. Rep. 7, 3419 (2017).

33. B. Weigelin, G.-J. Bakker, and P. Friedl, "Third harmonic generation microscopy of cells and tissue organization," J. Cell Sci. 129, 245-255 (2016).

34. N. Olivier, M. A. Luengo-Oroz, L. Duloquin, et al., "Cell lineage reconstruction of early zebrafish embryos using label-free nonlinear microscopy," Science 329, 967-971 (2010).

35. F. Aptel, N. Olivier, A. Deniset-Besseau, et al., "Multimodal nonlinear imaging of the human cornea," Invest. Ophthalmol. Vis. Sci. 51, 2459-2465 (2010).

36. S. You, H. Tu, E. J. Chaney, et al., "Intravital imaging by simultaneous label-free autofluorescence-multiharmonic microscopy," Nat. Commun. 9, 2125 (2018).

37. A. Rakhymzhan, A. F. Fiedler, R. Günther, et al., "Optimized intravital three-photon imaging of intact mouse tibia links plasma cell motility to functional states," iScience 27, 110985 (2024).

38. N. Akbari, R. Tatarsky, K. Kolkman, et al., "Label-free, whole-brain in vivo mapping in an adult vertebrate with third harmonic generation microscopy," J. Comp. Neurol. 532, e25614 (2024).

39. S. Witte, A. Negrean, J. C. Lodder, et al., "Label-free live brain imaging and targeted patching with third-harmonic generation microscopy," Proc. Natl. Acad. Sci. USA 108, 5970-5975 (2011).

40. H. Lim, D. Sharoukhov, I. Kassim, et al., "Label-free imaging of Schwann cell myelination by third harmonic generation microscopy," Proc. Natl. Acad. Sci. USA 111, 18025-18030 (2014).

41. M. Redlich and H. Lim, "A method to measure myeloarchitecture of the murine cerebral cortex in vivo and ex vivo by intrinsic third-harmonic generation," Front. Neuroanat. 13, 65 (2019).

42. P.-J. Wu, H.-C. Tseng, C.-C. Chao, et al., "Discontinuity third harmonic generation microscopy for label-free imaging and quantification of intraepidermal nerve fibers," Cell Rep. Methods 4, 100735 (2024).

43. H. Hajjar, H. Boukhaddaoui, A. Rizgui, et al., "Label-free nonlinear microscopy to measure myelin outcome in a rodent model of Charcot-Marie-Tooth diseases," J. Biophotonics 11, e201800186 (2018).

44. N. E. Chernavsky, N. Hong, M. Lamont, et al., "Label-free tracking of myelin dynamics in subcortical white matter of a mouse model of multiple sclerosis using third harmonic generation microscopy," in Optics and the Brain (Optica Publishing Group, 2024), paper BM1C-3.

45. N. R. Meijns, M. Blokker, S. Idema, et al., "Dynamic imaging of myelin pathology in physiologically preserved human brain tissue using third harmonic generation microscopy," PLoS One 20, e0310663 (2025).

46. J. Morizet, G. Sartorello, N. Dray, et al., "Modeling nonlinear microscopy near index-mismatched interfaces," Optics 8, 944-951 (2021).

47. J. Morizet, N. Olivier, P. Mahou, et al., "Third harmonic imaging contrast from tubular structures in the presence of index discontinuity," Sci. Rep. 13, 7850 (2023).

48. E. Klebs, "Die nerven der organischen muskelfasern," Virchows Arch. 32, 168-198 (1865).

49. G. F. Göthlin, Die doppelbrechenden Eigenschaften des Nervengewebes: ihre Ursachen und ihre biologischen Konsequenzen (Almqvist & Wiksells Boktryckeri, 1913).

50. W. J. Schmidt, Die Bausteine des Tierkörpers in polarisiertem Lichte (F. Cohen, 1924).

51. B. de Campos Vidal, M. L. S. Mello, A. C. Caseiro-Filho, et al., "Anisotropic properties of the myelin sheath," Acta Histochem. 66, 32-39 (1980).

52. F. O. Schmitt, R. S. Bear, and G. L. Clark, "X-ray diffraction studies on nerve," Radiology 25, 131-151 (1935).

53. W. Schmidt, "Doppelbrechung und feinbau der markscheide der nervenfasern," Z. Zellforsch. mikrosk. Anat. 23, 657-676 (1936).

54. D. Oron, E. Tal, and Y. Silberberg, "Depth-resolved multiphoton polarization microscopy by third-harmonic generation," Opt. Lett. 28, 2315-2317 (2003).

55. J. Morizet, G. Ducourthial, W. Supatto, et al., "High-speed polarization-resolved third-harmonic microscopy," Optica 6, 385-388 (2019).

56. C.-F. Chang, C.-H. Yu, and C.-K. Sun, "Multi-photon resonance enhancement of third harmonic generation in human oxyhemoglobin and deoxyhemoglobin," J. Biophotonics 3, 678-685 (2010).

57. J. Ferrer Ortas, P. Mahou, S. Escot, et al., "Label-free imaging of red blood cells and oxygenation with color third-order sum-frequency generation microscopy," Light 12, 29 (2023).

58. N. Spassky, C. Olivier, I. Cobos, et al., "The early steps of oligodendrogenesis: insights from the study of the plp lineage in the brain of chicks and rodents," Dev. Neurosci. 23, 318-326 (2001).

59. M. Gleeson, J. Morizet, P. Mahou, et al., "Kidney stone classification using multimodal multiphoton microscopy," ACS Photonics 10, 3594-3604 (2023).

60. C. Buckley, A. Marguerie, W. Alderton, et al., "Temporal dynamics of myelination in the zebrafish spinal cord," Glia 58, 802-812. (2010).

61. M. Preston and W. Macklin, "Zebrafish as a model to investigate CNS myelination," Glia 63, 177-193 (2015).

62. T. Czopka, "Insights into mechanisms of central nervous system myelination using zebrafish," Glia 64, 333-349. (2016).

63. N. Olivier and E. Beaurepaire, "Third-harmonic generation microscopy with focus-engineered beams: a numerical study," Opt. Express 16, 14703-14715 (2008).

64. L. Novotny, Principles of Nano-Optics (Cambridge University Press, 2003).

65. M. Menzel, M. Axer, H. De Raedt, et al., "Toward a high-resolution reconstruction of 3D nerve fiber architectures and crossings in the brain using light scattering measurements and finite-difference time-domain simulations," Phys. Rev. X 10, 021002 (2020).

66. D. Débarre and E. Beaurepaire, "Quantitative characterization of biological liquids for third-harmonic generation microscopy," Biophys. J. 92, 603-612 (2007).

67. D. Hildebrand, M. Cicconet, R. Torres, et al., "Whole-brain serial-section electron microscopy in larval zebrafish," Nature 545, 345-349 (2017).

68. R. Turcotte, D. J. Rutledge, E. Bélanger, et al., "Intravital assessment of myelin molecular order with polarimetric multiphoton microscopy," Sci. Rep. 6, 31685 (2016).

69. T. Xie, A. Novoseltseva, N. Blanke, et al., "Birefringence microscopy for quantifying birefringence of myelin in thin brain section," in Microscopy Histopathology and Analytics (Optica Publishing Group, 2024), paper JS4A-42.

70. M. Zimmerley, P. Mahou, D. Débarre, et al., "Probing ordered lipid assemblies with polarized third-harmonic-generation microscopy," Phys. Rev. X 3, 011002 (2013).

71. H. Werner, M. Jung, M. Klugmann, et al., "Mouse models of myelin diseases," Brain Pathol. 8, 771-793 (1998).

72. T. Nomura, Y. Bando, H. Bochimoto, et al., "Three-dimensional ultrastructures of myelin and the axons in the spinal cord: application of SEM with the osmium maceration method to the central nervous system in two mouse models," Neurosci. Res. 75, 190-197 (2013).

73. D. Débarre, W. Supatto, and E. Beaurepaire, "Structure sensitivity in third-harmonic generation microscopy," Opt. Lett. 30, 2134-2136 (2005).

74. N. G. Horton, K. Wang, D. Kobat, et al., "In vivo three-photon microscopy of subcortical structures within an intact mouse brain," Nat. Photonics 7, 205-209 (2013).

75. K. Guesmi, L. Abdeladim, S. Tozer, et al., "Dual-color deep-tissue three-photon microscopy with a multiband infrared laser," Light 7, 12 (2018).

76. T. Wang, D. Ouzounov, C. Wu, et al., "Three-photon imaging of mouse brain structure and function through the intact skull," Nat. Methods 15, 789-792 (2018).

77. L. Streich, J. Boffi, L. Wang, et al., "High-resolution structural and functional deep brain imaging using adaptive optics three-photon microscopy," Nat. Methods 18, 1253-1258 (2021).

78. M. Gleeson, S. Escot, X. Solinas, et al., "Multiplexed orthogonal polarizations harmonic imaging for probing dynamic biological processes," ACS Photonics 11, 3259-3266 (2024).
79. A. T. Mok, T. Wang, S. Zhao, et al., "A large field-of-view, single-cell-resolution two-and three-photon microscope for deep imaging," eLight 4, 20 (2024).

80. C. Rodríguez, A. Chen, J. A. Rivera, et al., "An adaptive optics module for deep tissue multiphoton imaging in vivo," Nat. Methods 18, 1259-1264 (2021).

81. M. Thetiot, R. Ronzano, M. Aigrot, et al., "Preparation and immunostaining of myelinating organotypic cerebellar slice cultures," J. Vis. Exp. 145, e59163 (2019).

82. J. Morizet, N. Olivier, M. S. Aigrot, et al. "Multiscale characterization of myelin distribution with polarized THG microscopy," GitHub (2025), https://github.com/LaboratoryOpticsBiosciences/FDTD_Neurons.