ML‐Augmented Bayesian Optimization of Pain Induced by Microneedles

Microneedles (MNs) have emerged as a promising solution for drug delivery and extraction of body fluids. Pain is an important physiological attribute to be examined when designing MNs. There is no known representation of pain with geometric features of a MN despite the focus on experimental work. This study focuses on optimizing MN designs with the aim of minimizing pain through means of machine learning, finite element analysis, and optimization tools. Three distinct approaches are proposed. The first approach involves training multiple regression models on data obtained through finite element analysis in COMSOL. The second approach uses COMSOL's built‐in nonlinear optimization solver. Finally, the third approach utilizes the LiveLink interface between COMSOL and MATLAB, combined with Bayesian optimization. Each approach presents unique strengths and challenges, with the third approach demonstrating significant promise due to its efficiency, practicality, and time‐saving.


Introduction
9][10][11] What makes MNs even more unique is that the risk of infection due to bloodborne pathogens is low or near zero because the insertion depth is within the dermis layer and in many cases blood is not involved. [3,12,13]Additionally, MNs eliminate the reuse of surgical needles in many poor and underdeveloped countries suffering from shortages in medical supplies. [7,14]][17] One study has reported that MNs can expedite the transfer of drug into the outermost layer, stratum corneum (SC), without alerting the pain receptors.On the other hand, blood extraction and processing from the subcutaneous layer can be laborious, hard to handle, expensive, and painful when extracted because of the presence of pain receptors on site of insertion. [1,7,18][21] Nevertheless, drawbacks such as expensive production costs and susceptibility to fragility necessitated the exploration of alternative materials, leading to the proposition of polymers. [21,22]Polymers offer numerous benefits that align well with microneedle production, including established biocompatibility, biodegradability, and suitable mechanical properties for microneedle applications. [21]][28] Lately, a simplified fabrication strategy for MNs has been achieved using Additive Manufacturing (AM) & 3D printing which enabled a wide range of applications.It allows the entire process to be done in a single step; desired dimensions and size can be altered using the Computer Aided Design (CAD) software during the fabrication process, which aids in flexibility and modification of the initial draft design. [19,29,30]ain is an important physiological attribute that should be taken into account for MN designs.Yet, in 2001, a single study was performed to measure the pain due to MNs. [12,31] It was demonstrated that, in contrast to 2 mm depth insertion of 26gauge hypodermic needles, negligible pain was experienced for 400 MNs of 150 mm height. [12,31]Another related study observed that scraping the skin with MNs ranging from 50 to 200 mm in length resulted in minimal discomfort. [12,32]The conducted study by Gill, Denson, Burris, and Prausnitz [12] found that significantly less pain was felt for the shortest MNs and were 20 times lower compared to the hypodermic needles.The pain scores varied between 5% and 40% among different MN geometries. [12]hey initially restricted the insertion of hypodermic needles to a distance of 5 mm to decrease any bleeding although the reported clinical insertion depth is 8-12 mm into the skin.Thus, in clinical practice, the expected pain experienced by hypodermic needles is of a greater value and the associated pain in MNs are even less significant. [12,33]Furthermore, the insertion was done without injecting any substance and in the case of an injection, pain will be induced for hypodermic needles whilst solid MNs are not entailed for fluid injection, hence, no pain is sensed for MN's insertion and the relative reduction in pain is expanded. [12]he relationship between pain and some geometric parameters of MNs has been studied and reported by Gill, Denson, Burris, and Prausnitz [12] and Gill and Prausnitz [34,21] Pain is induced by larger MNs, although it can deliver more drug.Different MN lengths (480-1450 μm), widths (260-465 μm), thickness (30-100 μm), and other variables have been examined.Their findings indicate an escalation in pain is felt when there is a slight increase in MN length, while an increase in the number of MNs shows little impact on pain.Additionally, pain is linked to the tip sharpness instead of the overall angle of the needle.Thereby, the pain experienced relies highly on insertion force at the tip.For sharper tips, patients seemingly feel relief.No or small correlation on needle diameter [18,21] while needle thickness induces pain. [35]Moreover, MNs with lengths ranging up to 1.5 mm or less scored than 5-40% less painful than hypodermic needles; as such, an MN length of less than 1000 mm is recommended to reduce pain. [18]isual analog scale (VAS) is a well-known scoring tool that is considered a valid and assured scale. [1]Although pain is difficult to measure since it varies from one person to another and is considered as a phenomenon experienced, an accurate assessment device called McGill Pain Questionnaire (MPQ) is used to assess pain.MPQ incorporates a VAS to gauge and measure pain intensity through specific descriptive words and an evaluation index such as mild, discomforting, and excruciating Melzack. [36]t grants a qualitative and quantitative data. [1,37]During the VAS test, the participants are asked to assess the pain at the time of insertion of MNs and the infusion of fluid.A measurement ruler is placed labeled "No Pain" at the far left and "Worst Imaginable Pain" at far right then participants are asked to point a mark on the ruler scale for the pain score in the course of assessment for different types of MNs. [9]n this study, pain will be the primary parameter for MN design.Machine learning (ML) and optimization tools will be uti-lized to achieve an optimal design of a microneedle to reduce pain for various age group.Mainly, the workflow will be composed of three approaches as shown in Figure 1A; First, a large dataset is created to evaluate pain using finite element analysis (FEA) through COMSOL and then using ML in MATLAB to train the datasets to acquire the best regression model.The best regression model is used to find at which MN design is the minimum pain value.Second, a built-in nonlinear optimizer in COMSOL that configure a gradient-free algorithm is used to minimize pain and find the MN design by accurately formulating the optimization problem.Third, Bayesian optimization solver, bayesopt, is utilized by benefiting from the LiveLink option between COMSOL model script and MATLAB codes to explore the best microneedle design to minimize pain for each age group.Then, these approaches will be compared in terms of accuracy, time, and storage space.Finally, a graphical user interface app has been developed to be interactively used and explored through all age cases.

