Probing elastic anisotropy of human skin in vivo with light using non-contact acoustic micro-tapping OCE and polarization sensitive OCT

Skin broadly protects the human body from undesired factors such as ultraviolet radiation and abrasion and helps conserve body temperature and hydration. Skin’s elasticity and its level of anisotropy are key to its aesthetics and function. Currently, however, treatment success is often speculative and subjective, and is rarely based on skin’s elastic properties because there is no fast and accurate non-contact method for imaging of skin’s elasticity. Here we report on a non-contact and non-invasive method to image and characterize skin’s elastic anisotropy. It combines acoustic micro-tapping optical coherence elastography (AμT-OCE) with a nearly incompressible transversely isotropic (NITI) model to quantify skin’s elastic moduli. In addition, skin sites were imaged with polarization sensitive optical coherence tomography (PS-OCT) to help define fiber orientation. Forearm skin areas were investigated in five volunteers. Results clearly demonstrate elastic anisotropy of skin in all subjects. AμT-OCE has distinct advantages over competitive techniques because it provides objective, quantitative characterization of skin’s elasticity without contact, which opens the door for broad translation into clinical use. Finally, we demonstrate that a combination of multiple OCT modalities (structural OCT, OCT angiography, PS-OCT and AμT-OCE) may provide rich information about skin and can be used to characterize scar.


Supplementary Note 1. Elastic constants of a linear elastic, nearly incompressible, transversely isotropic (NITI) material
The elastic modulus matrix of a linear-elastic, transversely isotropic (TI) material can be represented in a form 1 where 12 = 11 − 2 66 due to symmetry conditions 1,2 . Considering that two of three principal planes in a TI material are isotropic, we will use a modified notation of moduli based on the notation common of isotropic materials 3 : An additional modulus shows that shear deformation can be different if shear stress is applied along the symmetry axis compared to that applied across it. Parameters 1 and 2 indicate that differences between 33 and 11 , and between 12 and 13 may be nonzero 3 .
Using the elasticity matrix (S.2), the stress-strain relation in a TI material for small deformations is described by Hook's law and takes the form (in Voigt notation) 3 where and are the components of the stress and strain tensors, respectively.
For any arbitrary load applied to a TI medium, the strain field is computed by inverting the stressstrain relation, where the inverse of the elastic modulus matrix is the compliance matrix = −1 : where the components of the compliance matrix can be expressed in the form: where E and E are the transversal and longitudinal Young's moduli, respectively, and and are Poisson's ratios, respectively. These parameters can be expressed through the components of the compliance matrix as 3 : where We also note also that = .
The incompressibility condition appropriate for materials such as soft tissue constrains the strain such that For an accurate description of an incompressible TI material, this condition must be applied directly to the stress strain relations in the limit that θ approaches zero and the longitudinal modulus λ approaches infinity, i.e.

= lim
where is defined as the scalar pressure. As demonstrated in Ref 3 , a strict application of the incompressibility condition produces stress-strain relations in which the longitudinal terms include the scalar pressure P and inversion of this equation to compute the Young's moduli and Poisson's ratios of an incompressible TI material is not strictly defined. However, following an approach presented in Ref 3 for an isotropic incompressible material, the Young's moduli and Poisson's ratios can be approximated using the formulas presented above and taking the limit as λ approaches infinity. Under these conditions, the Young's moduli and Poisson's ratios for a nearly-incompressible TI material, i.e. NITI material, can be approximated as: where = 2 − 2 1 .
We will further differentiate a fast-axis NITI material (i.e. the material with a unidirectional fiber orientation and, thereby, the symmetry axis collinear with fibers) considered in this paper for skin, from a slow-axis NITI material (i.e. the material where fibers arbitrary distributed in a plane and, thereby, the symmetry axis is perpendicular to the fiber direction) previously introduced for cornea 4 . Supplementary Figure 1 shows how Young's moduli and Poisson's ratios can change with . Clearly, has a very narrow range, with a lower limit of 2 and upper limit of 4 as → ∞, and is a linear function of .
A lower limit for the argument / is defined by the restriction for Poisson's ratio to be positive.
The case when / < 0 corresponds to the slow-axis NITI medium, when 33 < 11 . This situation is realized in cornea, when lamellae are oriented randomly in-plane of the cornea, and, therefore, the symmetry Z-axis corresponds to the direction perpendicular to fibers.
The case when / > 0 corresponds to the fast-axis NITI medium, when 33 > 11 . This situation is typical for skin and muscle, where fibers can be considered locally unidirectional, and, therefore, the symmetry Z-axis corresponds to the direction along the fibers.

