A platform for automated and label-free monitoring of morphological features and kinetics of spheroid fusion

Spheroids are widely applied as building blocks for biofabrication of living tissues, where they exhibit spontaneous fusion toward an integrated structure upon contact. Tissue fusion is a fundamental biological process, but due to a lack of automated monitoring systems, the in-depth characterization of this process is still limited. Therefore, a quantitative high-throughput platform was developed to semi-automatically select doublet candidates and automatically monitor their fusion kinetics. Spheroids with varying degrees of chondrogenic maturation (days 1, 7, 14, and 21) were produced from two different cell pools, and their fusion kinetics were analyzed via the following steps: (1) by applying a novel spheroid seeding approach, the background noise was decreased due to the removal of cell debris while a sufficient number of doublets were still generated. (2) The doublet candidates were semi-automatically selected, thereby reducing the time and effort spent on manual selection. This was achieved by automatic detection of the microwells and building a random forest classifier, obtaining average accuracies, sensitivities, and precisions ranging from 95.0% to 97.4%, from 51.5% to 92.0%, and from 66.7% to 83.9%, respectively. (3) A software tool was developed to automatically extract morphological features such as the doublet area, roundness, contact length, and intersphere angle. For all data sets, the segmentation procedure obtained average sensitivities and precisions ranging from 96.8% to 98.1% and from 97.7% to 98.8%, respectively. Moreover, the average relative errors for the doublet area and contact length ranged from 1.23% to 2.26% and from 2.30% to 4.66%, respectively, while the average absolute errors for the doublet roundness and intersphere angle ranged from 0.0083 to 0.0135 and from 10.70 to 13.44°, respectively. (4) The data of both cell pools were analyzed, and an exponential model was used to extract kinetic parameters from the time-series data of the doublet roundness. For both cell pools, the technology was able to characterize the fusion rate and quality in an automated manner and allowed us to demonstrate that an increased chondrogenic maturity was linked with a decreased fusion rate. The platform is also applicable to other spheroid types, enabling an increased understanding of tissue fusion. Finally, our approach to study spheroid fusion over time will aid in the design of controlled fabrication of “assembloids” and bottom-up biofabrication of living tissues using spheroids.


Supplementary
Radial frequency (M) is 3 and angular frequency (N) is 12, resulting in 52 features. These features were extracted for the binary image of the doublet and both separated spheroid regions. Therefore, 156 (3x52) features were obtained per sample. The feature extraction was implemented using a Matlab function [2].

157:159 Roundness
Minor/major axis length of the ellipse that has the same normalized second central moments as the region. Computed for the whole and both separated regions. Based on the maximally inscribed circle method (described in Figure 3), multiple features relating to the overlap were calculated: -Intersection over Union (IoU).
-Number of overlap pixels between both circles respectively divided by the number of pixels in circle 1, circle 2 and the merged circles.
-For the whole and both separated regions, the number of region pixels not overlapping with their corresponding inscribed circles, divided by the number of pixels in that region.

175
Contact ratio The contact length (as described in Figure 3) over the averaged spheroid width.

176:178
Focus measure Diagonal Laplacian [3]. The measure was implemented in Matlab by [4], but adapted to only consider the specified region. Computed for the whole and both separated regions.

179:187 Correlation measure
The gray-level co-occurrence matrix was computed with the following offsets: [0 1; 0 2; 0 3]. The graycoprops (Matlab) function was used to extract the 'correlation' measure. Computed for the whole and both separated regions.  Figure 3.

Main features Description/extraction method [Unit]
File index Frame number of the image sequence.

Doublet area
The number of pixels in the doublet region [Pixel] (regionprops -'Area').

Doublet length
The average number of pixels around the center of mass of the doublet along the horizontal axis [Pixel].

Doublet width
The average number of pixels around the center of mass of the doublet along the vertical axis [Pixel]. The doublet width is only considered after a doublet roundness of 0.8 is reached.

Doublet roundness
The doublet roundness was computed according to the following formula: Minor/major axis length (regionprops -'MinorAxisLength'/'MajorAxisLength') of the ellipse that has the same normalized second central moments as the doublet region. The value ranges from 1 (circle) to approximately 0 (infinite line segment).