Microneedle Design
A concave profile MN was adopted to study design geometry.With the help of the Quadratic Bezier Curve (QBC) in COMSOL, smooth curves between points at the predefined starting and ending points produced the designated MN design.The control Bezier point (bz) helped defining the path of the curves through interpolating the points in a smooth manner which was what gave it the concave profile.Successively creating a QBC in 2D sketch could then be revolved around the vertical axis to construct a closed 3D model of the MN.
The MN consisted of four essential parts, as shown in Figure 2A; inlet diameter (in), outer diameter (out), length (len), and thickness (th).The inlet diameter defined where the target fluid to be extracted out from the skin while the outlet diameter acted as suction process.Both the inlet and outer diameter were controlled through the bz point.The wall thickness and length of the MN were also included in the geometry description.The lower and upper limit of the geometric values complied with the pain tolerance and fabrication limits.Additionally, Table S1 (Supporting Information) presents the geometric parameters thatare currently utilized and were previously used by Ref. [38,  39] Although the same geometric parameters were used, [39] the first improvement was transformation from 3D to 2D axisymmetric, which significantly reduced the simulation time and allowed to produce larger dataset of various MN designs, materials, and conditions (Section 2.4).The simulations were validated via a mesh convergence test in 2D and 3D models (Section 2.4).The second improvement compared to, [39] was introduction of Non-linear algorithm based on Nelder-Mead (or so-called Builtin COMSOL), which pathed a new way of running optimization problems for similar cases (Section 2.3.2).Another modification applied into the optimization problem was including the factor of safety (FOS) constraint between the allowable and applied compressive strength (Sections 2.3.2 and 2.3.3).Finally, an addition was made to approach 2 by including Monte Carlo algorithm, which turned out to work far better than the Nelder-Mead (Tables 6 and 7 and Sections 3.2-3.4).The algorithm could correctly model a constrained optimization problem.

