Dynamics and mechanisms for tear breakup (TBU) on the ocular surface

1.   Introduction

Rats are commonly used as experimental objects to study various mechanisms in the current medical research, such as the field of medicine ^[1,2,3,4,5]. To date, lots of biomedical laboratory animals are rats, and they play a vital role in the experiments of newly developed medical preparations ^[2,6,7,8]. The main reason is that their genetic composition and biological behavior characteristics are extremely similar to humans, and about 80% of them are genetically identical to humans. Moreover, rats modified by genetically modified technology ^[9,10] can even carry genes that cause illness in humans.

Due to the outbreak of the new coronavirus at the end of 2019 (COVID-19) ^[11,12,13], many researchers in the medical field are urgently concerned about respiratory diseases research. As a very important branch of respiratory research, chest images ^{[14,15,16,17]} obtained through imaging equipment can clearly show abnormal and normal lung lesions. Therefore, clinicians can diagnose most lung diseases through chest images. For most respiratory diseases, such as lung infections including COVID-19 ^{[12,13,17,18]}, bronchiectasis ^[19,20], chronic obstructive pulmonary diseases (COPD) ^[21,22], or idiopathic pulmonary fibrosis ^[23] etc., the computed tomography (CT) has unique advantages to achieve the purpose of diagnosis. Especially for a patient with COVID-19, although the nucleic acid test ^[24,25,26] is negative, as long as there is a lesion on the CT image, infection of COVID-19 can be suspected and treated as this disease.

Many researchers currently use rats for respiratory research ^[27,28,29], but few researchers can effectively measure lung parenchyma parameters of rats. The quantitative analysis of the rat lungs is largely limited. Currently, the gold standard for determining lung function in human is routine pulmonary function testing (PFT) ^[30,31]. The researches ^[32,33] have shown that the lung volume, density, pixel index, and other parameters of lung images acquired by CT technology have a good correlation with the pulmonary function data acquired by the PFT. Therefore, the corresponding lung parenchyma parameters from the lung images of rats may be useful for respiratory researches.

At present, most researchers are committed to using artificial intelligence methods to diagnose lung diseases from CT images ^{[14,18,34,35,36]}. However, after the diagnosis of the disease, corresponding medical measures need to be taken to treat the disease. For a newly developed drug, it is first necessary to conduct drug experiments on animals to observe the response of diseased animals to the drug. For respiratory diseases, inhalation therapy ^[37,38,39] is often a very effective solution. Therefore, rats are bred into models with corresponding diseases and use the models to experiment with newly developed inhaled pharmaceutical preparations. The easiest and most feasible way is to image rats through CT and measure the changes in lung parenchyma parameters from the images to understand the actual effects of the drugs.

Our paper proposes a method for measuring rats' lung parenchyma parameters of the right lung (RL), the left lung (LL), and the bilateral lung (BL) including the CT value, the density, the lung volume (LV), and the surface area (SA). After getting the CT value, the density of the BL, LL, and RL, the average CT value (ACT) and the average density (AD) which can reflect the situation of the BL, LL, and RL can be calculated.

The above lung parenchyma parameters all have their important clinical significance ^[40,41,42]. If the density in local regions is higher than that in other regions of lung tissue, and the ACT in local regions will be correspondingly higher than that in other regions. In clinical diagnosis, an increase or decrease in the ACT of lung tissue will indicate lung tissue disease. An increase in the ACT of lung tissue is generally seen in the lung parenchyma inflammation, masses, nodules, etc., and a decrease in the ACT is more common in emphysema and pneumothorax. If the AD indicating the amount of lung tissue per unit volume is abnormal, there may be some changes such as lung tissue proliferation and remodeling, which may indicate different degrees of disease in the clinic. The LV refers to the volume occupied by the lung tissue in the thoracic cavity. Clinically, the increase in the LV is commonly seen in pathologies such as emphysema and masses, and the decrease in the LV can be seen in pathologies such as lung consolidation and fibrosis. The SA often easily be overlooked is an important parameter for researching the increase or the decrease of lung tissue per unit area due to emphysema.

