Emergence of a geometric pattern of cell fates from tissue-scale mechanics in the Drosophila eye

Pattern formation of biological structures involves the arrangement of different types of cells in an ordered spatial configuration. In this study, we investigate the mechanism of patterning the Drosophila eye into a precise triangular grid of photoreceptor clusters called ommatidia. Previous studies had led to a long-standing biochemical model whereby a reaction-diffusion process is templated by recently formed ommatidia to propagate a molecular prepattern across the eye epithelium. Here, we find that the templating mechanism is instead, mechanical in origin; newly born columns of ommatidia serve as a template to spatially pattern cell flows that move the cells in the epithelium into position to form each new column of ommatidia. Cell flow is generated by a pressure gradient that is caused by a narrow zone of cell dilation precisely positioned behind the growing wavefront of ommatidia. The newly formed lattice grid of ommatidia cells are immobile, deflecting and focusing the flow of other cells. Thus, the self-organization of a regular pattern of cell fates in an epithelium is mechanically driven.

Methods
Samples were imaged with a laser confocal microscope system at room temperature with z-stacks collected at precise 5 min intervals. We used  Leica Sp5 (wildtype replicate 1), Zeiss LSM 800 with no airy scanner (wildtype replicate 2), and Leica Sp8 (scabrous strong and weak mutants) systems. Wildtype replicate 1 was imaged on a Leica Sp5 with a detector size of 1024 x 1024 using the 100x/1.44NA objective, 3x internal zoom, and the pinhole set to 1 airy unit, resulting in a voxel size of 0.0505 x 0.0505 x 0.5003 µm and a xy field-of-view of 51.72 x 51.72 µm. The sample was scanned at 8000 Hz using the resonance scanner with bidirectional scanning and 16x line averaging. A total depth of 29.71 µm was collected using 60 steps of 0.5 µm each. The E-cadherin-GFP reporter was excited using the 488 nm channel of the argon laser (set to 25%) at 15% and collected on a HyD-PMT (GaAsP) detector with 300% gain. Wildtype replicate 2 was was imaged on a Zeiss LSM 800 with a detector size of 634 x 634 using a 63x/1.2NA objective, 2x internal zoom, and the pinhole set to 1.29 airy units, resulting in a voxel size of 0.08 x 0.08 x 0.7 µm and a xy field-of-view of 50.71 x 50.71 µm. The sample was scanned at 1000 Hz using the point scanner with bidirectional scanning and 4x line averaging. A total depth of 28.7 µm was collected in 42 steps of 0.683 µm each. The E-cadherin-GFP reporter was excited using the 488nm laser at 0.8% and collected on a GaAsP-Pmt detector using 652 V of gain. Both scabrous datasets were collected on a Leica Sp8 with a detector size of 704 x 704 using a 63x/1.4NA oil immersion objective, 2x internal zoom, and the pinhole set to 1 airy unit, resulting in a voxel size of 0.0829 x 0.0829 x 0.4 µm and a xy field-of-view of 61.60 x 61.60 µm. A total depth of 9.2 µm was collected in 24 steps of 0.38 µm each. The E-cadherin-GFP reporter was excited using the Diode 488nm laser at 0.8% and collected on a HyD-PMT (GaAsP) detector with 200% gain.

Note that the different wildtype replicates were imaged using different microscope systems and yet had highly reproducible behaviors for a wide spectrum of tissue and cellular features. Moreover, the scabrous mutant discs also exhibited highly reproducible behaviors similar to wildtype. Thus, different imaging platforms had little or no effect on the biological behaviors that we analyzed. All samples experienced photobleaching, evidenced by diminishing GFP signal over the course of imaging, and most likely experienced phototoxicity. We did not systematically test the effect of phototoxicity on the tissue. Generally, bleaching of the GFP reporter leveled off after about 3 or 4 hours of imaging. Laser powers were selected to generate sufficiently bright signal for segmentation after the point when GFP bleaching stabilized. It was not necessary to excite the GFP reporter up to the point of saturating pixel brightness, as is often done with fixed samples. Instead, we found that exciting the GFP reporter such that the signal was easily visible by eye was sufficient for cell segmentation.