Pain Evaluation
Pain was a difficult phenomenon to be measured as a single number, yet it could be declared by an absolute interrelation between MN's geometric relation and other properties.The evaluation of pain could be reasonably associated with the length, inlet diameter, and height as geometric parameters correlated with the flow rate within the MN as suggested before.This relation can be approximately described as a proportionality given by Equation ( 1): where Q is the volumetric flow rate within the MN during extraction of ISF (m 3 s −1 ), and P is the skin force applied on the MN (N) [39] which were obtained using the skin elasticity reported in Ref. [40].Table S2 (Supporting Information) provides the different cases of applied force corresponding to different ages and genders.
The physics problem behind the MN extraction was an integration of solid and fluid mechanics also known as fluid-solid interaction (FSI).FSI problems were intensively complex and nonlinear in nature and therefore required an FEM to be solved.The struggle also lied in the coupling of both solid and fluid mechanics, where the fluid motion affected the deformation of solid and vice versa.COMSOL was a great tool in processing these types of problems, especially through the built-in Multiphysics components.Thus, the MN problem could be split into two modules, solid and fluid mechanics within the COMSOL environment.In solid mechanics, the material of choice for the MN was Polydimethylsiloxane (PDMS) with 0.49 and 970 kg m −3 poisson's ratio and density, respectively.Fixed boundary conditions at the top surface and applied load force at the outer body of MN were established.Laminar flow for the fluid mechanics was defined due to the micro-scale physics and low Reynolds number associated with ISF extraction.The fluid volume for the fluid domain had been extracted from the solid MN geometry.Pressures at the inlet diameter 530 Pa was assumed due to lymph flow and 4500 Pa extraction pressure developed at outer diameter based on the human pulse.Clearly, faster fluid flow was linked to larger outlet diameter's pressure, therefore pain was reduced in exchange.Body forces and suppress backflow were also included into the problem modeling in COMSOL.

Approach
In this section, a breakdown of this comprehensive approach to minimize pain induced by microneedles was developed.The methodology was designed to address, compare, and find optimal approaches for the key aspects of the problem while considering potential challenges.The common base of the different approaches was using optimization techniques to explore possible solutions.The three approaches included: 1) Utilizing parametric sweep in COMSOL simulation accompanied by ML regression models in MATLAB.2) Using a built-in COMSOL optimizer that operates on gradient-free algorithms.3) Implementing COMSOL simulations LiveLink with MATLAB codes interface.Figure 1A summarizes the workflow of the different ap-proaches aimed at minimizing pain induced by MNs.The next sections explain the detailed workflow developed by each approach.

Parametric Sweep with Regression Models
The framework of using COMSOL simulation and then feeding the results into MATLAB's available ML regression models could be effectively used.The parametric sweep preference in COM-SOL could generate various MN designs as per defined geometric datasets.Those datasets were obtained during the design step.Next, FEA for all combinations were solved using a stationary study solver followed by the pain evaluation results.Afterward, the results obtained through the evaluation were used as an input during MATLAB training phase.The dataset was preprocessed into features and a single outcome.The features that uniquely defined the MNs' geometry include in, out, len, bz, th, and force cases.While the only expected outcome was the pain value.The features and outcome were normalized between [−1, 1] and,[0 1] respectively.Finally, the normalized dataset was split into training and validating sets to run the machine learning algorithms.Different regression models were trained including but not limited to Linear, Support Vector Machine (SVM), Ensemble, Decision Tree (DT), Gaussian Process (GP), and Neural Network (NN) regressions with holdout validation of 25%.The models were evaluated via their mean squared error (MSE) and R2 measurements through minimizing the possible error arose during the computational and running the of models.Figure 1B also describes the whole process.

