Comprehensive Geometric Parameterization and Computationally Efficient 3D Shape Matching Optimization of Realistic Stents

Shape representation is critical in design optimization as it defines the search space of possible geometries. Vast application areas emerge, including the design and optimization of the most commonly used medical device, stents [1]. However, to generate findings for meaningful translation, representation schemes must be flexible and compact to generate a wide range of geometries and maximize memory utilization, resolution, and accuracy, especially when coupled with optimization methods. The early efforts on shape representation focused on defining 2D shapes, using either explicit or implicit methods. In explicit methods, using a points array, the coordinates of boundary points are either restricted to specific grid points/pixels or selected from a bounded interval within the design space. The explicit schemes are not compact because of the usually large number of pixels or points required to define the geometry. Moreover, customized point perturbation mechanisms are needed to derive shape variants, which significantly reduce the scheme’s flexibility. On the contrary, implicit methods are more flexible and compact but also have more complex representations. Different implicit methods have been used to represent 2D shape boundaries via, e.g., chain coding, polygons, or different splines. In chain coding, a curve is defined using a sequence of angles between the unit length line segments [2], and in polygon representation, the curve is defined using an ordered set of vertices [3]. The former can represent open and closed curves, while polygon representation can only describe closed curves. For complex curves, both chain coding and polygon representation schemes tend to be less compact as numerous segments or vertices need to be defined. To balance the needs of compactness and flexibility, the use of parametric spline forms such as Bezier curves [4], B-splines [5,6], and their variants, e.g., non-uniform-rational-B-splines (NURBS) [7] is suggested. B-spline-based schemes offer additional control over continuity and differentiability along the length of the curve. The number of control points, their coordinates, and the order of the curve allow for global (Bezier curves) or local (B-splines) shape control. This advantage makes them the most common schemes adopted in shape optimization. Domain-specific custom shape parameterization schemes exist, such as the PARSEC representation of airfoil shapes [8], where the variables defining the shape directly influence the performance (lift and drag in this case).

The representation of 3D shapes is far more involved than 2D, with both explicit and implicit methods reported in Refs. [9,10]. There are three forms of explicit representations, relying on either points, voxels, or the geometry mesh. In point-based representation, a point cloud describes the 3D shape, while voxel-based means assigning value to small 3D grids of cubes to delineate the shape’s volume and boundary. In mesh-based representation, the 3D shape is defined using edges and faces, in addition to the points that form the vertices. Each approach has shortfalls; point-based methods struggle to describe complex topologies, and voxel-based representations require large memory and complex screening steps to generate valid geometries. The management of vertex positions and continuity of faces is non-trivial in mesh-based methods. In addition to these explicit schemes, several implicit schemes have also been suggested in Refs. [10–12], including B-spline surfaces, summation of primitives, and the use of continuous distance functions. While using B-spline surfaces offers a good compromise between compactness and flexibility, manipulating the control points to generate valid geometries is challenging for complex shapes. Complex 3D shapes can also be generated via summation of simpler primitives, and the primitives themselves can be selected from a variety of shapes, e.g., implicit algebraic surfaces [12], superquadrics [13], etc. While such schemes offer a memory-efficient representation of 3D shapes, the range of topologies generated can be restricted due to its reliance on boolean operations of relatively simple primitives. Continuous distance functions can also be used to define shapes where the function returns the distance of a point from the surface, with a zero-level set lying on the surface of the 3D shape. While we can trace the origins of distance functions to the classification of points inside or outside an object using radial basis function kernel support vector machines [14], several signed distance function models [10] and extensions have been proposed in recent years that can also handle non-watertight models [15]. Distance functions can be learned using neural networks, and theoretically, one can achieve an infinite-resolution representation of a shape with minimal memory requirements.

Optimization methods must be coupled with the chosen shape representation scheme to uncover designs with excellent performance metrics. For practical applications, the performance metrics are typically derived using computationally expensive simulations such as finite element methods (FEM) or computational fluid dynamics (CFD) analysis. Thus, there is a limit to the number of designs that we can evaluate using such simulations during the optimization course. Surrogates/approximations are usually coupled with the underlying optimization schemes to work within a limited computing budget. Before embarking on an actual optimization exercise, the flexibility of the representation scheme and the efficiency of the optimization algorithm are usually assessed using a shape matching exercise. In shape matching, the aim is to identify a shape that resembles a target shape [16,17]. The vast majority of shape matching studies have focused on 2D shapes [17,18], or generic geometries such as that of a fish [19], a dolphin [17], etc. The 2D shapes themselves were parameterized using B-splines with either a fixed number of control points [17,19], doubling of control points until stagnation during optimization [20], or increasing the number of points by relying on adaptive approaches to control the fidelity of the representation [6,16]. Only a handful of reports described the attempts to use B-spline control polygons to match 3D shapes [21]. The quality of the match was assessed using Hausdorff distances [17,22] or a maximum of Euclidean and Hausdorff distances [19,21]. In all of the above studies, the optimization algorithm was directly coupled with the shape representation scheme, and the number of design evaluations ran from thousands to hundreds of thousands [17,18,21,23]. In real-life settings, however, one cannot afford such a large number of evaluations, and there is a need to couple approximations/surrogates within the optimization algorithm. To the best of the authors’ knowledge, no reports of 3D shape matching with a limited computing budget exist.