Supplementary Note 2. Bulk waves in a linear elastic, fast-axis NITI material
An isotropic solid supports two unique waves -one longitudinal with speed and one shear with speed . Biological tissue is nearly incompressible, meaning that the longitudinal wave speed (typically ≈ 1540 m/s) is much larger than the shear wave speed (typically ≈ 1-10 m/s).
In biological tissue such as the skin, the natural orientation of collagen fibers in the dermis are traditionally described using Langer's lines and correspond with estimates of mechanical anisotropy. Thus, in skin, a single shear modulus is not sufficient to describe elastic behavior. The simplest constitutive model that may account for skin's structure is a NITI model, recently introduced for describing elastic anisotropy in cornea 4 . In cornea, however, the orientation of fibers is in plane but quite random. This creates a symmetry axis perpendicular to the fiber structure. Thus, the symmetry axis for cornea is a slow axis. In contrast, the orientation of collagen fibers in skin is locally unidirectional, which creates a symmetry axis (which is the fast axis) in the plane of skin.
First, we consider propagation of bulk acoustic waves in a fast-axis NITI medium in the YZ plane (see Fig S2a below), parallel to the skin surface, but within the skin volume. Wave speeds can be found from the solution of the Christoffel equation 2 : where = is the Christoffel tensor; , ( , = 1, 2, 3) are the directional cosines of the propagation vector ; is the Kronecker operator; is the wave polarization unit vector, which is also the eigenvector of tensor ; and is the velocity of the propagating mechanical wave.
In the YZ plane, the unit vector for directional cosines takes the form , = [0, sin , cos ]. In this plane, the solution points to three bulk waves: a quasi-longitudinal, a quasi-shear ( ), and a shear wave, that propagate with speeds , and , respectively 5 : where Velocities or can be calculated using the positive or negative discriminants, respectively. Because ≫ , , in nearly incompressible media, variations of quasi-longitudinal wave speed are negligible over all angles and can be approximately set as = � / , and the speed of quasi-shear wave can be expressed as: The shear wave is always polarized along the X-axis. It can also be shown, that in the nearincompressible limit, the polarization of a quasi-shear wave is orthogonal to its propagation direction and lies in plane (YZ plane). Clearly, the shear wave speed (black dotted line) monotonically decreases from � / along the Z direction (e.g., along the fibers in skin) to � / along the Y direction (e.g., across the fibers).
Propagation of the quasi-shear wave (solid lines) is more complicated. Along the fibers, it is converted to a pure shear wave with a speed of � / . It is also converted to a pure shear wave propagating across the fibers (at = 90°). Because this wave is polarized in the YZ plane, the propagation across fibers corresponds to its polarization along fibers and vice versa. Both situations result in the same propagation speed � / . It is interesting that at = 45°, the propagation speed does not depend on and approaches � / as → 0.