COMSOL Built-In Optimizer
One feature that COMSOL had was the built-in solver for nonlinear optimization problems.This was typically performed by expressing the problem correctly within the software interface then the solver searched for the best solution.The optimization problem was formulated according to the available criteria, namely, objective function, linear and nonlinear constraints, control variables, and lower and upper bounds variables.The objective function was a mathematical expression that needed to be minimized or maximized.This could be a design goal of maximizing efficiency or minimizing a specific parameter in a system.The control variables, which were the sparking segments in the optimization problem, represented the elements that could be manipulated and tuned to optimize the objective function.Additionally, those variables had bounds or limits.While the set of constraints were linear and nonlinear depending on their structure, those were related to physical or design specification where the solution should lie and satisfy the constraints at the same time.By filling the above criteria, the optimizer could provide a powerful, flexible, and efficient approach to optimizing complex systems and processes.Clearly, the goal of the optimization problem was to minimize the pain induced by MNs that was formerly correlated in Equation (1).The control variables that were used were the geometric parameters in addition to the lower and upper limits.No constraints had been included because the MN's geometry parameters were within the specific range of pain and fabrication limits as shown in literature.Thus, the optimization problem can be simply expressed as: where f(x) is the pain equation representing the objective function to be minimized, and x is the set of control variables represented in R 5 .h(x)and g(x) are equality and inequality constraints, respectively.An inequality constraint is introduced that describes the ratio between the allowable and applied compressive strengths, namely, the factor of safety (FOS) condition to avoid failure of the material into the optimization problem.The modified optimization problem becomes for g(x): and the value 2 is the minimum ratio permitted for the FOS design criterion.Finally, the solver that was utilized was the Nelder-Mead method, also known as simplex which was a gradient-free algorithm that worked best for unconstraint optimization problem without using any derivatives.

Bayesian Optimization Using LiveLink
The ground of the third approach was based on the LiveLink between COMSOL and MATLAB.The LiveLink option worked as an intersection between both environments where geometry, simulations, and solvers in COMSOL were shared and interacted with scripts, functions, and codes of MATLAB and vice versa.This allowed to manipulate COMSOL results to be visualized and create custom numerical methods through MATLAB capabilities.
The goal of this approach was to call the COMSOL model directly from MATLAB and interactively develop scripts to control the simulations.First, the model file that contained the simulation information was ready in COMSOL and then saved as an m file extension to be read in MATLAB workspace.A single script was enough to handle and process all the codes needed to run the optimization.The Bayesian optimization command, bayesopt, was used with the goal of running the optimization problem by defining the correct optimizable variables and bounds for the model.Bayesopt used GP regression to approximate the objective function where it constructed a posterior distribution of functions that fit best in data supported and then used that posterior distribution to search into next point.Therefore, fewer functions were required to evaluate functions unlike the previous methods which involved iterative techniques for evaluation of the objective functions.The solver was able to evaluate different objective functions by combining various geometric parameters within the model to search for the optimal solution.Figure 1B also describes the whole process.