Efficient shape representation and optimization methods are highly sought-after as they have applications across multiple domains. Improved design of coronary stents is one such example which can have far-reaching effects since stents are the most commonly used medical device. In fact, coronary stenting is the preferred treatment for coronary artery disease, the leading cause of morbidity and mortality in the world [24]. A coronary stent is a cylindrical mesh structure that is placed at a disease-narrowed coronary artery before expanding it. A retrievable medical balloon is used to permanently scaffold the narrowed artery and restore sufficient blood in the coronary arteries to supply the heart muscle. Commercially available stent geometries today are either (a) independent ring (IR) or (b) helical stent (HS) designs, whereby IR are marked by connectors linking individual rings or elements, and HS are comprised of a single or double helix structure with connectors providing additional axial connections [1]. IR designs [1,25] have been available for some time, whereas HS designs are relatively recent additions [26,27], with double helix Orsiro (Biotronik, Inc.) receiving the United States Food and Drug Administration approval in 2019. All stent optimization studies reported in the literature to date are limited to either idealized IR stent designs [28–30] or outdated IR designs that are no longer in clinical use [31,32], and they concern a limited number of design variables (⁠$3−4$⁠) [33,34]. In most cases, the shape of the stent strut centerline has been simplified using, for example, straight lines and semicircular crowns [31,32,35,36] or represented using splines [30,37]. Commonly considered variables include strut width [34,36,38], strut thickness [30,39], or strut height [34,35,39]. Other relevant variables such as the strut cross-sectional type, number of connectors, and row spacing have only recently been considered in optimizing idealized stent designs of up to seven design variables [28]. However, important design variables such as the number of connectors, which are known to play a critical role in clinical stenting success by governing arterial apposition [40], for example, have not been considered in the parameterization of contemporary stent designs.

Interventions with stents continue to face challenges [41,42], and the latest improvements in drug-coating, delivery systems, hemodynamic optimization, and material developments can significantly benefit from the development of a shape parameterization scheme to accurately represent current and future market stent designs to optimize across multiple performance considerations of stents [1]. To maximize the effectiveness of our approach, we aim to allow a modification of the parameters such as the number of connectors and the stent shape, together with standard stent optimization variables, including strut length, width, and thickness—to account for all design characteristics affecting a stent’s performance. Moreover, with the HS designs being a new type of commercially available stent designs, we aim to develop a shape parameterization scheme for these and other possible classes of stent design to uncover superior design variations that could previously not be captured. Therefore, the key contributions of our work to the field are as follows:

1. Design space: We present the first complete 3D shape parameterization scheme of contemporary IR design with 17 relevant size and shape parameters (compared to the previous four or fewer variables for outdated designs, and seven or fewer for idealized contemporary designs in literature to date).

2. Design class: We present the first shape parameterization scheme of a contemporary double helix HS design, which has never been attempted before. The geometry is defined using 18 design variables via a combination of NURBS, PARSEC, quarter circle, and straight line segments.

3. Surrogate-assisted 3D shape matching: We perform a computationally efficient surrogate-assisted 3D shape matching optimization to demonstrate the utility of the developed parameterization schemes. Shape matching optimization has only been used to validate parameterization and optimization techniques in settings with large evaluation budgets, hindering its meaningful application to real-world optimization problems (such as coronary stent design optimization), which are accomplished within limited objective functions. Here, we overcome this by demonstrating 3D shape matching application for coronary stents with a limited function evaluation setting using surrogate-assisted optimization.

The paper is organized into five sections. Section 1 provides the background of the work and contributions. Section 2 describes the parameterization schemes for the IR and double HS designs. Section 3 provides an overview of the surrogate-assisted optimization algorithm and details of the 3D shape matching problem with a pre-defined objective function evaluation limit. Section 4 discusses and compares the results obtained from three different optimization schemes. Section 5 concludes the findings and discusses the limitations of the current work and its possible future extensions.

