Assessment of nonalcoholic fatty liver function by photoacoustic imaging

Abstract. Significance The assessment of liver function plays an important role in the diagnosis of nonalcoholic fatty liver disease (NAFLD). Current noninvasive imaging methods have limited applicability in this regard. Aim We report an application of multispectral photoacoustic imaging (PAI), an emerging modality, to visualize lipid accumulation and liver function in NAFLD. Approach We first demonstrated the liver function reserve with indocyanine green (ICG) to verify the organ’s dysfunction due to NAFLD in a rabbit model. We then noninvasively quantified lipid content in the liver using multispectral PAI. The in vivo PAI results were compared and verified with photoacoustic ex vivo images and liver biopsy. Results A significant difference in the lipidmean value was observed [lipidmean  =  0.081  ±  0.0161 arbitrary units (a.u.) control versus NAFLD 0.198  ±  0.048 a.u., P  =  0.003]. Similar to in vivo analysis, a significant difference in lipidmean was observed (lipidmean  =  0.0673  ±  0.0165 versus 0.486  ±  0.073 a.u., P  <  0.0001) between control and NAFLD group ex vivo. For liver function, the control group showed a rapid decrease after the peak point, whereas the elimination of ICG for the NAFLD group was slower. Conclusions Our study shows that PAI has the potential to provide a noninvasive biomarker for the assessment of liver function and lipid accumulation for NAFLD diagnosis and treatment.


Introduction
Nonalcoholic fatty liver disease (NAFLD) is characterized by excessive hepatic fat accumulation associated with an increasing rate of obesity, diabetes, and metabolic syndrome. The prevalence of NAFLD is ∼25% globally, afflicting about 1 billion individuals worldwide, which is becoming the predominant cause of chronic liver disease and a substantial clinical and economic burden around the world. [1][2][3] In the context of NAFLD, the first diagnostic requirement is to show the presence of fat in the liver accurately. The presence of fat >5% to 10% is considered abnormal. 1,4 Liver biopsy is the gold standard for the diagnosis and classification of NAFLD. However, it suffers from being invasive, expensive, associated with potential complications and sampling error, and plagued with interobserver variability. [5][6][7] Magnetic resonance spectroscopy or magnetic resonance imaging (MRI) quantifies the proton density fat fraction (PDFF) to measure hepatic triglyceride amount. However, both techniques are costly with restricted accessibility and high requirements for expertise. 2,[8][9][10] Compared to MRI-PDFF, the transient elastography-based controlled attenuation parameter is less accurate for the diagnosis of hepatic steatosis. 11 Optical method, such as near-infrared hyperspectral imaging, can be an alternative that is a noninvasive, noncontact, nonionizing, and label-free technique. However, it suffers from penetration depth limitation due to the scattering of light. 12 Although a variety of techniques are available for hepatic steatosis evaluation, little is known about the liver function reserve (LFR) with the progression of NAFLD. 13 LFR is the remnant functional capacity of the liver after liver injury, which is considered an important factor for treatment and prognosis in chronic liver diseases. Indocyanine green (ICG) clearance test is the acknowledged approach for LFR assessment, the nature of which is to accurately detect ICG concentration change over time in blood. ICG can be selectively up taken by liver cells, eliminated to bile and excreted with feces through the bowel in the original form. ICG was not uptake by other organs or tissue and is an ideal drug for LFR reflection. In hepatic surgery, ICG retention rate at 15 min ðICGR15Þ < 30% is an essential condition for safe hepatectomy. ICGR15 > 40% is the contraindication for hepatectomy because posthepatectomy liver failure can occur with a large probability. [14][15][16] Photoacoustic imaging (PAI) is an emerging noninvasive modality focusing on the measurement of optical absorption of endogenous or exogenous absorbers, which can afford high-resolution images in superficial or deep tissue. Its clinical applicability has been demonstrated in breast cancer, thyroid cancer, melanoma, arthritis, bowel disease, [17][18][19][20][21] and in the detection of endogenous molecular chromophores, such as melanin, hemoglobin, lipids, and collagen. [22][23][24][25] Steatosis in fatty liver mice diagnosed with PAI indicated a prominent difference from normal mice. [26][27][28] However, previous applications were focused on small animals. Here we want to test the feasibility in a larger animal model with similar imaging depth to check its potential for clinical translations. In addition, combining the LFR assessment by PAI with general contrast agents 29 and the multispectral method has not been demonstrated on large animal models. Based on this prior study, 29 the goal of the current study is to show that PAI is capable of not only assessing LFR but also visualizing lipid content in fatty liver rabbits and to demonstrate PAI as a new noninvasive imaging biomarker for the assessment of NAFLD.