Mesh Convergence
The accuracy of computational results was important, yet it came with price of time and workload.The mesh convergency test was a well-known test that used a computational software to check and validate the solution of a simulation.This could simply be done by iteratively decreasing the element size while monitoring the parameters of interest in the problem.To assess the accuracy and reliability of the results, parameters of interest no longer varied with an increase in mesh resolution (refinement of mesh size).Thus, no significant change in results was observed and it could be concluded that the mesh size was independent of the solution.Additionally, the mesh convergence test could ensure that the errors between the experiment and the computational works were minimal.During this study, mesh convergence was employed to compare and measure the accuracy of the simulations of a 3D and 2D model.A single case was picked accordingly, in 2D and 3D, to ensure a flawless replication of the two scenarios and attain error-free correspondence to the original case.Considering the reference case as a baseline, shown in Table 1 and 3D geometry in Figure 2A.The mesh convergence test for 3D model was carried out ranging from coarse mesh to a finer mesh type.Table 2 uncovers some of the important outcomes of the 3D model with the flow rate as the parameter of interest.With an increase in the number of elements, the equations of degrees of freedoms (DOFs) to be solved were significantly impacted, resulting in a corresponding increase on the computed workload and time.Tetrahedral elements were used as a default element type due to the nature of a fluid-solid interaction problem. Figure 2B represents the mesh size with respect to the parameter of interest (flow rate).As the mesh size decreased at different points, the flow rate tended to diverge to a constant value and/or relatively small changes as mesh size was reduced.According to the displayed values, the optimal mesh size element for the 3D model lied at 7.77 μm of size.The relative change between mesh sizes of 7.77 and 2.73 μm was only 0.35% with evidently states mesh size was independent of flow rate values.
The main objective of the mesh convergency was to quantify the use of the 2D axisymmetric model instead of the 3D model due to simulation complexity and time.Similarly, mesh dependency in 2D axisymmetric model was used.Figure 2C defines the geometry in 2D axisymmetric with the domains respectively.Due to the size of the 2D model, it was expected that the number of DOFs were remarkably less than 3D model owing to the assumption of representing the geometry in a 2D symmetry only.Table 3 contains the mesh convergency test specifically for the 2D model for the same referenced case defined earlier.Undoubtedly, the computational time was cutdown to seconds throughout the convergence test with less equations to be solved.Moreover, Figure 2D explains mesh size indecency to flow rate at 2.94 μm size.With the provided information, the relative change between 2.94 and 0.3 μm was fairly small at a value of 0.59%.It was important to compare and validate the use of 2D model instead of 3D model via mesh convergence.The following Figure 2E summarizes the two models as a function of normalized mesh sizes from,[0 1] "0" being finest mesh size and "1" being coarsest mesh size.Again, flow rate was the decisive parameter of interest in this comparison, each model has its own flow rate values at different mesh size points where it was computed.Nonetheless, it was apparent that the difference in flow rate values between the two models was tiny, within the limit of 2.6% at a converged mesh.One thing to notice was that the behavior of the two model lines was the same, which can easily be explained since the same problem definition was used.In other words, provided that the problem was well-defined within the software and proper attention was given to appropriate meshing and boundary conditions in 2D model, the outcome wouldn't deviate too far from the 3D solution.Besides, the consumed time in 3D model was continually increasing as the mesh size decreased during the test, hence, time factor played an important role in the solution for larger datasets.Additionally, the whole dataset was examined for both models.Table 4 summarizes these findings for memory storage, computational time, and mesh size.Henceforth, 2D axisymmetric models will be implemented to simulate incoming models.

Ethical Approval Statement
No human or animal materials were used in this study, all work was exclusively simulation based.Therefore, no ethical approval was required.

Parametric Sweep with Regression Models
A total of 9072 outcome simulation results were evaluated by the software via the parametric sweep defined earlier.These datapoints were extracted to MATLAB to begin the regression modeling process.During the fitting modeling process, the dataset was normalized to speed up the training and to improve the accuracy of the models.Only six of the regression models were trained on the dataset.Two of the regression models were badly trained and poorly performed in accuracy, while two had somewhat acceptable performance.Yet, two were evaluated closely and performed better than all other models, namely GP and NN regression.Linear and SVM regression scored the lowest in terms MSE and R2.This can be related to linearity, presence of outliers, and high dimensionality which are properties that have made them inappropriate models.On the other hand, decision tree and ensemble of tree are related to the tree structure.The ensemble of tree uses and combines multiple predictors of a decision of tree as the base model.Thus, the results of their regression models were acceptable in terms of R2 score but couldn't fit the model due to the noise in the large dataset and outliers.Finally, the GP had the best training results considering fewer assumptions made about the data where it can model complex data more accurately for non-parametric models.Additionally, it handles the noise coming from the dataset and uncertainties during computational result calculations.Figure 3 and Table 5 provide an informative comparison among all models.