In this section, we outline the methodology for parameterization of contemporary IR and double HS design.

The geometry of the IR stent is parameterized using 17 independent variables as presented in Fig. 1. The stent structure consists of two primary components—the end and main rows. The repeated main row is divided into three sub-components—connector, short strut, and tall strut. Among the 17 independent design variables, 10 describe the physical size and shape of the stent components, i.e., stent length (SL), crimped stent radius (CR), number of rows (NROW), end row strut height (HTR), strut centerline width (WS), number of connectors per row (NCON), connector lateral radius (RC), connector lateral length (LC), short strut height (HS), and tall strut height (HT). Three of the 17 independent variables (i.e., radius ratio (RR), stent width factor (⁠$βw$⁠), and stent thickness factor (⁠$βt$⁠)) are used to derive additional dependent variables using constraint equations. The remaining four variables represent the weights of NURBS [43] that control the shape of the strut centerline splines of the end and main rows. The baseline values of these variables used as the target for the shape matching are listed in Table 1.

The end row of the stent (Fig. 1—bottom left) comprises repeating crowns at the sides and struts in the middle. The crown strut centerline is modeled by a semicircle (radius RT). We define the strut centerline by a NURBS curve of degree 7 with eight control points (CPs) and four weights W1–W4 applied to CPs $1−4$ and $8−5$⁠, respectively. The control polygon’s WS and HTR define the spline geometry completely. A PARSEC [44] connection between the strut and crowns maintains $c0$ and $c1$ continuity. The main row (Fig. 1—bottom right), as described earlier, has three sub-components—connector, short strut, and tall strut. All three have a strut at the middle, defined using the same NURBS spline as the end row, but with different control polygon of heights HS, HT, and connector strut height (HC), respectively. Again, the strut ends are connected to other components through PARSEC to maintain $c0$ and $c1$ continuity. The connector top is composed of a lateral region defined using radius RC and straight length LC, whereas the crown at the connector bottom is a semicircle (radius RO). The short strut’s top crown has a larger radius (RS), while all other crowns for the short and tall struts have the same smaller radius (RO). The complete stent structure has constant width (SW) and thickness (ST).

The stent length, cylindrical rows, and the definition of radius ratio provide four equality constraints that must be satisfied to construct the stent geometry successfully. The constraints are listed as follows:

• Total length constraint

  $2RT+HTR+(NROW−1)×(RO+4RC+HC)=SL$

  (1)
• End row cylindrical constraint

  $4NCON×RT=2πCR$

  (2)
• Main row cylindrical constraint

  $2NCON×(4RO+RS)=2πCR$

  (3)
• Radius ratio constraint

  $RS=RR×RO$

  (4)

The geometry of the double HS design is parameterized using 18 independent variables as presented in Figs. 2 and 3. The stent structure consists of three primary components—the main helices with connectors, the end row struts with connectors, and the transition rows connecting the main helix to the end row. Among these 18 design variables, 7 of them, i.e., (SL, end row length (EL), CR, connector strut centerline width (CW), number of unit struts (NUS), ST, and SW) describes the size and shape of the main helix and its connectors. The remaining 11 variables represent the weights of NURBS [43] that control the shape of the strut centerline of the main helix and the connector. The baseline values of these variables based on contemporary ultra-thin stent designs are listed in Table 2. Since our focus is on the parameterization of the main helix and the connectors, the axial and circumferential dimensions of the end rows and the transition rows are set to be linearly dependent on the EL and the CR of the stent.

We model the helix geometry by repeating the unit strut, with each new unit shifted by distance AAX and AAY in the circumferential and axial directions, respectively, as presented in Fig. 3. The unit strut comprises four sub-segments: AB, BC, CD, and DA. The AB and CD sub-segments are identical and defined by a NURBS curve of degree 7 with eight CPs and four weights W1–W4 applied to CPs $1−4$ and $8−5$⁠, respectively. The control polygon’s width (WAB) and height (HAB) define the spline geometry completely. The DA sub-segment is subdivided into a middle spline and quarter-ellipse crowns with semi-axis lengths of $a$ and $b$⁠. The middle spline portion of DA is represented by a NURBS curve of degree 5 defined using six CPs with three weights W5–W7 applied to CPs $1−3$ and $6−4$⁠, respectively. Again, the width (WDAN) and height (HDAN) define the control polygon of the spline entirely. The middle spline of sub-segment BC is obtained by forcing a zero-width condition on the NURBS segment DA through a linear width correction. We introduce an intermediate PARSEC [44] correction to maintain $c0$ and $c1$ continuity at the intersection of the quarter-ellipse crowns and BC spline. Finally, the connector is represented by a NURBS with characteristics similar to the AB sub-segment but with a different control polygon shape. As with other NURBS defined earlier, CW and connector height (CH) define the CP’s position in the control polygon.

