**Large-Scale Image Analysis for Investigating Spatio-Temporal Changes in Nuclear DNA Damage Caused by Nitrogen Atmospheric Pressure Plasma Jets**

#### **Xu Han 1,2,**† **, James Kapaldo 1,2,**† **, Yueying Liu <sup>3</sup> , M. Sharon Stack 3,4, Elahe Alizadeh <sup>5</sup> and Sylwia Ptasinska 1,2,\***


Received: 9 May 2020; Accepted: 8 June 2020; Published: 10 June 2020

**Abstract:** The effective clinical application of atmospheric pressure plasma jet (APPJ) treatments requires a well-founded methodology that can describe the interactions between the plasma jet and a treated sample and the temporal and spatial changes that result from the treatment. In this study, we developed a large-scale image analysis method to identify the cell-cycle stage and quantify damage to nuclear DNA in single cells. The method was then tested and used to examine spatio-temporal distributions of nuclear DNA damage in two cell lines from the same anatomic location, namely the oral cavity, after treatment with a nitrogen APPJ. One cell line was malignant, and the other, nonmalignant. The results showed that DNA damage in cancer cells was maximized at the plasma jet treatment region, where the APPJ directly contacted the sample, and declined radially outward. As incubation continued, DNA damage in cancer cells decreased slightly over the first 4 h before rapidly decreasing by approximately 60% at 8 h post-treatment. In nonmalignant cells, no damage was observed within 1 h after treatment, but damage was detected 2 h after treatment. Notably, the damage was 5-fold less than that detected in irradiated cancer cells. Moreover, examining damage with respect to the cell cycle showed that S phase cells were more susceptible to DNA damage than either G1 or G2 phase cells. The proposed methodology for large-scale image analysis is not limited to APPJ post-treatment applications and can be utilized to evaluate biological samples affected by any type of radiation, and, more so, the cell-cycle classification can be used on any cell type with any nuclear DNA staining.

**Keywords:** atmospheric pressure plasma jets; large-scale imaging; machine learning; cancer treatment; cellular imaging

#### **1. Introduction**

In recent years, numerous in vitro studies have shown the considerable anticancer effects of nonthermal atmospheric pressure plasmas in approximately 20 types of malignant cell lines, including lung cancer [1], prostate cancer [2], ovarian cancer [3], osteosarcoma [4], and oral cancer [5]. Furthermore, several in vivo investigations using tumor models of pancreatic cancer [6], glioblastoma [7], melanoma [8,9], ovarian cancer [10], and breast cancer [11] have demonstrated the significant inhibition of cellular growth and tumor damage following atmospheric pressure plasma treatment. The ability of atmospheric pressure plasma jets (APPJs) to inactivate or kill malignant cells relies strongly on the production of a variety of plasma reactive species [12,13]. APPJs synergistically provide free electrons, positive ions, radicals, photons, and electromagnetic fields, which can damage biological targets without elevating the temperature of the treated area [14]. More importantly, plasma treatments in animal models have been reported to selectively damage targeted cancer cells, without affecting surrounding healthy tissues [15,16]. These features suggest that nonthermal atmospheric pressure plasmas may represent a promising alternative to conventional cancer treatments [14,17].

Although some primary clinical studies have previously been performed [18–20], the extensive clinical applications of APPJs require more detailed investigations to examine their effects on a variety of cancer cell lines, both in vitro and in vivo [21,22]. There is concern regarding the potential carcinogenic risk and side effects of prolonged clinical use due to the formation of free radicals. These can cause adverse and acute impacts that can present safety risks in long-term APPJ applications [14,23,24]. Also, technical issues, such as the optimal plasma dosage inside tissues, the penetration depth of reactive species, and the distribution of cellular damage, remain poorly understood and require further investigations. A variety of bioanalytical tools and imaging techniques have been used to quantify the induced damage and cellular responses following plasma irradiation, including fluorescence microscopy [25–27] and flow cytometry [28]. While these techniques can be utilized to perform routine cellular analyses, each possess both advantages and limitations, in terms of sample preparation requirements, sensitivity, measurable parameters, throughput, and costs. For example, fluorescence microscopy can capture images of small sample regions with high spatial resolution, facilitating the assessment of quantitative morphology [29]. In contrast, flow cytometry can facilitate the analysis of cellular kinetics and cell-cycle phases, but cannot provide spatial information; however, highly sensitive multicolor phenotypic data can be obtained from populations of different cells, within minutes [30].

In the current study, first we explored two dimensional (2D) spatial distributions of damage to deoxyribonucleic acid (DNA) induced by the APPJ treatment of cancer and nonmalignant cells. DNA damage was assessed by measuring double-strand break (DSB) formation in cell nuclei. In the cellular environment, DSBs trigger the phosphorylation of histone H2AX near the break site, resulting in the appearance of γH2AX foci and leading to local changes in the chromatin structure. These modifications are macroscopic structures that can be directly visualized with the assistance of antibody staining inside the cell nuclei.