COMSOL Built-In Optimizer
The results obtained from the COMSOL's built-in solver for nonlinear optimization problems demonstrated significant findings.The Nelder-Mead method effectively identified optimal solutions for the control variables within the specified bounds, minimizing the pain induced by the MNs as outlined in the objective function.
No violation of the specified range of pain and fabrication limits  occurred, indicating that the geometric parameters of the MNs successfully adhered to the predefined bounds.This compliance with the established parameters signifies that the optimization process functioned as intended, ensuring that the solutions provided fell within practical and acceptable limits.This was utilized for all cases and demonstrated in one case (male 20-30 years old).
The solver was able to find the minimum pain value by searching for specific combination of parameters.Table 6 provides M_20-30 case with the geometric parameters that yielded the minimum pain value.Consequently, the lack of constraints did not appear to adversely affect the outcome of the optimization process, suggesting that the control variables were sufficiently constrained by their lower and upper limits.In addition, the generated solutions were both feasible and applicable, by using simplex method, aligned with the available literature on MN geometry parameters for pain reduction.It's worth noting that the application of the Nelder-Mead method underlines its effectiveness in solving unconstrained optimization problems without the use of derivatives, making it a suitable choice for this particular optimization task in-short time.
Although the Nelder-Mead proved its effectiveness for unconstrained problems, clearly, the proposed optimization problem required an algorithm that can handle the constrained part properly.The Monte Carlo, which is a stochastic algorithm, that explores and extracts statistical information from the design space was utilized as a measure of improving and handling the full optimization problem.The process of using the Monte Carlo method is of no difference to the Nelder-Mead, simply, changing the optimizer solver type while introducing the constrained condition of the problem.The results are provided in Table 6 as a contrast to the earlier algorithm.

Bayesian Optimization Using LiveLink
Utilizing the LiveLink between COMSOL and MATLAB proved to be a highly effective approach in achieving the optimal solution.By directly calling the COMSOL model from MATLAB and interactively developing scripts to control the simulations, we were able to seamlessly manipulate and visualize the COMSOL results through MATLAB's extensive capabilities.It was successfully achieved by completing and managing the optimization process through a single MATLAB script by loading the model file containing the simulation information into the MATLAB workspace.This unified approach significantly straightened the simulation process, enhancing the overall efficiency of the task.The use of the bayesopt command for Bayesian optimization was instrumental in guiding the search for optimal solutions.By defining the appropriate optimizable variables and bounds, the "bayesopt" function effectively leveraged Gaussian Process regression to approximate the objective function.This resulted in the creation of a posterior distribution of functions that best fit the data, facilitating the search for the next point to evaluate.Unlike traditional methods that rely on iterative techniques for evaluating objective functions, the Bayesian optimization approach resulted in fewer required function evaluations.Similarly, MN design for male 20-30 years old was optimized using this approach.The algorithm computed the function evaluations for 150 iterations; however, it converged only after 30 iterations toward the optimal solution, as shown in Figure 4. Thus, the minimum pain value using this approach has been recorded in Table 7 along with the geometric parameters.Therefore, by combining various geometric parameters within the model, the solver was able to evaluate different objective functions and successfully search for the optimal solution.This demonstrates the considerable potential of utilizing the LiveLink between COM-SOL and MATLAB in conjunction with Bayesian optimization for solving complex optimization problems.

Comparisons
In order to compare the performance and degree of success among all three approaches, a single case had been chosen as a baseline for a 2D axisymmetric problem of the MN.For Approach 1, which involved running thousands of simulations then training multiple regression models on the data, GP and NN regression models had superior performance due to their ability to handle complex, nonlinear relationships and noisy data.It was concluded that they were the best in terms of fitting the data and accuracy.However, despite the thoroughness of this approach, it was tedious to handle such large dataset and the presence of outliers posed difficulties for this approach and made it computationally expensive in terms of storage space and time.
Approach 2 was managed through COMSOL's built-in solver for nonlinear optimization problems.This method allowed the identification of optimal solutions within specified bounds effectively.It was done by ensuring that the solution fell within  practical and acceptable limits.This approach also had an advantage in terms of maintaining accordance to pre-defined geometric parameters.The Nelder-Mead method used here, quickly found the optimal solutions for our control variables, all within our specified bounds.This means we were able to minimize the pain induced by the MNs.Moreover, our control variables were effectively contained by their upper and lower limits.So even without any additional constraints, our solutions were still practical and applicable, aligning perfectly with the existing research on MN geometry parameters.The algorithm demonstrated its strength and speed in dealing with unconstrained optimization problems without needing to use any derivates which was ideal for this case.As a result, we found ourselves with a solution that was not only optimal but also one that was achieved in a relatively short amount of time compared to approach 1.Furthermore, the appreciation of having multiple optimizers in the COMSOL groceries helped us to model the optimization problem correctly when introducing the constrained condition using the Monte Carlo algorithm.Its layout better performance compared to Nelder-Mead and regression approach.The simulation outcomes are presented in Figure 5 for the related optimized geometric parameters such as velocity, stress, and pressure contours.
Finally, approach 3 utilized the LiveLink between COMSOL and MATLAB.This method integrated the capabilities of both software environments.With the Bayesian optimization approach, this method required fewer function evaluations and facilitated an efficient search for the optimal solution.Additionally, it demonstrated considerable potential for solving complex optimization problems and highlighted the value of combining COMSOL and MATLAB with Bayesian optimization for other works.This approach had the shortest running time compared to others.