The helix length, helix angle (⁠$α$⁠), and the relative arrangement of the two helices provide four equality constraints that need to be satisfied to construct the stent geometry successfully. The constraints are listed as follows:

• Total length constraint

  $(3NUS+1)2×AAY+HAB=SL−2EL$

  (7)
• Helix angle constraint

  $tanα=AAYAAX$

  (8)
• Circumferential alignment constraint

  $3AAX+2WC−WDAN=2πCR$

  (9)
• Axial alignment constraint

  $3AAY=HC$

  (10)

Four dependent variables, HAB, $a$⁠, $b$⁠, and HC, were analytically obtained by solving Eqs. (7)–(10). In addition to these four, we still need the values of four other dependent variables, $α$⁠, WDAN, HDAN, and WAB to construct the stent geometry. These were obtained by solving a single-objective-constrained optimization problem using sequential quadratic programming (SQP), referred to as repair here. Based on the contemporary stent designs, we set the initial guess values for these variables as (⁠$38deg,0.125mm,$$0.675mm,0.3mm$⁠). The objective of the optimization problem was the minimization of the difference between HAB and HC, subject to constraints listed in Eqs. (11) and (12). The variables were allowed to vary between $±30%$ of the initial guess values.

Based on the above methods of geometry parameterization, a 3D stent geometry builder was developed by coupling matlab 2023b (Mathworks Inc.) programming interface with Solidworks 2024 API (Dassault Systems Inc.) through python 3.9 and pywin32. The 3D stent geometry builder would be fed with independent design variables (rounded to five decimal places) as inputs, and it would, in turn, generate a neutral Parasolid CAD file and a stereolithography (STL) file of the resultant stent geometry. Building a stent geometry takes about 136 s for the IR stent and 197 s for HS on a workstation with Intel(R) Xeon(R) Gold 6226R CPU @ 2.90 GHz, 2.89 GHz processors and 128 GB of installed RAM.

Shape matching is often performed as the first step before undertaking a real design optimization exercise involving FEM analysis or CFD tools. A shape matching optimization exercise aims to identify a geometry that resembles a given target geometry, and the process is depicted in Fig. 4. To successfully solve a shape matching optimization problem, one would need to couple a flexible shape representation scheme with a computationally efficient optimization algorithm. Since we need to operate within a limited computing budget scenario, we rely on a surrogate-assisted believer and expected improvement (EI) maximization formulations that are regularly used in applications with computationally expensive analysis. It is important to note that such formulations/schemes have never been attempted for 3D shape matching problems.

The minimization problem is constrained by the hidden feasible design space $Ωf$⁠, which represents the regions in the design space where the candidate design geometry can be generated using the geometry builder. A design that cannot be successfully constructed using the geometry builder is thus referred to as an infeasible design. Infeasible designs may be a result of various factors, e.g., the limitations of the computer-aided design software, such as singularities during geometry construction, inability to generate complex splines or high-curvature surfaces, etc. The objective function for shape matching is the norm of the distances between each point in the target design and its nearest point in the candidate design. Such a measure is similar to the inverted generational distance used in assessing the performance of multiobjective optimization algorithms.

In order to get an objective assessment of the objective function landscape, $20D$ samples were created using Latin hypercube sampling (LHS), where $D$ denotes the number of independent variables of the optimization problem (i.e., 17 for the IR design and 18 for the HS design). Among the 340 and 360 designs generated by LHS for the IR and the HS designs, the geometry builder was successful in constructing 267 and 195 feasible designs. We assigned an objective function value of $1×106$ to all infeasible solutions, i.e., ones that failed in the geometry construction phase.

With the information that $21%$ and $46%$ of the design space is likely infeasible based on results from LHS sampling, we investigated the performance of the following three approaches with a limited number of design evaluations (set as 200), (a) gradient-based optimization (hereafter referred as direct optimization (DO)) with the best LHS solution as the initial guess (starting solution), (b) a surrogate-assisted believer optimization (referred as predictor believer (PB)) approach which incorporated a classifier to predict feasibility of a solution and a regressor to predict its performance, and (c) a surrogate-assisted Bayesian optimization approach (referred as EI approach) wherein an acquisition function (EI) is maximized [45], taking the uncertainty of the predictor model into account to suggest a new candidate solution.

