Method for predicting euploidy based on blastocyst development kinetics and spherical harmonic decomposition fusion

Abstract

This invention discloses an euploidy prediction method based on the fusion of blastocyst developmental dynamics and spherical harmonic decomposition. The method acquires multifocal plane image sequences of blastocysts using a time-difference imaging system and extracts dynamic parameters. After semantic segmentation and depth estimation, a three-dimensional surface model is obtained through surface reconstruction. The trophoblast cell instances are segmented to obtain three-dimensional centroid coordinates. These centroid coordinates are radially projected onto a unit sphere, and a spherical density function is constructed after excluding the inner cell mass mask. Spherical harmonic decomposition is then normalized using Monte Carlo zero-model normalization to obtain normalized power spectra at each degree. The basic morphological features, dynamic parameters, and spherical harmonic features are fused and filtered before training an ensemble learning model to output prediction results. This invention is the first to introduce spherical harmonic analysis into blastocyst assessment, capturing multi-scale spatial distribution information that global statistics cannot obtain. The power spectrum exhibits rotational invariance, and zero-model normalization eliminates cell number confounding.

Description

Technical Field

[0001] This invention belongs to the field of biomedical engineering technology, specifically relating to an euploidy prediction method based on the fusion of blastocyst development dynamics and spherical harmonic decomposition. Background Technology

[0002] Embryo selection is crucial in in vitro fertilization-embryo transfer (IVF-ET) treatment. In clinical practice, embryologists need to select the individuals with the highest developmental potential from multiple transferable embryos to improve implantation and pregnancy success rates. Embryo chromosomal ploidy is one of the key factors affecting implantation outcomes, and aneuploid embryos are a major cause of implantation failure and early miscarriage. Preimplantation genetic testing-aneuploidy screening (PGT-A) is currently the gold standard method for assessing embryonic chromosomal ploidy status; however, this method requires invasive biopsy sampling of trophoblast cells from the blastocyst, posing a potential risk of embryo damage, and is also expensive. Whether it should be routinely used in IVF treatment remains controversial.

[0003] In search of alternatives, various non-invasive assessment methods have emerged in recent years. Deep learning methods based on two-dimensional time-lapse images utilize convolutional neural networks to directly extract features from single-frame or video images of blastocysts for ploidy classification. A representative work is the STORK-A system, which achieves a predictive performance of approximately 0.70 AUC by combining a ResNet18 architecture with maternal age, kinetic parameters, and morphological scores. Machine learning methods based on kinetic parameters construct predictive models using embryonic development time points recorded by time-lapse culture systems (such as cleavage time, morula formation time, and blastocyst expansion time). However, these methods all rely on two-dimensional images or temporal dimension information and fail to fully exploit the three-dimensional spatial structural features of the blastocyst.

[0004] The introduction of 3D reconstruction technology has opened up a new dimension for blastocyst evaluation. Existing studies have involved rotating the blastocyst to obtain multi-angle images, stitching them together on a sphere to form a 3D surface model, and then extracting 3D morphological parameters such as trophoblast (TE) cell count, TE density, and blastocyst cavity volume, establishing the correlation between these parameters and euploidy. Other studies have attempted to directly perform 3D reconstruction using multifocal plane images from time-lapse culture systems, avoiding the physical manipulation of rotating the embryo and demonstrating advantages in clinical compatibility.

[0005] However, the aforementioned three-dimensional feature extraction methods share a common limitation: they primarily focus on global statistics, such as mean, standard deviation, density, and volume, while neglecting the spatial distribution pattern of TE cells on the three-dimensional surface of the blastocyst. In fact, the distribution of TE cells on the blastocyst surface is not uniform. Blastocysts naturally exhibit polarity differentiation; polar TE cells facing the endometrium differ from wall-borne TE cells on the opposite side in cell density and functional state. In aneuploid embryos, TE cells often exhibit abnormal spatial distribution patterns such as local clustering, empty areas, and density mutations. For example, two blastocysts with the same number of TE cells and area variation coefficient may present different spatial distributions—one with a uniform TE distribution, and the other with obvious local clustering—which existing global statistical features cannot distinguish between.

[0006] In summary, the existing technology lacks a method to extract multi-scale features of TE cell spatial distribution from the three-dimensional structure of the blastocyst without damaging the embryo, and to fuse these features with traditional morphological parameters and developmental dynamics information to predict chromosome euploidy. Summary of the Invention

[0007] This invention proposes an euploidy prediction method based on the fusion of blastocyst developmental dynamics and spherical harmonic decomposition, to address the problem in existing technologies where blastocyst three-dimensional morphology assessment only uses global statistics and ignores the spatial distribution pattern of trophoblast cells.

[0008] To address the aforementioned technical problems, this invention provides an euploidy prediction method based on the fusion of blastocyst developmental dynamics and spherical harmonic decomposition, comprising the following steps: Step S1: Obtain image sequences by taking multi-focal plane images of the blastocyst using a time-difference imaging system, and extract developmental kinetic parameters from the image sequences; Step S2: Semantically segment and depth estimate the image sequence sequentially, then reconstruct the surface to obtain a three-dimensional surface model. Perform instance segmentation on the trophoblast cells on the three-dimensional surface model to obtain the coordinates of each three-dimensional centroid, mark the boundary of the inner cell mass region, and extract basic morphological features. Step S3: Project the coordinates of the three-dimensional centroid of the three-dimensional surface model radially to the unit sphere with the centroid of the three-dimensional surface model as the origin. After excluding the inner cell mass region as a mask, construct the spherical density function in the effective trophoblast region using kernel density estimation. Normalize the spherical density function by Monte Carlo zero model normalization after spherical harmonic decomposition of the spherical density function to obtain the normalized power spectrum of each degree. Step S4: Combine the basic morphological features, the developmental kinetic parameters, and the normalized power spectrum into an extended feature set and then filter them; Step S5: Train an ensemble learning model using the selected feature set, and use the ensemble learning model to output holoplasticity prediction results for the blastocyst under test.