Experimental Setup
Both in vivo and ex vivo PAI were performed in this study. Before photoacoustic (PA) scanning, the anatomical regions of interest are scanned using an ultrasound system with a linear array working at 8.5 MHz, 100 Hz frame rate (iNSIGHT 23R, Saset Healthcare, Inc. Chengdu, China) by a single professional operator. The setup for in vivo PA imaging [ Fig. 1(a)] is based on a Qswitched Nd:YAG-pumped optical parameter oscillator (Surelite, Continuum, California), a 128element concave ultrasound transducer array (Japan Probe Corporation, Yokohama, Japan), and a 64-channel DAQ system (PXIe5105, National Instrument, USA). The diameter of the array is 100 mm. The pitch size of each element is 2 mm. The laser provides tunable laser pulses from 700 to 960 nm at a repetition rate of up to 20 Hz. The output laser light is coupled into a custommade fiber bundle. The output shape of the fiber bundle is a rectangular region with a 10 mm × 5 mm area. The maximum light fluence from the fiber bundle at 800 nm is estimated to be 13.6 mJ∕cm 2 (within the safety limit of the American National Standards Institute regulations). The concave ultrasound transducer array has a center frequency of 5 MHz (∼90% two-way relative bandwidth) with about 2.1 mm in radius one-way Nyquist zone 30,31 and a spatial resolution of ∼150 μm for cross-sectional imaging. The array is enclosed with a resin shell covered by a membrane (100-μm thickness), and water is filled into the space between the membrane and transducer. In the animal study, the PAI probe was put on the animal body (ear or skin surface), and the fiber bundle output was adjusted so that the illumination area was right at the center of the transducer [about 25 deg illumination angle, 15 to 20 mm height as shown in Figs. 1(b) and 1(c)].
Detected PA signals are amplified through a custom-built low-noise preamplifier (55 dB) and transferred to a 64-channel analog-to-digital system (PXIe5105, National Instrument, Texas) after 2:1 multiplexing with 12-bit digital resolution and 50 MS∕s sampling rate. Thus the frame rate of the system is 10 Hz. Signals are saved in the onboard computer (PXIe8840, National Instruments, Texas, USA) which also works as a control panel for the PAI system.
The setup for the ex vivo PA imaging is adapted from the in vivo system, as illustrated in Fig. 2(a). The output laser light is coupled into a custom-made fiber bundle to illuminate the whole area of a sample from the above, achieving an averaged light fluence of 13.6 mJ∕cm 2 on  (2) fiber bundles. DAS, data acquisition system (preamplifiers, multiplexers, and analog-to-digital converters); Nd:YAG, Nd:YAG pumping laser; OPO, optical parameter oscillator system; PC, personal computer; and imaging probe, a 128-element transducer array, a fiber bundle, and a resin shell enclosed with a membrane. the sample surface at 800 nm (again within the safety limit of the American National Standards Institute regulations). A sample piece is cut from an intact lobe of a rabbit liver with a thickness of 1 to 2 mm. The same concave ultrasound transducer array is used to perform 2D tomographic imaging of an ex vivo liver. Both the tissue and transducer are put in a small water tank. The liver tissue is fixed in a cylindrical sample holder made from resin and put in the center of the transducer to ensure maximum ultrasound sensitivity and the laser beam from the fiber was 15 mm above the tissue to produce a rounded area of illumination [ Fig. 2(b), water tank and sample holder is not shown here]. After amplification by our custom-built preamplifier, the PA signal is recorded with the DAQ system.

