The study is based on a three-dimensional parametric model of a leaflet prosthesis for the mitral position, generated algorithmically using MATLAB. The geometry of the leaflet was determined by a combination of six parameters, as shown in Table 1 and Figure 2.

By varying the values of these parameters, within specified ranges, the MATLAB algorithm generated a set of lines, a surface, and finally an STL mesh composed of geometric primitives (triangles). All models are represented as shells, described by 1,109–10464 S3-type elements per leaflet, aligning well with the approach endorsed in modern scientific approaches (Johnson et al., 2021). A total of 11,565 leaflet designs were generated.

The models were analyzed using FEM to simulate the effect of pressure on the leaflet opening. Within the MATLAB environment a calculation file for “Abaqus/CAE” was created, specifying the following:

• Boundary conditions. The valve model was fixed to prevent displacement along the lower edge (“encastre” boundary condition) in all six degrees of freedom. This was done to mimic the method of attachment of the valve to the prosthetic valve.

The results of all 11,565 numerical experiments were then systematically analyzed using our Python algorithm. We note that this work represents a proof-of-concept approach to semi-automated optimization; therefore, at this stage, we deliberately limit our analysis to the most important parameters from our perspective: the opening area (LMN) and the maximum principal stress (STS). The former is a primary indicator of prosthetic performance and illustrates the effects of optimization. The latter serves as a key criterion for leaflet material failure under sustained loads (Claiborne et al., 2013; Chen et al., 2019), leading us to focus primarily on optimizing the mechanical stress distribution.

All valve design generation, FEM, and data collection processes were performed simultaneously on three high-performance computers, each equipped with a 4-core, 8-thread Intel Core i7-3770K processor with 3.4 GHz base clock frequency, 32 GB of RAM, and an Nvidia Quadro RTX 4000 graphics processor with 8 GB of video memory. The average wall clock time for each calculation and analysis was 766 s (12.8 min).

We use six parameters–HGT, DIA, THK, CVT, ANG, and ELM - as features to train our model with two parameters, LMN and STS, used as targets. Our choice of ML algorithms was motivated through an extensive research into applied techniques within the medical domain, including Decision Tree (Quinlan, 1986), Random Forest (Breiman, 2001), Extra Tree (Geurts et al., 2006), Neural Network (Hornik et al., 1989), eXtreme Gradient Boosting (XGBoost) (Chen and Guestrin, 2016), Light Gradient Boosting Machine (Light GBM) (Ke et al., 2017), CatBoost (Prokhorenkova et al., 2017), and Ensemble (Valentini and Masulli, 2002) with the addition of Stacking (Wolpert, 1992). These regressors are integrated into a single AutoML framework for training and validation. In addition, we used a “Baseline” regressor, which makes predictions by averaging the target from training data, and acts as a comparison point for the 8 regressors evaluated but is not used for PHV prototyping.

Normalization can significantly influence the performance of ML algorithms, prompting us to investigate how different normalization techniques influence input parameters and model output. We applied normalization separately on input (geometric parameters of a PHV design) and model output, focusing on 5 types of normalization: 1) Raw–no transformation is applied, 2) MinMax–features are transformed by scaling each feature to a given range (from 0 to 1), 3) Standard–features are standardized by removing the mean and scaling to unit variance, 4) Robust–features are scaled using statistics that are robust to outliers, removing the median and scaling according to the quantile range, 5) Power–a family of parametric, monotonic transformations that make data more Gaussian, using the Yeo-Johnson transformation (Yeo and Johnson, 2000) to optimize variance and minimize skewness.

Our study consist of 11,565 designs with 9,252 designs randomly selected for training and 2,313 designs for validation. Since the focus of this study is the prediction of numerical scores (regression), Root Mean Square Error (RMSE) was used as the primary loss function. We implemented and evaluated two single-output ML models, with the first model predicting LMN scores and the second model predicting STS scores.

We applied several high-quality techniques to optimize the performance of our ML model during training. These techniques include golden features, k-means features, feature selection, hill climbing, and boost on errors, which have demonstrated their effectiveness in various areas of ML.