[0009] Preferably, the developmental kinetic parameters in step S1 include the time of pronucleus emergence. Pronuclear disappearance time The start time of each stage from the two-cell stage to the nine-cell stage to 1. Start time of mulberry embryonic stage Time of blastocyst cavity appearance The cavity expands to half the size of the blastocyst in time. and the time for complete cavity expansion .

[0010] Preferably, the basic morphological features in step S2 include size and volume features, cytological features, morphological parameters, and spatial relationship features; wherein the size and volume features include the total volume of the blastocyst, the volume of the cavity, the volume of the inner cell mass, and the radius of the blastocyst; the cytological features include the number of trophoblast cells and the density of trophoblast cells; the morphological parameters include the mean area, standard deviation of area, mean perimeter, mean major axis, mean minor axis, standard deviation of the ratio of major to minor axis, and blastocyst sphericity; and the spatial relationship features include the eccentricity of the inner cell mass.

[0011] Preferably, in step S2, the depth estimation uses an improved Laplacian operator as the focal metric function. The focal metric curve formed by each pixel along the focal plane direction of the image sequence is fitted to obtain the depth value by sub-pixel-level peak localization. The instance segmentation adopts a two-dimensional segmentation and three-dimensional mapping strategy. The trophoblast cells are segmented into two-dimensional instances on the best focal plane with the highest focal metric response and its adjacent focal planes. After deduplication and merging according to the intersection-union ratio greater than or equal to 0.5, the segments are then mapped to the three-dimensional centroid coordinates through the depth map.

[0012] Preferably, the kernel function for kernel density estimation in step S3 is... The expression is: ; in Let be the unit direction vector of the query point. Let be the unit direction vector of the centroid of the trophoblast cell. This is a concentration parameter; the concentration parameter pass Confirmed, among which The number of trophoblast cells, It is an adjustment coefficient with a value range of 8 to 15.

[0013] Preferably, step S3, when excluding the inner cell mass region as a mask, further includes: A width of [width value] is set outside the mask boundary. A linearly gradient transition zone of curvature, within which the value of the spherical density function linearly decays to zero from the effective nourishment layer region side toward the mask side.

[0014] Preferably, the Monte Carlo null model normalization method in step S3 includes: performing spherical harmonic decomposition on the spherical density function to obtain the original power spectrum at each degree, wherein the degree... In spherical harmonic decomposition, a non-negative integer representing the spatial scale, with values ​​ranging from 0 to the maximum degree. ; Generate at least 100 groups randomly and evenly distributed on the effective nutrient layer area. A set of simulated points, of which The number of the trophoblast cells; Perform the same kernel density estimation and spherical harmonic decomposition as the observed data on each set of simulated points, and calculate the degree of all simulated sets. The average power spectrum on the desired power spectrum ; Observation data in degrees The original power spectrum Divide by the desired power spectrum Obtain the normalized power spectrum : ; in A value greater than 1 indicates a degree The distribution variation of the trophoblast cells at the corresponding spatial scale exceeds the statistical expectation of a random uniform distribution.

[0015] Preferably, in step S3, after obtaining the normalized power spectrum... This is followed by the calculation of the Spatial Organization Index (SOI) and its logarithmic transformation as supplementary features; the SOI is defined as a degree. Normalized power spectral components when the value is 0 with degrees From 1 to The ratio of the sum of the normalized power spectral components at time: ; The logarithmic transformation is Its calculation method is as follows , where ln is the natural logarithm; The higher the SOI value, the more uniform the spatial distribution of the trophoblast cells.

[0016] Preferably, the maximum degree of the spherical harmonic decomposition in step S3 is... The value range is 4 to 10. The spherical meshing adopts the HEALPix scheme of equal area hierarchical pixelation and the pixelation parameters are... Greater than or equal to 32.

[0017] Preferably, the screening in step S4 adopts a two-stage method: the first stage uses statistical tests to retain features with a significance level less than a preset significance threshold, and the second stage uses a gradient boosting model to rank the features retained in the first stage by importance, and determines the optimal number of features based on the inflection point of the cumulative importance curve. The ensemble learning model described in step S5 uses the LightGBM gradient boosting decision tree algorithm.

[0018] The beneficial effects of the present invention include at least the following: First, this invention decomposes the distribution of TE cells on the blastocyst surface into components of different spatial frequencies using spherical harmonic power spectroscopy, with each degree corresponding to a specific biological meaning: degree =1 reflects the hemispherical asymmetry of the blastocyst (the difference between polar TE and wall TE), degree =2 reflects the quadrupole deformation mode, degree. =3 to 6 reflect local clustering or empty zone patterns. This information cannot be provided by existing global statistical features, enabling multi-scale quantification of the spatial distribution patterns of trophoblast cells.

[0019] Second, the area covered by the inner cell mass (ICM) on the surface of the blastocyst naturally lacks TE cells. If left untreated, this area would be misjudged as a "gap" in TE distribution by spherical harmonic analysis. This invention excludes the ICM region before constructing the density function, ensuring that the spherical harmonic features reflect the abnormal distribution of TE cells themselves rather than the structural gaps caused by the presence of ICM, thus eliminating the interference of boundary effects.

[0020] Third, the number of TE cells in different blastocysts naturally varies (usually between 60 and 100). When the number of cells is low, the statistical noise increases, leading to a systematic increase in the high-frequency power spectrum. By dividing the observed power spectrum by the expected power spectrum that is randomly and uniformly distributed under the condition of the same number of cells, the normalized power spectrum eliminates this confounding factor. This invention effectively controls the systematic confounding effect of cell number on the power spectrum through Monte Carlo null model normalization, making blastocysts with different cell numbers comparable.