In this paper, the left lung and right lung of rats are segmented using a convolutional neural network segmentation model, and then the BL, LL, and RL of rats are extracted, and finally lung parenchyma parameters of the BL, LL, and RL can be measured. Based on the available original images, our purpose is not only to show the differences in the same lung parenchymal parameter of the two groups of rats, but to provide a method for measuring the lung parenchymal parameters of rats to promote current medical respiratory researches. To verify the correctness of the method proposed in our paper, the ACT and the AD are calculated and the relationship between the LV and the SA is analyzed. The results show that they are in line with the law of changes in the lung parenchyma. At the same time, it is concluded that the lung volume is negatively correlated with the CT value or the density.

2.   Materials and methods

2.2.   Materials

The dataset of rats' pulmonary window images with window level -600 and window width 1200 used in our study is from the department of radiology in the first affiliated hospital of Guangzhou medical university. Original images are 14 sets of datasets in DICOM format obtained by CT scanning for patients. The parameter of the CT made by SIEMENS is set to 80 KV, 50 mAs, Pith: 0.5, and FOV: 100, respectively.

Seven healthy SPF SD rats are named the control group before inhaling the nebulization. After continuously inhaling the nebulization, these seven rats were then named the treatment group. Weight of the seven rats is in the range of 90-100 g. During the 12 hours before the inhalation, the seven rats were fasted and could only drink water. A kind of X-ray contrast agent named the iodixanol injection or the Visipaque was used as the nebulized medicine. Basic information of the dataset is that the size of each pulmonary window images is 512 × 512 × N and the voxel spacing is 0.1953 * 0.1953 * 1 (mm). N represents the number of slices of each pulmonary window images, and its range is 37 to 52 in our dataset.

2.2.   Methods

To measure the lung parenchyma parameters of rats, the lung region of rats is segmented to obtain the lung segmentation masks, and then the BL, LL, and RL of rats with Hounsfield unit (HU) values are extracted. Finally, the lung parenchyma parameters of the BL, LL, and RL can measured or calculated, respectively. Figure 1 shows the overall implementation plan of our article to measure the lung parenchyma parameters of rats including a lung region segmentation model, a lung extraction model, and a lung parameter measurement model.

As shown in Figure 1(a), the lung region segmentation model is used to segment the lung region of rats to obtain the lung segmentation masks. Specifically, the original DICOM images are input to the lung region segmentation network named U-net(R231) ^[43], and the U-net (R231) outputs the lung segmentation masks of the LL and RL. A detailed description of the lung region segmentation model for rats is shown in Section 2.2.1.

After getting the lung segmentation masks in Figure 1(a), the lung parenchyma extraction model in Figure 1(b) is used to extract the lung parenchyma of the BL, LL, and RL with HU values. A detailed description of the lung parenchyma extraction model with HU value is shown in Section 2.2.2.

After obtaining the lung parenchyma through the lung extraction model in Figure 1(b), the lung parameter measurement model is used to measure the CT value and calculate the LV. Then, the HU value to density conversion unit is used to convert the CT value to the density. Before calculating the SA, the tracheal boundary line needs to be removed from the boundary of the lung segmentation masks in Figure 1(c). A detailed description of the parameter measurement model for rats is shown Section 2.2.3.

2.2.1.   Lung region segmentation model for rats

Since all automatic segmentation methods cannot achieve 100% accuracy, a semi-automatic segmentation method is used in our paper. In the lung region segmentation model, the deep learning network model is used for the coarse segmentation to obtain the mask of the LL and the mask of the RL. To ensure the accuracy of the segmentation, the fine segmentation is manually performed on the results of the coarse segmentation.

1) Coarse segmentation

There is not enough information of the pulmonary window images with 36 to 52 slices, so it is difficult to use traditional methods ^[44,45,46] to segment the BL, LL, and RL. Until the network model named U-net (R231) ^[43] is discovered, we try to input the original DICOM images of rats into the U-net (R231) which had trained by the human lung and the segmentation results is beyond our expectations. The U-net (R231) is actually a U network model with residual building block (ResU-Net) and the architecture of ResU-Net model is shown in Figure 2.

