Design factors for determining the radula shape of Euhadra Peliomphala

Biomimetics present useful ideas for various product designs. However, most biomimetics only mimic the features of living organisms. It has not been clarified how a given shape is attained through natural selection. This paper presents the design factors that optimize the radula shape of Euhadra peliomphala. Clarifying the important design factors would help designers in solving several problems simultaneously in order to adapt to complicated and multi-functionalized design mechanisms. We measured the radula of Euhadra peliomphala by using a microscope and modeled the grinding/cutting force using the finite element analysis (FEA). We reproduced the natural selection using multi-objective genetic algorithm (MOGA). We compared the solutions when optimizing the radula shape using objective functions of each combination of stress, cutting force, abrasion, or volume. The results show that the solution obtained through two-objective optimization with stress and cutting force was the closest to the actual radula shape.

Biomimetics mimic the processes and structures that living organisms adopt to survive in the environment for various product designs [1][2][3] . Otto Schmitt proposed "biomimetics" for the first time to mimic the signal processing of a neural network and develop an electric circuit for eliminating noise and converting it into a rectangular wave called "Schmitt trigger" 4,5 . There are several studies on biomimetics. Velcro was developed to mimic burdock seeds 6 . Adhesive tape was developed to mimic the Van der Waals forces on the skin of a Gecko's feet that support the body weight of the Gecko 7 . TEIJIN Inc. developed rain wear with ultra-hydrophobicity by mimicking lotus leaves 8,9 . LIXIL developed self-cleaning building materials to mimic the biological tissue of snails 10 . Furthermore, Jin studied the nanomechanical properties of dragonfly wings (Anax parthenope) 11 . Micro air vehicles were developed to implement biomimetic design 12 . However, these studies have only mimicked the features of shapes, and have not discussed the factors behind shape designs.
In fact, in optimization in manufacturing, there are requirements and restrictions in terms of various performances, manufacturing cost, and design variables within a certain range like dimensions. An optimization method is used to obtain the most appropriate solution among the given conditions. In optimization, the relationship between the design variables and the requirement item is shown as the objective function. The solution is derived to maximize or minimize the objective function within the range of the design variables. For actual product design, the satisfaction of multi-objective function is required. For example, quality, cost, delivery (QCD) is the main issue for designers [13][14][15] . In addition, a universal design is regarded important 16,17 . Furthermore, reduction of power consumption 18 and weight 19 is also indispensable for product design. Therefore, designers require a method to solve multiple problems simultaneously to adapt to complicated and multi-functionalized design mechanism [20][21][22] . Conventionally, designers have solved problems depending on extensive experience and intuition 23,24 . However, experience and intuition vary among individuals and this will lead to different solutions and will not necessarily result in an optimal solution 25 . It is necessary to clarify a method to satisfy multiple design indicators. Wolfram was the first to suggest that the shapes of many living organisms can be given as mathematically optimal solutions 26 . However, conventionally, there are few studies on the factors determining the shape of living organisms. Berthe studied the frictional properties of snake skins to analyze the detailed relationship between biological and frictional structure 27 , but did not clarify which factors helped in determining shape.
In this paper, we aim to clarify the important indicators for an optimal design by simulating natural selection based on biomimetics. Most biomimetics have heuristically mimicked shapes 28,29 . We examine which design indicators are used from the shapes of living organisms by simulating the natural selection with the multi-objective Scientific REPoRTS | (2019) 9:749 | DOI:10.1038/s41598-018-36397-x genetic algorithm (MOGA). We simulate the dynamic behavior and shape of living organisms by using CAE and a multi-objective optimization method to examine what factors optimize the shape. In this study, we selected a land snail, Euhadra peliomphala as a model case, as shown in Fig. 1. Euhadra peliomphala can mostly be found in Kanto district and eat plants and concrete 30,31 . The tongue of Euhadra peliomphala is called the "radula" as shown in Fig. 2, and it consists of regularly aligned projections 32 . Euhadra peliomphala uses the radula to scrap hard concrete and eat 33 . The lamellar microstructure on Euhadra peliomphala's radula is consisted of the basic building blocks like the Pectinidae, Aragonite and Red Abalone [34][35][36] . The old radula at the tip is replaced with new radula generated from the radula sac in a manner similar to that of human nails 37 . In natural selection, radula would have been optimized to satisfy the multi-objective. The candidates for objective function are the following four. One is minimizing the cutting force of the radula. Next is minimizing the stress on the radula. Three is minimizing the radula volume because calcium contained in radula is difficult to ingest. Four is minimizing the amount of abrasion per unit time during grinding of the radula. We hypothesized that the radula would have been optimized by any of these four objective functions, or a combination thereof. We simulate the natural selection  by MOGA and reveal the shape of the radula so that one or combination of these four objective functions minimizes. We clarified that the objective function having solution of radula shape closest to the actual radula shape is emphasized by natural selection. The factors that optimize the shape of the radula would be useful in designing a scraping mechanism such as a lathe.