Imaging resolution was selected to be sufficient for accurate cell segmentation within the morphogenetic furrow (MF). Cells in the MF are maximally constricted in their apical domains with the smallest cell diameters being around 0.5 µm. The imaging resolution of wildtype replicate 1 resulted in ~10 pixels along a 0.5 µm cell edge, whereas the imaging resolution of wildtype replicate 2 and both scabrous datasets led to ~6 pixels along a 0.5 µm edge. While we were able to segment wildtype replicate 2 and the scabrous datasets, wildtype replicate 1 was considerably easier to segment because of the higher resolution. We would advise future datasets to be collected at this resolution.

All samples were imaged about half-way between the margin of the eye disc and the equator. We did not specifically image dorsal only or ventral only regions. We found this region to be most mechanically stable. Often, the equatorial region of the disc would slowly sink basally over the course of imaging, requiring extremely large z-volumes to ensure the apical surface remained within the z-stack over the course of the movie. This resulted in excessive bleaching and diminishing signal intensity as the distance between the apical surface and cover slip increased. The average spacing of ommatidia changes as a function of dorsal-ventral (DV) position within the tissue, with tighter spacing near the equator and larger spacing near the margin. Imaging at the same relative DV position in discs from animals of approximately the same age was essential for keeping the spacing of ommatidia similar across datasets.

Image processing

We created an image processing pipeline to segment and track all cells within a chosen field of view (Figure S1). This pipeline involved the following steps in this order: 1) image collection, 2) surface projection, 3) image registration, 4) cell segmentation, 5) cell tracking. These steps are explained in more detail below.

Surface detection and projection

Cells were fluorescently labeled using an E-cadherin-GFP fusion protein that localizes to the adherens junction (AJ) surface of the eye disc epithelium. The E-cadherin-GFP protein also labels cells in the peripodial membrane, an epithelium of squamous cells located apical to the disc epithelium and separated from the disc epithelium by a fluid-filled lumen (Figure S1). In order to get a 2D projection of the AJ plane with signal-to-noise that permits downstream image analysis, and without inclusion of the peripodial membrane signal, we used the open-source MATLAB based software ImSAnE (Image Surface Analysis Environment) (Heemskerk and Streichan, 2015).

The first step in the surface detection workflow is the surfaceDetection module. We used the MIPDetector class of the planarDetector method of the surfaceDetection module. The MIPDetector and planarEdgeDetector classes are useful for finding approximately planar surfaces such as the AJ surface in imaginal discs. The MIPDetector relies on finding the brightest z-position for every xy coordinate in the gaussian derivative of the 3D image stack (i.e. after smoothing with a Gaussian filter). ImSAnE does this in a rescaled version of the 3D image stack where the pixels have been interpolated so that the physical voxels they represent have a unit aspect ratio of 1 (our image voxels were anisotropic with z-size ~10x larger than xy size, see section above). Because cells are more densely packed within the AJ surface than in the peripodial membrane surface, the AJ surface creates a more continuous and brighter structure in the Gaussian filtered 3D stack, which allows for detection via the MIPDetector method. This method additionally relies on exclusion of outliers (with regard to z-position) by looking at the distribution of z-positions within local xy neighborhoods. Representative parameters for the MIPDetector are sigma 5, channels 1, zdir -3, maxIthresh 0.05, summedIthresh 0.5, sigZoutliers 1, scaleZoutliers 50, and seedDistance 20. For some time points, it proved impossible to get the MIPDetector to specifically recognize the AJ surface of the disc epithelium without inclusion of portions of the peripodial membrane. This would happen at xy coordinates where there were dark interior of large cells (often dividing cells) in the AJ surface aligned with bright portions of the peripodial membrane signal. For these time points, we found that manually masking the peripodial membrane signal using ImageJ allowed us to isolate the AJ surface in ImSAnE.