Figure 2(a) shows the architecture of U-Net model. The input of the ResU-Net model is the original DICOM images with the size 512 × 512 × N, and the output of the ResU-Net model is the lung segmentation masks of the LL and RL with the size 512 × 512 × N. The U-Net model includes a contracting path (red dotted box) and an expansive path (bule dotted box). The contracting path includes 5 sub-encoder modules named UNetConvBlock (0) to UNetConvBlock (4). Each sub-decoder module followed by an average-pooling operation (red arrow) is as down-sampling modules, and 5 sub-encoder modules are connected in sequence on the left of Figure 2(a). Similarly, the expansive path includes 5 sub-decoder modules named UNetUpBlock (0) to UNetUpBlock (4). Each sub-decoder modules copies the feature map generated by the corresponding layer from sub-encoder modules in the contracting path and connects it with the feature map generated from up-sample operation (green arrow) on the right of Figure 2(b). The output of the ResU-Net model as segmentation result is generated by a 1 × 1 convolution with Softmax.

Figure 2(a), (b) shows the detailed architecture of the contracting path (red dotted box) and the expansive path (bule dotted box) in the U-Net model, respectively. All sub-encoders are residual building blocks with two feedforward branches in the contracting path. The outputs of the two feedforward branches are added together by a shortcut connection as an output of each sub-encoder. The original DICOM images as the input of UNetUpBlock (0) is decoded by a residual building block. A feedforward branch of the UNetUpBlock (0) includes a 3 × 3 convolution, sequentially followed by a rectified linear unit (ReLU), a batch normalization (BN) and a 3 × 3 convolution. The other feedforward branch of the UNetUpBlock (0) includes a 1 × 1 convolution followed by a BN. The UNetConvBlock (1) to UNetConvBlock (4) continue to decode the original DICOM images with the same residual building blocks as shown in Figure 2(b). The residual building blocks of the UNetConvBlock (1) to UNetConvBlock (4) also have two feedforward branches. A feedforward branch of the UNetConvBlock (1) to UNetConvBlock (4) includes a ReLU sequentially followed by a BN, a 3 × 3 convolution, a ReLU, a BN, and a 3 × 3 convolution. The architecture of the other feedforward branch of the UNetConvBlock (1) to UNetConvBlock (4) is the same as that of the UNetUpBlock (0).

Similarly, the UNetUpBlock (0) to UNetUpBlock (4) are the same residual building blocks as shown in Figure 2(c). There are a residual building sub-block (1) and a residual building sub-block (2) in each UNetUpBlock (i), and the residual building sub-block (1) sequentially followed by a 1 × 1 convolution and a BN as a feedforward branch of the residual building sub-block(2). The input of the residual building sub-block (1) is directly as the output the other feedforward branch of the residual building sub-block (1), and the outputs generated from two feedforward branches of the residual building sub-block (2) are added together by a shortcut connection as an output of each UNetUpBlock (i). Specifically, there are two feedforward branches in the residual building sub-block (1). A feedforward branch of the residual building sub-block (1) includes a 3 × 3 convolution, sequentially followed by a ReLU, a BN, a 3 × 3 convolution, a ReLU, and a BN. The other feedforward branch of the residual building sub-block (1) includes a 1 × 1 convolution followed by a BN. A output of the residual building sub-block (1) is generated by an addition of the two feedforward branches in the residual building sub-block (1).

Experienced radiologists can compare the original DICOM images with the segmentation masks of the corresponding slice to judge whether the segmentation result is correct. If two experienced radiologists out of three believe that a segmentation result is acceptable, then the segmentation result can be considered acceptable. Figure 3 shows the typical segmentation results of a set in the dataset. all slices except the 14th slice of the segmentation results is acceptable judged by experienced radiologists, and only the 14th slice of the segmentation results needs to be manually corrected by experienced radiologists.

2) Fine segmentation

To compensate for the segmentation errors of the lung region segmentation model in certain slices, a tool named ITK-SNAP is used to manually correct those incorrect segmentation results. The software named ITK-SNAP can be available on the website, http://www.itksnap.org/pmwiki/pmwiki.php?n=Downloads.SNAP3.

