Accurate Machine Learning Model for Human Embryo Morphokinetic Stage Detection

Automatic Morphokinetic stage detection

Machine learning methodologies have been employed in different stages of embryo selection processes. Most of these models focus on directly predicting the success rate of an embryo to reach the foetal heartbeat stage [14]. Some tools reported surpassing embryologists in the accuracy of identifying viable embryos[15]. There are three main approaches to utilise TLI for embryo selection, models that work directly with the video generated from TLI images and no frame selection is needed by user [15] [16], methods that only work on single image and user need to select specific frame of the TLI images and use that image for prediction [17][18], and there is a hybrid approach where timings are extracted manually and then fed into a ML model such as Kidscore [19] [9].

While approaches such as IDAScore V1 [15] and V2 [16] utilize embryo time-lapse videos, others including ERICA [20], The Life whisperer [17], and Stork [18] rely on a static single image of the blastocyst. The pipeline for selecting the blastocyst image is not entirely automated for ERICA [20], LIFE whisperer [17], and Stork [18], necessitating manual selection. However, the automation of morphokinetic event detection can streamline the entire process by identifying the most appropriate images and incorporating them into the pipeline as a hybrid approach, where timings are extracted and then fed into a ML model such as Kidscore [19] [9].

Early methods of automatic morphokinetic stage detection relied on manually designed features to identify morphokinetic events [21]. For example, Feyeux et al. used grey-level analysis of microscopic images to predict morphokinetic stages [10]. In recent years, however, modern approaches have primarily utilized deep learning techniques [14], focusing on the supervised training of convolutional neural networks to facilitate the automatic annotation of morphokinetic events. Models employing ResNet-50, Long Short-Term Memory (LSTM), and ResNet-3D type models are now commonly adopted [3]. While the models that predict Morphokinetic events show similar levels of accuracy and overall performance, they are developed on proprietary datasets that are not available. Two prominent examples are the commercial models for EmbryoScope [21] and Gerry Incubators [22]. Zabari et al. also proposed a DL methodology that leverages video frame-based initial predictions, which are further refined through Monotonic regression. This approach aims to mitigate prediction noise by limiting predictions to either the current or an immediately subsequent developmental class [21]. Recently, Gomez et al. have introduced a large, annotated dataset, comprising 704 videos of developing embryos, featuring 337,000 images across 16 developmental stages. This dataset provides a critical resource for model benchmarking and development. In this study two different DL models were developed and trained on Gomez et.al dataset [3]. a transformer [23] architecture will be used to fuse the input image features and information about time post-fertilisation when the image was captured. This novel, multimodal methodology incorporates both embryo images and the elapsed time since fertilisation as inputs to enhance the model’s performance.

Dataset

Time-lapse imaging of human embryos was obtained from Gomez et al. [3]. The dataset contains 704 Embryoscope videos recorded at 7 focal lengths and annotated for 16 morphokinetic events (Vitrolife ©). This dataset is labelled by one embryologist. To label the videos, the embryologist first identified the frame in which each event occurs and assigned the label to these frames. Then, all subsequent frames until the next morphokinetic event occurs are assigned the current label [3].

Each video has on average, 8 or more events. 499 of the videos are of viable embryos, with the remaining 205 non-viable, to attempt to capture the myriad features that can cause embryonic development failure, such as high levels of fragmentation, 3PN, necrosis, etc. Only the central focal plane was used and only files that had uncorrupted jpeg images were analysed, resulting in 273,438 images with labelled events.

In order to mitigate some of the subjectivity of morphokinetic transitions effect, two images forward and backward from the recorded transition time were removed from the training set. All images were evaluated in the test set and included in performance metrics.

During the dataset quality review, it was noted that, as is standard practice, the embryos are removed for freezing or transferring on Day 5; however, the label still reflects the last event, i.e. expanded blastocyst. This mislabelling of the wells introduces ground truth errors in the dataset. To identify and relabel these images, an empty/non-empty model was developed and applied to all images(Appendix 1). The results were verified visually, and 9,734 images were labelled as empty. Additionally, it was noted that there were only 41 hatched blastocyst images (Table 1). This small number of examples makes predicting this class difficult and subject to increased variability.

Table 1 presents the data annotations along with the number of samples corresponding to each class within the new fixed dataset.

Deep learning models architecture

For model 1, the backbone is an efficient net V2 large [24] fine-tuned to categorize 17 morphokinetic classes (Table 1). The input to the model was static greyscale JPG images of 380×380 resolution and the greyscale values were copied 3 times into the R,G,B channels. Weights were initialised with ImageNet dataset [24]. Figure 1 (a) illustrates the structure of this network.

• Download figure
• Open in new tab