Supplementary Note 3: Rayleigh wave Propagation in a linear elastic, fast-axis NITI material
A well-defined surface imposes additional constraints on propagating waves so that traction forces should be taken into account. A real root of the wave motion equation corresponds to a Rayleigh wave with amplitude decreasing exponentially with depth. In nearly incompressible isotropic materials, its speed is directly related to , ≈ 0.9553 . Propagation velocities of shear and Rayleigh waves can then be converted into the Young's modulus = 3 2 . Now consider wave propagation in the same YZ plane, but over the surface of a fast-axis NITI medium. In addition to the unit propagation vector = [0, sin θ, cos θ], a surface unit normal vector = [1, 0, 0] should be introduced. The solution for the Rayleigh wave speed can be found using the Stroh formalism with 3x3 matrices 6 : = , (S17) = . (S18) Matrices (S16) -(S18) can be combined to form the 6x6 Stroh eigenvalue problem: where is the wave number and the matrix N Eigenvectors contain polarization and traction vectors of harmonic waves that satisfy the free surface boundary conditions and whose displacement and traction are represented in the form: The amplitude of a Rayleigh wave must decay with depth, and therefore has a positive imaginary part. Six eigenvalues occur in complex conjugate pairs � �, but only three depthdecaying modes form Rayleigh waves. In addition, the free surface boundary condition implies that the linear combination of tractions must equal zero: which can be rewritten as: with a nontrivial solution if and only if the determinant of is zero.
To solve for the Rayleigh wave speed, an iterative algorithm is used. At each iteration, a trial wave speed is considered, and the Stroh eigenvalue problem is solved to obtain the relevant eigenvectors, as described above. The traction vectors are extracted and used to form the matrix and calculate its determinant. This process is repeated until the absolute value of the determinant is minimized. The resulting wave speed corresponds to the Rayleigh wave speed .
The simplest scenario occurs when a Rayleigh wave propagates across fibers in the fast-axis NITI medium. Its speed has the same relationship with the shear wave speed as in an isotropic material, and it is fully defined by modulus without any dependence on and .
Obtaining an analytic solution for the Rayleigh wave speed at an arbitrary angle to the symmetry axis Z (fiber direction for skin) and its dependence on / and / is complicated and possible only for a few cases. It can be shown (see Ref 7 ) that that the velocity of Rayleigh waves along fibers obeys the secular equation below, similar to that for the isotropic case 8 , but with different coefficients: where = 2 / , Equation (S25) can also be written in an alternative form: As noted in Ref 5 , the speed of the quasi-shear wave is not limited by � / and can change broadly with ; however, the speed of a Rayleigh wave in a NITI medium cannot exceed � / . This limit is reached when / → ∞.
Note that in the isotropic limit, the equation reduces to the well-known Rayleigh equation 7 : Supplementary Figure S3 shows the behavior of the Rayleigh wave speed in the fast-axis NITI medium (typical for skin and muscles) at different values of / . Clearly, the absolute limitation on the Rayleigh wave speed is the speed of bulk shear waves polarized perpendicular to the surface (along the X-axis).