Figures 4(a1)-6(c1) show two segmentation errors with blue circles and blue rectangular frames from the view of the coronal plane, the sagittal plane, and the axial plane with window setting value 3144 and level setting value 533. Radiologists use the ITK-SNAP to fill in or modify the segmentation errors to get accurate segmentation results for this slice. Specifically, after selecting the Paintbrush Mode the Paintbrush Inspector panel can be seen in the ITK-SNAP Toolbox to fill in or modify the segmentation errors. Due to the structural features of the lung, when modifying or filling in the lung boundary of the lung segmentation masks, we recommend choosing the third type of Brush Styles like the quarter moon and using a small Brush Size in the Paintbrush Inspector panel. In the Paintbrush Inspector panel, the Adaptive Algorithm is the default value with Granularity 20 and Smoothness 15.

After filling or modifying the segmentation errors, the mask of the LL and the mask of the RL is exported as the NIFTI format for each slice shown in Figures 4(a2)-6(c2). If the segmentation errors are modified in Figure 4(a1), the corresponding masks in Figure 4(b1), (c1) can be also automatically modified.

2.2.2.   Lung parenchyma extraction model

After the fine segmentation, the lung parenchyma extraction model is used to extract the lung parenchyma of the BL, the BL, LL, and RL with HU values. In the lung parenchyma extraction model, the original DICOM images multiply by the corresponding slice of the lung segmentation masks of the BL, LL, and RL to get the lung parenchyma with HU values. For details, please refer to our previous paper ^[47]. Note that when the original DICOM images multiply by the lung segmentation masks, make sure that each slice of the original DICOM images and the corresponding slice of the lung segmentation masks are on the same slice.

Figure 5 shows the process of typical lung parenchyma extraction of the BL, LL, and RL with HU values. Figure 5(a1)-(a3) show the original DICOM images from three perspectives respectively. Figure 5(b1)-(b3) show the three images of lung segmentation masks from the three perspectives respectively, and Figure 5(c1)-(c3) show the lung parenchyma extraction results of the lung extraction model.

2.2.3.   Lung parenchyma parameter measurement model for rats

Now that the lung parenchyma with HU values has been extracted, the lung parenchyma parameters of the BL, LL, and RL can be measured, respectively. the HU value corresponding to each pixel in the lung parenchyma can easily got, and then the average HU values to get the ACT are calculated. Before calculating the average HU values, because HU values in the lung parenchyma cannot be greater than or equal to 0 HU, these HU values have been deleted.

When calculating the density, the following formula named the HU value to density conversion is used as mentioned in Figure 1(c) to convert HU values to density values.

The conversion formula:

where $\alpha = \left({{D_{density\_air}} - {D_{density\_water}}} \right)/{H_{HU\_air}}; {D_{density\_water}}{\rm{ = 999}}.{\rm{8395mg/mL}}; {D_{density\_air}}{\rm{ = 1}}.{\rm{2922}}$ mg/mL; H_{HU__air} = -1000 HU.

The parameter H_{HU__lung} is the HU values corresponding to each pixel less than 0 HU. Through the above conversion formula, the density of the BL, LL, and RL can be got, respectively. The unit of density values obtained from the conversion formula is mg/mL, then converted to g/mL.

The LV can be obtained by counting the total number of pixels less than 0 HU and the volume of each pixel. Specifically, the total number of pixels of the BL, LL, and RL is multiplied by the corresponding pixel volume to get the LV of the BL, LL, and RL. The pixel volume is also known as the voxel spacing and can be obtained through the header file of the DICOM file. In our paper, the value of the pixel volume is 0.03814 mm³ (0.1953 mm × 0.1953 mm × 1 mm).

When measuring the SA of the BL, LL, and RL, only the boundary pixels of them are needed by the edge detection method with the first-order Canny operator. The SA of the BL, LL, and RL can be obtained by multiplying the total number of boundary pixels by the corresponding pixel volume. The result of the BL, LL, and RL after using the edge detection method is shown in Figure 6, and there are trachea edges seen in the red circle of Figure 6 (a1)-(c1) which will impact on our calculation result of the SA. To eliminate the impact, a connected domain strategy is adopted. The connected domain strategy can delete trachea edges whose area is smaller than the value P, using 8 neighborhoods. In our paper, the value P can be configured as 50. After using the connected domain strategy, the boundary of the BL, LL, and RL is without the tracheal boundary line shown in Figure 6(a2)-(c2). In our paper, we actually detect the lung region segmentation masks obtained from the lung region segmentation model and get the boundary of masks with the tracheal boundary line. Of course, the boundary of the BL, LL, and RL is multi-slice boundaries. When calculating the SA, the total number of the boundary pixels of all slices together can be counted, and then multiply the total number by the pixel volume with 0.03814 mm³.