Second, we developed a large-scale image analysis technique, using machine learning-based cell-cycle classifications, requiring only one staining dye. Generally, the cell cycle is divided into two major phases: interphase, including gap 1 (G1), DNA synthesis (S), and gap 2 (G2), and mitotic (M) phase. During the G1 phase, the cell grows in size at a high biosynthetic rate, producing proteins and copying organelles such as mitochondria and ribosomes to prepare for DNA synthesis (S phase). After DNA duplication, cells enter the second gap phase, G2, during which they grow rapidly and synthesize proteins and organelles in preparation for mitosis. As cells enter the M phase, they stop growing and synthesizing proteins to focus their energy on the complex and highly regulated cell division process. During standard analyses, different dyes are used to stain nuclear DNA in each cell phase. In this study, we used only 40 ,6-diamidino-2-phenylindole (DAPI) fluorescence stain to image and classify the cycle of every cell on a coverslip. Although we employed DAPI staining to build a correlation between nuclear DNA versus nucleus size for cell classification, our methodology can be applied to other types of nuclear DNA staining. However, some dyes can also stain mitochondrial DNA outside of the nucleus, and therefore will not be eligible for this method, because this method uses segmentation of nuclei due to their ellipse-like shape.

By combining microscopy images with a machine learning tool, we were able to study the spatial distribution of plasma-induced DNA damage in nuclei within each cell-cycle phase. By mapping the damage distribution over various post-treatment (incubation) times, in both malignant and nonmalignant oral cells, we revealed the spatio-temporal dependence of cellular responses to plasma-influenced regions. These results improve our understanding of APPJ effects on biological targets and their applications in plasma medicine. Furthermore, our proposed methodology for analyzing DNA damage in a large number of irradiated cells could facilitate the quantitative evaluation of the DNA damage caused by any other sources of cellular radiation, with widespread application in radiobiology research. *Int. J. Mol. Sci.* **2020**, *21*, x FOR PEER REVIEW 3 of 17 nonmalignant oral cells, we revealed the spatio-temporal dependence of cellular responses to plasmainfluenced regions. These results improve our understanding of APPJ effects on biological targets and their applications in plasma medicine. Furthermore, our proposed methodology for analyzing DNA damage in a large number of irradiated cells could facilitate the quantitative evaluation of the DNA damage caused by any other sources of cellular radiation, with widespread application in

#### **2. Results and Discussion 2. Results and Discussion**

radiobiology research.

surface was scanned.

#### *2.1. Computational Analysis 2.1. Computational Analysis*

Figure 1a shows representative fluorescence images of nuclear DNA (denoted as the DAPI channel) and DSBs in nuclear DNA (denoted as the γH2AX channel) on a large scale. The slide scanner's movements followed a zig-zag pattern, first scanning the entire length of the coverslip along one direction (called the fast axis), then taking a short step to the side (called the slow axis), before again scanning the entire length of the coverslip in the opposite direction, until the whole surface was scanned. Figure 1a shows representative fluorescence images of nuclear DNA (denoted as the DAPI channel) and DSBs in nuclear DNA (denoted as the γH2AX channel) on a large scale. The slide scanner's movements followed a zig-zag pattern, first scanning the entire length of the coverslip along one direction (called the fast axis), then taking a short step to the side (called the slow axis), before again scanning the entire length of the coverslip in the opposite direction, until the whole

**Figure 1.** (**a**) Representative fluorescence images of cells, after 2-min plasma treatment and 1-h incubation. Green and blue fluorescence signals correspond to γH2AX and DAPI, respectively. (**b**) Effects of background correction. Plots show the nuclear DAPI intensity along the slow axis of the fluorescence slide scanner. Colors encode the nuclear density (along the slide scanner's fast axis). A histogram of the nuclear DAPI intensity is indicated to the right of each plot. The top and bottom **Figure 1.** (**a**) Representative fluorescence images of cells, after 2-min plasma treatment and 1-h incubation. Green and blue fluorescence signals correspond to γH2AX and DAPI, respectively. (**b**) Effects of background correction. Plots show the nuclear DAPI intensity along the slow axis of the fluorescence slide scanner. Colors encode the nuclear density (along the slide scanner's fast axis). A histogram of the nuclear DAPI intensity is indicated to the right of each plot. The top and bottom plots show data before background correction and after flattening and stripe artifact removal, respectively.

The image analysis is comprised of four primary steps, including nuclei segmentation and background/foreground correction, feature extraction, machine learning to classify the cell phase, and quantification of DNA damage.

#### 2.1.1. Nuclei Segmentation and Image Correction

Nuclei segmentation was achieved by filtering the DAPI images with a Gaussian blur filter (σ = 1) and then using an adaptive, log-weighted Otsu threshold [31]. This method successfully segmented all nuclei that were not in contact with other nuclei; however, in our studies, the initial confluence of the cells was 90%, and confluence increased during the incubation time, resulting in many overlapping nuclei. These overlapping nuclei needed to be further segmented for additional computational analyses, which was accomplished by applying a short-range attraction and a long-range repulsion (SALR) clustering algorithm [32], in which we identified the center of each nucleus following a geometric partitioning algorithm to segment the overlapping nuclei [33].

In addition, intensity correction of the image background and foreground is crucial for successful analysis. The image foreground represents all regions (pixels) containing nuclei; these regions were determined by the segmentation in the previous step. Likewise, the image background represents all regions not containing nuclei. Three primary causes of uneven background/foreground intensities were identified across our images: uneven fluorescence staining, which results in slowly changing background/foreground intensities; differing amounts of microscope illumination/collection, which can be caused by improper focusing across the whole coverslip; microscope calibration errors, which led to the appearance of bright stripes along the fast axis of the scanned images.