Golden features are a set of features with high predictive power. They are generated by applying mathematical operators (e.g., summation) to pairs of original features. To identify the best features, we generated all possible unique pairs of original features and randomly subsampled them to 250,000 if there were more than 250,000 pairs. The predictive power of a newly created feature was estimated using the decision tree algorithm with a maximum depth of 3. Only the top 10 new features were added to the training data as golden features.

K-means features are created using K-means clustering, which adds features based on cluster separability. These features can include a sample’s distance from all cluster centroids and a sample’s cluster label.

The feature selection is done in two steps. First, a random feature with a uniform distribution between 0 and 1 is inserted into the data set. Subsequently, all irrelevant features are discarded. If an existing feature has an importance lower than the random feature added in the first step, it is dropped.

Hill climbing is a local search algorithm that continuously moves in the direction of increasing height to find the top of the mountain or the best solution to the problem (Hernando et al., 2018). It stops when it reaches a peak value where no neighbor has a higher value. Hill climbing is a variant of the generate and test method that generates feedback to decide which direction to move in the search space. Since the hill climbing algorithm is a greedy approach, its search moves in the direction that optimizes the cost (loss function). Note, however, that this algorithm does not backtrack the search space because it does not remember previous states.

Finally, boost on errors is an approach similar to boosting techniques where a future model is improved using past errors.

The models were trained on a computer system equipped with an Intel Core i9-10940X processor clocked at 3.30 GHz, an NVIDIA GeForce RTX 3090 graphics card, 256 GB of RAM, and the Windows 10 operating system.

To determine the optimal geometries of PHV, we implemented six state-of-the-art optimization algorithms for sampling geometric parameters and efficiently pruning suboptimal designs: Random Search (RS) (Bergstra and Bengio, 2012), Tree-structured Parzen Estimator (TPE) (Bergstra et al., 2011; Bergstra et al., 2013), CMA-ES-based Algorithm (CMA) (Auger and Hansen, 2005; Ros and Hansen, 2008; Hansen, 2016; Nomura et al., 2021), Nondominated Sorting Genetic Algorithm (NSGA) (Deb et al., 2002), Multiobjective Tree-structured Parzen Estimator (MOTPE) (Ozaki et al., 2022), and Quasi-Monte Carlo Algorithm (QMC) (Bergstra and Bengio, 2012). To efficiently perform the design search procedure, we implemented a pruning mechanism based on the Hyperband pruner (Li et al., 2017), which demonstrated its performance in various optimization experiments (Ozaki et al., 2020) and outperformed benchmarks with Kurobako (Imamura, 2020) for non-deep learning tasks. In our study, each optimization algorithm was limited to 2000 iterations.

To evaluate the geometric design parameters, we introduced our own objective function based on the weighted harmonic mean, allowing for direct estimation of lumen and stress scores for each PHV design. Each PHV geometry is evaluated by the design score D S N : D S N = 1 + α 2 ∗ L M N s ∗ S T S s α 2 ∗ S T S s + L M N s (1) L M N s = L M N = O R F C O R F C max (2) S T S s = 1 − S T S U T S , w h e n   S T S U T S < 1 0 , w h e n   S T S U T S ≥ 1 (3) where L M N and S T S represent the absolute lumen and stress values obtained for a given set of geometric parameters in a design. α is a positive factor that allows for a trade-off between L M N s and S T S s , with a value chosen such that L M N s is considered α times as important as S T S s . Five different values of α were estimated, including 0.2, 0.5, 1.0, 2.0, and 5.0 ¹ . Ultimate tensile strength, U T S , represents the maximum stress a material can withstand before breaking when stretched or pulled. During the optimization phase, twenty U T S values were used, ranging from 2.34 to 57.1, corresponding to common materials used in medical device prototyping. O R F C and O R F C max are the effective and geometric orifice areas, respectively. O R F C is calculated based on FEM modeling, while O R F C max is equal to π r 2 , where r is the radius of the PHV. We note that L M N is equal to L M N s , which is obtained directly from trained models. This is due to the fact that L M N does not require further scaling as it falls within the interval used for quantitative design comparison. To ensure that design comparisons are unambiguous, all design scores, D S N , L M N s , and S T S s , were normalized in the range of 0 , 1 .