Animal Model Preparation
New Zealand rabbits (n ¼ 11, male, weight: 2.5 to 2.7 kg) were used in this study. The NAFLD group (n ¼ 8) was fed with a standard diet for the first 8 weeks, followed by a high-fat diet (75% standard feed with 10% butter, 2% cholesterol, 5% sugar, and 8% yolk powder) for the 12 weeks thereafter. The control group (n ¼ 3) was fed with the standard diet (100% standard feed) for 20 weeks. 32 These diets were provided by Dashuo Laboratory Animal Co., LTD, Chengdu, China.

In Vivo LFR Assessment
For the in vivo LFR assessment, 11 rabbits in two groups underwent PAI scanning of the liver region and ICG injection. All the animals were first anesthetized by isoflurane gas inhalation and were put on a warming pad and fixed to a custom-made animal pad to set in the right lateral position. We first used traditional ultrasonic imaging (US) to find the left lateral lobe of the liver for reference. The PA imaging probe was then placed at the same position as the US probe to acquire PA images, as shown in Fig. 1(c). Posture and imaging probe placements (both US and PA) were standardized and marked to ensure scanning reproducibility.
For in vivo imaging, multispectral PA signals were acquired at 760, 800, 850, and 930 nm to obtain oxygenated hemoglobin (HbO 2 ), deoxygenated hemoglobin (HbR), and lipid according to previous studies. 21 Twenty images were continuously acquired at one wavelength without averaging. After the multispectral imaging, the PAT imaging probe was placed above the right central auricular artery (CAA) as shown in Fig. 1(b) and the ICG was injected through the left auricular vein. Corresponding US images were acquired at the same position.
Before the PAI scan, ICG was dissolved in distilled water to obtain a concentration of 5 mg∕ml. PAI scanning at 805 nm was performed for 20 min over the right CAA. The acquisition rate was 1 frame per second (time-averaged for 10 consecutive frames). One minute after a baseline PAI scan, a dose of 0.1 ml∕kg ICG solution was injected via the left auricular vein along with a 5-ml saline flush for 10 s. Then the PA signal was continuously acquired for additional 19 min.

Ex Vivo and Histological Examination
After the in vivo imaging examination, representative tissue specimens were taken from previously imaged anatomical regions. After harvesting, tissue specimens from the same lobe were divided into two parts, one for ex vivo PA imaging, and the other one was fixed in a 4% formaldehyde/PBS solution and then embedded in paraffin.
For the ex vivo experiments, multispectral PA signals were acquired at 760, 800, 850, and 930 nm. For histological examination, liver tissue sections were stained with haematoxylineosin (HE) and Masson. The histologic features were evaluated with the NAFLD activity score (NAS), which was calculated as the unweighted sum of the scores for steatosis (0 to 3), lobular inflammation (0 to 3), and ballooning (0 to 2). Based on the NAS score, NAFLD was classified as "NAFLD but not-NASH" (NAS < 3), "borderline-NASH" (NASH ¼ 3 to 4), and "definite-NASH" (NAS ¼ 5 to 8). The classification process was performed by a pathologist with 8 years of experience in liver histology, who was blinded to the PA results.
To estimate the lipid area in H&E staining images, we used a machine learning toolbox (WeKa3) 33 in Fiji (National Institutes of Health). One expert labeled H&E staining images at NAS ¼ 3 which have a large area of lipid cells to train the network. Then the trained network was applied to classify all slices then the area of lipid cells was segmented by setting a threshold. The results were shown using the lipid cells area divided by the total area of the image.