Using the background region, we computed the spatially variable background intensity for each channel and removed it from the image. After this background flattening, we computed the background stripe artifacts and removed them. A representative stripe artifact can be observed in the large γH2AX image in Figure 1a. DAPI channel foreground correction was performed before the foreground correction of other channels. Two bands of DAPI intensities are observed in Figure 1b, which likely correspond to the G1 and G2 phases, since cells in the G2 phase contain twice as much DNA as cells in the G1 phase and because cells spend the most time in the G1 and G2 phases. We fitted these two bands with locally weighted linear fits (lowess) to flatten the G1 and G2 bands and position them at DAPI intensity values of 1 and 2, respectively. Figure 1b shows the nuclear DAPI intensities across the slow axis of the scanner, before (top) and after (bottom) correction. After the correction, the G1 and G2 DAPI bands are flat and positioned correctly. The foreground correction for the γH2AX channel was performed with caution, as the intensity values should not be flat because plasma treatment causes localized DSBs in the nucleus. We first qualitatively selected all G1 cells, using the corrected DAPI intensities from above (cells with DAPI values in the range of 0.7–1.3). Using these cells, we fitted a surface to the locally varying data in the bottom 2%–4% of the γH2AX intensities and then subtracted this surface from the data. This process will correctly flatten the γH2AX channel as long as some of the cells in a large region of the image (~1 mm<sup>2</sup> ) are not influenced by the plasma, which is true in our images based on the fact that the maximum damage ratio of G1 cells is 50% (as described below) and therefore, the bottom 2%–4% of γH2AX intensities represent undamaged cells.

#### 2.1.2. Feature Extraction

From each nucleus, we extracted a set of features from the DAPI channel for classifying cell-cycle phases. The features describe the nucleus shape, intensity, radial (shell) intensity, Haralick texture [34], and granularity. The Haralick and granularity features were computed using the same method as the commonly used CellProfiler software (version 2.1.1) [35]. The nuclear shapes were described, in part, using Fourier descriptors [36]. A list of the extracted features can be found in Appendix A. From the γH2AX channel, we extracted intensity features only. The integrated intensity of the γH2AX channel for each nucleus is proportional to the quantity of DSBs in each nucleus.

The code for feature extraction was written to compute the features of multiple nuclei in parallel on a graphical processing unit (GPU), in contrast with CellProfiler, which processes images in parallel on central processing units (CPUs) and processes nuclei features in serial. Extracting the nuclei features on a GPU results in a significant decrease in processing time when a large number of CPUs are not available. When applying our code during some basic tests, the entire processing time for the images (including segmentation and correction, which also used GPU acceleration), including feature extraction, required approximately 40–60 min per image, for both the DAPI and γH2AX channels, based on approximately 3 GB per channel. When using four CPUs, the same processing in CellProfiler took 2–3 h, without the application of SALR clustering for the accurate partitioning of nuclei clumps [35]. In addition, when using a simple watershed-based partitioning of overlapping cells, which is more similar to the method used by CellProfiler, our processing time is shortened to 20–30 min per image. The code was written in MATLAB (The MathWorks, Inc., Natick, MA, USA), and the hardware used for comparisons were an Nvidia GTX770, with 2 GB of memory, and an Intel 4 core i7-4790, 3.6 GHz. of the γH2AX channel for each nucleus is proportional to the quantity of DSBs in each nucleus. The code for feature extraction was written to compute the features of multiple nuclei in parallel on a graphical processing unit (GPU), in contrast with CellProfiler, which processes images in parallel on central processing units (CPUs) and processes nuclei features in serial. Extracting the nuclei features on a GPU results in a significant decrease in processing time when a large number of CPUs are not available. When applying our code during some basic tests, the entire processing time for the images (including segmentation and correction, which also used GPU acceleration), including feature extraction, required approximately 40–60 min per image, for both the DAPI and γH2AX channels, based on approximately 3 GB per channel. When using four CPUs, the same processing in CellProfiler took 2–3 h, without the application of SALR clustering for the accurate partitioning of nuclei clumps [35]. In addition, when using a simple watershed-based partitioning of overlapping cells, which is more similar to the method used by CellProfiler, our processing time is shortened to 20–30 min per image. The code was written in MATLAB (The MathWorks, Inc., Natick, MA, USA), and the hardware used for comparisons were an Nvidia GTX770, with 2 GB of memory, and an Intel 4 core i7-4790, 3.6 GHz.

Appendix A. From the γH2AX channel, we extracted intensity features only. The integrated intensity

#### 2.1.3. Cell-Cycle Classification 2.1.3. Cell-Cycle Classification