The next step in the surface detection workflow is the surfaceFitting module. We used the tpsFitter method of the surfaceFitter class of the surfaceFitting module. Briefly, the tpsFitter fits thin plate splines to the surface produced by the MIPDetector, which is stored as a 3D point cloud. This method only works for surfaces found using the MIPDetector or planarEdgeDetector classes of the planarDetector method, which have a single z point for every xy coordinate. We set the smoothing parameter to 1000 and used a grid size of 80 x 80 for all datasets. The typical image size ranged from 600 x 600 to 1024 x 1024 pixels. This smaller grid size was chosen because it significantly sped up the spline fitting process and did not affect the quality of the fit. ImSAnE also allows the surface to be shifted in the z direction using the fitOptions shift method. We shifted the surface by ~10 pixels in the interpolated z-space, which corresponds to ~0.5 µm in physical space, in order to intersect the denser region of signal. The final step of the surfaceFitting module is creating a surface pullback using the pullbackStack method. We additionally used the onionOpts class of the pullbackStack method, which creates a ‘z-stack’ around the tpsFitter surface with a specified number of slices at a specified spacing in the interpolated z-stack space; we used 21 slices with a spacing of 1 for all datasets. The onionOpts class generates a max intensity projection of this ‘z-stack’ around the surface. This MIP was used as the final 2D projection of the image for segmentation and further analysis. Using onionOpts was integral to the surface projection workflow, as there are always regions where the surface isn’t perfectly fit to the brightest portions of AJ signal and onionOpts helps capture the entirety of the signal in these regions.

ImSAnE also includes a surfaceAnalysis module with a Manifold2D class that provides a map between the 3D surface and 2D projection and tools for making distortion free measurements of area and paths along the surface. However, we did not use this module of ImSAnE for our analysis because the image registration we used to handle tissue drift (see below) created a new reference frame for the segmented images that is distinct from the reference frame of surface projection. We tested the error resulting from distance measurements in the 2D projection without regard for the 3D curvature of the surface vs. distance measurements made using the surfaceAnalysis module. Curvature is greatest along the AP axis, with little curvature along the dorsal-ventral (DV) axis. For distance measurements around the length scale separating cell centroids (~1-2 µm) along the anterior-posterior (AP) axis, which is also about the average velocity cells move within the morphogenetic furrow (MF) over 1 hour, we found that disregarding 3D curvature introduced a ~1.5% error compared to making the same measurements with ImSAnE’s surfaceAnalysis module (data not shown). Measurements made in the DV direction, where there is very little curvature, had only a ~0.1% error. The curvature of the surface is greatest in the AP direction around the location of the MF, where we observe periodic cell flows in the anterior direction. Therefore, this error is affecting our measurements of the velocity of all cells in the MF region and the magnitude of these periodic cell flows. However, this measurement error is not affecting our analysis of the DV organization of these cell flows because there is very little measurement error along that axis.

Image registration and alignment of the AP/DV axes

After creating 2D surface projections using ImSAnE, we next eliminated tissue drift using image registration. The ex vivo tissue had a tendency to drift in the anterior direction while imaging. This was often most severe towards the start of the ev vivo culture. The drift was sometimes so severe that it required adjustment of the imaging field-of-view over the course of imaging, such as with wildtype replicate 1. This was to prevent the region-of-interest near the morphogenetic furrow (MF) from drifting out the field-of-view. In order to eliminate the effect of tissue drift from our analysis of cell movement, we tested a number of different intensity and feature based registration algorithms in MATLAB and ImageJ and found that the Linear stack alignment with SIFT plugin for FIJI (Schindelin et al., 2012) produced the most stable movie - i.e. the smoothest movie with the least amount of random jittering between consecutive time points. The Linear stack alignment with SIFT plugin in FIJI is a JAVA implementation of the Scale-Invariant Feature Transform method (Lowe, 2004). Prior to registration, we padded the images with zeros in order to create a buffer region so that no portion of the image would be translated out of the field-of-view during the registration process. We registered together consecutive time points chronologically in time using a Rigid transformation and the following parameters. For the Scale Invariant Interest Point Detector, we used a 1.6 pixel kernel for the gaussian blur, 3 steps per octave, 64 pixel minimum image size, and 1024 pixel maximum image size. For the Feature Descriptor, we used a feature descriptor size of 4, 8 feature descriptor orientation bins, and a closest/next closest ratio of 0.92. For the Geometric Consensus Filter, we used a maximum alignment error of 10 pixels and an inlier ratio of 0.05.

