The liver is one of the most important organs and is capable of regenerating itself  . Many researchers have become interested in development of methods for reconstruction of the liver based on tissue engineering techniques  . The sinusoid is defined as a group of small blood vessels and the bile canaliculus is a thin tube that collects bile secreted from hepatocytes; these two tissues are arranged along cords of hepatocytes, without crossing each other, as shown in Figure 1. From a mathematical perspective, these types of periodic networks cannot be fully understood in two dimensions   . However, to date, imaging analyses in hepatology have largely relied on two-dimensional (2D) images because appropriate methods for three-dimensional (3D) reconstruction and quantification have not been developed.
Confocal microscopy has been used to analyze 3D structures of cells and tissues after immunofluorescence staining allowing the examination of relationships between cell arrangement and metabolic function   . To analyze the 3D sinusoidal network, it is necessary to explicitly segment the sinusoid network because the sinusoids are formed by sinusoid endothelial cells, as shown in Figure 1(a), and these cells can be imaged by immunostaining using confocal microscopy.
The segmentation process can be very time consuming; therefore, it is fundamental to choose the right techniques for filtering the image properly. An automated method for segmentation that saves time and human labor is always desirable.
A number of automated methods have been proposed   . It is known that methods using a binary method and a region extraction processing using edge extraction as methods were proposed by the 1990s  . The methods using statistical analysis and deformable shape model have been mainly used in recent researches from the latter half of the 1990s  .
Here, we propose a new method for segmentation of 3D sinusoidal networks utilizing the information performed with the reaction-diffusion (RD) model. Segmentation and edge detection of 2D images based on this model have been reported   . In these previous studies, the excitable media-type RD model
Figure 1. (a) A schematic structure of microstructures in the liver; (b) 3D fluorescence image reconstruction of the rat liver. Red and yellow-green points represent sinusoidal endothelial cells and bile canaliculi, respectively. By image analysis, the positions of hepatocytes were determined (a representative is indicated in blue, although hepatocytes were positioned in all empty spaces in the reconstruction).
was used for segmentation one or two materials in the 2D image. In contrast, in this study, a Turing-type of RD system was used for segmentation of 3D periodic networks, in which various types of network patterns self-organized    .
Liver disease is often associated with alterations in the microarchitecture of the liver  . Therefore, using the proposed method, the morphology of the 3D network patterns was evaluated during the changes in morphology observed in the context of nonalcoholic steatohepatitis (NASH), an advanced stage of nonalcoholic fatty liver disease (NAFLD), which is related to metabolic syndrome. Rats fed a high-fat/high-cholesterol (HFC) diet, which causes pathological features similar to those of human patients with NASH  , were used to evaluate the 3D pattern index of 3D segmentation sinusoidal network patterns.
2. Material and Methods
2.1. Animal Experiments
All animal experiments were performed in accordance with the Guideline for Animal Experiments of Kwansei Gakuin University and with approval of the Committee for Animal Research at Kwansei Gakuin University. Six-week-old male Wistar rats were purchased from Shimizu Laboratory Supplies Co. Japan, and randomly divided into eight groups. For the diets, stroke-prone control chow diet (20.8% crud protein, 4.8% crude liquid, 3.2% crude fiber, 5.0% crude ash, 8.0% moisture, and 58.2% carbohydrate) was used as a control diet, and the HFC diet was a mixture of 68% control diet, 25% palm oil, 5% cholesterol, and 2% cholic acid. Control groups (Cont) were fed the control diet for 3, 6, 9, or 12 weeks. HFC groups (HFC) were fed HFC diet for 3, 6, 9, or 12 weeks. Both diets were obtained from Funabashi Form (Chiba, Japan). The rats were housed in groups of three in standard breeding cages (27 × 22 × 12 cm) with freely available food and water under a 12-h light/12-h dark cycle (light on at 08:00).
After 18 - 20-h fast from the last feeding, all rats were sacrificed under pentobarbital (70 mg/kg)-induced anesthesia, and the livers were removed. A part of each liver was fixed in 4% buffered paraformaldehyde for histological analysis.
2.2. Immunofluorescent Staining and Observations by Confocal Microscopy
An immunofluorescence technique was applied to 40-μm-thick frozen sections of liver using a monoclonal antibody specific for hepatic sinusoidal endothelial cells (Mouse anito-SE1; Immuno-Biological Laboratories Co. Ltd., Japan)   . Alexa 594-conjugated rabbit secondary antibodies were purchased from Molecular Probes (Eugene, OR, USA).
Primary antibodies were diluted in phosphate-buffered saline (PBS) with 1% bovine serum albumin (BSA) and 0.2% Triton X-100 (at 1:50) and applied overnight at 4˚C. Fluorescence-conjugated antibodies were diluted in PBS with 1% BSA and 0.2% Triton X-100 (at 1:150), and the sections were incubated for 1 day at 4˚C. Sections were washed three times with Tris buffer (0.1 M Tris-HCl, pH 7.6, and 0.15 M NaCl) and twice with PBS with 0.2% Triton X-100, followed by overnight incubation with primary and fluorescence-tagged antibodies, respectively. After sections were mounted in fluorescent mounting medium (Vector Laboratories, Inc., Burlingame, CA, USA), confocal Z-stack images were obtained using an Olympus FV 1000 confocal microscope running FluoView version 2.0 c software (Olympus, Tokyo, Japan). For each 3D fluorescence image, 50 frames (640 × 640 pixel images) were obtained with a length of 0.50 μm between pixels and frames.
2.3. 3D Sinusoidal Segmentation
3D sinusoidal segmentation was based on a Turing RD model. To derive the sinusoid area in 3D reconstructed images, information processing was performed with a Turing RD model, in which temporal and spatial patterns are self-organized   . In this system, the interesting phenomena of Turing pattern formation had been reported    . The method for image segmentation was based on the RD mode (FitzHugh-Nagumo equation), with modifications as follows:
where, , , , , , and are positive constants; is the control parameter; the variable and are local concentrations of the activator and inhibitor, respectively; and indicates the intensity of 3D fluorescence images of sinusoid endothelial cells. Figure 2 shows the schematic diagram of the segmentation process. We first scale the [0, 255] scale image into [−0.5, 0.5] range linearly. The initial conditions of and were given equilibrium value with white noise without any spatial correlations.
Setting and , Turing conditions have been known as follows. Equation (1) includes a time-independent and uniform solution , which is defined through and . The steady state is stable in the ordinary differential equations given by f(u, v), and g(u, v), we have
, and ,
where the partial derivatives are evaluated at equilibrium . Since the steady-state solution becomes unstable in the partial differential equations given by Equation (1),
The model Equation (1) satisfying these conditions is called Turing system
Figure 2. A schematic showing of the diagram of numerical calculations.
 . Parameters have been selected as satisfying these conditions.
