2.1.

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 $Q$-switched Nd:YAG-pumped optical parameter oscillator (Surelite, Continuum, California), a 128-element 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 custom-made fiber bundle. The output shape of the fiber bundle is a rectangular region with a $10\mathrm{mm}\times 5\mathrm{mm}$ area. The maximum light fluence from the fiber bundle at 800 nm is estimated to be $13.6\mathrm{mJ}/{\mathrm{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 ($\sim 90\%$ two-way relative bandwidth) with about 2.1 mm in radius one-way Nyquist zone^30,31 and a spatial resolution of $\sim 150\mu \mathrm{m}$ for cross-sectional imaging. The array is enclosed with a resin shell covered by a membrane ($100\text{-}\mu \mathrm{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)].

Fig. 1

A schematic of the in vivo PAI experimental system. (a) The hardware components employed in the PAI system. (b), (c) Schematics showing two in vivo imaging positions at the central artery on the ear and liver, respectively. (d) Photograph of the PA probe: (1) ultrasound transducer and (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.

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\mathrm{MS}/\mathrm{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\mathrm{mJ}/{\mathrm{cm}}^2$ on 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.

Fig. 2

(a) The hardware components employed in ex vivo experiments. (b) Picture of the PA probe. DAS, data acquisition system (preamplifiers, multiplexers, and analog-to-digital converters); OPO, optical parameter oscillator system; PC, personal computer; and Nd:YAG, Nd:YAG bumping laser.

2.2.

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.³² These diets were provided by Dashuo Laboratory Animal Co., LTD, Chengdu, China.

2.3.

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 (${\mathrm{HbO}}_2$), deoxygenated hemoglobin (HbR), and lipid according to previous studies.²¹ 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\mathrm{mg}/\mathrm{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\mathrm{ml}/\mathrm{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.

2.4.

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 haematoxylin–eosin (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” ($\mathrm{NAS}<3$), “borderline-NASH” ($\mathrm{NASH}=3$ to 4), and “definite-NASH” ($\mathrm{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)³³ in Fiji (National Institutes of Health). One expert labeled H&E staining images at $\mathrm{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.

2.5.

Image Reconstruction and Spectral Unmixing

PA images were first reconstructed in Labview (National Instruments, USA) in real time using a back-projection algorithm³⁴ 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.³⁵

To remove the breathing-related artifact, we built a MATLAB code as described previously.³⁶ 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, ${\mathrm{HbO}}_2$, and lipid were obtained. Lipid unmixing was based on all wavelengths (760, 800, 850, and 930 nm), whereas the HbR and ${\mathrm{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.³⁸ The total hemoglobin ($\mathrm{HbT}=\mathrm{HbR}+{\mathrm{HbO}}_2$) and lipid values were calculated from $4\mathrm{mm}\times 2\mathrm{mm}$ (long × width) regions of interest close to the surface of the skin to minimize error due to the linear unmixing.³⁹ The averaged HbT signal, ${\mathrm{HbT}}_{\text{mean}}$, lipid signal, and ${\text{lipid}}_{\text{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⁴⁰ 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.

2.6.

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,²⁹ averaged PA signals in the CAA were normalized as the relative PA signal intensity (${\mathrm{PSI}}_{\mathrm{re}}$) to reflect the ICG concentration change as follows:

A mean total hemoglobin (${\mathrm{HbT}}_{\text{mean}}$) and lipid signal (${\text{lipid}}_{\text{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 $t$-tests. 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.

3.1.

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.

Fig. 3

ICG clearance studies between NAFLD and control groups of rabbits (NAFLD, $n=8$ and control group, $n=3$). (a) Transverse PA (top) and US (bottom) images of the rabbit ear with the CAA (yellow dashed circle). ROI (yellow dashed circle) was selected to determine the PA signal intensity. (b) The elimination rate constant $k$ (left) and elimination time $t_0$ (right). (c) Temporal PA signal trace from 11 rabbits in two groups. For each curve, data were normalized to the highest intensity acquired.

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\pm 0.0029$ versus $0.00422\pm 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\pm 33.50\mathrm{s}$ versus $661.5\pm 131.2\mathrm{s}$, $P=0.0002$).

3.2.

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 ${\mathrm{HbT}}_{\text{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 ${\mathrm{HbT}}_{\text{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.

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 ${\mathrm{HbT}}_{\text{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) ${\text{Lipid}}_{\text{mean}}$ as well as ${\text{lipid}}_{\text{mean}}/{\mathrm{HbT}}_{\text{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.

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 ${\mathrm{HbT}}_{\text{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,³⁹ a rectangle region was selected close to the top of the liver to calculate the mean lipid values (${\text{lipid}}_{\text{mean}}$) and lipid/HbT ratio (${\text{lipid}}_{\text{mean}}/{\mathrm{HbT}}_{\text{mean}}$) in this region.

A significant difference in the ${\text{lipid}}_{\text{mean}}$ value was observed (${\text{lipid}}_{\text{mean}}=0.081\pm 0.0161$ arbitrary units (a.u.) control versus NAFLD $0.198\pm 0.048$ a.u., $P=0.003$). Second, when comparing total hemoglobin between control and NAFLD groups (${\mathrm{HbT}}_{\text{mean}}=0.857\pm 0.104$ a.u. versus $0.95\pm 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 (${\text{lipid}}_{\text{mean}}/{\mathrm{HbT}}_{\text{mean}}=0.101\pm 0.018$ a.u. versus $0.218\pm 0.072$ a.u., $P=0.0255$) and a showed significant difference.

3.3.

Ex Vivo Fatty Liver Imaging

To verify the in vivo findings, corresponding tissue specimens were taken from both groups. Representative PA images of ${\mathrm{HbT}}_{\text{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 (${\text{lipid}}_{\text{mean}}$) and lipid/HbT ratio (${\text{lipid}}_{\text{mean}}/{\mathrm{HbT}}_{\text{mean}}$) were compared between two groups. First, a significant difference in ${\text{lipid}}_{\text{mean}}$ was observed (${\text{lipid}}_{\text{mean}}=0.0673\pm 0.0165$ versus $0.486\pm 0.073$ a.u., $P<0.0001$). Second, a significant difference was also seen in total hemoglobin between the two groups (${\mathrm{HbT}}_{\text{mean}}=0.163\pm 0.0194$ versus $0.0396\pm 0.0201$ a.u., $P<0.0001$). Finally, lipid/HbT ratio showed a significant difference between the two groups (${\text{lipid}}_{\text{mean}}=0.273\pm 0.179$ versus $8.82\pm 3.40$ a.u., $P=0.002$).

Fig. 5

Ex vivo PA imaging from liver samples (control = 3 and NAFLD = 8). (a) Representative images of liver samples for healthy (upper row) and NAFLD (bottom row) rabbits. Averaged PA signals, ${\mathrm{HbT}}_{\text{mean}}$, and lipid distribution were calculated based on a small sample from the whole liver (red dashed line), respectively. Differences in spectrally unmixed PA images between healthy and NAFLD rabbits are shown. (b) Representative H&E staining for the control group (upper row) and NAFLD group (bottom row) from corresponding liver tissue. (c) ${\text{Lipid}}_{\text{mean}}$ as well as ${\text{lipid}}_{\text{mean}}/{\mathrm{HbT}}_{\text{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.

3.4.

Relationship Between PA Signal and Histology

We also compared the PA signal (i.e., ${\text{lipid}}_{\text{mean}}$) with histopathological classification. The details of the histopathological scores on steatosis (0 to 3), lobular inflammation (0 to 3), and ballooning (0 to 2) of each rabbit are shown in Table 1. Based on the NAS score, NAFLD group rabbits were classified as “NAFLD but not-NASH” ($\mathrm{NAS}<3$), “borderline-NASH” ($\mathrm{NASH}=3$ to 4), and “definite-NASH” ($\mathrm{NAS}=5$ to 8).

Table 1

Detailed histopathological scores and PA values NAFLD rabbits and control group (C1 to C3).

As shown in Fig. 6(a), there was a progressive elevation in the in vivo ${\text{lipid}}_{\text{mean}}$ with increased severity from not-NASH to borderline-NASH. There was a significant difference in ${\text{lipid}}_{\text{mean}}$ between control and not-NASH ($0.0802\pm 0.0159$ versus $0.160\pm 0.0237$ a.u. $P=0.04$). The ${\text{lipid}}_{\text{mean}}$ increased significantly in borderline-NASH groups and showed significant differences between not-NASH and borderline-NASH ($0.160\pm 0.0237$ versus $0.236\pm 0.033$ a.u. $P=0.009$). Significant differences were also shown between control and borderline-NASH ($0.0802\pm 0.0159$ versus $0.236\pm 0.033$ a.u., $P=0.0007$).

Fig. 6

Scatter plots of the comparisons between PA signal and histopathological results. (a) In vivo quantitative assessment of steatosis by PAI in comparison to histological grading. Control: $n=3$; not-NASH: $n=4$; and borderline: $n=4$. (b). Ex vivo quantitative assessment of steatosis by PA in comparison to histological grading. (c). Linear regression between ex vivo and in vivo lipid data. (d) In vivo quantitative assessment of lipid by PAI in comparison to histological estimates. (e) Ex vivo quantitative assessment of lipid by PAI in comparison to histological estimates. (f) Representative H&E staining for the control group (left) and NAFLD group (right) from corresponding liver tissue.

In the ex vivo results, significant differences in ${\text{lipid}}_{\text{mean}}$ were shown between control and not-NASH ($0.0673\pm 0.0166$ versus $0.395\pm 0.0691$ a.u. $P=0.0001$). As the ${\text{lipid}}_{\text{mean}}$ increased, significant differences were shown between control and borderline-NASH ($0.06733\pm 0.01656$ versus $0.544\pm 0.0518$ a.u. $P=0.0109$). The ${\text{lipid}}_{\text{mean}}$ between not-NASH and borderline-NASH also showed significant differences ($0.395\pm 0.0691$ versus $0.544\pm 0.0518$ a.u. $P=0.0109$).

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 (${\text{lipid}}_{\text{mean}}$) had a good correlation with histology estimates ($r^2=0.76$, $P<0.001$), and ex vivo ${\text{lipid}}_{\text{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.