After removing tissue drift with image registration, the movie was translated such that the MF was parallel to the y-axis of the image. Since the MF is also parallel to the dorsal-ventral (DV) axis, this ensure that the DV and y axes were aligned, and the anterior-posterior and x axes were also aligned.

After creating 2D surface projections using ImSAnE and eliminating tissue drift with image registration, the final pre-processing steps before analysis were cell segmentation and tracking. Cell segmentation was achieved using a combination of a convolutional neural network (CNN) for pixel classification and MATLAB scripts for cell detection and tracking. Cell tracking relied on a MATLAB implementation of the Munkres assignment algorithm (Cao, 2021; Kuhn, 1955). Our code is publicly available (https://github.com/K-D-Gallagher/eye-patterning).

We trained a CNN to classify pixels as either cell edges or background/cell interior (https://github.com/K-D-Gallagher/CNN-pixel-classification) (Figure S1). We did this using a Pytorch package that provides transfer learning via pre-trained encoders paired with a variety of untrained decoder architectures that can be learned towards your specific pixel classification task (Yakubovskiy, 2020). This allowed us easily explore a variety of CNN architectures and find the most accurate one for our data.  We used a watershed transform to detect cells from the CNN pixel classification output and determined that our CNN was capable of accurately segmenting ~99.5% of cells relative to manually curate grown truth dataset (see below). However, when we attempted to track cells in datasets with ~0.5% error in cell segmentation, we were only able to accurately track ~80% of cells relative to our growth truth data. Errors in cell segmentation compound geometrically into errors in cell tracking. Therefore, we developed a custom MATLAB software to manually correct errors in segmentation. This software uses errors in cell tracking to discover the underlying errors in segmentation. Errors in cell segmentation were corrected using this software until no further errors in cell tracking could be detected.  After manually correcting cell segmentation, ~35% of segmented cells could be tracked through the entire course of the movie. The remainder of tracked cells either appeared/disappeared off the boundary of the field-of-view (FOV) (~55% of segmented cells) or appeared/disappeared within the interior of the FOV (~10% of segmented cells). This latter category we surmised to be cell birth and death events (see below). In sum, we could account for the behavior of 100% of segmented cells as either 1) being tracked throughout the entire duration of the movie (~35% of all segmented cells), 2) biological birth/death events (~6% / ~4% of all segmented cells), or 3) appearing/disappearing over a boundary of the FOV (~55% of all segmented cells).

There are no publicly available datasets for training CNNs to segment fluorescently labeled epithelial cells. Therefore, to train our CNN, we generated our own ground truth training set. This was done using Ilastik, an open-source machine learning (random forest classifier) based image processing software that offers segmentation, classification, and tracking workflows (Berg et al., 2019). Mistakes remaining after pixel classification in Ilastik were hand corrected using the above-mentioned MATLAB software. Relative to this hand corrected ground truth dataset, Ilastik was capable of accurately segmenting ~93% of cells (compared to ~99.5% with our CNN). The initial manually corrected ground truth training dataset for our CNN was approximately 120 images (one wildtype replicate) but progressively grew as new datasets were segmented and manually corrected. We also obtained confocal data of the Drosophila wing imaginal disc and thorax generously shared by Ken Irvine and Yohanns Bellaïche, respectively. Including this data in our training data not only increased the training library volume, which increases CNN model accuracy, but also increased the variability in the training dataset, which confers greater generalizability to new datasets not represented in the training data. Our trained CNN model is publicly available (https://github.com/K-D-Gallagher/CNN-pixel-classification).

The entire FOV was not segmented in order to reduce the necessary amount of manual correction of segmentation errors, as this was the most rate limiting step in our pipeline. We only segmented ~20 µm +/- the location of the morphogenetic furrow (MF), as this was our region-of-interest for analysis. Because the MF moves towards the anterior edge of the FOV, the number of segmented cells anterior to the MF progressively diminishes over the course of the movie. For example, the segmented FOV of wildtype replicate 1 begins with 732 cells (compared to ~1000 cells in the total image, including non-segmented cells) and ends with 458 segmented cells at the last frame of the movie; the reduced number of segmented cells is primarily to the anterior of the MF.

生物结构的模式形成涉及不同类型细胞在有序空间构型中的排布。本研究聚焦于果蝇复眼如何形成精确的光感受器簇（ommatidia，小眼）三角网格的模式形成机制。此前的研究提出了一种长期沿用的生化模型，认为反应扩散过程由新近形成的小眼作为模板，在眼上皮中传播分子预模式。而本研究发现，该模板机制实则起源于力学过程：新生的小眼柱作为模板，在空间上调控细胞流动，使上皮细胞迁移至对应位置以形成新的小眼柱。细胞流动由压力梯度驱动，该压力梯度源于恰好位于小眼生长波前后方的狭窄细胞扩张区域。已形成的小眼晶格细胞处于静止状态，偏转并汇聚其他细胞的流动。因此，上皮细胞命运的规则模式自组织过程是由力学机制驱动的。

### 方法
样本于室温下通过激光共聚焦显微镜系统（laser confocal microscope system）成像，以精准的5分钟间隔收集z-stack图像。我们使用了徕卡Sp5（野生型重复实验1）、蔡司LSM 800（无airy扫描器，野生型重复实验2）以及徕卡Sp8（scabrous强突变体和弱突变体）系统。
野生型重复实验1采用徕卡Sp5成像，探测器分辨率为1024×1024，使用100×/1.44NA物镜，3倍内部变焦，针孔设置为1艾里单位（airy unit），所得体素尺寸为0.0505×0.0505×0.5003 µm，xy视场范围为51.72×51.72 µm。样本以8000 Hz的共振扫描器进行双向扫描，采用16次线平均。共收集29.71 µm的总成像深度，分为60层，每层厚度0.5 µm。E-钙粘蛋白-GFP（E-cadherin-GFP）报告基因通过氩离子激光器的488 nm通道（功率设置为25%，实际激发功率15%）激发，使用HyD-PMT（GaAsP）探测器采集信号，增益设置为300%。
野生型重复实验2采用蔡司LSM 800成像，探测器分辨率为634×634，使用63×/1.2NA物镜，2倍内部变焦，针孔设置为1.29艾里单位，所得体素尺寸为0.08×0.08×0.7 µm，xy视场范围为50.71×50.71 µm。样本以1000 Hz的点扫描器进行双向扫描，采用4次线平均。共收集28.7 µm的总成像深度，分为42层，每层厚度0.683 µm。E-钙粘蛋白-GFP报告基因通过488 nm激光器激发（功率0.8%），使用GaAsP-PMT探测器采集信号，增益设置为652 V。
两组scabrous突变体数据集均采用徕卡Sp8成像，探测器分辨率为704×704，使用63×/1.4NA油浸物镜，2倍内部变焦，针孔设置为1艾里单位，所得体素尺寸为0.0829×0.0829×0.4 µm，xy视场范围为61.60×61.60 µm。共收集9.2 µm的总成像深度，分为24层，每层厚度0.38 µm。E-钙粘蛋白-GFP报告基因通过二极管488 nm激光器激发（功率0.8%），使用HyD-PMT（GaAsP）探测器采集信号，增益设置为200%。

需要注意的是，尽管不同野生型重复实验使用了不同的显微镜系统，但其在多种组织和细胞特征上均表现出高度可重复的行为。此外，scabrous突变体的眼盘也表现出与野生型高度相似的可重复行为。因此，不同的成像平台对我们分析的生物学行为几乎没有影响。
所有样本均发生了光漂白（photobleaching），表现为成像过程中GFP信号逐渐减弱，且大概率伴随光毒性（phototoxicity）。我们未系统测试光毒性对组织的影响。通常，GFP报告基因的漂白会在成像约3~4小时后趋于平缓。我们选择的激光功率足以在GFP漂白稳定后获得足够明亮的信号以用于细胞分割，无需如固定样本通常所做的那样将GFP激发至像素亮度饱和状态。实际上，我们发现将GFP激发至肉眼可轻松识别的信号强度便足以完成细胞分割。

成像分辨率的选择需确保形态发生沟（morphogenetic furrow, MF）内的细胞分割精度。MF内的细胞顶端区域处于最大收缩状态，最小细胞直径约为0.5 µm。野生型重复实验1的成像分辨率可在0.5 µm的细胞边缘上获得约10个像素，而野生型重复实验2及两组scabrous突变体数据集的成像分辨率仅可在0.5 µm的细胞边缘上获得约6个像素。尽管我们仍可对野生型重复实验2和scabrous突变体数据集进行分割，但野生型重复实验1因分辨率更高，分割难度显著更低。我们建议未来的数据集采用该分辨率进行采集。

所有样本均位于眼盘边缘与赤道区域的中点附近成像，未专门成像背侧或腹侧区域。我们发现该区域的力学稳定性最佳。通常，眼盘的赤道区域会在成像过程中逐渐向基底沉降，这需要极大的z轴成像体积才能确保顶端表面始终处于z-stack范围内，进而导致过度的光漂白和信号强度衰减，因为顶端表面与盖玻片之间的距离会逐渐增大。小眼的平均间距随组织内背腹轴（dorsal-ventral, DV）位置变化，赤道附近间距更紧密，边缘附近间距更大。在年龄相近的动物眼盘中，于相同相对背腹轴位置成像，是确保不同数据集中小眼间距保持一致的关键。

#### 图像处理
我们构建了一套图像处理流程，用于对选定视场内的所有细胞进行分割与追踪（图S1）。该流程按顺序包含以下步骤：1）图像采集；2）表面投影；3）图像配准；4）细胞分割；5）细胞追踪。下文将对各步骤进行详细说明。

##### 表面检测与投影
细胞通过定位在眼盘上皮黏着连接（adherens junction, AJ）表面的E-钙粘蛋白-GFP融合蛋白进行荧光标记。E-钙粘蛋白-GFP蛋白同时也标记了围腹膜（peripodial membrane）的细胞，该上皮为扁平细胞层，位于眼盘上皮顶端，通过充满液体的管腔与眼盘上皮分隔开（图S1）。为获得信噪比满足后续图像分析要求的黏着连接平面二维投影，同时排除围腹膜信号的干扰，我们使用了基于开源MATLAB的软件ImSAnE（Image Surface Analysis Environment）（Heemskerk和Streichan，2015）。

表面检测工作流的第一步为surfaceDetection模块。我们使用了该模块中planarDetector方法的MIPDetector类。MIPDetector与planarEdgeDetector类可用于检测近似平面的表面，如果蝇成虫盘的黏着连接表面。MIPDetector的工作原理是：对3D图像堆栈进行高斯滤波后，在其高斯导数图像中，为每个xy坐标找到最亮的z位置。ImSAnE会在缩放后的3D图像堆栈中执行该操作，其中像素已被插值，使所代表的物理体素的长宽比为1（我们的图像体素各向异性，z轴尺寸约为xy轴的10倍，详见前文）。由于黏着连接表面的细胞密度高于围腹膜表面，因此在高斯滤波后的3D堆栈中，黏着连接表面会形成更连续、更明亮的结构，这使得MIPDetector方法可实现有效检测。该方法还会通过分析局部xy邻域内的z位置分布，排除离群值（相对于z位置）。MIPDetector的代表性参数为：sigma 5、channels 1、zdir -3、maxIthresh 0.05、summedIthresh 0.5、sigZoutliers 1、scaleZoutliers 50、seedDistance 20。
在部分时间点，若不包含围腹膜的部分区域，MIPDetector无法特异性识别眼盘上皮的黏着连接表面。这种情况常发生在黏着连接表面的大细胞（通常为分裂细胞）的暗内部与围腹膜信号的明亮区域对齐的xy坐标处。对于这些时间点，我们使用ImageJ手动掩膜围腹膜信号，从而在ImSAnE中分离出黏着连接表面。

表面检测工作流的下一步为surfaceFitting模块。我们使用了surfaceFitter类的tpsFitter方法。简言之，tpsFitter将薄板样条（thin plate splines）拟合至MIPDetector生成的、以3D点云形式存储的表面。该方法仅适用于使用planarDetector方法的MIPDetector或planarEdgeDetector类生成的表面，这类表面对每个xy坐标仅有一个z点。我们将平滑参数设置为1000，并为所有数据集设置80×80的网格尺寸。典型图像尺寸范围为600×600至1024×1024像素。选择较小的网格尺寸可显著加快样条拟合速度，且不会影响拟合质量。ImSAnE还允许通过fitOptions shift方法在z方向偏移表面。我们将表面在插值后的z空间中偏移约10个像素，对应物理空间中的~0.5 µm，以确保与信号更密集的区域相交。
surfaceFitting模块的最后一步为使用pullbackStack方法创建表面拉回（surface pullback）。我们还使用了pullbackStack方法的onionOpts类，该类可在tpsFitter生成的表面周围创建一个“z堆栈”，包含指定数量的切片和指定间距的插值z堆栈空间；我们为所有数据集设置了21个切片，间距为1。onionOpts类会对该表面周围的“z堆栈”进行最大强度投影，该最大强度投影图像将用作最终的二维投影图像，用于细胞分割和后续分析。使用onionOpts是表面投影工作流的关键，因为总会存在表面无法完美匹配黏着连接信号最亮区域的情况，而onionOpts可帮助捕获这些区域的全部信号。

ImSAnE还包含一个surfaceAnalysis模块，其中的Manifold2D类提供了3D表面与二维投影之间的映射，以及用于对表面区域和路径进行无畸变测量的工具。但我们并未将该模块用于本研究的分析，因为我们用于处理组织漂移的图像配准（详见下文）为分割后的图像创建了新的参考框架，该框架与表面投影的参考框架不同。我们测试了不考虑表面3D曲率的二维投影距离测量与使用ImSAnE的surfaceAnalysis模块进行的距离测量之间的误差。表面曲率沿前后轴（anterior-posterior, AP）最大，沿背腹轴几乎无曲率。对于沿前后轴、长度尺度约为细胞中心间距（~1~2 µm）的距离测量（这也约等于细胞在形态发生沟内1小时内的平均移动速度），我们发现忽略3D曲率相较于使用ImSAnE的surfaceAnalysis模块进行测量，仅引入了约1.5%的误差（未展示数据）。而沿背腹轴的距离测量误差仅约0.1%，因为该轴方向几乎无曲率。表面曲率在形态发生沟附近的前后轴方向最大，此处我们观察到沿前方的周期性细胞流动。因此，该误差会影响我们对形态发生沟区域内所有细胞的速度及这些周期性细胞流动幅度的测量，但不会影响我们对细胞流动背腹轴组织模式的分析，因为该轴方向的测量误差极小。

##### 图像配准与前后/背腹轴对齐
通过ImSAnE完成二维表面投影后，我们接下来使用图像配准消除组织漂移。离体组织在成像过程中倾向于向前方漂移，这种现象在离体培养初期通常最为严重。漂移有时会非常显著，以至于需要在成像过程中调整成像视场，如野生型重复实验1，以防止形态发生沟附近的感兴趣区域漂移出视场。为了消除组织漂移对细胞运动分析的影响，我们在MATLAB和ImageJ中测试了多种基于强度和特征的配准算法，最终发现FIJI的Linear stack alignment with SIFT插件（Schindelin等，2012）可生成最稳定的成像序列——即连续时间点间抖动最少的最平滑影片。该插件是尺度不变特征变换（Scale-Invariant Feature Transform, SIFT）方法的JAVA实现（Lowe，2004）。
在配准前，我们使用零值对图像进行填充以创建缓冲区域，确保在配准过程中不会有图像区域被平移出视场。我们按时间顺序对连续时间点进行配准，使用刚性变换及以下参数：对于尺度不变兴趣点检测器，我们使用1.6像素内核的高斯模糊、每个八度3个步骤、最小图像尺寸64像素、最大图像尺寸1024像素；对于特征描述符，我们使用尺寸为4的特征描述符、8个特征描述符方向箱、0.92的最近/次近邻比率；对于几何一致性滤波器，我们使用最大对齐误差为10像素、内点比率为0.05。

完成图像配准消除组织漂移后，成像序列被平移，使形态发生沟与图像的y轴平行。由于形态发生沟同时也与背腹轴平行，这可确保背腹轴与y轴对齐，前后轴与x轴对齐。

在通过ImSAnE完成二维表面投影并通过图像配准消除组织漂移后，分析前的最终预处理步骤为细胞分割与追踪。细胞分割结合了用于像素分类的卷积神经网络（Convolutional Neural Network, CNN）与用于细胞检测和追踪的MATLAB脚本。细胞追踪基于MATLAB实现的Munkres分配算法（Cao，2021；Kuhn，1955）。我们的代码已开源（https://github.com/K-D-Gallagher/eye-patterning）。

我们训练了一个CNN以将像素分类为细胞边缘或背景/细胞内部（https://github.com/K-D-Gallagher/CNN-pixel-classification）（图S1）。我们使用了一个Pytorch工具包，该工具包通过预训练编码器结合多种未训练的解码器架构提供迁移学习，可针对特定的像素分类任务进行学习（Yakubovskiy，2020）。这使我们能够轻松探索多种CNN架构，并为我们的数据找到最准确的模型。我们使用分水岭变换从CNN的像素分类输出中检测细胞，并确定我们的CNN相对于手动整理的地面真值（ground truth）数据集（详见下文）可准确分割约99.5%的细胞。但当我们尝试对细胞分割误差约为0.5%的数据集进行细胞追踪时，相对于地面真值数据，我们仅能准确追踪约80%的细胞。细胞分割的误差会以几何方式累积为细胞追踪的误差。因此，我们开发了一款自定义MATLAB软件用于手动校正分割误差，该软件利用细胞追踪中的误差来识别潜在的分割错误。我们使用该软件校正细胞分割错误，直至无法再检测到细胞追踪误差。手动校正细胞分割后，约35%的分割细胞可在整个成像序列中被追踪。其余被追踪的细胞要么从视场（FOV）边界出现/消失（约占分割细胞的55%），要么在视场内部出现/消失（约占分割细胞的10%）。我们推测后一类情况属于细胞的诞生与死亡事件（详见下文）。综上，我们可以解释100%分割细胞的行为：1）在整个成像过程中均可被追踪（约占所有分割细胞的35%）；2）生物学意义上的诞生/死亡事件（约占所有分割细胞的6%/4%）；3）在视场边界出现/消失（约占所有分割细胞的55%）。

目前尚无公开可用的数据集可用于训练CNN以分割荧光标记的上皮细胞。因此，为了训练我们的CNN，我们生成了自己的地面真值训练集。我们使用Ilastik完成该操作，这是一款基于开源机器学习（随机森林分类器）的图像处理软件，提供分割、分类和追踪工作流（Berg等，2019）。使用上述MATLAB软件手动校正Ilastik像素分类后残留的错误。相对于该手动校正的地面真值数据集，Ilastik可准确分割约93%的细胞（而我们的CNN可达约99.5%）。我们用于训练CNN的初始手动校正地面真值数据集约包含120张图像（一个野生型重复实验），随着新数据集的分割和手动校正，数据集规模逐渐扩大。我们还获得了由Ken Irvine和Yohanns Bellaïche分别慷慨共享的果蝇翅成虫盘和胸部共聚焦数据。将该数据加入训练数据集不仅增加了训练库的规模，提升了CNN模型的准确性，还增加了训练数据集的多样性，使模型对训练数据中未涵盖的新数据集具有更好的泛化能力。我们训练的CNN模型已开源（https://github.com/K-D-Gallagher/CNN-pixel-classification）。

为减少所需的手动分割误差校正工作量（这是我们流程中最耗时的步骤），我们并未对整个视场进行分割。我们仅对形态发生沟周围约±20 µm的区域进行了分割，因为这是我们分析的感兴趣区域。由于形态发生沟向视场的前缘移动，在成像过程中，MF前方的分割细胞数量会逐渐减少。例如，野生型重复实验1的分割视场在初始时包含732个细胞（总图像中约有1000个细胞，包括未分割的细胞），在成像序列的最后一帧时仅剩458个分割细胞；分割细胞数量减少主要集中在MF前方区域。