3.   Results

In our paper, the lung parenchyma parameter of the BL, the LL and RL from the pulmonary window images of 7 healthy rats before and after inhaling the nebulization is measured, and mainly compare the CT value and the density before and after inhaling the nebulization and analyze the relationship of the lung parenchyma parameter. Healthy rats before inhaling the nebulization are named the control group, and the rats after inhaling the nebulization are named the treatment group. The following analysis of the lung parenchyma parameter is carried out by statistical software named GraphPadPrism (8.0.1).

3.1.   Comparison of the CT value and the density

To verify the method proposed in our paper, the changes of the CT value and the density of the BL, LL, and RL from the two groups are compared by observing whether the changes conform to the law of the inhaling the nebulization.

The ACT of the LL and RL from the control group in Figure 6(a1) are -578.3 and -560.3 HU, and that of the LL and RL from the treatment group in Figure 6(a2) are -526.3 and -518.8 HU. Difference between mean (DM) ± standard error of mean (SEM) between the ACT of the LL and the ACT of the RL from the control group is about 17.96 ± 0.50, and that from the treatment group is about 7.52 ± 0.48. Figures 6(a1) and 7(a2) show that the ACT of the LL and RL from their own group are no obvious difference. However, Figure 7(b1-c) clearly show that the ACT of the BL, LL, and RL from the treatment group are significantly greater than that from the control group, respectively. The significant ACT difference between the control group and the treatment group is determined by both t tests (P < 0.0001) and nonparametric tests (P < 0.0001). Specifically, the ACT of the LL from the two groups in Figure 6(b1) are -578.3 and -526.3HU, and the ACT of the RL from the two groups in Figure 6(b2) are -560.3 and -518.8 HU, and the ACT of the BL from the two groups in Figure 6(c) are -567.9 and -521.9 HU. DM ± SEM of the BL, LL, and RL from the two groups is about 51.98 ± 0.52, 41.54 ± 0.45, and 46.04 ± 0.34, respectively. Above all, it is found that the ACT in the treated group have significantly increased.

After analyzing the CT value, we continue to analyze the changes in the density in the BL, LL, and RL under the conditions of the control group and the treatment group. The density values can be calculated by the conversion formula. The relationship between density values and HU value is essentially a linear relationship. The AD of the RL and the LL from the control group is about 0.4404 and 0.4224 g/mL, and that from the treatment group is about 0.4818 and 0.4743 g/mL. DM ± SEM of the RL and the LL from the control group is about 0.0179 ± 0.0005, and that of the RL and the LL from the treatment group is about 0.0075 ± 0.0005. The AD difference of the RL and the LL from their own group is very small. However, the AD of the RL, the LL and the BL in the treatment group is significantly greater than that in the control group. The significant AD difference between the control group and the treatment group is determined by both t tests (P < 0.0001) and nonparametric tests (P < 0.0001). Specifically, the AD of the LL from the two groups are 0.4224 and 0.4743 g/mL, and that of the RL from the two groups are 0.4404 and 0.4818 g/mL, and that of the BL from the two groups in Figure are 0.4327 and 0.4787 g/mL. DM ± SEM of the BL, LL, and RL from the two groups is about 0.0519 ± 0.0005, 0.0414 ± 0.0005, and 0.0460 ± 0.0003, respectively. Above all, it is also found that the AD in the treatment group have significantly increased.

3.2.   Relationship between the LV and the SA

To further verify the accuracy of the method proposed in our paper, the relationship between the LV and the SA is tested.

Figure 8 shows that the relationship between the LV and the SA of the BL, LL, and RL in the control group and the treatment group, respectively. When the SA increases, the LV will increase with the SA, which is consistent with the law of lung changes during breathing. Although the regression relationship changes slowly in the treatment group as shown in Figure 8(b1)-(b3), it can also reflect the positive proportional relationship between the LV and the SA. But the value of both the LV and the SA in the treatment group is significantly less than that in the control group. Specifically, Figure 8 shows that the value of both the LV and the SA in the LL is less than that in the RL.