Image Reconstruction and Spectral Unmixing
PA images were first reconstructed in Labview (National Instruments, USA) in real time using a back-projection algorithm 34 with real-time display. The raw PA data were also stored for offline processing by MATLAB (R2016b, Mathworks, Inc., MA, USA). In the offline processing, 0.2 to 10 MHz filtering was applied to remove high-frequency noise, and then PA images were reconstructed offline using back-projection schemes. 35 To remove the breathing-related artifact, we built a MATLAB code as described previously. 36 Briefly, after we got 20 frames (each frame obtained without averaging) from one wavelength, a correlation coefficient (CC) was calculated among these frames. Frames with movement showing significantly low CC (slight movement: 0.9 to 1, breathing: 0.5 to 0.7) were rejected. In most cases, about 12 to 14 frames were effective and averaged as one wavelength image. We then repeated the same process above for different wavelengths and registered all averaged images to that for the first wavelength for spectral unmixing.
A linear unmixing was applied to obtain the signal from three endogenous absorbers following previous studies, 21,37 and then values of HbR, HbO 2 , and lipid were obtained. Lipid unmixing was based on all wavelengths (760, 800, 850, and 930 nm), whereas the HbR and HbO 2 signals were calculated from a subrange (760, 800, and 850 nm), which is more accurate in unmixing due to increased water absorptivity at higher wavelengths. 38 The total hemoglobin (HbT ¼ HbR þ HbO 2 ) and lipid values were calculated from 4 mm × 2 mm (long × width) regions of interest close to the surface of the skin to minimize error due to the linear unmixing. 39 The averaged HbT signal, HbT mean , lipid signal, and lipid mean were calculated for statistical analysis.
To minimize the error from unknown light fluence, we first measured the laser energy below the imaging probe and normalized the energy among different wavelengths. Then we divided our image into two parts, skin and liver based on the clear structural differences in PA images and referred to corresponding US images. A depth compensation method based on recently published research was applied to each wavelength 40 for in vivo data. The relevant optical parameters were obtained from two other papers. 41,42 For ex vivo data, we measured the laser energy on the tissue surface and normalize the energy. Since only the liver tissues were imaged in 2D form. No depth compensation was applied.

Statistical Analysis
Statistical analysis and graphical display of data were performed using GraphPad software (version 7.00; GraphPad Software, San Diego, CA, USA). In the ICG results, continuous variables were given as means and SEM in the manually selected ROI region at the CAA. Continuous PA images were registered to the first frame. Similar to the previous study, 29 averaged PA signals in the CAA were normalized as the relative PA signal intensity (PSI re ) to reflect the ICG concentration change as follows: PSI base is the base PSI before ICG injection which is generated from hemoglobin. PSIðtÞ is the PSI at time t. To quantify the kinetics of PSI re change after the ICG injection, an exponential decay model 43 was used as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 2 . 6 ; 1 1 6 ; 1 3 9 PSI re ðtÞ ¼ PSI re ð0Þe −kt ; where the elimination time (t 0 ) indicates the elapsed time from PSI re ð0Þ to the baseline. The rate constant (k) gives the decay rate. A mean total hemoglobin (HbT mean ) and lipid signal (lipid mean ) were calculated over anatomical regions. In vivo and ex vivo statistical analyses between two groups (NAFLD versus control) of PA parameters and histological results were conducted using independent samples ttests. A value of p < 0.05 was considered statistically significant. The correlation between in vivo and ex vivo lipid signals, between lipid signal and lipid/unit area, were calculated using linear regression.

In Vivo ICG Clearance Studies
ICG clearance studies were conducted first by real-time PAI for both groups (NAFLD and control group). PAI and US sections of the CAA in vivo were displayed as the region of interest (ROI) in Fig. 3(a) (top and bottom panels, respectively) (yellow dashed circle and arrow). Figure 3(c) shows the in vivo ICG measurements from the CAA. The CCA can be found easily in both the PA and US images. After the injection of ICG at the 60-s time point, a high-intensity peak was observed and gradually decreased over time [ Fig. 3(c)]. Statistical results for k display a high value in the control group and a lower value in the lipid group. The time to recovery also showed significant differences between the two groups.
The control group showed a rapid decrease after the peak point, whereas the elimination of ICG for the NAFLD group was slower [ Fig. 3(c)]. Two major parameters were analyzed here. We calculated k to show a difference in elimination rate (k 1 versus k 2 , 0.0192 AE 0.0029 versus 0.00422 AE 0.0023, P < 0.0001). In addition, the elimination time (t 0 ), which indicates the elapsed time after the summit of the injection to the baseline, was also analyzed (t 01 versus t 02 , 174.7 AE 33.50 s versus 661.5 AE 131.2 s, P ¼ 0.0002).

