Intelligent Fluorescence Image Analysis of Giant Unilamellar Vesicles and Protein Liquid Phase Droplets by Whole Z-stack Analysis

Fluorescence image analysis in biochemical science often involves the complex tasks of identifying samples for analysis and calculating the desired information from the intensity traces. Analyzing giant unilamellar vesicles is one of these tasks. Researchers need to identify many vesicles to statistically analyze the degree of molecular interaction or state of molecular organization on the membranes. This analysis is complicated, requiring a careful manual examination by researchers, so automating the analysis can signicantly aid in improving its eciency and reliability. We developed an intelligent analysis routine based on the 3D information of whole z-stack images. The program identies the valid vesicles to analyze and calculates the desired data automatically. The program can examine the amount of protein binding on the membranes and determine the state of domain phase separation by calculating the uorescence intensity trace along the membranes. We also show that the method can easily be applied to similar analyses, such as intensity analysis of phase-separated protein droplets. A deep learning-based classication approach enables the identication of vesicles even from relatively complex samples. We demonstrate that the proposed articial intelligence-assisted classication can further enhance the accuracy of the analysis close to the performance of manual examination. using a denition for the statistical GUVs. contour intensity of each z-section image of vesicles and uorescence intensity of that used the partition coecient the reporter


Introduction
A giant unilamellar vesicle (GUV) is an arti cial lipid membrane system characterized by its micrometerwide globular shape. Because GUVs have the advantage of deformability and a relatively large size, they have been used in a variety of biophysical and biochemical experiments. 1 Example studies include the study of lipid membrane phase separation, 2-4 phase modulation by membrane anchored proteins, 5,6 reconstitution of membrane remodeling processes, 7,8 the viral assembly process, 9 measurement of mechanical tension, 10,11 effect of osmotic pressure, 12 membrane protein reconstitution, 13 and monitoring protein binding to the membrane. 14 Fluorescence imaging and uorescence image analysis are extensively used in these types of experiments. The most common mode of analysis is manually studying each individual vesicle to examine enough vesicles for statistical analysis from a population of GUVs. 5 Although this may be suitable for experiments with a relatively small number of GUVs, for larger-scale sample numbers, manual analysis is very time consuming and often di cult to implement, particularly when researchers are not very experienced with the nature of the GUV system and uorescence imaging. Automated analysis can provide an effective solution to these problems by standardizing the analysis procedure. GUV image analysis is complicated enough that it requires well-designed algorithms to automate the analysis to the comparable level that can be performed by a trained researcher, and there has been some effort to address similar problems by other groups. [15][16][17][18] In 2014, Hermann et al. introduced the circular Hough transformation (CHT) algorithm for the automated detection of GUVs in uorescence images. 15 CHT is an algorithm originally developed to detect circular objects from complex images in computer science. 19 Brie y, CHT rst converts the original image into an edge image to clearly de ne the edge of objects. Basing on the edge image, it starts the voting process to score the likelihood of having a circle with radius r at any given pixel position. The pixels at the center of circular objects appear as high-score pixels in the transformed image space; thus, the position of circles with radius r can be decided based on likelihood. If radius r is to be exibly decided, radius r also becomes a parameter of the likelihood voting process, which makes it computationally more expensive.
Fluorescence images of GUVs usually show very clear boundaries of GUV contours, so CHT is an excellent method to automatically detect individual GUVs or any other similar circular objects from the uorescence image. CHT has been widely used and improved, so it is also readily available to apply in commonly used mathematical software and programing languages.
In this report, we extended the CHT detection algorithm by developing a deep learning-based approach for the selection and intensity analysis of full GUV image stacks, which are a collection of image sections along the z-axis, including the information of full 3D GUVs. An intensity analysis algorithm was developed to study protein binding and phase separation state decisions. We applied automated analysis methods to examine several typical sets of data. In the binding analysis, we show how his-tagged protein bound to Ni-DGS lipid can be detached by competitive binding of imidazole solution. In the phase state analysis, we present how the existence of a separated lipid phase can be detected by contour intensity analysis. In addition to the GUV analysis, we applied the CHT detection algorithm to automatically analyze uorescence images of liquid phase-separated protein droplets. Proteins with multivalent binding interaction or intrinsically disordered domain interaction at a high concentration spontaneously form protein droplet domains that may serve as independent organelles in our cells. [20][21][22][23][24][25][26][27][28][29][30] Given the appreciation of thermodynamic phase separation as a general principle of protein sorting in the eld, in vitro reconstitution experiments of protein droplets have become very common these days. Protein droplets have sizes and shapes comparable to those of GUVs, so a similar automated analysis may be introduced for the high-throughput analysis of protein droplet domains. We show how the CHT algorithm could be successfully applied to protein domain analysis, and we also present an example analysis of differences in protein cargo concentration at different conditions.
Lastly, we show how a deep learning approach can be leveraged for the problem as a complementary and potentially more capable solution. Deep learning is a process in which an arti cial network is trained with sample images, so it can recognize different images without speci cally designing and extracting features. It has been proven to be powerful in descriptor extraction and classi cation for visual recognition tasks. 31 Because of recent progress in deep learning, its capability and e ciency reach even beyond human recognition, and it has been applied to a wide range of related elds in leading-edge science and engineering. 32-40 As our work involves the classi cation of images, we integrated this approach into our algorithm for enhanced performance. We show an e cient implementation of deep learning in vesicle classi cation for choosing the right vesicles.
Method CHT-based GUV stack analysis The rst step of the analysis is to identify the GUVs and apply selection lters in order to choose GUVs with the right condition. In uorescence imaging, the images of 3D objects are collected as z-stacks or a series of image sections at different heights. In some simple analyses, examining single section images may be su cient, but whole stack analysis ensures that we are studying the entire 3D information of GUVs, and some ltering strategy and analysis are only possible when we have access to the whole 3D information.
Z-stacks are analyzed from the rst to the last images one by one, as shown in Figure 1a. In each image, the CHT algorithm is applied to detect GUV uorescence circles. The sensitivity parameter of CHT detection can be adjusted to change the strictness of circle de nition, and the radius range can also be set to only detect GUVs with a reasonable size. Once all GUV uorescence circles are detected from each section of the image, 3D GUVs are identi ed by grouping uorescence circles by their center positions. All uorescence circles with center (x, y) positions within the set range are considered as one GUV and are grouped as an individual GUV entity. Each grouped GUV is recorded as a matrix containing the position, radius, and intensity of each image section that can be used for further analysis.
GUV selection ltering strategies are applied in two separate stages, as shown in Figure 1b. The rst ltering is applied after detecting GUV uorescence circles from each image CHT. In this stage, the contrast between the circular edge and the internal space is used to lter out ill-de ned circles. The intensity per pixel is calculated for pixels near edges and pixels inside the circle, and any circles with insu cient relative intensity near edges are ltered out. This way, the images of multi-lamellar vesicles or vesicles with too much uorescence content can be ltered out. GUV circles that are not clearly analyzable using the intensity contour because of ambiguous positions can be removed in this stage, as well. The second ltering strategy is done after grouping individual GUVs. In this stage, any GUVs with insu cient information are ltered out. Speci cally, GUV groups with only one or two image sections are ltered out because they do not have enough images to be considered as an intact GUV.
The result of the initial CHT detection and ltering is displayed and recorded visually, so users can adjust the selection parameters to use until the outcome shows the most reliable detection result for the purpose of analysis.