To verify the effectiveness of the algorithms in determining the optimal geometry of a designed device, numerical modeling was utilized for performance evaluation. This process involved selecting 500 random combinations of geometric parameters (100 geometries per each α value) from the datasets obtained during optimization and using ML models to predict LMN and STS values. These combinations were then analyzed by numerical modeling using the method outlined in “Data Collection”. The criterion for accuracy, in this case the predicted versus modeled results, was determined by calculating RMSE.

When examining differences between samples, such as predicted vs. modeled, a paired t-test was employed as the samples were considered dependent and used a common set of geometric parameters. Differences were considered reliable at a significance level of p < 0.05. The analysis was performed using the Python library “SciPy”, which specializes in statistical processing.

The dataset utilized for training consists of combinations of geometric parameters within predefined ranges, in addition to their corresponding results obtained through the Abaqus/CAE simulation environment. The distribution of the features, including HGT, DIA, ANG, CVT, THK, and ELM, was determined to be uniform due to the random nature of their generation, as depicted in Figure 3. Comparatively, the targets in the dataset were represented by shifted ranges. The median for LMN was 0.55, with an interquartile range (IQR) of 0.05–0.73 and a minimum-maximum range of 0.002–0.999. For STS, the median was 1.13, with an IQR of 0.23–2.67 and a minimum-maximum range of 0.04–12.81. Analysis of the data set revealed the existence of highly effective designs, with 175 of them having LMN values greater than 0.95. It is worth noting that 25% of the total dataset (2,892) failed to exceed 0.05 in terms of LMN. A similar pattern was observed for STS, with the majority of the sample represented by small values; specifically, 5,496 results did not exceed 1.0 MPa. However, there were also instances of high STS values, with 512 geometries exceeding the threshold of 5.0 MPa.

The AutoML framework was successfully implemented and used to train 341 models for LMN prediction and 347 models for STS prediction. The results indicate that the lowest RMSE was achieved by ensembles, with an RMSE of 0.126 for LMN and 0.121 for STS. Comparison with other models showed significant inferiority in terms of error, as illustrated in Figure 4. Further analysis and optimization efforts will focus on the use of ensemble-based models.

The predictive performance of the obtained ensembles on the validation subset was evaluated by comparing the predicted to actual values, as shown in Figure 5A, and by examining the standardized residuals, as shown in Figure 5B. Ideally, the results should be closely aligned with a regressed diagonal line (Figure 5A) or horizontal zero axis (Figure 5B). However, further analysis revealed a relatively large error for 64 out of 2,313 designs (2.8% of the total), exceeding the limits of three standard deviations. Despite this, statistical analysis showed no significant difference between actual and predicted samples ( p = 0.151 ).

The evaluation of the STS score (Figure 6) showed similar results to the LMN ensemble. However, the number of geometries exceeding the three standard deviation threshold was lower, with 35 out of 2,313 designs (1.5% of the total) observed. Statistical analysis revealed no significant difference between the actual and predicted samples ( p = 0.092 ).

To formally evaluate the performance of the ensembles, various metrics were estimated on the training and validation subsets. Table 2 presents a summary of the results, including five error-based metrics that provide an average measure of the proximity of the predictions to the expected values. Furthermore, Supplementary Appendix SB provides a comprehensive set of regression metrics that fully describe the ensembles obtained and the distribution of RMSE across different ML algorithms used during training. Detailed descriptions and equations for these metrics can be found in the publicly available project repository ² .

In recent years, the use of ML models has become increasingly widespread, but the interpretability of these models remains a significant challenge. One method that has been proposed to increase the transparency and interpretability of ML models is the SHapley Additive exPlanations (SHAP) method (Lundberg et al., 2017). SHAP scores are based on cooperative game theory and have been widely used to understand the main features that affect the performance of a given model.

In this study, we applied SHAP to evaluate feature importance of two ensembles using a validation subset. To provide a clear visualization of the results, we employed a bee swarm plot, which is designed to display an information-dense summary of how top features in a dataset affect the model’s output. In the bee swarm plot, each PHV geometry in the given explanation is represented by a single dot on each feature row. The X position of the dot is determined by the SHAP value of that feature, and dots “pile up” along each feature row to show density. In addition, a blue-purple-red color scale was used to display the original feature values. It is important to note that while the SHAP values provide insight into the contribution or importance of each feature to the model’s prediction, they do not evaluate the quality of the prediction itself. The results of our analysis are shown in Figure 7.