Lipid Detection in a NAFLD Model
ICG clearance can serve as a biomarker for liver function but cannot provide lipid content, critical for accessing NAFLD. Hence, following the ICG clearance studies, multispectral PA studies were conducted.
US and PA imaging was performed in NAFLD rabbits and the control group. Spectral unmixing was applied to obtain HbT mean and lipid images [ Figs. 4(a) and 4(b)]. From the PA images, distinguishable features are present, indicating different tissue types around the liver. High HbT mean value regions reveal vascular-rich tissue areas, the dermal layer, and liver vessels. The high lipid signal regions were the subcutaneous lipid and liver. Due to the spectral unmixing, we can find significant differences in lipid concentration across the liver.
To select the ROI in the liver region for analysis, ultrasound images were used as a guide to ensure the liver's location. Then the HbT mean images showed distinguishable features which helped us to find the liver region below the skin layer and superficial vessels. Moreover, due to the limitation of our linear unmixing method, 39 a rectangle region was selected close to the top of the liver to calculate the mean lipid values (lipid mean ) and lipid/HbT ratio (lipid mean ∕HbT mean ) in this region.
A significant difference in the lipid mean value was observed (lipid mean ¼ 0.081 AE 0.0161 arbitrary units (a.u.) control versus NAFLD 0.198 AE 0.048 a.u., P ¼ 0.003). Second, when comparing total hemoglobin between control and NAFLD groups (HbT mean ¼ 0.857 AE 0.104 a.u. versus 0.95 AE 0.40 a.u., P ¼ 0.410), no significant difference was found. The lipid to hemoglobin ratio per independent rabbit between groups was also computed (lipid mean ∕HbT mean ¼ 0.101 AE 0.018 a.u. versus 0.218 AE 0.072 a.u., P ¼ 0.0255) and a showed significant difference. Fig. 4 In vivo PA imaging: representative images of liver from (a) healthy and (b) NAFLD rabbits. ROIs (yellow boxes) are determined in the US images for analysis. The PA/US images show HbT mean and lipid distributions within the liver and ROI, respectively. Differences in spectrally unmixed PA images between healthy and NAFLD rabbits are shown. (c), (e) Lipid mean as well as lipid mean ∕HbT mean for healthy and NAFLD rabbits. Each blue-filled circle and red-filled cube represents the mean PA signal per independent liver region for the healthy and NAFLD rabbits, respectively.

Ex Vivo Fatty Liver Imaging
To verify the in vivo findings, corresponding tissue specimens were taken from both groups. Representative PA images of HbT mean and lipid were shown in Fig. 5(a) and corresponding histological results were shown in Fig. 5(b). Similar to in vivo analysis, mean lipid values (lipid mean ) and lipid/HbT ratio (lipid mean ∕HbT mean ) were compared between two groups. First, a significant difference in lipid mean was observed (lipid mean ¼ 0.0673 AE 0.0165 versus 0.486 AE 0.073 a.u., P < 0.0001). Second, a significant difference was also seen in total hemoglobin between the two groups (HbT mean ¼ 0.163 AE 0.0194 versus 0.0396 AE 0.0201 a.u., P < 0.0001). Finally, lipid/HbT ratio showed a significant difference between the two groups (lipid mean ¼ 0.273 AE 0.179 versus 8.82 AE 3.40 a.u., P ¼ 0.002).
As shown in Fig. 6(a), there was a progressive elevation in the in vivo lipid mean with increased severity from not-NASH to borderline-NASH. There was a significant difference in lipid mean and NAFLD group (bottom row) from corresponding liver tissue. (c) Lipid mean as well as lipid mean ∕HbT mean for the healthy and NAFLD rabbits. Each blue-filled circle and red-filled cube represents the mean PA signal per independent liver region for the healthy and NAFLD rabbits, respectively.
In the ex vivo results, significant differences in lipid mean were shown between control and not-NASH (0.0673 AE 0.0166 versus 0.395 AE 0.0691 a.u. P ¼ 0.0001). As the lipid mean increased,  According to data in Table 1, the linear regression between ex vivo and in vivo lipid data is shown in Fig. 6(c) (r 2 ¼ 0.80, P < 0.001).
To evaluate the correlation between PA lipid estimation and the histology estimates, we calculated the relative area of lipid/unit area in each staining slice. Figures 6(d) and 6(e) showed that lipid value from PA images has a relatively high correlation with the lipid area. In vivo mean lipid values (lipid mean ) had a good correlation with histology estimates (r 2 ¼ 0.76, P < 0.001), and ex vivo lipid mean were well correlated with histology estimates (r 2 ¼ 0.67, P < 0.001). The H&E staining images showed normal hepatic tissue and steatosis pathology.