Figures 3(a)-(c) show the time-evolution of the distributions of u and v in one dimension. Previous studies have shown that static periodic patterns are self-organized   . Moreover, 3D Turing patterns have been studied previously   . In cases in which , the self-organized patterns were entrained to the distribution of . Figures 3(d)-(f) show the time-evolution of distributions of u and v with external data U in one dimension. Considering the situation in which local differences in intensities of fluorescence images existed and the distributions were of the kink type, with dents and different periodicities, the prepared distribution of U was utilized, as shown in Figure 3(d). Figure 3(f) shows the obtained distribution. Although the prepared distribution U was bumpy and exhibited spatially different amplitudes, the amplitude of the obtained distribution u became almost the same throughout the area, and the local periodicities of u and U became identical after the numerical calculation process.
To extend this method in 3D spaces, numerical simulations of Equation (1) were carried out in 3D. The space was divided into rectangle cells of sizes , , and , which were the same size as the pixel sizes of the images obtained by confocal microscopy. The Neumann boundary conditions were imposed at the system boundaries. The simple Euler algorithm was used with a time step , and the stop time to calculate was set to t = 10.0 in the case with pixel information . As shown in Figure 2, after computation, we obtained as the output image, which was scaled back into the [0, 255] scale image to produce the final output.
Figure 3. Time evolution of (thick red line), and (thick blue line). (a)-(c) were obtained numerically from Equation (1) with = 0 for Du = 0.15, Dv = 15.0, = 0.50. = 0.02, = 26.0, and = 0.5; and (d)-(f) were self-organized under the pixel information for the same with = 0.1. The dotted line indicates U(x). Since the patterns in (d)-(f) were generated much faster than that for (a)-(c), the times of the patterns are different between (b) and (e), and between (c) and (f).
One of the most crucial problems is to determine the parameter choice of . By changing , the spatial period of the self-organized pattern would changes, verifying segmentation performed precisely with eyes. When was comparable to the spatial period of fluorescence images, the patterns obtained using Equation (1) were almost identical. The simulations were repeated by changing the parameter . To identify the suitable for fluorescence images, the pattern index was calculated as follows
where , , , and are the mean and the variance of u and U, V is the total volume of the region, and is the correlation between u and U. To each 3D fluorescence image obtained by confocal microscopy, the with largest was selected as the chosen parameter for segmentation of 3D sinusoidal network patterns (Figure 4).
2.4. Statistical Analysis
The results obtained for δ*, and volume ratios of 3D sinusoidal networks, were compared using two-tailed Mann-Whitney U-test (MW test) between Cont and HFC groups. Moreover, the results were evaluated by the Kruskal-Wallis (KW) test followed with the two-tailed MW test among HFC groups. The statistical analyses were performed using SPSS 16.0 for Macintosh (SPSS Inc., Chicago, IL, USA).
3.1. 3D Segmentation Images of Sinusoidal Networks
Figure 5 show representative results for 3D segmentation of sinusoid networks
Figure 4. An example of 3D segmentation of sinusoidal network patterns of HFC6w. (a), (b): 3D segmentation patterns using raw pixel data of the fluorescence image of sinusoidal endothelial cells, and (c), (d): 3D segmentation patterns obtained by RD processing of Equation (1) with Du = 0.15, Dv = 15.0, = 0.50. = 0.02, = 26.0, = 0.1, and δ* = 1.35. (a) and (c) indicate the 3D segmentation patterns of sinusoidal network (red tubes), and (b) and (d) show the slices at the middle position of the z-axis from (a) and (c), where the white area indicates the positions inside the sinusoids.
Figure 5. Examples of 3D segmentation of sinusoidal networks. (a)-(d) and (e)-(h) show representative images from the Cont and HPC, respectively. (a) and (e): 3 weeks, (b) and (f): 6 weeks, (c) and (g): 9 weeks, and (d) and (h): 12 weeks. The red area indicates the sinusoidal network. The illustration of optical transmittance of 50% has been used.
from fluorescence pixel information utilizing the RD algorisms. We note that the black rectangle domains close to the boundaries are not closely seen there since the illustration of optical transmittance of 50% has been used. Changing the parameter in Equation (1), we calculated = 1.00, 1.05, 1.10, 1.10, 1.45, 1.65, and 1.70 for Cont at 3, 6, 9, and 12 weeks, and the HFC at 3, 6, 9, and 12weeks, respectively. We independently calculated for four segmentations of 3D sinusoidal networks for each individual. Since each experimental group has three individuals, we calculated twelve for each experimental group. Significant differences were observed between Cont and HFC at 3, 6, 9, and 12 weeks.
Moreover, we found that were increasing depending on the weeks of feeding of HFC diets, as shown in Figure 6. Significant differences were also observed among the of HFC groups. The results of statistical analyses are presented in Figure 6.
Figure 6. Change in δ* with largest I(δ). The number of stars indicates the statistical level of significance (★: p < 0.05. ★★: p < 0.01). MW test between Cont at 3 weeks (Cont3w) vs. HFC at 3 weeks (HFC3w), p = , Cont6w vs. HFC6w, p = , Cont9w vs. HFC9w, p = : Cont12w vs. HFC12w, p = . KW test among HFC3w, HFC6w, HFC9w, and HFC12w, p = ; MW test between HFC3w vs. HFC6w, p = , HFC6w vs. HFC9w, p = ; HFC 9w vs. HFC12w, p = .
3.2. Volume Ratios of Sinusoidal Network Volumes
In the previous section, we demonstrated that were capable of differentiation response the weeks of feeding, that is the degree of progress of NASH disease (shown in Figure 6). To evaluate the index proposed, the ratios of sinusoidal volumes were compared. Figure 7 shows the summary of the results. Significant differences were found between Cont and HFC at 9 and 12 weeks, and only between 6 and 9 weeks among HFC groups. The results of statistical analyses are presented in Figure 7.
In this study, a method for segmentation of 3D sinusoidal networks using the Turing RD model was proposed. The model could generate the distributions with the same degree of amplitude in all areas, although the pixel information obtained by confocal microscopy included the locality of the amplitude. Moreover, local periodicities of the distribution obtained from Equation (1) were the same as the pixel information, as shown in Figure 3(f). This was one of the advantages of the proposed method utilizing the Turing RD model.
Two other advantages of this method were the cost and time for performing segmentation of 3D network patterns. Although some of commercial software for 3D reconstruction of cell tissue images is available, such software is expensive and requires high-performance computers with huge amount of memory to perform 3D structure analysis. However, the proposed system could be applied even with laptop computers using the simple code for numerical simulation shown in Equation (1).
In this study, 3D structure analysis of rat livers was carried out using proposed
Figure 7. Change of volume ratios of 3D segmentation of sinusoidal networks. The number of stars indicates the statistical level of significance (★: p < 0.05. ★★: p < 0.01). MW test between Cont9w vs. HFC9w, p = : Cont12w vs. HFC12w, p = ; but Cont3w vs. HFC3w, p = , Cont6w vs. HFC6w, p = . KW test among HFC3w, HFC6w, HFC9w, and HFC12w, p = ; MW test between HFC6w vs. HFC9w, p = , but HFC3w vs. HFC6w, p = , HFC 9w vs. HFC12w, p = .
algorithm. The parameter captured the variations in feeding patterns for rats fed the HFC diet: these patterns were related to the degree of progress of NASH. This is another advantage of the proposed method.
Since the parameter is related to the period of the patterns obtained using Equation (1), it is possible that the periodicity of the 3D network pattern should be essential for detecting differences among the obtained patterns. However, strong localities in periodicities of 3D sinusoidal network pattern were observed. Therefore, it was not possible to detect clear periodicities in 3D sinusoidal patterns using calculations such as Fourier analysis or spatial correlation analysis.
Here, it is noted that the simulation time of Equation (1) is very short, 10 times less than native self-organized mode generated as shown in Figure 3(c) and Figure 3(f). If we performed much longer time simulation to segmentation sinusoidal patterns, the periodic patterns with same periodicity with naive self-organized mode kc would be obtained.
Although the index about the volume-ratios of sinusoid was good ones, was better one. This may be cause by the fact that is related to diffusion coefficients that are squarely affected to the periodicities of pattern, whilst volumes are affected lineally.
We also examined the 3D patterns of HFC at 16 weeks (n = 1) in the same way as described above, as precautionary measure. The index was not as large as ones of of HFC group at 12 weeks. The significant difference between HFC at 12 weeks and HFC at 16weeks was not detected (MW test HFC at 12 weeks vs. HFC at 16weeks, p = 0.743). It is known that the severe fibrosis would occur as the disease progresses on behalf of the morphological changes in those stages  . Hence, it was likely that the increasing in might reach maximum between HFC at 9 weeks and HFC at 16weeks. Therefore, we did not test HFC at 16weeks in details.
Finally, the periodic patterns appearing in the biological and medical field often include more localities than those appearing in physical or electric processes. Thus, the method proposed in this study would be useful for detecting periodicities with some localities. Therefore, the proposed method is expected facilitate further analysis of the 3D network structures encountered in biology or medicine.
We would like to thank K. Ono for his technical assistance and H. Hontani for his helpful comments. This research was supported by JSPS KAKENHI Grant Number 17H05302.