We used SQP implementation of matlab optimization toolbox for the gradient-based DO approach with step size set to $1×10−4$⁠. The surrogate-assisted optimization framework is presented in Algorithm 1. In brief, the classifier was trained using all LHS samples, while the regressor was trained using only feasible solutions. The optimization approach followed a steady-state paradigm where a single solution was evaluated in each step, followed by retraining both the classifier and the regressor models. We employed six different classification models—logistic regression, random forest classifier, Gaussian Naive Bayes, nearest neighbors classifier, C-support vector classifications, and gradient boosting classifier from the Scikit-Learn python module [46]. The classifier models were trained with 80% data using default hyperparameters. The classifier with the highest accuracy score on unseen test data (⁠$20%$⁠) was selected as the predictor. Thereafter, the chosen classifier was retrained on the complete data. For the regressor, we used eight different models—linear regression, Bayesian ridge, $ϵ$ support vector regression, stochastic gradient descent, Gaussian process regressor with radial basis function kernel, random forest regressor, K-nearest neighbors regressor, and gradient boosting regressor from the Scikit-Learn python module [46]. The training data for the regression models only included feasible solutions. The regressor training also followed the $80−20$ train-test split, and the model with the least mean squared error on unseen test data was selected as the predictor in each iteration. Again, the selected regressor was retrained on the complete training and test data. Finally, we evaluated the standard deviation for the EI approach using the bootstrapping technique, wherein 100 regressor models with the same characteristics as the selected regressor were trained on resampled data with replacement. The splitting of the data using the train-test split may affect the choice of the surrogate model in an iteration. However, the train-test split is repeated for each iteration of optimization. In every iteration, the test-train split method performs a random shuffling of the data to generate a training and test set. This shuffling is expected to minimize the effect of the splitting technique on the overall optimization results after 200 iterations. In the past, the selection of surrogate models for EI-based formulations was limited to the regressors that inherently provide uncertainty estimates, such as the Gaussian process regressors. The introduction of the bootstrapping technique in our framework allows error estimation for any regressor, thereby expanding our choice of regression models for EI-based optimization. For the optimization exercise, differential evolution, i.e., DE/rand/1/bin [47] scheme, was used with an initial population size of $20D$ and a crossover constant of $0.3$⁠. We used the Pymoo open-source python module [48] for the optimization task with default termination criteria. For the PB approach, the classifier output was added as a constraint, and the predicted performance from the selected regressor was minimized. We selected the candidate solution with the best-predicted performance for actual evaluation and re-trained the surrogate models after every new design evaluation. The process continued until 200 additional designs were evaluated. For the EI approach, the expected improvement metric was obtained using the mean performance from the selected regressor model, and the standard deviation was evaluated from the bootstrapped regressor models. The expected improvement metric was multiplied by the probability of feasibility from the selected classifier model, and the resulting metric was maximized by the optimization algorithm to obtain the suggested design solution.

Among 340 designs generated by LHS, $79%$ resulted in feasible designs (successful in constructing 3D geometry), with the best design having an objective function value of $771.4$ (Table 3, column 2). The sampled feasible solutions are marked with hollow circles in Fig. 5. To visualize the difference between the target and the generated candidate designs, a normalized distance metric, i.e., $norm(x→candidate−x→targetx→target)$ was used in the figure. A clear separation is visible in the normalized distance of three-connector designs from two- and four-connector designs. Both three- and four-connector designs show low objective function values, while the two-connector designs have a higher objective function value. The distribution of the feasible LHS samples is presented in Fig. 6(a).

To observe the performance of the gradient-based DO approach, the best LHS design was used as the initial guess/starting point. The search resulted in a design with an objective function value of $454.4$ (Table 3, column 3) within 200 function evaluations. During the course of the search, all sampled solutions were feasible three-connector designs and are marked with dots in Fig. 5. The distribution of the obtained solutions (presented in Fig. 6(b)) shows the concentration of the design points around objective function values of 600. This indicates that the optimization process started from the best LHS point and reached a local minimum. The gradient-based search from the best LHS design offered a performance improvement of $41.1%$ within a limited budget of 200 for additional design evaluations.

Next, we assess the performance of the surrogate-assisted optimization approaches PB and EI. Since a stochastic population-based algorithm (DE) is coupled with the surrogates to perform optimization, multiple independent runs with different seeds are necessary. In our case, we performed 11 runs with different seeds, and the median, best, and worst objective function values are listed in Table 3. All objective function values reported in the paper are true values obtained by evaluating the shape matching objective.