The cell-cycle phase was classified using two steps. During the first step, supervised machine learning was used to classify five visually distinct classes. Class 1 included interphase (G1, S, and G2), and the other four classes were derived from the different mitosis stages, including prometaphase (class 2), metaphase (class 3), early anaphase (class 4), and late anaphase, telophase, and early G1 (class 5). Examples of these classes are presented in Figure 2a. We manually classified approximately 500 nuclei in each class to create a balanced dataset, which was divided into training, validation, and test sets at a 0.7/0.15/0.15 ratio. The employed learning method was a shallow neural network, with two hidden layers (sizes 30 and 10), and a softmax layer to output the final classifications. The training used early stopping with the validation data to prevent overtraining. The confusion matrix on the test dataset is shown in Figure 2c, indicating that the overall classification accuracy was 97%. The network creation and training were all performed using MATLAB built-in routines, namely the "patternnet" and "train" functions. The cell-cycle phase was classified using two steps. During the first step, supervised machine learning was used to classify five visually distinct classes. Class 1 included interphase (G1, S, and G2), and the other four classes were derived from the different mitosis stages, including prometaphase (class 2), metaphase (class 3), early anaphase (class 4), and late anaphase, telophase, and early G1 (class 5). Examples of these classes are presented in Figure 2a. We manually classified approximately 500 nuclei in each class to create a balanced dataset, which was divided into training, validation, and test sets at a 0.7/0.15/0.15 ratio. The employed learning method was a shallow neural network, with two hidden layers (sizes 30 and 10), and a softmax layer to output the final classifications. The training used early stopping with the validation data to prevent overtraining. The confusion matrix on the test dataset is shown in Figure 2c, indicating that the overall classification accuracy was 97%. The network creation and training were all performed using MATLAB built-in routines, namely the "patternnet" and "train" functions.

**Figure 2.** (**a**) Examples of nuclei from each classification group; 1: interphase, 2: prometaphase, 3: metaphase, 4: early anaphase, 5: late anaphase, telophase, and early G1. The images are 34 µm × 34 µm, and they are modified to show a white background, for ease of viewing. (**b**) Nuclear area versus DAPI intensity (color gives nuclear density), with classification contours that enclose regions that contain primarily G1 (solid), S (dashed), G2 (dot-dashed), and M (black/white dash) phase cells. The right M contour primarily contains classes 3 and 4, whereas the left M contour primarily contains class 5. (**c**) Confusion matrix of the test set, based on the classification of different cell classes. **Figure 2.** (**a**) Examples of nuclei from each classification group; 1: interphase, 2: prometaphase, 3: metaphase, 4: early anaphase, 5: late anaphase, telophase, and early G1. The images are 34 µm × 34 µm, and they are modified to show a white background, for ease of viewing. (**b**) Nuclear area versus DAPI intensity (color gives nuclear density), with classification contours that enclose regions that contain primarily G1 (solid), S (dashed), G2 (dot-dashed), and M (black/white dash) phase cells. The right M contour primarily contains classes 3 and 4, whereas the left M contour primarily contains class 5. (**c**) Confusion matrix of the test set, based on the classification of different cell classes.

During the second step, the interphase cells were classified using a mixture-of-Gaussians model, with a uniform background [37], using the nuclear area and DAPI intensity (Figure 2b). The classification used seven Gaussians, equally spaced along the line connecting the center of the G1 and G2 peaks (Figure 2b), and an eighth Gaussian, located on the G1 peak. This classification of During the second step, the interphase cells were classified using a mixture-of-Gaussians model, with a uniform background [37], using the nuclear area and DAPI intensity (Figure 2b). The classification used seven Gaussians, equally spaced along the line connecting the center of the G1 and G2 peaks (Figure 2b), and an eighth Gaussian, located on the G1 peak. This classification of interphases

was confirmed by experiments using an anti- chromatin licensing and DNA replication factor 1 (CDT1) antibody to label G1 phase cells and 5-ethynyl-2'-deoxyuridine (EdU) to label S phase cells (results not shown). Figure 2b shows the distribution of the nuclear area versus the DAPI intensity, along with contours that denote regions containing cells in the G1, S, G2, and M phases. From this plot, we visualized the cell progression through the entire cell cycle. After the G1 phase (IDAPI = 1), cells enter the S phase, during which they duplicate their DNA (IDAPI = 1 → 2). Then, they enter the G2 phase and prepare for mitosis (IDAPI = 2). As they enter mitosis, they begin to condense, and their areas become smaller (IDAPI approximately 2.3, area approximately 250 µm<sup>2</sup> ). During anaphase, the image segmentation splits the two halves into distinct objects, resulting in an area and a DAPI intensity equal to half of the previous values (IDAPI approximately 1.1, area approximately 120 µm<sup>2</sup> ). The nuclei then start to grow and reenter the G1 phase.

#### 2.1.4. Damage Quantification

During the final step, the damage to each nucleus is defined as a fraction of the total damaged DNA, which can be computed as the ratio between the γH2AX and DAPI intensities, IγH2AX/IDAPI. To quantify the susceptibility of each cell phase to damage induced by the plasma jet, we used the damage ratio (DR), which indicates the number of cells with damage above a specified threshold divided by the total number of cells.

#### *2.2. APPJ Irradiation*

The above-described computational procedure for large-scale image analysis, with machine learning, was used to acquire the spatio-temporal distributions of DNA damage in malignant and nonmalignant cells after treatment with a nitrogen APPJ. The obtained images for two cell lines are presented, and the biological implications of plasma treatment are briefly discussed in the following subsections. These results are presented to illustrate the usefulness of the procedure for assessing the biological effects after plasma treatments, and certain conclusions have been derived; however, exploring the biological mechanisms induced by APPJs is not the focus of this study.