Discussion and Conclusion
In our previous work, we demonstrated the capability of PAI-based ICG clearance for LFR assessment in rabbits, 29 the accuracy of which was confirmed by the gold standard method spectrophotometry. In this work, we applied PAI/ICG to assess LFR in rabbits with NAFLD. We have shown that time-resolved PA signal can differentiate the kinetics of ICG uptake and excretion in control and NAFLD groups in vivo, indicating LFR impairment occurred during the development process of NAFLD.
Although hepatic ICG clearance serves as a biomarker for liver function in the assessment of NAFLD rabbit, NAFLD cannot be quantified by it. The pathology of NAFLD is a spectrum of lesions ranging from pure steatosis to a complex pattern with significant necroinflammatory injury, fibrosis, or even cirrhosis. 44 Therefore, LFR assessment was necessary but not enough for the diagnosis and treatment of the pathological change of NAFLD.
Multispectral PAI is a target-specific, quantitative, noninvasive, imaging modality 38 when compared with established imaging modalities, including ultrasound imaging and MRI. To obtain more specific information about NAFLD liver, in vivo and ex vivo multispectral PAI were performed and the results from our study demonstrated a significant increase in lipid content in the liver of the NAFLD group compared to the control group, which agreed well with the pathologic results, indicating lipid metabolism malfunction in the NAFLD liver tissue. 26,45 Both in vivo and ex vivo mean lipid values displayed a high correlation with histology estimates, which indicated the accuracy of using PAT for steatosis assessment of NAFLD. This was in line with previous studies in mice of NAFLD. 27,46 In addition, the linear correlation between lipid in vivo and ex vivo mean lipid is also relatively high.
The divergent outcomes for the in vivo and ex vivo HbT results may come from several reasons. First, the liver has a strong HbT mean signal as it is a vascular-rich tissue and HbT mean stands for strong endogenous absorbers in vivo. In an ex vivo setting, HbT mean signal decreased drastically due to blood loss in the tissue. In general, pathological changes associated with steatosis could be unequivocally detected using multispectral PAI both in vivo and ex vivo, consistent with similar findings from recent studies in mice. 46 In addition, the pathological changes in those two situations (in vivo and ex vivo), we note that there are some limitations in our current study. The light absorption of lipid in the 760-to 930-nm range is much weaker than that at 1200 and 1720 nm which might be a better wavelength to detect lipid regions. 47 As the limitation of our laser system, output energy at 1200 and 1720 nm is too low for PA imaging. As an alternative method, we refer to previous research 21,38 to use 930 nm and other wavelengths (e.g., 760, 800, 850, and 930 nm). Those wavelengths with some water absorption may induce inconsistent lipid areas as shown in Fig. 4(b) and show some "patchy" in the NAFLD liver.
Although we have a good correlation between PA lipid value and histology lipid area, we should notice that the r 2 is below 0.9 and more studies are needed. For example, we can briefly calculate the lipid area in H&E staining and compare it with our PA image results. Red oil O staining will be a better choice to estimate the lipid area. In addition, we need to optimize our system settings including the illumination geometry and fiber bundle parameters in both in vivo and ex vivo setups.
Except for the hardware, since we applied zero-thresholding to remove the negative signals in figures, some artifacts in reconstructed images may affect the quantitative results. 48 A recently developed method based on Hilbert transformation can be applied in the future to remove such noise and give more strict quantitative measurements. 31,48 In addition, the linear unmixing method may be another reason which affects our results. Advanced spectral unmixing methods 21,49 need to be developed to achieve the capability of refined quantification of PA signal.
In sum, this study showed significantly injured LFR in the NAFLD group compared to the control group. A future study with a larger sample size is needed to investigate the LFR change during various NAFLD pathologies. Due to the limitation of the sample size used in this study, refined PA classification of the NAFLD is not achieved. We are developing an advanced NAFLD animal model with various pathology states and upgrading our PAI system for this challenging study.

Disclosures
The authors declare no conflicts of interest.