3.3.   Relationship between the LV and the CT value or the density

After verifying the accuracy of the method proposed in our paper, the relationship between the LV and the CT value or the density is analyzed. We make the regression analysis on the CT value and the LV of the BL, LL, and RL in the control group and the treatment group, respectively. The regression equations which show that the CT value and the LV are negatively correlated.

Figure 9 shows that the LV of the BL, LL, and RL decreases as the CT value increases in both the control group and the treatment group. From the overall distribution of CT value in Figure 9, most CT value of the BL, LL, and RL in the treatment group is smaller than that in the control group, and the corresponding LV in the treatment group reduces. As the previous conversion formula, CT value and the density are linear. Similarly, the density and the LV of the BL, LL, and RL are also negatively correlated in both two groups.

4.   Discussion

The changes in the lung parenchyma parameters often reflect the development of the disease in medical respiratory researches. We propose a method to measure the CT value, the density, the LV, and the SA of the lung parenchyma parameters in our paper. Although the lung region segmentation model is a semi-automatic segmentation model, it has greatly reduced the workload of lung segmentation compared to the manual segmentation, and the semi-automatic segmentation can ensure the accuracy of lung parameters measurement. The dataset of seven healthy SPF SD rats obtained from CT for lung parenchyma segmentation in our paper, and we believe that if the dataset obtained from the micro-computed tomography for small animals has a higher resolution, there will be better results of the lung parenchyma segmentation. Using the lung region segmentation model and lung parenchyma extraction model, the lung parenchyma of the LL and RL is obtained, so quantitative analysis on lung parenchyma parameters of the LL and RL can be performed.

To verify the effectiveness of the method proposed in our paper, the ATC and the AD of the BL, LL, and RL under the conditions of the control group and the treatment group are analyzed, respectively. Theoretically, due to the continuous inhalation of the nebulization, the absorption of X-rays by the lung will increase and the CT value and the density of the treatment group would correspondingly increase. The ATC and the AD have significantly increased in the treatment group, which is consistent with the theoretical changes of the CT value and/or the density. the relationship between the LV and the SA of the BL, LL, and RL in the control group and the treatment group is also analyzed, and it is concluded that they are the positive proportional relationship. Because when the SA increases, the lung including the BL, LL, and RL is in a process of inhalation, and the LV will inevitably increase. Similarly, when the SA decreases, the lung is in a process of exhalation, and the LV will inevitably decrease. The change of the LV with the SA in both two groups conforms to the law of lung changes during breathing which means that the LV and the SA are positively related. We can see from the original DICOM images and Figure 8 that the LV of the LL from rats is also smaller than that of the RL which also conforms to morphological characteristics of the lung.

The relationship between the LV and the CT value or the density of the BL, LL, and RL is also preliminarily analyzed, and then it is concluded that the LV is negatively correlated with the CT value or the density. Due to the inhalation of the nebulization, the breathing movement of the lungs is restricted, so the LV of the BL, LL, and RL is significantly reduced. We think that if the relationship between lung parenchyma parameters can be clearly explained, the amount of dataset in the two groups need to be increased. At the same time, if the breathing process of rats can be controlled, there will be other important discoveries using the method proposed in our paper. Finally, what we have done above is to help promote the quantitative analysis of lung parenchyma parameters of rats in medical respiratory researches.

5.   Conclusions

We propose an effective method to measure the lung parenchyma parameters of rats from pulmonary window computed tomography images based on ResU-Net model. The above results show that the CT value, the density, the LV, and the SA of the lung parenchyma parameters conform to the law changes of rats' lung. After comparing and analyzing the lung parenchyma parameter CT value and the density of the two groups, it is discovered that the lung parenchyma parameter CT value and the density significantly increase in the treatment group. The relationship between the lung volume and the CT value or the density is also analyzed and get the conclusion that the LV is negatively correlated with the CT value and/or the density.

Acknowledgments

Thanks to the Department of Radiology, the First Affiliated Hospital of Guangzhou Medical University for providing the dataset of rats' pulmonary window images and the funding support from the Natural Science Foundation of Guangdong Province, China (2019A1515011382) and the National Natural Science Foundation of China (62071311).

Conflict of interest

The authors declare no conflict of interest in this paper.