#### 2.2.1. APPJ-Irradiated Malignant Cells

The first row of Figure 3a shows the 2D distributions of cancer cell damage after different incubation times. These distributions were obtained after subtracting the median damage value of the flow control cells, which did not show any effects following nitrogen flow treatment alone. After 1-h incubation, cells were observed to be damaged near the plasma jet treatment region, and the damage decreased radially outward. Damage that extends beyond the plasma jet treatment region may be caused by the following: (1) the diffusion of the plasma species above the liquid surface; (2) the diffusion of the reactive species in the liquid, induced by the plasma treatment; (3) cell−cell communications, which may contribute to the damage of bystander cells near the treated cells [38]. Similarly, an enlarged affected area has been observed in previous studies [5,15,39].

More importantly, the maximum damage was observed at locations approximately 1 mm away from the treatment center, as shown in cancer cells after 2 h of incubation, resulting in a ring-shape pattern centered on the treatment location. The radius of this ring was similar to the inner radius of the quartz tube orifice, which was 1 mm (Figure 4), suggesting that the most prominent damaging effects were caused by species located at the interface of the plasma jet and the surrounding air. At this interface, highly reactive species have formed that interact with cells, causing DSBs. For example, one highly reactive species produced at the interface is nitric oxide (NO) [5,40], which can cause DNA strand breaks via the production of other reactive nitrogen species (RNS), such as ONOO−, HNO2, and N2O<sup>3</sup> [41]. NO can also trigger the production of intracellular reactive oxygen species (ROS), which can initiate various pathways, including apoptosis. Similar observations of ring-shaped regions have also been reported in previous studies, such as the inactivation patterns of bacteria, the distributions of reactive species, theoretical modeling [42], and for a gas-shield, helium-based APPJ,

**3. Materials and Methods** 

during interactions with cancer cells [43]. After 8 h of incubation, the measured damage decreased, which was likely due to cell detachment. Those cells that suffered from severe damage could not recover through repair processes and, consequently, undergo cell death. As a result, the dead cells detach from the coverslip and are removed during the washing steps of immunofluorescence procedures. Similar cell detachment after plasma treatment has been observed in other studies [44,45].

Furthermore, cancer cells in different cell-cycle phases have been found to respond with different sensitivities to plasma treatments. As shown in Figure 3a, after 1 h of incubation, almost 100% of S phase cells were damaged at the treatment region, whereas a much smaller proportion of G1 and G2 phase cells showed damage. In addition, the damaged area for S phase cells was found to be much larger than those for cells in the G1 and G2 phases. These results imply that S phase cells are more sensitive (susceptible) to damage induced by plasma treatment than cells in the G1 and G2 phases. During the S phase, DNA replication requires the exposure of a single-stranded portion of DNA near the replication fork, which results in an increased vulnerability to external attacks compared with other phases. Such attacks can be due to reactive species that are solvated/generated in the liquid above the cells during the treatment [46], which generally extend radially outwards beyond the location where the plasma jet makes direct contact with the liquid, through diffusion in the liquid phase or the transport of the gas phase plasma species across the liquid surface by the nitrogen gas flow. The damage caused during DNA replication may also be associated with the production of intracellular ROS/RNS, as a cellular response to plasma treatment [47]. *Int. J. Mol. Sci.* **2020**, *21*, x FOR PEER REVIEW 8 of 17

**Figure 3.** 2D distributions of the damaged DNA fraction in each nucleus (top row) and the damage ratios (DRs) for cells in G1, S, and G2 phases (rows 2–4), based on γH2AX staining in malignant cells (**a**) and nonmalignant cells (**b**) grown on coverslips. The damage threshold used to determine the DR was 75 (noting the damage range is 0–150). Each circular distribution, with a diameter of 18 mm, is centered on the treatment locations. Locations with white patches indicate regions without cells or blurry regions of the image, which were discarded. The circles for the cancer cells are cut off on the left side, due to a limitation of the slide scanner. **Figure 3.** 2D distributions of the damaged DNA fraction in each nucleus (top row) and the damage ratios (DRs) for cells in G1, S, and G2 phases (rows 2–4), based on γH2AX staining in malignant cells (**a**) and nonmalignant cells (**b**) grown on coverslips. The damage threshold used to determine the DR was 75 (noting the damage range is 0–150). Each circular distribution, with a diameter of 18 mm, is centered on the treatment locations. Locations with white patches indicate regions without cells or blurry regions of the image, which were discarded. The circles for the cancer cells are cut off on the left side, due to a limitation of the slide scanner.

