This paper presents the development process for a simulation software with an application in designing biohybrid machines (BHMs). The software's main objective is to determine optimal values of the physical properties of materials used to manufacture BHMs by utilizing a genetic search algorithm. Consequently, a huge number of potential material configurations have been analyzed in an eficient way to provide practically employable results. The primary analysis of the simulation software pointed out that the initially developed simulator has eficiency issues within the physics engine-oriented component. To improve the eficiency of the software, a regression-based machine learning approach is used as a surrogate model to predict outcomes of voxel-based simulations in the search process.
Bio-hybrid machines (BHMs) combine living cell actuators with artificial materials in order to achieve
greater autonomy, flexibility and energy eficiency compared to standard robots [
A prerequisite for designing BHMs suitable for real-world applications is to identify and consider all
potential BHM morphologies in a desired environment. To achieve this goal, we use the Voxelyze [
Our BioMeld-Voxelyze simulator analyzes BHM’s behavior in an environment in which BHM (catheter) may collide with external objects (vascular walls). It also contains a genetic search feature that selects the optimal configuration of material properties of actuating and non-actuating voxels leading to the most eficient BHM catheter displacement, while at the same time minimizing its pressure on vascular walls.
BioMeld-Voxelyze genetic search functionality may be extremely slow since the full simulation has to be performed for each chromosome encoding values of physical properties. Consequently, one may ask if performing the full simulation can be replaced by a predictive model, and, in turn, significantly accelerate the execution of the genetic search algorithm. Therefore, in this paper we investigate whether the genetic search functionality can be replaced by eficient machine learning surrogate models based on regression algorithms.
The rest of the paper is structured as follows. The BHM modeling approach is presented in Section 2. In the following Section 3 we describe the BioMeld-Voxelyze simulator. Improvements achieved by developing an alternative solution for the time-consuming module of the simulator are explained in 4. The last section concludes the paper.
The first task in the modeling and simulation of BHMs is the design and analysis of their morphologies
(voxel-based geometries). In our solution, the Voxelyze physics engine [
Voxelyze is a highly sophisticated physics engine designed to model and simulate the mechanical behavior of soft and flexible materials. The system uses a finite element method, where a material is divided into a mesh of interconnected elements, to model the deformation and mechanical properties of the material. The elements in the mesh represent the behavior of small volumes of the material, and their interactions with each other are used to describe the overall mechanical response of the material.
The main modeling approach in Voxelyze is based on voxels and mass-spring lattice. Each voxel in the lattice is modeled with six degrees of freedom: three translational and three rotational degrees of freedom. Furthermore, each voxel contains its mass and rotational inertia to support more complex simulations with lateral shearing and rotation. Finally, each voxel is connected to 6 neighboring voxels with connection elements, where every connection contains its own transnational and rotational stifness (Figure1). Such a setup provides a very realistic deformation simulation under diferent forces and moments.
The task of the BioMeld-Voxelyze component in the whole pipeline is to consider and select parameters related to the physical properties of materials used to build BHM and its morphology (geometry). More precisely, in Voxelyze engine we apply force to the BHM (e.g., catheter) to move it through an environment that contains other objects (e.g., vascular walls). The Voxelyze takes care of various physical aspects including: forces, collision detection, gravity, material properties and temperature.
Biohybrid machines (BHMs) morphologies are created and tested in the soft-body physics engine
Voxelyze on a 3D Cartesian lattice workspace, where voxels are connected together to form a functional
machine. The simulator models the distance between adjacent voxels as Euler-Bernoulli beams [
At the lowest level, close to the base Voxelzye engine, we have the bhmsim package. The idea behind this part of the BHM simulator is to provide basic structures for building diferent simulation scenarios in terms of morphologies, hyperparameter values and simulator goals.
Morphology classes contained in the morph package should be considered as an auxiliary tool for
easier setting of the BHM simulator and simulated objects. It holds plain, low-level information about
the physical properties of simulation objects. It is possible to visualize assembled morphologies through
instances of classes inside the camera package. Classes defined within this package are designed to
enable two sights over morphologies and their motions during simulation in two and three-dimensional
space. For this purpose, The OpenGL Utility Toolkit (GLUT) [
Finally, the last module of our BHM simulator is intended to run a configurable genetic algorithm through numerous simulations to identify the best-featured combination of materials. The component name is ga and it contains class definitions for basic concepts of the genetic algorithm.
Generally speaking, genetic algorithms are heuristic algorithms applicable to a wide range of optimization problems with the goal of achieving a solution by simulating the mechanism of biological evolution. In preparation for building a genetic algorithm, the initial problem must be amenable to redesign in the GA framework. This framework expects to have a set of problem instances (valid problem solutions) as an initial population which is repeatedly modified until the desired condition is fulfilled. One instance of the problem has to encode properties whose values impact the optimal solution of the problem. In the case of BHM simulation, the goal is to identify properties of the materials that potentially would be suitable for the production of BHMs intended for use in real-world scenarios. In our BHM simulator, the goal is to identify the stifness and density of the material from which BHMs are supposed to be made.
At the highest level of abstraction, the optimization problem is solved by the following steps: • (1) generation of the initial population which produces a set of diferent BHMProperties assets, • (2) creation of Voxelyze physical construction through import files with model data, • (3) assignment of the generated properties over the physical presentation and evaluation incorporated with GA optimization process.
GA optimization is a cyclic process that involves GA operators over previously evaluated individuals. Typically used GA operators are selection, crossover, and mutation. The usage of these operators ensures the quality of the genetic evaluation process by maintaining the quality of the population in each iteration. The algorithm is performed with the evaluation of every individual genome of the population, after which the best N individuals are selected and marked for the next round of the genetic algorithm. At this point, population size is decreased and must be compensated (by creating new individuals) before the next evaluation round to match the initial population size. New individuals are created from the existing ones, through the selection of two to serve as a source of genetic material. Parents are selected randomly in the selection process over which the crossover operator is applied. Crossover creates a child individual through a combination of parents’ genetic material. For this GA operator, we consider that the child genome contains an arithmetic value of the corresponding parents’ attributes. The production of new individuals is finished after restocking the number of genomes in the population to the initially defined size.
Preparation of the new ofspring for the next algorithm generation requires mutation over some of the individuals in the population. The mutation is performed over every individual with a certain probability and it considers a probable slight adaptation of every observed property. The function looks for an old value from the genome of the current property and changes it by adding a randomly generated noise. As it is already pointed out, not every individual has to be mutated in the process of new ofspring preparation, just has to be defined as a part of evolution. The goal of the GA implementation is to achieve a certain number of iterations and the best configuration in the latest iteration is considered as a champion (optimal solution).
In terms of the quality assessment mechanism, the fitness function is defined to maximize the BHM displacement and minimize the pressure on external objects in the environment (e.g., vascular walls). In other words, the goal is to find individuals that achieve the largest possible displacement and at the same time make the lowest pressure on the vascular wall. After every individual of the population is evaluated by computing its fitness value, the population is sorted by the obtained fitness values in order to keep the top best configurations for the generation of new ofspring.
Figure2 presents the champion morphologies produced within BioMeld modeling and simulation
framework that deals with the morpohology optimizations based on compositional pattern producing
networks [
For every considered morphology, GA simulation examines the space for the selection of a winning configuration. The examined space refers to the possible range of values for monitored properties of interest crucial for understanding the goal features. One considered configuration consists of stifness and density material properties whose correlative behavior is studied through Voxelyze simulation. GA component here has the role of origination of the configurations whose analysis should be considered within the physics engine simulator. The idea was to try to have as many as possible provided combinations of values for the BHMs component building features to consider so that the defined goal is accomplished with greater certainty.
GA parameters, such as the size of the population, number of iterations, and threshold value for new generation genesis were selected by using a probabilistic method. Considering the number of possible values for the properties of BHM components, and the possible combinations, the idea was to try to achieve the highest percentage of examinations of parameter space but also to identify simulation eficiency for diferent sizes of problems selected to solve.
Table 1 provides perspective over executed experiments. Every row in the table corresponds to the specific morphology from Figure 2 and a corresponding simulation executed. One GA simulation is determined with GA parameters used within, the number of iterations , population size , value (GA threshold – how many instances are kept for the next GA iteration), and the longest and shortest iteration execution time. Iteration time duration is selected within the first 50 iterations of every simulation.
As it can be observed executed simulations had a very long execution time, and most of them were interrupted before the end. The analysis of the completed GA iterations identified that none of the considered experiments did not evolved its fitness function to the point of saturation. We will focus on the results of the experiment performed for configuration A. The scope for the simulation was a relatively small dimension, counting the 200 iterations of GA for the population size 80, where almost one-quarter of the population was retained for the next iteration. The time required for the longest iteration to be performed was more than six hours, whereas the shortest iteration needed only two hours less than that. Considering that, we can set 5 hours as the average time execution of one iteration, which makes this simulation have a total duration of approximately 42 hours. The most time-consuming component of the software was the calculation that a soft body physic machine Voxelyze had.
To answer the above-stated research question we have examined several regression models trained to
predict the outcomes of simulations participating in the cost function of the genetic search algorithms:
displacement of BHM and pressure on non-BHM objects. The examined regression models are selected
to have fast training and inference algorithms [
A B C D E F G H I J
• Elastic net regression (ELNET), • Support vector machine regression (SVM), • Random forest regression (RF), • K-nearest neighbors regression (KNN), and • Adaboost regression (ADAB).
The linear regression algorithm determines the coeficients of a linear model by minimizing the residual sum of squares between real values of the outcome variable and predictions obtained by the linear approximation. Lasso and ElasticNet find a linear model by minimizing the residual sum of squares with incorporated regularisation penalties: Lasso uses the L1 regularisation penalty, while ElasticNet uses both L1 and L2 regularisation penalties. Regression models trained over BHM Simulator data had default values for L1 and L2 regularisation penalties, defined within Scikit-learn implementation.
SVM, RF, and KNN are adaptations of the corresponding classification algorithms for regression tasks. Support vector machine classifiers are based on the idea of using linear models to identify non-linear boundaries of categories. This is achieved by transforming data instances into a new higherdimensional space using non-linear mapping. Quadratic programming algorithms are then employed in the higher-dimensional space to determine maximum margin hyperplanes separating instances from diferent categories. KNN is an instance-based (lazy) classifier: the category for a given data instance is determined from categories of the closest instances from the training dataset determined using some distance function (e.g., Euclidean or Manhattan distance). In our experiments, we have set to 20. A random forest is an ensemble of decision trees learned from bootstrapped samples of the training data. The random forest algorithm employs so-called feature bagging to determine a random subset of features for learning individual decision trees. The category for a given input instance is then the most frequent category derived from decision trees in the ensemble. In our experiments, we have used random forest ensembles containing 50 trees.
AdaBoost is a meta-learning algorithm boosting an arbitrary regression method (decision-tree-based regression in our case) towards data instances having high prediction errors. The algorithm starts by assigning equal weights to all instances in the training dataset. Those weights indicate how hard it is to predict the outcome variable for a particular instance, i.e., instances with higher prediction errors have higher weights. Some base regression method is employed to form an initial regression model. Instances in the training datasets are reweighted according to errors obtained by the trained regression model. Then, a new regression model is learned using a loss function that is pondered with weights. The two steps are repeated either for a fixed number of iterations or until the error of the model becomes acceptable. Predictions by AdaBoost are made considering all sequentially trained regression models.
As the baseline for evaluating the above-mentioned regression models, we use the so-called DUMMY regression model. The DUMMY model always predicts the same value: the mean of the outcome variable computed from the training dataset.
For the evaluation of models’ performance, we used 10-fold cross-validation. The following metrics are used to assess the quality of predictive models: • mean absolute error (MAE) • mean squared error (MSE) • coeficient of determination (R2) • Pearson’s correlation coeficient (PC): ∑︀=1 | − |
, ∑︀=1( − )2 ,
1 − ∑︀=1( − ˆ)2 , ∑︀=1( − ¯) 2 ∑︀=1( − ¯)( − ¯) √︀∑︀ =1( − ¯) 2 ∑︀=1( − ¯) 2 .
All the measures are formally defined above so that is the number of dataset instances, is the target attribute value of -th instance, is the predicted value of the target attribute of the -th instance, is the mean of the target attribute’s values, and is the mean value of all predicted values for the target attribute. Lower MAE and MSE values indicate better-performing models (for perfect regressors MAE and MSE are equal to 0), whereas higher values of R2 and PC indicate accurate predictors (the maximal value of R2 and PC indicating perfect predictors is equal to 1).
The evaluation of regression models used to predict the outcomes (pressure and displacement) of
simulations is conducted on 10 champion morphologies [
The result of 10-fold cross-validation for examined morphology configurations A, B, and C are presented in Table 2- 4. Tables show the estimate of mean absolute error (MAE), mean squared error (MSE), coeficient of determination (R2), and Pearson’s correlation coeficient (PC) for the target variables (PRESSURE and DISPLACEMENT). It is important to mention that the PC for DUMMY cannot be computed since DUMMY models always predict the same value. Resulted values are averaged using values calculated after each iteration of the 10-fold cross-validation performed.
A regression model can be considered as an accurate predictor for a target variable if the following conditions are satisfied: • MAE() ≪ MAE(DUMMY), • MSE() ≪ MSE(DUMMY), • R2() > 0.9, and • PC() > 0.9.
The symbol ≪ denotes ”much smaller”. According to the above-stated criteria, the list of accurate predictors for PRESSURE and DISPLACEMENT is given in Table 5. As it can be observed:
SVM RF KNN
SVM RF KNN
MAE
• RF is always an accurate predictive model for both PRESSURE and DISPLACEMENT. • Besides RF, PRESSURE can be accurately predicted by KNN for 8 morphology configurations and by ADAB for 4 configurations. • DISPLACEMENT can be accurately predicted by RF, KNN, and ADAB in all examined configurations. • LINEAR and SVM models have a minor accurate power (for some configurations they can accurately predict DISPLACEMENT), while ELNET and LASSO are always bad predictors of considered target variables.
The main focus of this paper is on the design of the BioMeld-Voxelyze Simulator software, which is part of a framework for designing bio-hybrid machines (BHMs). The main purpose of the simulator is to determine optimal values of parameters related to the physical properties of materials (stifness and density of actuating and non-actuating voxels) used to build BHM. To achieve this goal, a huge space of potential BHM configurations has to be examined.
Initially, the combination of a genetic algorithm and simulations based on the Voxelyze engine was used to identify optimal BHM configurations. However, performing BHM simulations required a lot of time to process and simulate the mechanical behavior of the materials assigned to BHM components within the simulation.
To overcome the time usage problem and accelerate the execution of the genetic search algorithm, we have devised a diferent approach. By employing the results of the initially proposed solutions, regression models were built to predict the outcomes (pressure and displacement) of simulations conducted on 10 champion morphologies. For the evaluation of models’ performance, we used 10-fold cross-validation and a set of metrics to assess the quality of predictive models.
Our experimental results indicate that random-forest regression models can accurately predict the outcomes of simulations. This result has an important practical consequence: the genetic search algorithm of the BioMeld Voxelyze simulator looking for optimal physical properties of BHM voxels can be significantly accelerated by using regression predictive models instead of performing simulations for each candidate chromosome encoding physical properties. This project has received funding from the European Union’s Horizon Europe research and innovation program under grant agreement No.101070328.