[0021] Fourth, this invention combines spherical spatial distribution features with three-dimensional basic morphological features and developmental dynamic parameters into an extended feature set, which is then input into the ensemble learning model after two-stage screening. This results in improved prediction performance compared to using any single category of features alone.

[0022] Fifth, the entire method utilizes only multifocal plane images acquired during routine culture using a time-lapse culture system, eliminating the need for additional physical manipulation or fluorescent labeling of the embryos and completely avoiding the risk of biopsy damage. From image analysis, 3D reconstruction, feature extraction to model prediction, the entire process can be automated, reducing inconsistencies in embryologists' subjective assessments. Attached Figure Description

[0023] Figure 1 This is a schematic diagram of the overall process of the method in an embodiment of the present invention; Figure 2 This is a schematic diagram of the blastocyst three-dimensional reconstruction process according to an embodiment of the present invention; Figure 3 Box plot of statistical analysis of TE cell morphological parameters in an embodiment of the present invention; Figure 4 This is a schematic diagram of the spherical parameterization of the blastocyst surface according to an embodiment of the present invention; Figure 5 This is a schematic diagram of the ICM mask and transition zone processing according to an embodiment of the present invention; Figure 6 This is a schematic diagram comparing the spherical density function and the normalized power spectrum in an embodiment of the present invention; Figure 7 This is a schematic diagram illustrating the LightGBM feature importance analysis according to an embodiment of the present invention; Figure 8 This is a schematic diagram of SHAP feature attribution analysis according to an embodiment of the present invention; Figure 9 This is a schematic diagram comparing the ROC curves of six ensemble learning models according to an embodiment of the present invention; Figure 10 This is a schematic diagram showing a comprehensive comparison of the ROC curves of various models in the same coordinate system according to the embodiments of the present invention. Detailed Implementation

[0024] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the protection scope of the present invention.

[0025] like Figure 1As shown, the overall process of the method of the present invention includes: Step S1: acquiring image sequences and extracting developmental dynamics parameters through a time-of-flight imaging system; Step S2: semantic segmentation and depth estimation followed by surface reconstruction and segmentation of trophoblast cell instances to obtain three-dimensional centroid coordinates and basic morphological features; Step S3: radially projecting the three-dimensional centroid coordinates onto a unit sphere to construct a spherical density function and obtaining normalized power spectra at each degree through spherical harmonic decomposition and Monte Carlo null model normalization; Step S4: combining the basic morphological features, developmental dynamics parameters, and normalized power spectra into an extended feature set and filtering it; Step S5: using the filtered feature set to train an ensemble learning model to output euploidy prediction results for the blastocysts under test. The following examples are based on 325 human blastocyst samples with complete time-of-flight imaging data and PGT-A detection results, including 158 euploid cases and 167 aneuploid cases.

[0026] Step S1: Obtain image sequences by taking multi-focal plane images of the blastocyst using a time-difference imaging system, and extract developmental dynamic parameters from the image sequences.

[0027] The purpose of this step is to acquire imaging data of the entire blastocyst development process and extract dynamic parameters that reflect the developmental rhythm.

[0028] During data acquisition, a culture device equipped with a time-lapse imaging system is used to continuously monitor embryos during blastocyst culture (days 1 to 5 or 6 post-fertilization). The imaging system takes multi-focal plane images of each embryo at fixed time intervals (usually every 5-20 minutes). Each image captures 5 to 11 focal plane images with an interlayer spacing of 10 to 25 μm, covering cross-sectional information at different depths of the embryo. The accumulated time-series image sequence throughout the entire culture cycle provides a complete record of the embryo's development from fertilized egg to blastocyst. The time-lapse imaging system employs the Huffman modulation contrast (HMC) optical principle, and the safety of this imaging method has been widely validated in assisted reproductive clinical practice, showing no phototoxicity or other physical damage to the embryo.

[0029] In the kinetic parameter extraction stage, key time nodes in the embryonic development process are identified from the aforementioned time-series image sequences to form a set of kinetic parameters. This invention uses a total of 14 kinetic parameters, specifically defined as follows: pronuclear appearance time. The time when the pronucleus first becomes visible after fertilization; the time when the pronucleus disappears. This refers to the moment when prokaryotic fusion disappears; the start time of each stage from the two-cell stage to the nine-cell stage. to The time when the embryo reaches the stage of 2 to 9 blastomeres was recorded; the start time of the morula stage was recorded. That is, the moment when each blastomeres are tightly compacted to form the morula morphology; the time when the blastocyst cavity appears. This refers to the moment when a identifiable cavity first appears in the center of the embryo; the time when the cavity expands to half the size of the blastocyst. That is, the moment when the cavity area expands to about half of the embryo's cross-section; the time when the cavity fully expands. This refers to the moment when the blastocoel is fully expanded and the zona pellucida is significantly thinned. These parameters can be obtained from time-varying images through manual annotation or automatic recognition algorithms.