Results
Grinding force/cutting force. We observed the eating trace of Euhadra peliomphala to decide grinding and cutting forces, as listed in Tables 1-3. Similar to the determination of Cp, we compared the movement of the radula and the grinding/cutting force model, as summarized in Tables 4 and 5.   We made the response surfaces between each objective function I-IV and each radula parameters such as the rake angle, round, tooth angle, and tooth height. From these response surfaces, we clarified the error between the measured value and the theoretical value by the response surface on each Objective function, shown in Fig. 3. In Fig. 3(c,d), there is a spike so we analyzed the measurement value as the outlier based on the statistical analysis.
Optimized solution search using MOGA. We validate the objective function that derives the nearest optimized solution to the actual radula. Single objective, two-objective, three-objective and four-objective optimization solutions each comprise 4 C 1 = 4 ways, 4 C 2 = 6 ways, 4 C 3 = 4 ways and 4 C 4 = 1 way. Sum of solutions is     Tables 6 and 7.  From Tables 6 and 7, the value of each single objective function was conflicted; thus, the radula was optimized for multi-objective optimization in Appendix 3.

Comparison of analysis and estimated values.
We analyzed the radula model by using finite element analysis. In this study, the radula would perform local grinding, one of the radula would perform cutting work, and grinding processing would be realized entirely by the entire radula. Therefore, we observed the local radula motion microscopically and estimated the total radula motion macroscopically. We used finite element analysis to compare the estimated and analysis values, and the values are listed in Table 8. The estimated vertical force on the radula was 1.33 N; the vertical force obtained in the radula shape analysis was 1.06 N. Therefore, the estimated and analysis values showed good correspondence. However, the tangential force obtained in the analysis was 5.81 N, but the estimated value was 3.31 N. The tangential force obtained in the analysis has some error. On the other hand, a comparison of the estimated value and that obtained in the cone shape analysis indicated that the cone shape matches the estimated shape. In the cone shape analysis, the tangential force was 3.27 N and the estimated value was 3.31 N. With respect to the vertical force, the analyzed value was 1.46 N and the estimated value was 1.33 N. These results show that the cone shape model would be suitable for analysis microscopically.

Comparison of optimized solution and radula shape.
We investigated what design factors determined the radula shape. Table 9 shows a comparison of single objective optimization results. Table 9 lists the contradictory design variables among the objective functions. It indicates that the radula would have an optimized shape involving a few objective functions. Next, we studied the range of each objective function on all objective (four objectives) optimization. Table 10 lists the results of all objective function optimization. Table 11 shows that the cutting force was minimized by 36.3%, stress was minimized by 84.2%, abrasion was minimized by 56.6%, and volume was minimized by 7.31%. Especially, minimizing the stress by 84.2% considerably affected the shape of the radula. The results suggested that the shape of the radula would be optimized if stress and other objective functions are minimized. On the contrary, volume minimization would hardly influence the shape of the radula.
Objective of optimization on radula. The range of pareto optimization is different for each objective function. To compare each solution, we normalized the solution so that the maximum value is 1 and the minimum value is 0. Next, to evaluate the optimization of each objective function, we calculated the distance L min between the actual radula and the nearest normalized pareto solution. Figure 4 shows the distance between each   Fig. 4, the lowest distance L min is 0.329 for the stress-cutting force optimization. Thus, the solution of the stress-cutting force optimization is nearest to the performance of the actual radula. It indicates that the radula should be designed by optimizing the stress and cutting force.
In future works, we will try to the radula observation more precisely and model the radula using Hertzian-Mindlin solution. Furthermore, we will adapt the stress-cutting force objective optimization method to one with drilling functions such as the lathe. In addition, we will study not only radula, but also shell of Euhadra peliomphala. Since the shell of Euhadra peliomphala protect the soft bod from attacks like as Nacre 38-40 , we need clarify how the combination of shape and microstructure contribute to the mechanical robustness. Moreover, we would compare the mechanical properties and the objective function obtained from the radula shape of snail and shellfish with different food properties.