gas speed of approximately 7.8 m/s. The gas flowed through a quartz tube, with an inner orifice of 2 mm. Oral cancer cells originally derived from squamous cell carcinoma of the tongue (SCC25) was obtained from American Type Culture Collection (ATCC, Rockville, MD, USA) and grown to approximately 90% confluence (approximately 105 cells·cm−2), on coverslips, in p35 dishes. Prior to plasma exposure, the cell culture medium in the p35 dish was replaced with 2.4 mL phosphatebuffered saline (PBS, 1 X, Mediatech. Inc., Manassas, VA, USA), forming a 3-mm-thick liquid layer above the cells. The detailed recipes for the cell culture medium used with cells in this study are listed in Appendix B. During the treatment, the dish was placed 2 cm below the quartz tube orifice (Figure 4). A plasma treatment time of 2 min was selected based on our previous study as the optimal treatment time for the detection of DNA damage without causing excessive buffer evaporation [5]. After 2 min of treatment, the PBS was replaced with 2 mL fresh culture medium, and the dishes were transferred to an incubator for 1, 2, 4, and 8 h. The samples were incubated at 37 °C, in an atmosphere of > 95% humidity and 5% CO2. Before image analysis, cells were fixed, permeabilized, and blocked, including PBS washes between each step, as described previously [5]. For imaging purposes, antiphospho-histone H2AX (Ser139) antibody (Mouse, EMD Millipore Corp.) and goat anti-mouse IgG (H+L) (Alexa Fluor 488 from Thermo Fisher Scientific Inc., Waltham, MA, USA) were applied, to evaluate the DSBs in nuclear DNA. Fluorescent DAPI Mounting Solution (Vector Laboratories Inc.,

A schematic diagram of the APPJ source and the experimental setup is illustrated in Figure 4. A detailed description of the plasma source used in this study has been previously reported [5]. The Appendix E.

provided in Appendices C and D, respectively.

and quantified.

label S phase cells.

Burlingame, CA, USA) was used to stain nuclear DNA, and then the DNA contents were visualized

We developed this large-scale image analysis method to facilitate cell-cycle classifications using only DAPI fluorescence. Cell-cycle classifications were confirmed by experiments using anti- CDT1 (Abcam plc.) antibody, which label G1 phase cells, and EdU (Thermo Fisher Scientific Inc.), which

To compare the effects of plasma treatments between two cell lines, the same preparation and treatment procedures were performed for oral nonmalignant cells. Telomerase reverse transcriptaseimmortalized normal oral keratiocytes (OKF6/T) cell were the generous gift of James Rhinewald, Brigham and Women's Hospital, Harvard Institutes of Medicine (Boston, MA, USA). Additionally, in the case of SCC25 and OKF6/T cells, two groups of samples were prepared using the same incubation times, including one group treated with nitrogen flow but no plasma irradiation and one

When the cells were not actively in culture, the cells were frozen and preserved for long-term storage. The procedures for cell culture and for freezing and thawing cells conducted in this work are

To obtain cellular data, the coverslips were scanned with an Aperio fluorescence slide scanner (Leica Microsystems, Wetzlar, Germany), using the DAPI and Alexa Fluor 488 (denoted as γH2AX)

group that received no treatment, which were used as control groups.

**Figure 4.** A schematic diagram of the nitrogen APPJ source used to treat cultured cells on coverslips. All dimensions are in mm. The inset shows the high voltage (HV) and current of alternating current (AC) waveforms of the nitrogen APPJ discharge. **Figure 4.** A schematic diagram of the nitrogen APPJ source used to treat cultured cells on coverslips. All dimensions are in mm. The inset shows the high voltage (HV) and current of alternating current (AC) waveforms of the nitrogen APPJ discharge.

**4. Conclusions**  We introduced a large-scale image analysis method, combined with machine learning, to study the spatial distributions of plasma jet-induced nuclear DNA damage over time in two cell lines from the oral cavity. The analysis included nuclear segmentation utilizing SALR clustering, the Additionally, the communication between treated cells and adjacent cells (i.e., bystander cells) can initiate cellular damage pathways in the bystander cells, similar to those initiated in treated cells [38]. Therefore, DNA damage in a cell layer can propagate radially outwards from the treatment region. Notably, G1 phase cells also showed circular damage patterns, with a damage ratio of approximately 50% after 1 h of incubation, whereas the damage ratio observed for G2 phase cells was much lower and not prominently localized at the plasma treatment region. These results are different from the observations reported previously [28], in which no significant γH2AX signals were found for either the G1 or G2 phases. These differences may be due to the different plasma sources used (i.e., helium plasma jet in the previous study [28], compared with the nitrogen plasma jet in our study), as the composition and distribution of plasma species may differ, resulting in different levels of damaging effects.

#### 2.2.2. APPJ-Irradiated Nonmalignant Cells

In contrast with the malignant cell line, plasma-treated nonmalignant cells did not show damage after 1 h of incubation (Figure 3b). After 2 h of incubation, the distribution of total damage in nonmalignant cells showed a circular pattern, with a much lower intensity than that observed in cancer cells (<1/5 the maximum value). This result suggested that under our experimental conditions, nonmalignant cells were minimally influenced by the plasma treatment. There are several factors: biological (e.g., gene expression, member structure, tolerance to oxidative stress) and experimental (e.g., dose and type of delivered radicals produced in APPJ from different plasma sources) factors reported that could contribute to different responses in various cell lines [14]. Because experimental conditions for both cell lines in this study are the same, the massive difference in damage susceptibility between two cell lines following plasma treatment is due to their biological differences. The observation that plasma-treated nonmalignant cells displayed mild damage at a later time (2 h of incubation

compared with 1 h of incubation for cancer cells) suggested that the two types of cells respond at different rates. Therefore, to accurately compare the plasma-induced effects on different types of cells, the cellular responses must be monitored over time, instead of using a single time point after treatment. However, more studies are needed to assess the APPJ effectiveness in biological targets and to determine the mechanisms through which plasma interacts with cells, which are outside of the scope of this study.