Graphical User Interface
A practical way to calculate and interactively display the optimal MN design for pain reduction has been develop in a Graphical User Interface (GUI) in MATLAB environment.This app can handle different cases of age and gender as previously mentioned in Table S2 (Supporting Information).The user-friendly app contains two drop-down menus for gender type and age ranges to test different combination of MNs.The framework approach that has been used in the app was Bayesian Optimization Using LiveLink as it provided a comprehensive details of MN design values and achieved the lowest pain score compared to other approaches.The app, has been named "PainBayesOpti", can additionally display the time needed to evaluate the chosen case.Furthermore, the iterations for the bayesopt function have been set to 30 objective evaluations as tested its effectiveness and functionality in the previous 150 iterations (converged nearly at 30 iterations).The optimization was once more run for 30 iterations through the app while obtaining similar geometric parameters and pain score shown previously in Table 7. Figure 6 provides a preview of the app design and available information.

Limitations and Concerns
One major difficulty is to understand the mechanical elasticity of skin structure.As reported, the human skin displays a viscoelastic behavior which provides an idea of the elasticity of human's skin which can be hypothesized as a function of the thickness of its layer, [41] where the stiffness decreases with aging due to the loss of the dermal layer as a result of increase in fat infiltration into the skin. [42,43]Pain evaluation for targeted sites slightly changes in various oral cavities, [44] sensitive to pain in forehead,  8.
forearm, and less in abdomen skin with a dependence on physical parameters of MN [45] and when used at common injection location, especially during facial procedures where patients pain threshold, depth and differences in skin structures at those locations are taken into account. [35]ommon MNs have been fabricated using molding methods such X-ray lithography and UV lithography but those are time consuming methods and require costly cleanrooms, other better fabrication methods have been explored using laser ablation through CO2 laser cutter where sharp conical shape is possible to be achieved. [46]Additive manufacturing techniques such as twophoton polymerization (TPP) and 3D printing combining with UV light via Stereolithography (SLA), [47] are being developed for complex geometries, fast fabrication, nontechnical expertise and requires no expensive cleanrooms where layer by layer fabrications are processed. [48]On the other hand, fabrication technique varies with the MN and material types that being used.Typically PDMS platforms have been widely studied for short, pain-free MNs that are used in drug delivery usage where micromolding, TPP with hot embossing, and casting are most feasible to be used. [48]