Supplementary Note 4: Influence of a thin isotropic layer (epidermis) on Rayleigh wave propagation in a fast-axis NITI material
Skin is a multilayered structure consisting of epidermis, dermis, and hypodermis layers. Here we investigate skin's elastic properties assuming that they are dominated by the dermis and that the dermis represents a NITI medium. The epidermis in skin is much thinner than the dermis and, therefore, cannot strongly affect the skin's elastic behavior and its deformation. However, adding a layer on top of a NITI material might affect the propagation speed of a Rayleigh wave depending on the relationship between its wavelength and the layer thickness, i.e. on the parameter ℎ (ℎ is the layer thickness and is 2π times the wavelength). Indeed, when the layer thickness is much smaller than the wavelength, ℎ ≪ 1, the speed of a Rayleigh wave is defined by the speed in the substrate, i.e., the layer can be ignored 8 . In contrast for a very thick layer, ℎ ≫ 1, the Rayleigh wave speed is solely determined by the properties of the layer 8, 10 .
In our AµT-OCE measurements in skin, the maximum frequency of Rayleigh waves is usually limited to about 3 kHz, which corresponds to about 1 mm wavelength in dermis. It is evident that the epidermal thickness around 50 µm to 100 µm is much smaller that the wavelength. Thus, intuitively, the epidermis should not affect the Rayleigh wave speed, but considering that the difference in Young's moduli between epidermal and dermal layers can be dramatic, it is worth checking this effect quantitatively.
The influence of medium loading with a thin layer on the propagation characteristics of Rayleigh waves was investigated in multiple past studies, mostly for solids. In the majority of cases, a softer coating was considered 8 . Studies for a stiffer coating than the substrate are limited 10 . However, Rayleigh waves in a layered soft media, to our knowledge, were not quantitatively described.
For skin the substrate (dermis) is an anisotropic NITI medium, whereas the epidermis was shown to be isotropic (at least optically) in multiple studies (see Ref 11 , for example). Considering this structure analytically is not a simple task; therefore, numerical simulations were performed in OnScale instead.
We have used OnScale to build numerical AµT-OCE models in several numerical studies before, including simulating nearly incompressible anisotropic media 4 . Details on the simulation parameters for nearly incompressible media and numerical model testing can be found in Refs 4,11 . The simulation script for the model used here is shared in Supplementary Software Library. Supplementary Figures S4a,b present diagrams of wave propagation along and across fibers. Examples of Rayleigh wavefields, and corresponding 2-D spectra, are shown in Supplementary  Figures S4c-f.
Because epidermis provides mainly protection, we set its shear modulus to be = 500 kPa, which is much larger than that in dermis. Note that the Young's modulus in epidermis may be different in different body sites, may vary with age, gender, and an individual's lifestyle and environment, and many other conditions. As seen in Suppl. Figure S4c-f, even if the modulus in the epidermis is 50 time larger than that for dermis, it does not change the Rayleigh wave speed appreciably. Thus, we can conclude that ignoring the epidermis in reconstructing elastic moduli in skin is justified for the case of ℎ ≪ 1. Figure S4: Diagram of Rayleigh wave propagation over the surface of a 2-layer (epidermis/dermis) system along (a) and across (b) the symmetry axis of the second layer. Epidermis is modeled as a 50 m thick isotropic material with = 500 kPa; dermis is modeled as a fast-axis NITI hemispace with = = 10 kPa, = 30 kPa. Simulated wavefields of vertical vibration velocity and 2D Fourier spectra of Rayleigh waves are shown accordingly for propagation along (c), (e) and across (d), (f) the symmetry axis of the NITI medium. Solid lines in panels (c)-(f) correspond to Rayleigh wave speeds in a bulk NITI medium, and dashed lines correspond to fitting. The color bar is linear in panels (c) and (d) and shown in arbitrary units. The color bar for spectral amplitudes in panels (e) and (f) is shown in dB relative to the maximum signal.