Intensity binding analysis
The purpose of intensity binding analysis is to quantify the change in uorescence intensity on the membrane in a uorescence channel for any membrane binding species. Reference channel A, in which the uorescence signal comes from the lipid membrane, is used to detect and group GUVs to analyze; analysis channel B, in which the uorescence signal comes from interacting species, such as membrane binding proteins, is analyzed for changes in the uorescence intensity signal on the membrane.
For all uorescence images that passed selection strategies, the matching uorescence intensity of channel B is analyzed using the position and radius information obtained from channel A. The uorescence intensity of all pixels adjacent to the circle de ned by the center position and radius is summed up, as shown in Figure 2a. The total intensity values and intensity per pixel values are recorded for all images. The CHT algorithm used in this study detects each circle in such a way that the radius is slightly greater than the actual intensity maxima of the uorescence circle, so counting pixels that are n pixels, by user de nition, inward successfully sums up all intensity values. As each GUV is grouped, the intensity per pixel from each GUV can be calculated.
To account for the background signal, the net intensity is de ned as follows in Eq. 1: This means that the general background signal per pixel was estimated and subtracted from the raw intensity signal. The background signal can be estimated by measuring the average signal coming from the dark part of the images. To make this process more systematic, we introduced automated methods to the program. In the automated estimation, a background signal is estimated based on the analysis of the intensity histogram of each image. Users may choose to use the most common intensity value or the triangular thresholding estimation value. 41 Both worked ne for most analyses we tried, but triangular thresholding is recommended, which uses the most probable intensity values from both background and net uorescence signals.
Binary masking for high-sensitivity analysis of low signal-to-noise ratio images Background intensity estimation was not only used to correctly calculate the net intensity signal but also to create a binary mask of the image to analyze. We learned that creating a binary mask, which is an image version in which the original image is converted into 1 or 0 image based on the estimated threshold intensity, signi cantly enhanced the performance of circular detection in CHT, especially when the signal-to-noise ratio of the image is relatively low. To make the CHT method generally applicable to images with any signal-to-noise ratio, binary masking was used for all the methods described. Successful examples of circular detection using binary masking are shown in Figure 2b, c.
Phase separation determination using intensity contour analysis The purpose of phase separation determination analysis is to systematically identify the binary phase state of GUVs, homogeneous or phase separated, using a well-de ned de nition for the statistical analysis of a population of GUVs. In addition to the phase determination, the analysis also calculates and reports the contour intensity of each z-section image of vesicles and the uorescence intensity of different phases that may be used to calculate the partition coe cient of the reporter molecule.
The analysis rst detects individual GUVs using the general CHT circular detection and ltering strategies to nd qualifying GUVs to analyze. Example GUV detection results are shown in Figure 3a, b for a typical uniform GUV sample with most of GUVs in a uniform phase state and a phase-separated GUV sample with most of GUVs in a phase-separated state. Binary masking CHT detection followed by the GUV stack analysis successfully detected GUVs to analyze in typical phase separation samples. Intensity contour analysis uses the intensity around the circular perimeter of each section image to make a decision about phase separation. If the GUV is phase separated, the intensity around the perimeter will show discrete discontinuity at the phase boundary, as shown in Figure 3c. If the GUV is uniform without phase separation, the intensity plot will not have any noticeable discontinuity around the perimeter, as shown in Figure 3d. As this decision is made for each section of a GUV, the decision parameter can be set to make a nal state decision for each GUV. For example, a GUV may be determined as a phase-separated vesicle when at least 60% of the z-section images show clear phase separation.
For contour intensity analysis, the intensity around the perimeter was divided into N segments. The average net intensity per pixel was calculated for each segment by using the pixels within (r-∆r, r) and (θ, θ + 360̊ / N), where r and θ indicate the radius and the angle in a polar coordinate of a circle image, respectively. For each contour intensity segment, histogram analysis is performed by aligning segments by average intensity from the highest intensity to the lowest one. The segment intensity of the 20% from the highest is de ned as high intensity, and the segment intensity of the 20% from the lowest or 80% from the highest is de ned as low intensity. These are estimated intensities of the uorescence intensity of two different phase domains. The mid-intensity between the high intensity and the low one is calculated, and discontinuity analysis is performed along the direction of the increasing degree θ to determine whether there is any directional decrease or increase in intensity spanning ±p % from the mid-intensity (Figure 3c, d). The p value may be adjusted based on the partition coe cient of the uorescent report used. For a high partition coe cient reporter, which strongly prefers one domain over the other, a higher p may be used. For a low partition coe cient reporter with a relatively even distribution but with a clear preference for one domain over the other, a smaller p value should be used. High and low intensity estimation from the 20% percentile may be adjusted if necessary, but choosing a too small percentile is not recommended, as it may select a false high or low intensity coming from an outlier intensity. Choosing an overly large percentile may fail to correctly select high and low intensity in the histogram distribution when the relative proportion of one domain is small. Through an analysis of the number of discontinuities in the intensity contour, the number of domains may be estimated if needed, when nely modulated phase behavior with many domains per vesicle is expected, as opposed to a binary segregation with only two domains.