The PB approach resulted in objective function values varying between $379.4$ and $786.8$⁠, although its median performance (⁠$774.3$⁠) was worse than the best LHS value of 771.4 (Table 3). The feasible solutions sampled by PB for the median run are marked with squares in Fig. 5(a). During this process, 154 feasible solutions were evaluated. Within the total 200 evaluations, the PB approach mostly sampled the design space with four-connector designs (175 designs) and only sampled 25 three-connector designs (see Supplemental Fig. 1 available in the Supplemental Materials on the ASME Digital Collection). No two-connector designs were evaluated by this approach. The result shows that the believer model was fixated on a local minimum of four-connector designs. The distribution of PB solutions in Fig. 6(c) reveals that most evaluations led to objective function values of around 900 and repeated evaluation of four-connector designs.

Both DO and PB approaches are expected to face difficulty if the function landscape is highly multi-modal. The EI approach is likely to be slow but is capable of dealing with such landscapes by virtue of its uncertainty considerations. The median EI result corresponds to an objective function value of $416.6$⁠, which translates to $8.3%$ improvement over DO (Fig. 5(b)). In contrast to the DO and PB techniques, the EI approach sampled 2,170, and 28 evaluations of two-, three-, and four-connector designs (refer to Supplemental Fig. 1 available in the Supplemental Materials). One can also notice that the distribution of designs shows a much wider variation in the objective function value (refer Fig 6(d)). Additionally, the EI approach had a higher number of feasible design evaluations at 183 compared to 154 evaluations by the PB method in their median run.

One can also observe the rate of convergence of DO and those of PB and EI for median runs (refer to Supplemental Fig. 2 available in the Supplemental Materials). We observe that the PB approach does not show any improvement over the LHS results, whereas the DO takes the complete 200 iteration budget to reach its minimum objective function value of $454.4$⁠. The EI approach shows the best convergence behavior by achieving its minimum in 141 iterations. A comparison of computational costs for different components of the DO, PB, and EI techniques in a typical run is provided (refer to Supplemental Table 1 available in the Supplemental Materials). We observe that the total objective function evaluation time (28, 984 s, 34, 678 s, 29, 395 s) is significantly higher than the optimization time (⁠$195s$⁠, $1,081s$⁠, $5,831s$⁠) for all three algorithms (DO, PB, EI). The optimization refers to only optimization time in DO, surrogate prediction and optimization time in PB, and surrogate prediction with bootstrapping and optimization time in EI. The optimization time increase between PB and EI is due to the bootstrapping technique used to estimate the standard deviation for expected improvement calculation. The total surrogate creation cost in PB was $991s$⁠. The addition of bootstrapping in EI increased the cost to $4,449s$⁠. In terms of objective function evaluation times, we observe a significant difference in the average evaluation time for feasible IR stent designs across different techniques: DO (⁠$144.9s$⁠), PB (⁠$202.3s$⁠), and EI (⁠$153.5s$⁠). We attribute the difference to the fact that the DO and EI algorithms mostly generated three-connector stents, while the PB primarily produced four-connector stents that require a longer generation time due to more features in the stent design.

Regarding the choices of the classifiers in the PB approach, one can notice that all classifiers except the K-nearest neighbor were selected with different frequencies (refer to Fig. 7(a)). Gradient boosting classifier was most frequently selected (86 out of 200), and Gaussian Naive Bayes was the least favored with ten selections. The EI approach followed a similar selection with all six classifiers chosen during optimization. Again, the gradient boosting classifier was most frequently selected in 76 function evaluations (refer to Fig. 7(b)), while the Gaussian Naive Bayes was least selected, being used in only 9 out of 200 evaluations. As for the regressors, only one regression model, i.e., gradient boosting regressor, was selected in all iterations of both PB and EI maximization approaches. Regarding the surrogate metrics, the classifier accuracy score for PB and EI in the first iteration was $0.9$⁠, as both surrogates were trained on the same LHS data. By the last iteration, the accuracy score decreased to $0.86$ for PB and increased to $0.91$ for EI. The regressor root mean squared error (RMSE) was $131.97$ for both approaches in the first iteration. At the last iteration, the RMSE decreased to $100.72$ for PB and increased to $141.2$ for EI. It should be noted that surrogate metrics assess the overall quality of the model. These overall metrics may not improve during surrogate-assisted optimization, although the accuracy would be higher around the promising regions due to intense sampling.