Methods
Modeling of grinding/cutting in Euhadra peliomphala. We observed the radula of Euhadra peliomphala. The radula was obtained using the following procedure.
(2) Boil the frozen Euhadra peliomphala. We obtained the 3D shape of the radula by using a microscope (VHX-5000, KEYANCE) to determine the height data of each pixel. We developed a 3D CAD model from these height data, as shown in Fig. 5 and Table 11. Round 1 and Round 2 in Table 11 are defined by the rho dimension of a conic section.
Round 1 and 2 are defined by the rho dimension of the conic curve. We describe rho in Fig. 6. The shape of the conic curve is defined by the rho dimension of the conical arc PQ. The straight-line segment PR is tangent to the ellipse at the point P and the straight-line segment QR is tangent to the ellipse at the point Q. The straight-line segment RD intersects the straight-line segment PQ at the point D, which is the midpoint of the straight segment PQ. The "rho" dimension is a value that specifies the ratio up to the point C to the vector up to the vertex  In Table 11, tooth angle 2α is 60°. The coefficient of friction μ is 0.2, which is similar to the coefficient of friction between wood and wood because the Euhadra peliomphala has mucus in damp condition. We carried out an experiment and calculated the grinding energy ratio C p [N/mm 2 ]. It is defined as the energy required for removal per unit volume of food, in the Appendix. We investigated a food trace to decide the cutting depth Δ [mm] and grinding diameter of the radula D [μm] and observed the motion of Euhadra peliomphala by using a microscope to determine the effective tooth number N, the feeding speed ν [mm/s], the radula moving speed V [mm/s], and grinding width b [mm], as shown in Table 12.
In this paper, we assume that a part of the radula is cut locally, but the overall radula grinds globally. Hertzian-Mindlin solution is a useful method which can know the outline of contact force separately from finite element analysis due to shape deformation at the time of contact 41 . On the other hand, obtaining the deformed shape is based on an analysis like the finite element method, and it is necessary to utilize the finite element method itself. For this reason, it is important to clarify whether Hertzian-Mindlin solution is applicable to the radula, but it will be a subject for the future. Based on the shape of radula, we develop the grinding and cutting model, as shown in Fig. 7 and Table 12. Thus, grinding force F t , F n is determined by the effective tooth number N  multiplied by the tangential cutting force per unit f t , f n . We determine the tangential grinding force F t and vertical direction grinding force F n by using equations (1) and (2).
Experiments were performed to determine the grinding energy ratio Cp, and the results are described in Appendix 1. The results of the Young's modulus and Poisson's ratio are described in Appendix 2.
Optimization condition. We optimized the shape of the radula based on the L32 orthogonal table in Table 13 after 32 trials and developed the response surface for each objective function. We searched the optimized solution from the response surface by using multi-objective genetic algorithm, MOGA.
In this research, the reason why we used the MOGA method is because this method was performed based on multiple initial solutions, and it was used from the viewpoint that it is hard to arrive at a local solution due to occurrence of a mutation. The SOBOL sequence is used to obtain the initial generation (first generation). The SOBOL sequence is an experimental design method based on the quasi-random number method (SOBOL),  which is suitable when the number of input variables is 2 to 20, and since the design is arranged more evenly than the uniform random number sequence, the multi-purpose GA is said to be suitable for generating the initial state of the first generation or MOGA. Genetic algorithm (GA) is a learning algorithm that imitates the evolution of living things, which is very broad in scope. In other words, it genetically models genetic laws that have evolved over thousands of years, hundreds of millions of years, and yields a learning method that is useful for engineering design. GA flowchart is as follows.
(1) Initialization: Generate N individuals with random chromosomes and set the initial generation population (2) Reproduction: Computation of each individual is calculated, and individuals are regenerated according to a certain rule depending on the goodness of fit. Here, some individuals with low fitness are culled out, and individuals with a high degree of adaptability will multiply as many as that number. (3) Crossover: Crossover is performed by the set crossover probability and crossover method to generate new individuals. (4) Mutation: Mutation is performed according to the set mutation probability and mutation method, and a new individual is generated. As a result, a population of a new generation is generated. (5) End judgment: If the termination condition is satisfied, the best individual obtained at that time is taken as the suboptimal solution of the problem. Otherwise return to step (2).  Table 14. Single objective optimization parameters.

Number of individuals per generation 24
Number of generations to evolve 50 Application probability of "directional crossover" 0.5 Application probability of "selection" 0.05 Application probability of "mutation" 0.1 Inversion rate of DNA by mutation 0.05

Elite strategy Effectiveness
Handling of constraint condition Adding penalty function to objective function Algorithm type MOGA-generation alternating function evolution Seed for uniform random number generation 1 Table 15. Two objectives optimization parameters.

Number of individuals per generation 24
Number of generations to evolve 50 Application probability of "directional crossover" 0.5 Application probability of "selection" 0.05 Application probability of "mutation" 0.1 Inversion rate of DNA by mutation 0.05

Elite strategy Effectiveness
Handling of constraint condition Adding penalty function to objective function Algorithm type MOGA-generation alternating function evolution Seed for uniform random number generation 1 The basic operation of such a GA is represented by a flow chart as shown in Fig. 8. The end judgment condition in GA's procedure (5) is one of the following three. One is that the maximum degree of conformity in the population exceeds the set threshold. Next is that the average fitness of the whole population exceeds the set threshold. Three is that the number of times of generation change exceeds a preset number of times. MPGA parameters are shown in Tables 14-17.
Number of generations to evolve 100 Application probability of "directional crossover" 0.5 Application probability of "selection" 0.05 Application probability of "mutation" 0.1 Inversion rate of DNA by mutation 0.05

Elite strategy Effectiveness
Handling of constraint condition Adding penalty function to objective function Algorithm type MOGA-generation alternating function evolution Seed for uniform random number generation 1 Table 17. All objectives optimization parameters.