Supplementary Note 5: Guided waves in a bounded (with hypodermis) NITI layer (epidermis).
As demonstrated in Supplementary Note 4 above, a thin epidermal layer has little influence on the propagation of Rayleigh waves in dermis due to its small thickness compared to the wavelength, i.e. under the condition ℎ ≪ 1. However, the thickness of dermis itself is not infinite and usually on the order of 1 mm, representing a ℎ ≈ 1. For this case, guided wave behavior can be expected. Below, we investigate how a mismatch between elastic properties of dermis and hypodermis influence mechanical wave propagation over the skin surface.
Because the hypodermis (or subcutaneous tissue) is mainly connective tissue with a large fraction of fat and almost no fibers, we assume subcutaneous tissue to be isotropic and having a lower Young's modulus than that in dermis due to the very low elastic modulus of fat. In simulations, we considered it to be = 2 kPa. The hypodermis was simulated as a hemispace bounded with the dermis on top (see Suppl. Figure S5a,b).
Since the dermal layer is anisotropic (assumed as a fast-axis NITI medium), we explored wave propagation along and across fibers. The epidermis was assumed to be 50 m thick and the dermis was assumed to be 1 mm thick. The elastic properties of epidermis and dermis are listed in the caption to Suppl. Figure S5 and are the same as those used in Suppl. Figure S4.
As seen in Suppl. Figure S5c,d, bounding of the dermal layer on its bottom modifies the wavefields for Rayleigh wave propagation both along and across the symmetry axis (fiber direction). The wavefield also changes for wave propagation along other directions in the YZ plane. Fourier analysis shows evidence of guided waves when propagation over multiple wavelengths is used for analysis. However, at small distances from the AµT source, as shown by the fit lines in Figure  S5c,d, the wavefields closely resemble the unbounded case, strongly suggesting that group velocity measurements represent Rayleigh wave speeds. (epidermis/dermis/subcutaneous tissue) system along (a) and across (b) the symmetry axis of the second layer. Epidermis is modeled as a 50 m thick isotropic material with = 500 kPa; dermis is modeled as a 1 mm thick fast-axis NITI medium with = = 10 kPa, = 30 kPa; and subcutaneous tissue is modeled as an isotropic hemispace with = 2 kPa. Simulated OnScale wavefields and 2D Fourier spectra of Rayleigh waves are shown accordingly for propagation along (c), (e) and across (d), (f) the symmetry axis of the NITI medium. Solid lines in panels (c)-(f) correspond to Rayleigh wave speeds in a bulk NITI medium, and dashed lines correspond to fitting. The color bar is linear in panels (c) and (d) and shown in arbitrary units. The color bar for spectral amplitudes in panels (e) and (f) is shown in dB relative to the maximum signal.

Supplementary Note 6: Experimental study of guided waves in ex vivo chicken skin.
AµT-OCE experiments were performed in ex vivo chicken samples to investigate the influence of subcutaneous tissue on the propagation of mechanical waves in skin. Chicken drumsticks were bought in a grocery store.
In the first set of experiments, whole chicken drumsticks (Suppl. Fig. S6a) were explored. Waveforms were recorded at different propagation directions for several positions in the drumstick. Ten different waveforms of surface-propagating mechanical waves were recorded. A typical wavefield for the whole drumstick is represented in Suppl. Fig. S6d. Clearly, guided behavior is not seen in this case, and the wavefield can be fit with a linear function, indicating nondispersive wave propagation. The slope of the wavefield corresponds to the group velocity of the propagating wave, which equals the Rayleigh wave speed when dispersion is not observed. The mean propagating wave speed calculated by the ensemble of 10 measurements was = (7 ± 1) m/s corresponding to the shear modulus = (54 ± 15) kPa, assuming chicken skin isotropy. The calculated standard deviation corresponds to the statistical variation of the propagation wave speed that can be due to sample inhomogeneity or some anisotropy.
In the second set of experiments, skin was removed from muscle (see Suppl. Fig. S6b) and placed on top of a container filled with water (see Suppl. Fig. S6c). In this case, skin had an interface with a medium having zero shear modulus, which should produce guided waves in skin with (see Suppl. Fig. S6e) a strong geometrical dispersion of wave velocity (see Suppl. Fig. S6f) typical for Lamb/Scholte waves 8 . Fitting the dispersion curve with the analytical expression for a Scholte wave propagating in the layer bounded with air on one side and with water on the other 12, 13 , gives a modulus estimate = (53 ± 2) kPa. Because estimation of the mean is performed by fitting over a wide frequency range, its accuracy is much higher than that obtained from the mean group velocity.
The results of this study strongly support the assumption that when skin is bounded with a tissue having a similar shear modulus, guided behavior is not observed. Hence, wave propagation over the skin surface may be considered as a Rayleigh wave propagating over the surface of a bulk material with moduli corresponding to that of skin. A detailed investigation of the effect of subcutaneous tissue on wave propagation in skin is outside of the scope of this study. Definitely, when skin is thin and located on top of bones, or in close proximity to blood vessels, guided propagation is much more possible.