The design with the worst (two-, three-, and four-connector stents) and the best objective function value obtained via LHS is presented in Fig. 8 along with the best designs obtained via DO, PB, EI approaches and the target design. The difference in the stent geometries with different numbers of connectors is visible from both the Stent top and the bottom view. The crown radii significantly decreased from two-connector stents to four-connector stents. The worst LHS designs have lower lengths for all cases and fewer rows for two- and three-connector designs. Within the best-matched designs, the best LHS and EI designs have wider struts than the target. The PB design significantly differs with four-connectors, while DO and EI approaches best match with the target. The best PB design is a four-connector stent that has significantly tighter packing of stent struts and higher metal volume than the target design. The overall match of the best PB design to the target is even worse than the best LHS design obtained through plain sampling of the design space. As for the best EI, it is the overall best match to the target with the closest stent length, number of connectors, strut shape, and strut dimensions.

Among 360 designs sampled by LHS, $54%$ of them resulted in feasible designs (successful in constructing 3D geometry), with the best design having an objective function value of $707.2$ (Table 4, column 2). The sampled feasible solutions are marked with hollow circles in Fig. 9. The objective function values of the LHS samples are distributed between 700 to 3800 with most designs around 1300 (refer to Supplemental Fig. 3 available in the Supplemental Materials).

To observe the performance of gradient-based DO, the best LHS design was used as the initial guess. The search resulted in a design with an objective value of 522 (Table 4, column 3) within 200 function evaluations. During the search, 168 feasible designs were evaluated, marked with dots in Fig. 9.

While the gradient-based search from the best LHS design offered a performance improvement of $26.2%$⁠, it was important to observe the benefits of the surrogate-assisted optimization approaches. The PB approach resulted in an objective function value of $286.1$⁠, providing a further $45.1%$ improvement over the DO approach. The objective function values varied between $174.5$ and 406.9 across 11 runs (Table 4, column 4). 113 feasible solutions were sampled by PB for the median run, marked with squares in Fig. 9(a).

Finally, the EI maximization approach showed a median objective function value of $347.1$⁠, providing a $33.5%$ improvement over the DO approach (Fig. 9(b)). The EI approach sampled 109 feasible solutions for the median run, and the distribution of solutions sampled is presented (see Supplemental Fig. 3 available in the Supplemental Materials). It is interesting to observe that, in this instance, the median performance of EI is worse than PB. This can be mainly attributed to two factors. First, the EI maximization method maintains a balance between performance and uncertainty, i.e., it attempts to sample in regions of high uncertainty, unlike the PB model, which samples in areas of best-predicted performance. Thus, if the EI algorithm samples a location close to the global minimum, the exploration component of the method will force the algorithm to sample the regions with higher surrogate uncertainty in future iterations. However, the PB approach quickly converges to the global optimum if the surrogate has identified the global minimum basin. This explains the higher median and standard deviation of the EI in this case. Second, the expected improvement value is multiplied by the probability of feasibility of the classifier model to discourage sampling solutions that may have low chances of being feasible. The accuracy of the classification model is also expected to impact the overall performance of the EI approach.

One can also observe the rate of convergence of DO and those of PB and EI for median runs (refer to Supplemental Fig. 4 available in the Supplemental Materials). We observe that the DO, PB, and EI approaches reach their minimum objective function values in 165, 188, and 153 iterations respectively, again highlighting the superior observed convergence of the EI optimization. A comparison of computational costs for different components of the DO, PB, and EI techniques in a typical run is provided (refer to Supplemental Table 2 available in the Supplemental Materials). We observe that the total objective function evaluation time (⁠$42,328s$⁠, $37,273s$⁠, $35,783s$⁠) is significantly higher than the optimization time (⁠$124s$⁠, $2,061s$⁠, $6,242s$⁠) for all three algorithms (DO, PB, EI). Again, we attribute the optimization time increase between PB and EI to the bootstrapping technique used to estimate the standard deviation for expected improvement calculation. The total surrogate creation cost in PB was $865s$⁠. The addition of bootstrapping in EI increased the cost to $4,120s$⁠.

Regarding the chosen estimators in the PB approach, all classifiers were selected with different frequencies (see Supplemental Fig. 5 available in the Supplemental Materials). Support vector classifier was most frequently selected (71 out of 200), and gradient boosting classifier was the least favored with two selections. The EI approach showed a similar selection with five of six classifiers chosen during optimization. In this case, the logistic regression classifier was most frequently selected in 115 function evaluations (refer to Supplemental Fig. 5 available in the Supplemental Materials), while the gradient boosting classifier was not selected in any iteration. As for the regressors, only two regression models, i.e., gradient boosting regressor (PB—176, EI—168) and random forest regressor (PB—24, EI—32), were selected in all 200 iterations of both PB and EI maximization approaches. The classifier accuracy score in the first iteration was $0.71$ for both PB and EI. By the last iteration, it decreased to $0.67$ for PB and increased to $0.75$ for EI. The regressor RMSE was $112.32$ for both approaches in the first iteration. It improved to $108.59$ for PB but worsened to $149.79$ for EI in the last iteration. Once again, these overall metrics may not improve during surrogate-assisted optimization, although the accuracy would be higher around the promising regions due to intense sampling.