#### **3. Materials and Methods**

A schematic diagram of the APPJ source and the experimental setup is illustrated in Figure 4. A detailed description of the plasma source used in this study has been previously reported [5]. The APPJ was ignited with a 22.4-kV peak-to-peak voltage and a 59-mA peak-to-peak leaking current, at 28 kHz (Figure 4). Ultrahigh, pure, 5.0-grade nitrogen (purity of 99.999%, Airgas, Radnor, PA, USA) was used as the feed gas, with a flow rate of 1.5 standard liters per minute (slm), corresponding to a gas speed of approximately 7.8 m/s. The gas flowed through a quartz tube, with an inner orifice of 2 mm. Oral cancer cells originally derived from squamous cell carcinoma of the tongue (SCC25) was obtained from American Type Culture Collection (ATCC, Rockville, MD, USA) and grown to approximately 90% confluence (approximately 10<sup>5</sup> cells·cm−<sup>2</sup> ), on coverslips, in p35 dishes. Prior to plasma exposure, the cell culture medium in the p35 dish was replaced with 2.4 mL phosphate-buffered saline (PBS, 1 X, Mediatech. Inc., Manassas, VA, USA), forming a 3-mm-thick liquid layer above the cells. The detailed recipes for the cell culture medium used with cells in this study are listed in Appendix B. During the treatment, the dish was placed 2 cm below the quartz tube orifice (Figure 4). A plasma treatment time of 2 min was selected based on our previous study as the optimal treatment time for the detection of DNA damage without causing excessive buffer evaporation [5]. After 2 min of treatment, the PBS was replaced with 2 mL fresh culture medium, and the dishes were transferred to an incubator for 1, 2, 4, and 8 h. The samples were incubated at 37 ◦C, in an atmosphere of > 95% humidity and 5% CO2. Before image analysis, cells were fixed, permeabilized, and blocked, including PBS washes between each step, as described previously [5]. For imaging purposes, anti-phospho-histone H2AX (Ser139) antibody (Mouse, EMD Millipore Corp.) and goat anti-mouse IgG (H+L) (Alexa Fluor 488 from Thermo Fisher Scientific Inc., Waltham, MA, USA) were applied, to evaluate the DSBs in nuclear DNA. Fluorescent DAPI Mounting Solution (Vector Laboratories Inc., Burlingame, CA, USA) was used to stain nuclear DNA, and then the DNA contents were visualized and quantified.

We developed this large-scale image analysis method to facilitate cell-cycle classifications using only DAPI fluorescence. Cell-cycle classifications were confirmed by experiments using anti- CDT1 (Abcam plc.) antibody, which label G1 phase cells, and EdU (Thermo Fisher Scientific Inc.), which label S phase cells.

To compare the effects of plasma treatments between two cell lines, the same preparation and treatment procedures were performed for oral nonmalignant cells. Telomerase reverse transcriptase-immortalized normal oral keratiocytes (OKF6/T) cell were the generous gift of James Rhinewald, Brigham and Women's Hospital, Harvard Institutes of Medicine (Boston, MA, USA). Additionally, in the case of SCC25 and OKF6/T cells, two groups of samples were prepared using the same incubation times, including one group treated with nitrogen flow but no plasma irradiation and one group that received no treatment, which were used as control groups.

When the cells were not actively in culture, the cells were frozen and preserved for long-term storage. The procedures for cell culture and for freezing and thawing cells conducted in this work are provided in Appendices C and D, respectively.

To obtain cellular data, the coverslips were scanned with an Aperio fluorescence slide scanner (Leica Microsystems, Wetzlar, Germany), using the DAPI and Alexa Fluor 488 (denoted as γH2AX) channels, at 20× magnification. The immunofluorescence staining procedures are described in Appendix E.

#### **4. Conclusions** *Int. J. Mol. Sci.* **2020**, *21*, x FOR PEER REVIEW 10 of 17

We introduced a large-scale image analysis method, combined with machine learning, to study the spatial distributions of plasma jet-induced nuclear DNA damage over time in two cell lines from the oral cavity. The analysis included nuclear segmentation utilizing SALR clustering, the background and foreground correction of images, nuclei feature extraction, and cell-cycle classification. We demonstrated that the application of fluorescence microscopy using only DAPI staining (or any other dye for nuclear DNA staining) could achieve successful cell-cycle classification. This classification system not only preserves spatial information but also allows the discrimination of M phase cells from interphase cells, which flow cytometry often fails to accomplish. background and foreground correction of images, nuclei feature extraction, and cell-cycle classification. We demonstrated that the application of fluorescence microscopy using only DAPI staining (or any other dye for nuclear DNA staining) could achieve successful cell-cycle classification. This classification system not only preserves spatial information but also allows the discrimination of M phase cells from interphase cells, which flow cytometry often fails to accomplish. By using this analysis method, we were able to visualize that the 2D DNA damage distributions depend on the different cell-cycle stages at various incubation times. Therefore, to realize the efficacy