Contact length
The length of the contact region, i.e. interface between the two spheroids [Pixel]. * Since every doublet is composed of two spheroids, this feature is calculated for both spheroids. ** Can only be determined up to 180 degrees. A small time interval for imaging will limit the rotation of the doublet between consecutive frames.   Figure 2 (A) The original image was cropped around the detected microwell, resulting in the raw image with the doublet candidate in the center of the image. (B) The doublet was segmented from the raw image by combining (OR operation) Roberts gradient operator (threshold = 0.045) [5] and intensity thresholding (derived from Otsu's threshold) [6]. The gradient output was not normalised. The extracted region was overlaid (AND operation) with the binary micro-well mask (radius of 79 pixels), holes smaller than 600 pixels (4-connectivity) were filled and the mask was eroded with a disk of 10 pixels radius. In order to remove small regions, regions smaller than 300 pixels were filtered out and the mask was dilated with a disk of 10 pixels radius. Finally, the largest object was retained.

(G)-(H)
The transformed mask was refined with respect to the previous (frame t) binary reference mask. Matching and small non-matching (up to 20 pixels) holes were filled, while larger non-matching holes were excluded. Moreover, 'missed' regions (with respect to the reference mask) larger than 25 pixels were extracted, 'smoothed' with the doublet border, and it was attempted to correct for under-segmentation based on the connectivity of these 'missed' regions with the doublet mask. Regions that were substantially (larger than 80%) connected to the doublet mask were retained and added to the doublet mask. Next, the previous operations (B) on erosion, removal of small regions and dilation were performed. Only the largest object (i.e. the doublet) was retained. The final doublet mask was re-aligned with the horizontal axis (based on the fitted ellipse) and replaces the binary and grayscale reference mask of the previous time-point in (C). The cycle was repeated for the next time point. Figure 3 Stage (1):

Detailed information on
The area of the doublet mask was extracted (regionprops -'Area'). An ellipse was fitted with the same normalized second central moments as the doublet region, and its minor and major axis length (regionprops -'MinorAxisLength'/'MajorAxisLength') were used to compute the doublet roundness. Next, the orientation (regionprops -'Orientation') of the ellipse was used for alignment with the horizontal axis, and the doublet length was computed by averaging the lengths (along the horizontal axis) within a 2-pixel neighbourhood of the center of mass.
(C) The doublet was separated in an upper and lower half (along the horizontal axis) and all minima were detected for both profiles. Their location and prominence was stored.

(D)
The aligned doublet region was vertically splitted in two separate (spheroid) regions and by maximising the minimum roundness of the separated regions, a rough separation was obtained. Next, a maximally inscribed circle was computed for each region, maximising the overlap with the respective region. The points where both circles intersected were also included as profile points in (C). In case the circles initially did not intersect, they were both grown until intersection. The prominence of these points was set at 1.

(E)
In order to obtain the final separation of the spheroids, the anchor (contact) point on the upper and lower half of the doublet was selected. All possible pairs of anchor points were considered and the pair with the lowest cost, according to a cost function based on the spheroids their circularity, prominence of the anchor points and shift between the upper and lower anchor point, was selected. The contact length was computed as the Euclidean distance between both anchor points. After separation, both spheroid areas (regionprops -'Area') and widths were computed. The spheroid widths were computed by averaging the lengths (along the vertical axis) within a 2-pixel neighbourhood of the center of mass of the spheroid region.
(F) The upper and lower intersphere angle were measured through extraction of the local spheroid border (4 in total, 2 up and 2 down), starting from the anchor points of the contact length. The number of pixels considered for each border was 40% of the horizontal distance from the anchor point to the center of mass of the respective spheroid, with a minimum number of 10 pixels. Next, a linear function was fitted to each local border, with the anchor point as a fixed point. In this way, the linear functions were used to determine the upper and lower intersphere angle.

Stage (2):
In order to determine the anchor points in the subsequent frames, the doublet was separated in an upper and lower half (along the horizontal axis, as in (C)) and all minima (with a minimum prominence of 2) were detected for both profiles. Their location and prominence was stored. Based on their prominence and horizontal distance to the averaged anchor point coordinate (based on previous frames), certain minima were retained. In addition to this, the averaged anchor points (upper and lower half) were also retained. All possible pairs (upper with lower half) were matched and the pair with the lowest Euclidian distance (minimal contact length) was selected as the new anchor points. The other features were extracted as described in (1).