Figure 1

Architecture of the proposed models (a) Model 1, utilises efficient-net-V2-large as backbone for feature extraction and it has a fully connected layer as classifier head. (b) Model 2, has efficient net-large-V2 as feature extractor, this model has two classifier heads one after the feature extraction step and one after the transformer. The transformer in model 2 fuses the time after fertilization and extracted features from the efficient-net. The two head structure of this network ensures the proper gradient flow to backbone.

• Download figure
• Open in new tab

Figure 2

Example of the effect of postprocessing algorithm on predictions of Model 2. Model predictions are uncertain in some moments of the video, specifically near the moment that a morphokinetic stage change is recorded. The postprocessing algorithm ensures a clean set of predictions by assessing and ignoring noisy predictions

Model 2 also utilises an Efficient-net V2 L as the backbone for processing images, followed by a transformer [23]. The transformer has only one encoder layer with hidden size of 512 with 4 self-attention heads. The second input to the network is the number of minutes that have elapsed since the time of fertilization. As transformers are designed to work on sequences of vectors, the time since fertilisation is converted into a binary vector representing two-hour windows of time. This approach allows the self-attention mechanism to interact with the image features extracted earlier in the pipeline. Specifically, we assume a maximum of seven days incubation time and encode every two hour window with a one-hot vector or length resulting in a vector of length 84. Where all v elements are 0 except of v_index. For example, if the time is 150 minutes, the index will be calculated as the time vector will be [0, 1, 0, 0, …, 0], the time vector will be the same for 160 minute, but for 240 minutes the index will be and the time vector will be [0, 0, 1, 0, …, 0].

The time vector and the image features generated by the backbone are passed to the transformer layer (Figure 1b) layer before a final Multilayer Perceptron (MLP) layer to make class predictions for each image. By design, this model has two classification heads, one immediately after the backbone and one after the transformer layer. During the training, two cross-entropy losses are summed. This ensures the proper gradient for backbone layers and prevents the model from overfitting only on the time vector inputs.

Model 3 has the same architecture as Model 1, but was trained on the original Gomez dataset that did not include the reclassified ‘empty’ well labels. Model 4 was trained with Resnet-50 backbone to recreate the previous work by Gomez et al [3] to serve as a baseline model for comparison.

All models were trained using the Adam optimizer, as provided in the PyTorch 2.0.0 library [25], with an initial learning rate of 0.001 with Cross-entropy loss function. The learning rate was dynamically reduced upon observing a plateau decrease in the validation set performance. The networks underwent training for 50 epochs on the training subset of the dataset and were subsequently evaluated on a test subset, which remained unseen during the training phase.

The issue of class imbalance represents a significant challenge within this dataset, as illustrated in Table 1. Such imbalance can induce bias in the neural network towards classes with a higher proportion of samples. To mitigate this issue inverse class frequency method was used. This technique enables the model to generate regulated gradients for the classes that are represented by fewer samples.

Embryo images often exhibit significant variations in brightness, with some regions appearing overly bright and others notably dark, while the embryo is typically positioned centrally therefore proper image normalization is crucial. The image augmentation techniques employed during the training included rotating and shifting the image with a 30% probability, flipping the image with a 50% probability, and applying noise or blur with a 50% probability, followed by a final normalization step to ensure consistent image quality across the dataset. The augmentation step was conducted using the Albumentations library in Python [26]. The normalization step was performed using image contrast enhancement (CLAHE) and the normalize function, with average pixel values across RGB channels set as (0.485, 0.456, 0.406) and the standard deviation as (0.229, 0.224, 0.225) to align with the ImageNet weights.

Postprocessing algorithm

The morphokinetic annotation of videos critically depends on the consistency of predictions on all video frames. Given the dynamic nature of embryos, there are instances when the classification of a frame is ambiguous, resulting in “noisy” predictions. It is essential to accurately identify the precise moment when the morphokinetic stage of the embryo changes because these exact time values are utilized in subsequent methods that assess the embryo’s viability, such as KidScore [12]. A heuristic method was developed to address this issue. This method is applied to the predictions made by networks for each frame of the embryo’s timelapse video. The algorithm accommodates shifts to more advanced morphokinetic classes as well as instances of reverse cleavage. It identifies a trend within each time interval and corrects mispredictions that disrupt this trend. Mispredictions often occur due to embryo movement and the movement of cells within it, where some cells may temporarily be obscured in the image, rendering any detection of morphokinetic state changes during these moments inaccurate. Consequently, the algorithm initially amends the prediction in our models by substituting low-confidence predictions with the last morphokinetic state prediction with over 80% confidence. Subsequently, it detects changes by comparing the predicted class with the prediction of the preceding image. Let P_I represent the value of the most probable predicted class at the index i, changes are defined as Δ : By summing the values of changes within a sequence, various consecutive groups can be delineated as follows: The values of G_I can be utilized for grouping purposes, indicating that all images possessing identical G_I values belong to a contiguous group and share the same prediction values without interruption. Through the comparison of the length of each continuous group and its neighbouring groups, a determination is made regarding whether the current group should be regarded as noise or retained as a correct morphokinetic stage change. The algorithm is outlined in the following:

• Let C = {c1, c2, …, cn} represent the set of unique predicted classes in the video.

• For each class c in C, create a subset Dc of all the predictions D such that all elements in Dc have prediction class equal to c. or each subset Dc, identify unique group values as Gc = {g1, g2, …, gm}.

• Calculate the length of each consecutive group g in Gc, denoted as Lg, where Lg is the number of elements in Dc that belong to group g.

• Determine the minimum and maximum values within G_c, denoted as g_min and g_max, respectively.

• For each potential group identifier i_g in the range , identify if i_g is an interruptive group by checking if ig ∉ Gc. For each interruptive group, calculate its length L_ig

• All interruptive groups are discarded and only the main groups are kept as final labels

Single Image Processing

During the testing step for the Models, standard classification metrics were calculated. All the images in the test dataset were processed by the models and the outputs compared with ground truth in the dataset. Confusion matrix, accuracy, F1-score, precision, and recall were calculated for each model. Figure 3 presents the confusion matrix of Model 2 applied to the test set of the dataset. It is evident that certain stages, such as t2 and t4, are easier to identify, whereas other stages, including tB and t5, are more challenging for the model to classify accurately. Misclassifications most commonly occur between consecutive stages, such as when images of class t5 are misclassified as either t4 or t6. It is highly unlikely for the model to confuse an image with a class significantly ahead or behind; for instance, an image labeled as t5 has less than a 1% chance of being classified as anything other than t4 or t6.

• Download figure
• Open in new tab

Figure 3

Model 2 confusion matrix on the test set in percentages. The blue gradient bar on the right of the matrix indicates the percentage graphically with 0% - white and 100% - dark blue. Values of 0.0 indicate that not a single miss classification occurred and perfect model prediction would be 100% down the lower right diagonal.

Standard classification metrics for the proposed models are shown in Table 2. Both exhibit high accuracy levels, with Model 2 demonstrating a slightly superior performance on all the metrics.

Table 3 presents the accuracy of the models trained by Gomez et al. alongside those proposed in this study. Model 1 and 2 demonstrate an approximate 17% improvement in accuracy compared to the models in the original dataset paper [3]. Models 1 and 2 have the best performance, and models 3 and 4, which did not have empty labels in their dataset, have less accuracy. It can also be seen that Model 4 which was trained to represent the Gomez work is performing only slightly better than resnet-50 that was reported by Gomez et al [3].

To visually assess the root causes of misclassification of the models, a random sample of images that were mis-predicted by Model 2 for the classes that had highest error rate are shown in Figure 7 to Figure 10.

Postprocessing

After applying the postprocessing method the exact predicted times for each morphokinetic event is extracted. For each detected morphokinetic stage change one single time is extracted. In the analysis here only the resulting morphokinetic state is considered, for example, changes from t2 to t4, and t3 to t4 are considered as transition to t4.

While the model results were accurate, a plus or minus one developmental stage classifications was observed in a sequential set of images (Figure 3). This is to be expected due to the embryos moving around the incubation dish i.e., not fixed, exposing or obscuring different features. To combat this, a heuristic post-processing was implemented.

Upon applying the model and the subsequent postprocessing algorithm to the entire video in the test dataset, the timings for changes in the MK stage were extracted. Error is calculated as time_prediction − time_{ground truth} the average and standard deviation and 25%, 50%, 75%, 80%, 85%, 90%, 95%, and 99% percentile of the errors belonging to each class is calculated and presented in Table 4. It is evident that errors are small for the majority of cases and large errors are seen in less than 5% of cases. The largest errors belong to tB class.

Time extraction had less average error in some classes such as tPNf with -0.01 hours average error while the error is more significant in t9+ with -2.94 average hour error. It is important to note that images are 20 minutes apart and only one frame error contributes as 0.33 error.

Due to the significant impact of labeling errors on the accuracy of stage change time extraction, we conducted an in-depth assessment of label subjectivity and identified the primary causes of model misclassification. Figures 4 and 5 provide examples of the largest errors observed, showcasing a random sample of frames from classes with the greatest time detection discrepancies. These visual examples help highlight the inconsistencies in labeling and the challenges they pose for model accuracy.