Protein droplet intensity analysis
Protein droplet phase separation is another phase separation behavior that is commonly studied with an in vitro reconstitution experiment. Because of the comparable size and globular shape of the protein droplets to the GUVs, the same analysis method can be applied to automatically calculate the uorescence intensity of each droplet. CHT analysis can detect circular images from z-section images, and it can be reconstituted to nd individual droplets. The uorescence intensity of each droplet is analyzed by calculating the net average intensity per pixel of all the pixels within the circles that can be used to quantify the amount of protein in the droplet. 21 Background intensity is subtracted to correctly calculate the net intensity of the droplets, as speci ed in Eq. 1. Some droplets that are at the stage of two droplets coalescing with each other may be detected as two droplets that are very close to each other. Such cases can be avoided by ruling out all overlapping droplets, or they may just be analyzed as individual droplets.
Deep learning-based image lter The purpose of a deep learning-based image lter is to enhance the capability of the GUV image selection strategy so that the program can recognize even more complicated cases of images that a simple calculation algorithm might miss. The convolutional neural network (CNN), a powerful and reliable network of deep learning, was used in this study to achieve this objective. Using the training images and de ned classes, the CNN is trained by minimizing the mean squared error between the inputs and outputs. The CNN is mainly structured with the normalization layer, the convolutional layer, the recti ed linear unit (ReLu) layer, the pooling layer, and the fully connected layer. During the training process, the sample images will be inputted to the normalization layer, which performs image normalization. The convolutional layer is applied for feature extraction. The ReLu layer will introduce nonlinearity to the CNN model. The pooling layer will reduce the image spatial dimension. The fully connected layer is used for high-level reasoning in the neural network. A total of 5,853 training images of individual vesicles were used. The images were converted into 50 × 50 pixels as inputs for the CNN. Each image represented a possible z-section image of a GUV. The images were classi ed into four classes. In the training process, these four classes were de ned as C1-valid vesicles, C2-multi-lamellar vesicles, C3-vesicles overlapped with other vesicles, and C4-hazy or unclear images. A total of 1,918 training images of C1, 877 of C2, 1,129 of C3, and 1,929 of C4 were used. The neural network was trained and validated using the training set. Then, the trained CNN model takes unlabeled images as the input and generates the corresponding classes. In this new version of the phase determination program, instead of the vesicles being ltered based on intensity calculation, each detected vesicle image was converted into 50 × 50 images for classi cation by the trained neural network. Only images of C1-valid vesicles passed the enhanced ltering strategy. The rest of the calculation was identical to the previously explained algorithm for vesicle state determination.