The design with the worst and the best performance obtained via LHS sampling is presented in Fig. 10 along with the best designs obtained via DO, PB, EI approaches, and the target design. The 3D stent geometry builder generates the stent design in a left-to-right manner with the origin at the left end. Therefore, the stent appears similar at the left end. However, the shapes become visibly dissimilar toward the right end. The worst LHS design has significantly shorter length and wider unit struts. The stent length increases as we move from the worst LHS toward the best LHS and the optimized designs obtained by the three approaches. The metric shows a significant difference between the PB and EI approaches. The best PB design has a lower objective function value than the best EI design. There are two major reasons for this difference. The first is that the best EI design is slightly shorter than the target and the length has a disproportionate impact on the overall metric value. The second reason is that a slight variation in the weights of the spline for the best EI design affects the shape of the struts, moving them away from the target geometry and increasing the objective function value. The best PB design has the closest length and strut shape to the target design.

The results demonstrate the benefits of the EI-based approach over PB and the direct gradient-based optimization approach. Although the PB approach resulted in a better median objective function value than EI for the HS design, it performed poorly for the IR design, indicating the high possibility of the PB model getting trapped in a local minima basin. We observe that the exploration component of the Bayesian Optimization using the popular EI acquisition function allowed the optimizer to search for global optima in the three-connector IR stents, while the complete exploitation focus of the believer model limited the ability of the optimization algorithm to search for the global minimum, demonstrated by the excessive sampling of the four-connector IR designs. Additionally, we successfully used the bootstrapping technique to extend the choice of regressors from traditional Gaussian models to any regressor, demonstrated here by the use of linear regression, Bayesian ridge, $ϵ$ support vector regression, stochastic gradient descent, random forest regressor, K-nearest neighbors regressor, and gradient boosting regressors for EI maximization.

We introduced two detailed shape parameterization schemes to represent contemporary stent geometries, i.e., one each for IR and HS designs. Since the representations are compact and flexible, they can generate 3D stent geometries that can be optimized. For the IR design, we developed complete parameterizations via 17 variables (strut width, thickness, length, number of connectors and rows, and strut shapes controlled by NURBS parameters) that are known to affect stent performance. We presented the first shape parameterization scheme to represent double helix HS designs derived via a combination of NURBS, PARSEC, quarter circle, and straight line segments. This recent stent design has never been optimized before, and we hope that the present geometry parameterization will pave the way for optimized future designs.

The above shape representation schemes were further coupled with surrogate-assisted optimization methods to solve 3D shape matching problems, where the objective is to match the target geometry with a limited number of design assessments. It is important to note that while there are studies on shape matching of 2D and 3D geometries, very few use approximations/surrogates, which is necessary to solve such problems with a limited computing budget. Two commonly adopted formulations of surrogate-assisted optimization were used in the study, i.e., the believer model and the other based on expected improvement and probability of feasibility. The above performance was also objectively compared with a gradient-based optimization method to highlight their differences and strengths.

These methods pave the way for full-fledged optimization studies of stent designs with subsequent domain-specific objectives (CFD, FEM). We have individually optimized the designs, i.e., conducted shape matching with a target stent individually, and there is the opportunity to develop efficient surrogate-assisted optimization algorithms that can search across multiple stent design concepts in future works. The underlying surrogates in such models can be further improved using hyper-parameter tuning. Recent developments in machine learning can also eliminate the notion of design concept altogether, wherein neural network models can generate arbitrary stent designs using a common representation framework based on signed distance functions. The authors are currently pursuing some of these directions.

A.K. gratefully acknowledges the support from the Commonwealth Government through the Australian Government Research Training Program Scholarship. S.B. would like to acknowledge the support from the National Heart Foundation Vanguard grant.

N.J. declares grant and honoraria support from Abbott Vascular. The other authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

The stent 3D geometry builder codes will be made publicly available upon publication (DOI: https://doi.org/10.5281/zenodo.11368913) and GitHub repository available online.²

CFD =

computational fluid dynamics

CPs =

control points

D =

number of variables

DO =

gradient-based direct optimization

EI =

expected improvement maximization based surrogate-assisted optimization

FEM =

finite element method

HS =

helical stent

IR =

independent ring

LHS =

Latin hypercube sampling

NURBS =

non-uniform rational B-spline

PB =

predictor believer based surrogate-assisted optimization