Conclusion
In conclusion, this study presented three different approaches to minimize pain associated with microneedle designs with the help of machine learning and optimization tools.Each of the approaches presented strengths and weaknesses.The first approach involved the application of multiple regression models trained on data generated by finite element analysis which provided a comprehensive analysis of various regression models.Despite the challenge of handling a large dataset, the Gaussian Process and Neural Network models outperformed others due to their adept handling of complex, non-linear relationships and noisy data.The second approach showcased the power of leveraging built-in software tools and the use of nonlinear algorithms like the Nelder-Mead for solving optimization problems.While approach 2 can be improved by using other available built-in nonlinear solvers, Nelder-Mead and Monte Carlo were the choice of use.The third approach presented the benefits of integrating software environments and the efficiency of Bayesian optimization.Since each approach had its merits, the third approach offered a combination of efficiency, practicality, and timesaving, thus making it the most promising for complex optimization problems.Table 8 presents an overall results for all three approaches related to M_20-30 case.The difference could be explained that each approach was able to find a local minima.Our results suggest that the third approach, using the LiveLink between COMSOL and MATLAB coupled with Bayesian optimization, offers a promising strategy for future research due to its efficiency, practicality, and timesaving.However, the decision of which approach to use ultimately depends on the specific problem at hand and the

Figure 1 .
Figure 1.A) Workflow of three different approaches for pain optimization.Each approach exhibits unique features.B) Comparison between Approach 1 (Regression) and Approach 3 (LiveLink).First, defining the geometric definition and modeling the boundary condition of the fluid-solid interaction to define parameter sets.At regression branch evaluation of pain values through FEA.Then, datasets are created and preprocessing through normalization of the dataset.Training the machine learning algorithms for regression models.At LiveLink brach, the COMSOL mode file (.mph) is converted into MATLAB script file (.m).Next, Bayesian Optimization using bayesopt function is activated to search for best design variables.Finally, a comparison of the two approaches is evaluated based on pain index.

Figure 2 .
Figure 2. A) Geometric Definition of MN Design; Inlet, Outlet, Length, Thickness.B) 3D Mesh Convergence Test Plot.Mesh size on horizontal axis while flow rate values on vertical.C) 2D axisymmetric geometry design of MN.D) 2D Mesh Convergence Test Plot.Mesh size on horizontal axis while flow rate values on vertical.E) Comparison of 2D & 3D model mesh convergence.Normalized mesh on horizontal axis and flow rate values on vertical axis.

Figure 3 .
Figure 3. Regression models obtained from approach 1. A) Gaussian process regression, B) Neural Network regression, C) Decision tree regression, D) ensemble of trees regression E) Support vector machine regression, F) Linear regression.

Figure 4 .
Figure 4. Minimum observed objective function that has been actually evaluated during the optimization and estimated minimum objective function that has been obtained from probabilistic model versus function evaluation for the Bayesian Optimization using third approach in Livelink environment.

Figure 5 .
Figure 5. COMSOL simulation results (for approach 2) from left to right; Velocity flow achieved within the MN during the ISF extraction where the maximum velocity is 3.5 m s −1 .The Von Mises stress due to the applied force on the body of the MN caused a deformation that threatens the failure of the material, yet the maximum stress occurred is around at 0.15 MPa where the allowable compressive force is 2 MPa (FOS = 14.2).Pressure distribution contour provided that high pressure exists at the entrance of flow supported the pressure-driven flow against the body forces.

Figure 6 .
Figure 6.Design App of "PainBayesOpti" content showing the evaluation time of 434 s (≈7 min) for 30 iterations and the obtained optimal geometric parameters as consistent with Table8.

Table 2 .
Evaluation of flow rate, number of elements, degree of freedom (DOF), and computed time for 3D mesh convergence test.

Table 3 .
Evaluation of flow rate, number of elements, degree of freedom (DOF), and computed time for 2D mesh convergence test.

Table 4 .
Comparing 2D and 3D models by storage space, computational time, and mesh size type.

Table 5 .
Model type used for each regression models.Accuracy results in terms of MSE and R2 validation of the training models.

Table 6 .
COMSOL built-in optimizer results for objective function (pain) and control variables (geometric parameters).

Table 8 .
Summary for the three approaches for nine different cases.Each resulted in different local minimum for induced pain in MN.Comparing in terms of computed time and storage memory.resourcesavailability.As we continue to advance in this area, it is important to keep evaluating and improving these methodologies, thus exploring new methods, and pushing the bounds of what we can achieve in minimally invasive microneedle designs.