Sample preparation
GUV was prepared with the electroformation method. 42,43 Brie y, a lipid mixture of a desired composition was deposited on an indium tin oxide (ITO)-coated glass (Delta Technologies) at 55°C. The lipid was dried by nitrogen gas and further dried in a vacuum chamber. The glass, along with another ITO glass, formed a sealed chamber separated by a silicon spacer (Grace Bio-Labs), and sucrose solution was introduced. Sinusoidal voltage of 2 V/mm space at 5 Hz was applied for about 2 h at 55°C to induce GUVs. The GUVs were collected and used in one day. The lipids used in this study were from Avanti Polar Lipids. The proteins utilized in this research were puri ed by E. coli overexpression, followed by Ni a nity puri cation and gel ltration. The plasmids were given by Michael Rosen (Addgene plasmid # 127093, #126946, and #127093). 21 The imaging samples were prepared similar to the method described in 5 . Brie y, an AttoFluor cell chamber (Invitrogen) and a cover glass cleaned by bath sonication in isopropyl alcohol:water = 1:1 were assembled. The glass surface was blocked by incubating with 5 mg/mL BSA solution for 30 min. After incubation, the chamber was washed ve times with the buffer for the addition of the GUV solution for imaging. To analyze protein interaction, the desired concentrations of proteins were added by pipette injection and mixed gently. To ensure a quick and homogeneous interaction, the protein solution was added at a volume comparable to that of the solution in the chamber (20%-50% by volume). Protein droplet samples were rst mixed in a test tube to incubate for at least 30 min, and they were introduced to the chamber for imaging. To decrease the volume of the chamber, a silicon O-ring was introduced.

Imaging condition
All the images shown in this report were collected using confocal uorescence laser scanning microscopy unless speci ed otherwise. Brie y, a Nikon Ti-E-based C2 confocal microscope was used. Excitation laser lights of 488 and 561 nm were used with matching emission lters to collect signals from the uorescent molecules. A Nikon Plan Apo 100× NA 1.45 oil immersion objective was used without further magnifying the lens in the optical path. The typical mode of scanning was to collect data as 1,024 × 1,024 pixels spanning a 127.3 µm × 127.3 μm area, whereas motorized z-axis movement allowed the automated acquisition of z-stack images.
Software platform used GNU Octave, a freely available scienti c programming language, was used to test all the non-deep learning algorithms shown in this report. The algorithm and the strategy we tested in this report are generally applicable in any language, but we chose this platform because it is freely available to researchers, the language is intuitive and approachable for students in training to learn, and exact procedures of the subroutines used can be easily examined. The version we used in this project was completed with GNU Octave 5.1.0 and its image package. A deep learning-based algorithm was developed and implemented in MATLAB R2020b (MathWorks) with image processing and deep learning toolboxes. All raw codes in this article are available as supplementary materials with example images.

Results
Detecting the unbinding of his-tagged uorescent proteins from the membrane by imidazole inhibition We tested multi-channel intensity-based binding analysis by examining data from relatively simple binding and unbinding experiments of uorescent proteins on the GUV membrane. In this experiment, histagged green uorescence protein (GFP) was bound to the GUV with functionalized Ni-DGS lipids by incubating the sample with the protein. The protein was then detached by introducing a high concentration (300 mM) of imidazole, without altering the osmolality of the buffer, to detach the proteins from the membrane via inhibition of Ni-his tag interaction. We analyzed GFP uorescence intensity on the membrane before and after introducing the inhibitor. CHT detection found circular uorescence using the TR-DHPE uorescent lipid signal, and qualifying GUVs after ltering strategies were used to quantify the uorescence signal from the GFP protein on the membrane in the protein uorescence channel. Multiple GUVs under the same conditions were quanti ed, and the average value was calculated. As shown in Figure 4a, the intensity-based binding analysis successfully quanti ed the difference before and after introducing the imidazole inhibitor to the sample. The amount of his-tagged GFP protein on the membrane after introducing the inhibitor was less than 5% of the intensity before introducing the inhibitor, suggesting that the majority of proteins that were originally bound to the membranes were removed from the membrane as a result of inhibition.

Statistical analysis of phase state
We tested the contour intensity-based phase state determination program by analyzing samples with relatively clear differences in phase behavior. GUVs with uniform phase were prepared as a 1,2-dioleoylsn-glycero intensity coming from a double-layered GUV with two bilayers very closely spaced to each other. Furthermore, contour analysis could not analyze the information of the image at the very top and at the bottom, as it shows phase state information as phase behavior happening on the plane instead of at the perimeter. This may be valuable information for human researchers, but this version of contour analysis will not be able to use this information. Nevertheless, the performance of the program was generally reliable and unbiased.

Protein droplet intensity analysis at different cargo concentrations
We tested the protein droplet intensity analysis program by analyzing protein droplet samples prepared at different concentrations of uorescent cargo proteins. As shown in Figure 5a, the same approach of the automated analysis of CHT-based circle detection, followed by the whole stack analysis, selected individual protein phase droplets very effectively. The example analysis was performed on protein droplets formed between two multivalent binding partners, SUMOx10 repeats and SIMx10 repeats, at about 10 µM concentrations each. GFP-SUMO3, a uorescent protein with three repeated SUMO domains, was introduced as a cargo protein that reported the protein droplet clearly. The cargo protein concentration was varied at a relatively low concentration range from 10 to 100 nM. At all concentrations, droplet detection was e cient, and uorescence intensity showed a monotonic increase from 20 to 100 nM concentration regions, although it did not show a noticeable difference in intensity from 10 to 20 nM. (Figure 5b) Deep learning-based GUV selection strategy for enhanced vesicle identi cation The deep learning-based approach uses CNN trained by the sample images. The sample images consist of pre-classi ed classes of images as the input. The CNN is trained until it can successfully classify images that were not part of the training sample images. This approach uses arti cial intuition by pattern recognition of the objects, so it is potentially more capable when the objects to classify are di cult to de ne with a simple calculation-based strategy. Figure 6a shows the typical images in each class of images we used to train the neural network. Class 1 represents the valid vesicles to analyze; other classes represent vesicles that should be excluded for various reasons. Class 2 represents mutilamellar vesicles or vesicles that contain unexpected lipid membrane structures inside; class 3 represents vesicles that are too crowded, which makes the analysis of an individual vesicle ambiguous; and class 4 represents vesicles with a too low image quality to analyze mostly because the contour of the lipid membrane is not fully in focus.
When the trained network was used to detect vesicles images, the ltering strategy very strictly excluded vesicles that had any features represented by classes 2-4, as shown in Figure 6b-c. Some vesicle images had multiple features of classes 2-4, such as a multi-lamellar vesicle with a hazy intensity trace. Therefore, some vesicles could not be classi ed uniquely, but the vesicles to be ltered almost always fell under one of classes 2-4, allowing only valid vesicle images of class 1 to survive the selection strategy.
Combined with the whole z-stack approach, the likelihood of invalid vesicles surviving the automated selection strategy was very small, and valid GUVs were successfully selected. This approach was especially bene cial to exclude some cases that were ill de ned by simple calculation. Multi-lamellar vesicles that have multilayers only near the outermost vesicle membranes and crowded vesicles in which a part of the membrane is touching other lipid membranes are examples.

Discussion
We conducted an intelligent uorescence image analysis based on whole z-stack images of GUVs. Our approach is unique in that it can effectively select the valid vesicles to analyze from a population of vesicles and perform uorescence intensity-based calculation on various pieces of information, such as the amount of lipid-protein interaction and the phase state of vesicles. The proposed automated selection and standardized calculation can also facilitate numerous large-scale and complex image analyses and statistics in both science and engineering areas. Our method can detect vesicles from relatively low signal-to-noise ratio samples, and it can easily be expanded to various applications of multi-channel intensity analysis. We also showed that the method of automated detection can be used for protein droplet analysis, which we believe will be especially useful when studying the interaction between lipid vesicles and protein phase droplets. Deep learning-based classi cation could successfully recognize vesicle types for the purpose of excluding unwanted vesicles. To the best of our knowledge, this is the rst implementation of deep learning to the problem. Deep learning is potentially a very powerful approach because it can classify very complicated classes as long as an arti cial neural network is well trained. However, its performance strongly depends on the quality and amount of training image sets, so a reasonable de nition of classes with image samples spanning many variations of each class will be important for the training and validation accuracy. As typical image analysis in the eld requires a process in which researchers identify each individual entity and make a decision on the class of each, the application of the deep learning approach can greatly improve the e ciency and reliability of the analysis. A deep learning-based decision is also fast, and in the long term, it can output class decision information in real time as the data are being collected, providing a dynamic feedback platform to enhance the effectiveness and quality of the experiment.