After conducting an analysis of the swarm plots for the LMN and STS ensemble models, our findings indicate that both models have relatively similar feature importance. Specifically, DIA and ANG were found to have a positive contribution to the prediction output when their values were high, and a negative contribution when their values were low. Conversely, THK and ELM were found to have a positive contribution to the prediction output when their values were low, and a small negative contribution when their values were high. Furthermore, CVT was found to have a weak prediction contribution whilst HGT had almost no contribution in both ensemble models. The global feature importances are depicted in the feature rows in Figure 7. For further analysis, the distribution of SHAP values for each feature is presented in Supplementary Appendix SC.

We used a Hyperband-based pruning mechanism in conjunction with six optimization algorithms to identify optimal PHV geometries. Here, the search space of explored parameters was strictly defined to match modern manufacturing standards, to ensure applicability in medical device manufacturing, as described in Table 1, and omit parameter values that are unfeasible in a real setting.

Figure 8 showcases the results from optimal design exploration. Each optimizer performed 2000 iterations per alpha value, using five different alpha values: 0.2 and 0.5 (optimization prioritizing STS); 1.0 (balanced optimization between LMN and STS); 2.0 and 5.0 (optimization prioritizing LMN). In total, each optimization algorithm evaluated 10,000 designs with different characteristics (2000 iterations × 5 alpha values). As shown in Figure 8, designs prioritizing STS (alpha values of 0.2 and 0.5) converged faster to higher design scores. In contrast, designs prioritizing LMN converged more slowly. Balanced designs (alpha value of 1.0) showed intermediate convergence dynamics whilst the best designs studied converged to a similar range of scores, with DSN above 0.90.

As part of our optimal design exploration, we evaluated two methods for assessing the importance of hyperparameters: Mean Decrease Impurity (MDI) and Functional ANOVA (fANOVA) (Hutter et al., 2014). Both methods use a random forest regression model to predict the objective value based on a given parameter configuration. The accuracy of this model directly affects the reliability of the importance scores provided by the evaluators. For both evaluators, the number of trees in the forest was set to 64, and the maximum depth of the trees was set to 64.

The results of the feature importances, as depicted in Figure 9, indicate that MDI and fANOVA results are similar in nature to those of SHAP. All optimization algorithms consider Young’s modulus (ELM) and PHV thickness (THK) as the most influential parameters with a significant impact on the optimization process. In contrast, four other geometry parameters including ANG, CVT, DIA, and GHT were found to have little effect.

The quantitative evaluation on the test subset is of utmost importance during the in-silico verification of the optimal designs. The results, as presented in Table 3, indicate that the metrics for the geometries obtained during the test phase are comparable to the values obtained during the training and validation phases. A comprehensive set of metrics designed to assess the accuracy of LMS and STS predictions for the three sample groups (training, validation, and test) is detailed in Supplementary Appendix SB.

The optimization algorithms themselves have demonstrated remarkable performance by selecting geometric parameters of the leaflet that result in high DSN scores. Despite a shift in the balance of optimality towards either LMN or STS, all algorithms were able to identify geometries with high integral DSN scores in the range of 0.86–0.97.

A qualitative examination of the optimization results confirms the efficiency of the algorithms. All algorithms effectively selected parameter combinations that resulted in optimal leaflet opening while minimizing stress, as demonstrated in Figure 10 and Supplementary Appendix SD. The optimal leaflet models typically exhibit high LMN opening areas while maintaining moderate STS values, which in most cases do not exceed 2.0 MPa. This is well below the strength limit of the material models. It is noteworthy that the optimization algorithms show similar trends in selecting the best geometric indicators. The CMA, NSGA, QMC, and TPE algorithms selected larger diameters (31.0–32.2 mm) as optimal, minimized the leaflet lift angle for a better opening, and showed a similar thickness (0.30–0.34 mm). However, the parameters of HGT and CVT varied significantly among the optimizers and no general trend was observed.