[0030] In the 3D image stack acquisition stage, the target evaluation time point was selected as the blastocyst time during the expansion phase on day 5 (i.e. (Nearby), extract the multi-focal plane image stack at that moment from the time difference system for subsequent 3D reconstruction processing.

[0031] Step S2: Semantically segment and depth estimate the image sequence sequentially, then reconstruct the surface to obtain a three-dimensional surface model. Perform instance segmentation on the trophoblast cells on the three-dimensional surface model to obtain the coordinates of each three-dimensional centroid, mark the boundary of the inner cell mass region, and extract basic morphological features.

[0032] The purpose of this step is to reconstruct a three-dimensional surface model of the blastocyst from multifocal plane images, perform single-cell-level segmentation and localization of trophoblast cells on the model, and calculate multi-dimensional basic morphological features.

[0033] The core process of 3D reconstruction includes three stages: semantic segmentation, depth estimation, and surface reconstruction. Figure 2 As shown.

[0034] First, semantic segmentation is performed on each layer of the multi-focal-plane image stack. In this embodiment, a semantic segmentation network based on the U-Net architecture is used, with an input size of 512×512 pixels. Pixels in each frame are categorized into three semantic classes: blastocyst cavity, inner cell mass (ICM), and trophoblast (TE). The segmentation network is pre-trained on a manually labeled blastocyst image dataset, and the segmentation accuracy for the three regions is evaluated using the Dice coefficient. Contour extraction is then performed on the segmentation results to obtain two-dimensional contour lines of the blastocyst outer contour, ICM boundary, and TE cell region on each focal plane.

[0035] Subsequently, depth estimation based on focus metric is performed on the segmented image. This is done for each pixel location in the multi-focal plane image stack. In each focal plane The local focus metric is calculated. The focus metric function employs the Sum-Modified-Laplacian (SML) operator, which calculates the sum of the absolute values ​​of the second-order differences in the horizontal and vertical directions within a local window centered on the target pixel (window size ranging from 5×5 to 11×11 pixels) to quantify the sharpness of that pixel on the corresponding focal plane. The SML operator is highly sensitive to defocus blur and reaches its maximum value at the optimal focus focal plane. A focus metric curve is generated along the focal plane direction for each pixel location. Gaussian or quadratic polynomial fitting is used to perform sub-pixel-level localization of the peak region of the curve to obtain the optimal focus depth value for that pixel. This depth value is combined with the known focal plane spacing to generate a depth map of the entire image. Finally, the image coordinates of the two-dimensional contour points are combined with the corresponding depth values ​​to form three-dimensional spatial points, which are then transformed using camera calibration parameters to form a three-dimensional point cloud.

[0036] Finally, a Poisson surface reconstruction algorithm was applied to the 3D point cloud to generate a closed 3D surface mesh model. The reconstructed model can clearly distinguish the spatial structures of the three regions: blastocoel, ICM, and TE, providing a geometric basis for subsequent single-cell segmentation and feature calculation.

[0037] After reconstructing the 3D surface model, it is necessary to perform single-cell-level instance segmentation on the TE cells on the model surface. This embodiment of the invention employs a strategy of "2D segmentation and 3D mapping" to achieve this goal.

[0038] The specific method is as follows: On the optimal focal plane with the highest global response in focus metric and its two to three adjacent focal planes, a deep learning-based instance segmentation network (such as Mask R-CNN or a watershed segmentation variant based on U-Net) is used to perform two-dimensional single-cell segmentation of TE cells, obtaining the two-dimensional mask and two-dimensional centroid coordinates of each TE cell. The segmentation result on the optimal focal plane is used as the main result, while the segmentation results on adjacent focal planes are used to supplement edge region cells that cannot be identified on the optimal focal plane due to defocus blur. For the same cell detected on multiple focal planes, deduplication and merging are performed according to the criterion that the Intersection over Union (IoU) is greater than or equal to 0.5, retaining the segmentation result on the focal plane with the highest focus metric value. Then, using the depth map obtained in the depth estimation step, the two-dimensional centroid of each TE cell is... Mapped to three-dimensional centroid coordinates Simultaneously, the pixels within the 2D mask are projected onto the 3D surface using the same depth mapping, thereby calculating the projected area and perimeter of the cell on the blastocyst surface. Furthermore, the boundaries of the ICM regions marked in the semantic segmentation step are also recorded for masking processing in step S3.

[0039] From the reconstructed 3D model and single-cell segmentation results, 58 basic morphological features were automatically calculated and summarized by category as follows: Size and volume characteristics: total blastocyst volume, cavity volume, ICM volume, total TE volume, blastocyst radius.

[0040] Cellular characteristics: number of TE cells and TE density, defined as the number of TE cells divided by the surface area of ​​the blastocyst.

[0041] Morphological parameters: mean TE cell area, standard deviation of area, median area, mean perimeter, standard deviation of perimeter, median perimeter, mean major axis, median major axis, mean minor axis, median minor axis, standard deviation of minor-to-major axis ratio, blastocyst sphericity, ICM shape factor, and TE surface area.

[0042] Spatial relationship characteristics: ICM eccentricity (the degree of deviation of ICM relative to the centroid of blastocyst), relative positional relationship between ICM and TE, and regularity indicators of TE cell spatial arrangement.

[0043] To verify the association between the aforementioned three-dimensional morphological features and embryonic chromosome ploidy, a statistical analysis was performed on the basic morphological features of 325 samples. The Mann-Whitney U test was used to assess the statistical significance of the differences in each continuous feature between the euploid and aneuploid groups, and the results are shown in Table 1.

[0044] Table 1: Statistical analysis results of three-dimensional morphological characteristics of blastocysts The above analysis revealed significant spatial morphological differences between aneuploid and euploid embryos. For example... Figure 3 As shown in the box plot, taking the median short axis of TE cells as an example, this parameter showed a significant difference between the euploid and aneuploid groups (p=7.3×10). -5 In embryo samples with abnormal PGT-A results, the overall spatial dimensions of the blastocyst were significantly reduced (smaller radius and volume), and the number and density of TE cells were significantly lower than in normal embryos. Simultaneously, the individual size of TE cells in aneuploid embryos was increased (larger mean area and perimeter), and they exhibited higher cellular morphological heterogeneity—the standard deviations of the perimeter, area, and ratio of the short and long axes of TE cells were all significantly elevated, indicating that the differences in TE cell size and shape are more pronounced in aneuploid embryos during development.

[0045] Step S3: Project the coordinates of the three-dimensional centroid of the three-dimensional surface model radially to the unit sphere, using the centroid of the three-dimensional surface model as the origin. After excluding the inner cell cluster region as a mask, construct the spherical density function in the effective trophoblast region using kernel density estimation. Normalize the spherical density function by Monte Carlo zero model normalization after spherical harmonic decomposition to obtain the normalized power spectrum of each degree.

[0046] The purpose of this step is to extract multi-scale information from the spatial distribution of TE cells on the three-dimensional surface of the blastocyst that cannot be obtained by global statistics.

[0047] Using the centroid O of the blastocyst's three-dimensional surface model as the origin, the least squares method is used to fit the distances from all vertices of the model surface to the origin, obtaining the optimal fitted spherical radius R. The blastocyst in the expansion stage is approximately a spherical structure, and the fitting residuals... (Defined as the standard deviation of the difference between the distances to each vertex and the fitted radius divided by the fitted radius) is typically in the range of 0.03 to 0.10, resulting in high accuracy for spherical approximation. When At that time, it was considered that the spherical surface was approximately valid.

[0048] The three-dimensional centroid of each TE cell obtained in step S2 Mapped to a unit sphere by radial projection, such as Figure 4 As shown, convert to spherical coordinates ,in Polar angle, It is the azimuth angle. This is the distance from the centroid to the origin.

[0049] Figure 8 The left side shows the spatial distribution of TE cell centroids (red scattered dots) on the surface of the blastocyst in a 3D surface model. The center of the sphere is marked as centroid O. The ICM region (blue arc region) at the top of the sphere shows a natural lack of TE cells. Figure 4 The right side shows the unit sphere as unfolded by the Mollweide equal-area projection, with the vertical axis corresponding to the polar angle. The horizontal axis corresponds to the azimuth angle. Each TE cell centroid in the left-hand 3D space has a unique corresponding projection point on the right-hand spherical unfolded diagram. The light blue area at the top is the spherical projection of the ICM mask region. The middle arrow indicates the radial projection transformation process.

[0050] The effect of radial projection on point position is limited to angular deviation, when Maximum angular deviation The angular standard deviation (approximately 19° to 22°) of the kernel function in subsequent kernel density estimation is relatively small, and its impact on the density function is negligible.

[0051] The area covered by the intracellular matrix (ICM) on the blastocyst surface naturally lacks tetracellular cells (TE cells). This biological characteristic is not due to an abnormal distribution of TE cells, but rather an inherent structural feature of the blastocyst. Without intervention, the area where the ICM is located will be misidentified as a "vacancy" for TE distribution in spherical harmonic analysis, interfering with the identification of the TE cell distribution pattern itself.

[0052] Therefore, this invention excludes the ICM region before constructing the spherical density function. Based on the semantic segmentation result of step S2, points on the boundary of the ICM region are mapped to the unit sphere through the same radial projection, forming a mask boundary, and the sphere is divided into effective TE regions. and ICM mask area ,like Figure 5 As shown, Figure 5 The spatial relationship of the three regions is shown using a top-down cross-section of a unit sphere: the light blue area represents the effective nutrient layer region. TE cells (red scattered dots) are only distributed within this area; the gray fan-shaped area is the ICM mask area. The density value in this region is set to zero; the narrow orange band between them is a linear gradient transition band, and its width is... From concentration parameter Decide. Figure 5 The density decay profile embedded in the lower right corner illustrates how the density function changes across the mask boundary: the complete estimate from the effective nutrient layer region side linearly decays through the transition zone to zero on the mask side, forming a trapezoidal transition profile. Since the ICM typically occupies only 10% to 15% of the total spherical area, the mask area is a limited proportion.

[0053] To further eliminate the Gibbs effect caused by the abrupt transition of the density function from its effective value to zero at the ICM boundary—that is, the artificial elevation of high-frequency components in the spherical harmonic decomposition due to discontinuous abrupt transitions—this invention sets a linearly gradient transition band on the outer side of the mask boundary. The transition band width... The density function value (approximately half the kernel function bandwidth) linearly decays from the full estimate on the effective TE region side to zero on the mask side within the transition band. This treatment effectively suppresses high-frequency spectral leakage caused by boundary jumps, and with the total mask area increasing by approximately 2% to 3% after the transition band, the impact on the overall power distribution of the spherical harmonic decomposition is limited.

[0054] In this embodiment of the invention, in the effective TE region Within this model, the spherical density function of the spatial distribution of TE cells was constructed using the von Mises-Fisher (vMF) nuclear density estimation method. The vMF kernel function is defined as follows: ; in Let be the unit direction vector of the query point on the sphere. It is the unit direction vector of the nuclear center (i.e., the projection direction of the centroid of the TE cell). For concentration parameters, It is the inner product of two direction vectors (equivalent to the cosine of the spherical angular distance). It is a hyperbolic sine function. When When the value is large, the kernel function has high concentration and low smoothing; conversely, the kernel function is diffuse and the smoothing effect is strong.

[0055] The TE cell spatial density function is defined as the normalized superposition of all TE cell kernel functions: ; in Total number of TE cells, normalization factor Ensure that the magnitude of the density function is not directly affected by the total number of cells.

[0056] Concentration parameter The value of directly affects the smoothness of density estimation. This invention adopts a determination method that is adaptive to the number of cells: ,in This is an adjustment coefficient, ranging from 8 to 15, with a default value of 10. The theoretical basis of this method is: assuming... There are points evenly distributed on a unit sphere, with each point occupying an average area of... The corresponding angular radius is approximately Radius; the ideal nuclear bandwidth should cover 2 to 3 average cell regions, with an angular standard deviation of approximately Since the angular standard deviation of the vMF distribution is approximately The combined equations yield the following results: Approximately Therefore, in this embodiment, we take... For further refinement, optimization can be performed on the training set using leave-one-out cross-validation. The value of .

[0057] The spherical density function requires discretization and meshing of the sphere during calculation. This invention employs the HEALPix (equal area hierarchical pixelation) scheme, setting pixelation parameters. Corresponding to the total number of pixels Each pixel has an equal area, with an angular resolution of approximately 1.8°. The HEALPix scheme ensures that each pixel has a strictly identical area, which is crucial for the numerical accuracy of subsequent spherical harmonic decomposition.

[0058] For the spherical density function constructed above Perform spherical harmonic expansion: ; in For spherical harmonic functions ( In degrees, (for order) is the spherical harmonic coefficient (complex number). The maximum degree. Spherical harmonic decomposition decomposes the density distribution on the spherical surface into a series of orthogonal components according to spatial frequency: low degree components ( Smaller (smaller) describes large-scale distribution patterns, height number component ( Larger (larger) describes local changes at a small scale.

[0059] The range is 4 to 10, with a default value of 6. This value is determined based on: for The spherical distribution at discrete points can reliably recover the maximum degree of spherical harmonics, which is approximately... The typical range of TE cell counts in the samples of this invention is 68 to 92, corresponding to a maximum degree of reliable recovery of approximately 7 to 9. It can capture meaningful spatial distribution patterns across the entire range without being affected by sampling noise.

[0060] Extracting the power spectrum at each degree from the spherical harmonic coefficients: ; For rotational invariants—regardless of the orientation of the blastocyst in the incubator, the same blastocyst's... The values ​​remain constant. The biological meaning of the power spectrum at each degree is as follows: Corresponding to the global average density component; Corresponding to hemispherical asymmetry, it reflects the density difference between polar TE and wall TE; Corresponding to the four-pole deformation mode; The number 6 corresponds to mid-to-high frequency spatial variations, which can capture local clustering or empty zone patterns.

[0061] The number of TE cells varies among different blastocysts. When the number of cells is low, the statistical noise of the discrete point distribution on the sphere is relatively large, which leads to a systematic increase in the high-frequency power spectrum. To eliminate this systematic bias of the power spectrum caused by the number of cells, this embodiment of the invention uses Monte Carlo simulation to construct a null model for normalization.

[0062] The specific method is as follows: for a given and Values ​​are generated, with at least 100 sets (preferably 200 sets) randomly and evenly distributed within the effective TE region. On A set of simulated points. For each set of simulated points, the same processing procedure as for the observed data is performed: including vMF kernel density estimation (including ICM masking and transition band processing) and spherical harmonic decomposition; calculating the power spectrum of each set at each degree, and taking the average of all simulated sets as the expected power spectrum at that degree. Record the standard deviation at the same time. The normalized power spectrum is defined as: ; The meaning is: if This indicates that the observed data is in degrees. If the TE distribution variation intensity at the corresponding spatial scale exceeds the statistical expectation generated by the same number of random uniform points, it suggests the existence of a real spatial distribution anomaly (such as local clustering or empty areas) at that scale; if This indicates that the variation at that scale is consistent with random noise, and it is impossible to determine the existence of a meaningful spatial structure. For example... Figure 6 As shown, euploid and aneuploid embryos exhibit identifiable differences in spherical density function and normalized power spectrum. Figure 6 The top row shows the analysis results of euploid embryos: the density function color distribution in the Mollweide projection thermogram on the left is relatively uniform, indicating that TE cells do not have obvious clusters or empty areas on the sphere; the various degrees in the bar chart on the right... The values ​​are all close to those indicated by the red dashed line. baseline, corresponding . Figure 6 The bottom row shows the analysis results of aneuploid embryos: the left-hand heatmap shows obvious color inhomogeneity, with localized high-density clusters and low-density sparse areas; the right-hand bar chart... and Significantly higher than The baseline indicates that the variation in the TE distribution at the spatial scales corresponding to quadrupole deformation and local clustering far exceeds the random level, while... A value below the baseline indicates insufficient global uniformity, corresponding to The comparison between the top and bottom rows intuitively reflects the ability of the spherical harmonic power spectrum to distinguish different spatial distribution modes of TE.

[0063] The null model can also be used to assess the statistical significance of power spectrum deviations at various degrees: calculating standardized scores. ,when This indicates that the power deviation from zero at this degree level reaches a 95% confidence level.

[0064] Special notes are required Meaning: Due to the attenuation or zeroing of the density in the ICM mask region, the observed data... The component (global average density) is lower than the uniform distribution on a perfect sphere. However, the null model also... Random points are generated on the region and the same masking process is applied, therefore ICM mask pairs are already included. The systemic impact, In practice, it reflects the degree of deviation of the uniformity of TE within the effective region from the random baseline.

[0065] In the embodiments of the present invention, the normalized power spectrum Based on this, the Spatial Organization Index (SOI) is defined: ; ; SOI molecules To measure global uniformity, the denominator is the degree. The sum of the normalized power spectra measures the total spatial variation at each scale. A higher SOI value indicates a more uniform spatial distribution of TE cells and less anomalous spatial variation at each scale. Since both the numerator and denominator of the SOI are normalized under the same ICM mask conditions using a null model, differences in ICM area have little impact on the overall spatial variation. The systematic effects have been eliminated, and SOI is comparable across samples among embryos with different ICM areas. Logarithmic transformation was employed. It can compress the dynamic range of numerical values, improving the distribution characteristics of features in subsequent machine learning models.

[0066] The final spherical harmonic feature set is defined as ,when The set contains a total of 8 features. These 8 spherical harmonic features, together with the 58 basic morphological features extracted in step S2 and the 14 developmental dynamic parameters extracted in step S1, constitute a total of 80 extended feature sets, which are then used in the screening process in step S4.

[0067] Step S4: Combine basic morphological features, developmental dynamics parameters, and normalized power spectra into an extended feature set and then filter it.

[0068] The extended feature set contains 80 features, which inevitably include redundant and noisy features. To extract the most predictive subset of features, avoid model overfitting, and improve interpretability, this embodiment of the invention employs a two-stage screening strategy.

[0069] The first stage is statistical screening. For each feature in the expanded feature set, it is grouped and compared according to euploid and aneuploid groups. The Mann-Whitney U test is used for continuous variables, and the chi-square test is used for categorical variables. Features with a significance level p-value less than a preset significance threshold are retained. The purpose of this step is to eliminate noisy features that are clearly irrelevant to the ploidy outcome. In the 325 samples of this invention, this step retained approximately 19 statistically relevant features from 80 features.

[0070] The statistical analysis results are shown in Table 2. All three characteristic categories contained indicators significantly associated with PGT-A outcomes. Among the kinetic parameters, the timing of late developmental events was particularly closely related to ploidy status. , , , , , Among the three-dimensional morphological features, the short-axis parameters of TE cells showed the most significant differences: , Secondly, the average area of ​​TE. and cavity volume In spherical harmonic characteristics, and The significant correlation with PGT-A outcome indicates a detectable difference in the mid-to-high frequency spatial distribution pattern of TE cells between euploid and aneuploid embryos.

[0071] Table 2: Statistical Analysis Results of Dynamic Characteristics The second stage is model-based importance ranking. A LightGBM gradient boosting decision tree model is used to train the features retained from the first stage, outputting an importance score for each feature to the prediction target. Cumulative importance curves are then plotted after ranking by importance from highest to lowest. Figure 7 As shown on the right, the curve exhibits a distinct "elbow" inflection point around the 15th feature. Further increasing the number of features yields very limited predictive power gains, while significantly increasing model complexity and the risk of overfitting. Therefore, the top 15 features are selected to form the final core feature subset, such as... Figure 7 The left side shows the ranking of the top 15 features by importance after filtering.

[0072] To further verify the contribution direction and pattern of each feature to the prediction model, this invention uses the SHAP (SHapley Additive exPlanations) method for attribution analysis. For example... Figure 8 As shown, the SHAP analysis results indicate that the top-ranked features include both kinetic parameters (such as...) It also includes three-dimensional morphological parameters (such as...) , ) and spherical harmonic space characteristics (such as This verifies the effectiveness of the multimodal feature fusion strategy—the model does not rely solely on information from a single category, but rather integrates features from both temporal and spatial dimensions. For example, SHAP values ​​exhibit a clear color gradient distribution: larger values... Values ​​(slower development) correspond to a negative SHAP contribution (leaning towards "aneuploid" prediction), with smaller values... The value corresponds to a positive contribution, which is consistent with the clinical understanding that embryos with developmental delays have a higher risk of chromosomal abnormalities. For example, a larger blastocoel volume makes a positive contribution to the prediction of "euploidy," providing an intuitive biological explanation.

[0073] Step S5: Train an ensemble learning model using the selected feature set, and use the ensemble learning model to output the euploidy prediction result of the blastocyst to be tested.

[0074] The core feature subset (approximately 15 fusion features) selected in step S4 is used as input, and the known PGT-A detection results (euploidy / aneuploidy) are used as training labels to construct an embryo euploidy prediction model. This embodiment of the invention employs the LightGBM gradient boosting decision tree algorithm as the prediction model. This algorithm, based on a histogram-based discretization strategy, can efficiently handle continuous numerical features, exhibiting excellent performance in both training speed and prediction accuracy, and is suitable for processing the primarily numerical tabular feature data in this invention.

[0075] Model training employs stratified K-fold cross-validation (K ranging from 5 to 10) for performance evaluation and hyperparameter optimization. Stratified sampling ensures that the ratio of euploid to aneuploid samples in each fold is consistent with the overall dataset. Hyperparameters are tuned using grid search, with the search space including key parameters such as learning rate, maximum tree depth, minimum number of samples per leaf node, and regularization strength, which will not be elaborated upon in this embodiment.

[0076] To determine the optimal prediction algorithm, this invention systematically trained and compared several mainstream ensemble learning algorithms, all using the same Top15 fusion feature set and cross-validation scheme. Evaluation metrics included AUC (Area Under the Receiver Operating Characteristic), accuracy, precision, recall, and F1 score. The evaluation results are shown in Table 3.

[0077] Table 3: Evaluation Results of Different Ensemble Learning Models The above results demonstrate that LightGBM achieves optimal performance across all evaluation metrics. For example... Figure 9 As shown in the subplots, LightGBM's ROC curve is closest to the top left corner, with an AUC of 0.7095, significantly outperforming random guessing (AUC=0.5) and other comparative models. Figure 10 As shown in the comparative chart, the LightGBM curve consistently ranks above the other models on the same coordinate system. LightGBM also achieves the highest accuracy, precision, recall, and F1 score, with similar and relatively balanced values, indicating no significant predictive bias in identifying euploid and aneuploid embryos. The VotingEnsemble model, with an AUC of 0.6866, ranks second, but is still lower than the standalone LightGBM model, suggesting that a finely tuned single strong learner may outperform a simple model stacking strategy in this task.

[0078] To quantify the independent contribution of spherical harmonic spatial distribution features to prediction performance, the following ablation experiments were designed in this embodiment of the invention to compare the prediction effects of different feature combinations. All experiments used the same LightGBM algorithm and hierarchical 5-fold cross-validation. The ablation experiment results are shown in Table 4.

[0079] Table 4: Ablation Experiment Results Ablation results showed that the AUC using the three-dimensional basic morphological features alone (C1) was 0.643, and the AUC using the kinetic parameters alone (C2) was 0.628. Combining both (C3) improved the AUC to 0.689, indicating complementarity between morphological and temporal progression information. Further adding eight spherical harmonic spatial distribution features to C3 (C4) resulted in an AUC of 0.710, a 3.0% improvement over C3. The DeLong test was used to assess the AUC difference between C3 and C4, with p = 0.031 < 0.05, indicating a statistically significant difference. This confirms that the spherical harmonic features provide independent predictive information beyond the basic morphological features and kinetic parameters. The 95% confidence interval was calculated using 2000 bootstrap resampling iterations.

[0080] The technical features of the above embodiments can be combined arbitrarily. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described; only preferred embodiments of the present invention are illustrated. The descriptions are relatively specific and detailed, but they should not be construed as limiting the scope of the present invention. As long as the combination of these technical features does not contradict each other, it should be considered within the scope of this specification.

[0081] It should be noted that those skilled in the art can make various modifications and improvements without departing from the inventive concept, and these all fall within the scope of protection of this invention. Therefore, the scope of protection of this invention should be determined by the appended claims.

Claims

1. A method for predicting euploidy based on the fusion of blastocyst developmental dynamics and spherical harmonic decomposition, characterized in that, Includes the following steps: Step S1: Obtain image sequences by taking multi-focal plane images of the blastocyst using a time-difference imaging system, and extract developmental kinetic parameters from the image sequences; Step S2: Semantically segment and depth estimate the image sequence sequentially, then reconstruct the surface to obtain a three-dimensional surface model. Perform instance segmentation on the trophoblast cells on the three-dimensional surface model to obtain the coordinates of each three-dimensional centroid, mark the boundary of the inner cell mass region, and extract basic morphological features. Step S3: Project the coordinates of the three-dimensional centroid of the three-dimensional surface model radially to the unit sphere with the centroid of the three-dimensional surface model as the origin. After excluding the inner cell mass region as a mask, construct the spherical density function in the effective trophoblast region using kernel density estimation. Normalize the spherical density function by Monte Carlo zero model normalization after spherical harmonic decomposition of the spherical density function to obtain the normalized power spectrum of each degree. Step S4: Combine the basic morphological features, the developmental kinetic parameters, and the normalized power spectrum into an extended feature set and then filter them; Step S5: Train an ensemble learning model using the selected feature set, and use the ensemble learning model to output holoplasticity prediction results for the blastocyst under test.

2. The method according to claim 1, characterized in that, The developmental kinetic parameters mentioned in step S1 include the time of pronucleus emergence. Pronuclear disappearance time The start time of each stage from the two-cell stage to the nine-cell stage to 1. Start time of mulberry embryonic stage Time of blastocyst cavity appearance The cavity expands to half the size of the blastocyst in time. and the time for complete cavity expansion .

3. The method according to claim 1 or 2, characterized in that, The basic morphological features mentioned in step S2 include size and volume features, cytological features, morphological parameters, and spatial relationship features; The size and volume characteristics include the total volume of the blastocyst, the cavity volume, the inner cell mass volume, and the blastocyst radius; the cytological characteristics include the number and density of trophoblast cells; the morphological parameters include the mean area, standard deviation, mean perimeter, mean major axis, mean minor axis, standard deviation of the ratio of major to minor axis, and blastocyst sphericity; and the spatial relationship characteristics include the inner cell mass eccentricity.

4. The method according to claim 1, characterized in that, In step S2, the depth estimation uses an improved Laplacian operator as the focal metric function. The focal metric curve formed by each pixel along the focal plane direction of the image sequence is fitted to obtain the depth value by sub-pixel-level peak localization. The instance segmentation adopts a two-dimensional segmentation and three-dimensional mapping strategy. The trophoblast cells are segmented into two-dimensional instances on the best focal plane with the highest focal metric response and its adjacent focal planes. After deduplication and merging according to the intersection-union ratio greater than or equal to 0.5, the segments are then mapped to the three-dimensional centroid coordinates through the depth map.

5. The method according to claim 1, characterized in that, The kernel function for kernel density estimation in step S3 The expression is: ； in Let be the unit direction vector of the query point. Let be the unit direction vector of the centroid of the trophoblast cell. This is a concentration parameter; the concentration parameter pass Confirmed, among which The number of trophoblast cells, It is an adjustment coefficient with a value range of 8 to 15.

6. The method according to claim 5, characterized in that, Step S3, which excludes the inner cell mass region as a mask, also includes: A width of [width value] is set outside the mask boundary. A linearly gradient transition zone of curvature, within which the value of the spherical density function linearly decays to zero from the effective nourishment layer region side toward the mask side.

7. The method according to claim 1, characterized in that, The Monte Carlo null model normalization method in step S3 includes: performing spherical harmonic decomposition on the spherical density function to obtain the original power spectrum at each degree, wherein the degree... In spherical harmonic decomposition, a non-negative integer representing the spatial scale, with values ​​ranging from 0 to the maximum degree. ; Generate at least 100 groups randomly and evenly distributed on the effective nutrient layer area. A set of simulated points, of which The number of the trophoblast cells; Perform the same kernel density estimation and spherical harmonic decomposition as the observed data on each set of simulated points, and calculate the degree of all simulated sets. The average power spectrum on the desired power spectrum ; Observation data in degrees The original power spectrum Divide by the desired power spectrum Obtain the normalized power spectrum : ； in A value greater than 1 indicates a degree The distribution variation of the trophoblast cells at the corresponding spatial scale exceeds the statistical expectation of a random uniform distribution.

8. The method according to claim 7, characterized in that, In step S3, the normalized power spectrum is obtained. This is followed by the calculation of spatial organization index and its logarithmic transformation as supplementary features; The Spatial Organization Index (SOI) is defined as a degree. Normalized power spectral components when the value is 0 with degrees From 1 to The ratio of the sum of the normalized power spectral components at time: ； The logarithmic transformation is Its calculation method is as follows , where ln is the natural logarithm; The higher the SOI value, the more uniform the spatial distribution of the trophoblast cells.

9. The method according to claim 1, characterized in that, The maximum degree of the spherical harmonic decomposition described in step S3 The value range is 4 to 10. The spherical meshing adopts the HEALPix scheme of equal area hierarchical pixelation and the pixelation parameters are... Greater than or equal to 32.

10. The method according to claim 1, characterized in that, The screening in step S4 adopts a two-stage method: the first stage uses statistical tests to retain features with a significance level less than a preset significance threshold; the second stage uses a gradient boosting model to rank the features retained in the first stage by importance, and determines the optimal number of features based on the inflection point of the cumulative importance curve. The ensemble learning model described in step S5 uses the LightGBM gradient boosting decision tree algorithm.

Images