By using this analysis method, we were able to visualize that the 2D DNA damage distributions depend on the different cell-cycle stages at various incubation times. Therefore, to realize the efficacy of plasma treatments, post-treatment time-dependent assessments are necessary to monitor the cellular responses in various irradiated cancerous cell lines and their nonmalignant counterparts. Determining the spatial and temporal distributions of nuclear DNA damage in plasma-irradiated cells can have significant impacts on revealing the molecular mechanisms of plasma effects. The robustness and versatility of this method will be beneficial for a more systematic and rigorous assessment of the strengths, such as selective targeting features and weaknesses, such as side effects of plasma clinical applications. Furthermore, our proposed methodology can be used to quantitatively examine the effects of other types of radiation on biological samples. of plasma treatments, post-treatment time-dependent assessments are necessary to monitor the cellular responses in various irradiated cancerous cell lines and their nonmalignant counterparts. Determining the spatial and temporal distributions of nuclear DNA damage in plasma-irradiated cells can have significant impacts on revealing the molecular mechanisms of plasma effects. The robustness and versatility of this method will be beneficial for a more systematic and rigorous assessment of the strengths, such as selective targeting features and weaknesses, such as side effects of plasma clinical applications. Furthermore, our proposed methodology can be used to quantitatively examine the effects of other types of radiation on biological samples. **Author Contributions:** Conceptualization, methodology, investigation, validation, software, visualization, writing—original draft preparation, J.K. and X.H; resources, M.S.S and S.P; review and editing, E.A. and S.P;

**Author Contributions:** Conceptualization, methodology, investigation, validation, software, visualization, writing—original draft preparation, J.K. and X.H.; resources, M.S.S. and S.P; review and editing, E.A. and S.P; supervision, M.S.S., Y.L. and S.P. All authors have read and agreed to the published version of the manuscript. **Funding:** This material is based upon work supported by the U.S. Department of Energy Office of Science, Office of Basic Energy Science under Award Number DE-FC02-04ER15533. This is contribution number NDRL 5282

supervision, M.S.S, Y.L. and S.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This material is based upon work supported by the U.S. Department of Energy Office of Science, Office of Basic Energy Science under Award Number DE-FC02-04ER15533. This is contribution number NDRL 5282 from the Notre Dame Radiation Laboratory. from the Notre Dame Radiation Laboratory. **Acknowledgments:** The authors would like to acknowledge Guillermina Garcia from Sanford-Burnham Medical Research Institute (www.sbpdiscovery.org) for providing the whole slide fluorescence imaging.

**Acknowledgments:** The authors would like to acknowledge Guillermina Garcia from Sanford-Burnham Medical Research Institute (www.sbpdiscovery.org) for providing the whole slide fluorescence imaging. **Conflicts of Interest:** The authors declare no conflict of interest.

**Conflicts of Interest:** The authors declare no conflict of interest. **Appendices** 

#### **Appendix A. Extracted Features** *Appendix A: Extracted Features*

A list of the extracted features is given in Table A1, following a few notes. The number of Fourier descriptors for description nuclei shape was 20. Haralick textures [34] were computed using periods of 3 and 7 pixels and the features were averaged over the four computation directions (0◦ , 45◦ , 90◦ , and 135◦ ). The granularity spectrum, length of 7 pixels, was computed after removing the background (morphological opening using a disk with radius of 10 pixels) using a background structuring element with radius of 10 pixels. The radial intensity features used the distance transform of a nucleus segmentation mask to extract intensity features from rings around the nucleus with different widths and radii. We computed the radial intensity values for three rings with a width of 4 pixels, as well as the center region of the nuclei as shown in Figure A1. The only features extracted from the γH2AX channels were Intensity features. A list of the extracted features is given in Table A1, following a few notes. The number of Fourier descriptors for description nuclei shape was 20. Haralick textures [34] were computed using periods of 3 and 7 pixels and the features were averaged over the four computation directions (0°, 45°, 90°, and 135°). The granularity spectrum, length of 7 pixels, was computed after removing the background (morphological opening using a disk with radius of 10 pixels) using a background structuring element with radius of 10 pixels. The radial intensity features used the distance transform of a nucleus segmentation mask to extract intensity features from rings around the nucleus with different widths and radii. We computed the radial intensity values for three rings with a width of 4 pixels, as well as the center region of the nuclei as shown in Figure A1. The only features extracted from the γH2AX channels were Intensity features.

**Figure A1.** Radial intensity features. Example of the 3 rings with a width of 4 pixels and the center region. In each of these four regions, the mean and standard deviation of the intensity are computed. **Figure A1.** Radial intensity features. Example of the 3 rings with a width of 4 pixels and the center region. In each of these four regions, the mean and standard deviation of the intensity are computed.


**Table A1.** List of Extracted Features.

#### **Appendix B. Cell Culture Medium Recipes**

The preparation of cell culture medium is required to be carried out in a laminar flow hood in a sterile condition. All the components of the medium are required to be filtered. The filter used in this study is with a polyethersulfone microporous membrane which has a pore size of 0.2 µm.

*Appendix B.1. SCC25 Cell Culture Medium: ~500 mL*


## *Appendix B.2. OKF6*/*T Cell Culture Medium: ~500 mL*


#### **Appendix C. Cell Culture Procedure**

The cell culture procedures are required to be carried out in a laminar flow hood in a sterile condition.

#### *Int. J. Mol. Sci.* **2020**, *21*, 4127
