FACTA UNIVERSITATIS Series: Mechanical Engineering Vol. 19, No 4, 2021, pp. 751 - 765 https://doi.org/10.22190/FUME201029023F © 2021 by University of Niš, Serbia | Creative Commons License: CC BY-NC-ND Original scientific paper CALIBRATION OF MATERIAL MODELS FOR THE HUMAN CERVICAL SPINE LIGAMENT BEHAVIOUR USING A GENETIC ALGORITHM Marina Franulović1, Kristina Marković1, Ana Trajkovski2 1Faculty of Engineering, University of Rijeka, Rijeka, Croatia 2Faculty of Mechanical Engineering, University of Ljubljana, Ljubljana, Slovenia Abstract. Research of biomaterials in loading conditions has become a significant field in the material science nowadays. In order to provide better understanding of the loading effects on material structures, complex material models are usually chosen, depending on their applicability to the material under consideration. In order to provide as accurate as possible the material behavior modeling of the human cervical spine ligaments, the procedure for calibration of two material models has been evaluated. The calibration of material models was based on the genetic algorithm procedure in order to make possible optimization of material parameters identification for the chosen models. The influence of genetic algorithm operators upon the results in evaluated procedure has been tested and discussed here and the simulated behavior of the material has been compared to the experimentally recorded stress stretch relationship of the material under consideration. Since various influential factors contribute to the genetic algorithm performance in calibration of complex material models and identification of material parameters, additional possible improvements have been suggested for further research. Key Words: Material Model, Behavior Modeling, Biomaterial, Genetic Algorithm, Inverse Analysis 1. INTRODUCTION Research and modeling of the biomaterials behavior have both secured their place in a broader field of material science research. In order to accurately model the behavior of any such material, an adequate constitutive material model must be implemented, regardless of whether using the known one or the one adapted for the investigated material or, in some cases, a developed new material model. The common feature in the biomaterials behavior modeling lies in the calibration of an adequate material model Received October 29, 2020 / Accepted February 08, 2021 Corresponding author: Kristina Marković Faculty of Engineering, University of Rijeka, Vukovarska 58, 51000 Rijeka, Croatia E-mail: kritina@riteh.hr 752 M. FRANULOVIĆ, K. MARKOVIĆ, A. TRAJKOVSKI through its validation. The validation process here is based on the material behavior simulation and its comparison to the material behavior acquired from the experimental procedures where test samples are loaded, usually with constant strain rate, while stresses are recorded until the sample separation in two parts. The data sets recorded through these experiments serve as a basis for acquiring the knowledge of the tissue properties. The basic material models used for modeling of soft tissues were the ones originally developed for rubber-like materials displaying hyperelastic behavior, such as Mooney’s model from 1940 [1] and improved version of the same model (Mooney-Rivlin’s material model) from 1948 [2]. With the development of materials research as well as the other fields where similar models are usually used, the number of available models has grown [3–6]. Some of the hyperelastic material models are relatively simple, while most of them consist of complex expressions and, consequently, a large number of unknowns, which are material parameters. The identification of these models’ parameters requires the use of a numerical procedure which is expected to result in an accurate set of parameters while capturing the nonlinearity of the material behavior during loading process. One of the methods that is widely used relies on the principle of evolutionary algorithms, such as the genetic algorithm, whose usefulness in multidimensional search problems has been confirmed [7] as well as in the cases where the objective function is “discontinuous, non- differentiable, stochastic, or highly nonlinear” [8]. While the beginning of genetic algorithms applications in search of variables domain dates back to 1975 [9] (changed edition 1992) and in some form to an older date, a wider use in the field of optimization was seen in the 1990’s, with the development of computer technology and software. Some of the early works in 1900’s were general parameter optimization[10] and application to computation [11], dynamic parameter encoding [12], physical motion analysis [13], identification of cloth animation models [14] and overview of possible uses in optimization and statistics [15], identifying parameters for dynamic simulation [16]. More recent use includes topology optimization [17], solving IVBV problems [18] and optimization of perforated elastic plate properties [19]. In addition to the genetic algorithm, many other evolutionary methods have found their applications in medicine and healthcare. This has provided quicker and easier diagnostics and improved the healthcare systems [20]. Examples for the use of the genetic algorithms in tissue mechanics and/or parameter determination include: breast tissues [21, 22], anterior cruciate ligament [23], mass-spring models for soft tissues [24], QLV theory [7, 23], ventricular myocardium [25], liver [26], porcine coronary media and adventitia [27], sheep tendons [28] and human placenta [29]. In this paper the influence of the genetic algorithm operators on the results inassesed procedure have been tested and discussed. Also, the simulated behavior of material has been compared to the experimentally recorded stress stretch relationship of the material under consideration, which is the material behavior modeling of the human cervical spine ligaments. 2. MECHANICAL PRINCIPLES In order to model and simulate the behavior of soft tissues, a material model suitable for modeling its hyperelastic behaviour must be implemented. Experimental observations of biological material, in presented case, the ligaments of the human cervical spine were conducted. It was found that the maximum load capacity of the Calibration of Material Models for Human Cervical Spine Ligament Behaviour Using a Genetic... 753 material is in the direction of their fibers [30]. Ligaments are composed of collagen fibers organized uniformly in one direction and integrated together by intercellular material. Consequently, resulting mechanical behavior of ligaments depends on the behavior of individual ligaments' constituents as well as on their correlation and relative position. Mathematical modeling of biological materials, especially with the presumption of modeling the said ligament material, is based upon the knowledge of the basic mechanical properties of materials. It is also necessary to know the type of interdependence of the stresses and strains. These characteristics can be determined by the experimental procedures that include subjecting samples to tensile loading until their fracture. Cca. eighty experiments were performed and Fig. 1 shows the load – elongation curve of the anterior longitudinal ligament (sampled at the cervical spine, level C4–C5, from the human donor, age 73) which was loaded in tension until fracture. Fig. 1 Experimentally determined load – elongation curve of the anterior longitudinal ligament Governing equations for calculating stress in hyperelastic material models are based on the strain energy density function, which contains the relationship between stress and strain through the materials life in loading conditions. Strain energy density function can be defined as a linking function between the strain energy density and deformation gradient. As stated in [2] and [31], it is possible to express strain energy density function using three strain invariants (I1, I2 and I3). It should be noted that not all invariants are used in every material model. For example, the simplest version of the Neo-Hookean material model uses only the first invariant and the Mooney-Rivlin material model uses only the second invariant [32], while the Ogden model uses stretches instead of invariants [33]. General expression of a strain energy density function for an isotropic hyperelastic material model is [34, 35]: 1 2 3 1 2 3 ( , , ) ( , , )W W I I I W   = = (1) where W represents strain energy density function, I1, I2 and I3 represent invariants of the Cauchy-Green deformation tensors and λ1, λ2 and λ3 represent principal stretches. Invariants are related to the principal stretches via the following expressions: 754 M. FRANULOVIĆ, K. MARKOVIĆ, A. TRAJKOVSKI 2 2 2 1 1 2 3 2 2 2 2 2 2 2 1 2 2 3 3 1 2 2 2 3 1 2 3 I I I             = + + = + + = (2) The relationship between stress and stretch is derived from the aforementioned strain density function. For this particular case, a complex material model was implemented in the material behavior analysis. The complexity here refers to more than one expression for a relationship between stress and stretch, based on a current value of stretch. In the presented analysis, only the first two parts of the curve (toe region and linear region) were modeled. Two complex material models were implemented, the first one (3) from [36] and the second one (4) which slightly differs. The stress-strain relationships in each of them are as follows: 4 ( 1) 3 5 6 [ 0 1 1] c c e c c         − = −      + (3) 𝜎 = { 0 𝐶3[𝑒 𝐶4(𝜆−1) − 1] 𝐶5𝜆 + 𝐶6 𝜆 < 1 1 < 𝜆 < 𝜆∗ 𝜆 ≥ 𝜆∗ 4 1)2 1 3 2 1 ( 2 5 6 1 2 1 1 2 0 [ ] 1 c c c e c c c                −       = − + −         − + +             (4) 𝜎 = { 0 𝐶1 (𝜆 2 − 1 𝜆 )+𝐶3[𝑒 𝐶4(𝜆−1) − 1] 𝐶1 (𝜆 2 − 1 𝜆 ) + 𝐶5𝜆 + 𝐶6 𝜆 < 1 1 < 𝜆 < 𝜆∗ 𝜆 ≥ 𝜆∗ Values c1 to c6 represent unknown variables - material parameters which have to be identified by the numerical procedures in order to calibrate material model, σ is stress value, λ represents the current stretch value and λ* represents the stretch value on the transition point between toe and linear regions (marked as A in Fig. 1). Calibration of Material Models for Human Cervical Spine Ligament Behaviour Using a Genetic... 755 3. PARAMETER DETERMINATION The material model calibration consists of determining the material parameters values, which in the cases of the models consisting of a larger number of parameters means determination of the parameter sets as ensemble. The number of parameters in the parameter set also correlates to the number of the phenomena under consideration that occur within the material structure and, therefore, higher complexity in models often offers requirements fulfillment. As the number and meaning of material parameters differ from one model to another, a short overview is given in this section. From Eqs. (3) and (4), it can be seen that the number of material parameters differs depending on the used material model. Parameters c1-c5 represent material parameters, while value c6 actually is not a common material parameter but a value that can be expressed in correlation to the other parameters. The expression for value c6 in the case of material model from (3) is the following [36]: ( )4 * 1 6 3 5 1 c c c e c   −  = − −   (5) On the other hand, when model (4) is used, the expression is: ( )4 * 12 6 3 5 1 c c c e c    −  = − −   (6) There is also another value that will be treated as a parameter in the genetic algorithm procedure, but which is not a material parameter in a more narrow sense. That value is λ*, which, as aforementioned, represents stretch value in transition point between toe region and linear region, as can be seen in Fig. 1. 4. INVERSE PROBLEM Inverse problem can be described as an assumption of an unknown parameter set and a reconstruction of the object’s mechanical properties using the known data [13]. In it, the response is observed and the system parameters are calibrated to reach that response [18]. For the identification of the values of the parameters, inverse analysis was used. Inverse analysis itself consists of three parts: system characterization, forward modeling and inverse or backward modeling. The first part consists of the definition of model parameters and was done in the previous section for both of the material models. Forward modeling deals with foreseeing the system behavior in correlation with the test results and is based on the mechanical principles of the tested material. The third part is founded upon the test results, which are used to calibrate the material model parameters as accurately as possible. 5. MATERIAL PARAMETERS In order to carry out the “inverse modeling” part of the inverse analysis, an appropriate objective function should be set. The basis for the definition of the objective function is the following expression [37]: ( ; ) i a=   (7) 756 M. FRANULOVIĆ, K. MARKOVIĆ, A. TRAJKOVSKI where  is the mapping function that connects experimental and theoretical stress- stretch relationship through the loading process, in the presented case for the first two characteristic regions of the materials life, with the material parameters, where ai represents the material parameters set. Parameters set ai for models (3) and (4), respectively, is: ai = [c3,c4,c5, *] (8) ai = [c1,c3,c4,c5, *] (9) 6. OBJECTIVE FUNCTION The choice of objective function in inverse modeling influences both the accuracy of the results and the velocity for acquiring them. Since evolutionary method, genetic algorithm is chosen to combine with the inverse modeling in this particular problem, even simple objective functions are expected to be efficient and to reach global minima in any given material behavior. Two different objective functions, which represent fitness functions in the genetic algorithm procedure, were tested here. The first one is a simple least squares method [13]: 2 1 [ ( ; )] n j j i j f a = = −    (10) where tilde (···) refers to the experimental data and n is number of recorded stress – stretch points in data set. The second function is an upgraded least squares method with additional influence of weighting factor, and is based on the one proposed in [37, 38]: 2 1 ( ; )n j j i j j a f =  − =           (11) 7. OPERATORS IN GENETIC ALGORITHM In order to obtain the optimal parameter sets for previously described material models, a numerical procedure should be applied. In this case, an evolutionary genetic algorithm is chosen because of its advantageous characteristics, such as applicability in nonlinear systems with a large number of unknowns, robustness considering the choice of objective function, ability to take into account a large number of data, ability to spread the possible results domain in order to avoid local minima and fast convergence to the results. Genetic algorithm can be defined as a heuristic search method with the main purpose of optimization of the described system [39], based on the principles of natural evolution. Thus the genetic operators in numerical procedure are inspired by the natural evolution principles, such as selection, reproduction, mutation and crossover [11, 19, 23, 39]. Together with the definition of its main properties, such as population initialization, definition of results domains, convergence rule and generations progress rule it is a powerful tool for the calibration of material models of nonlinear systems, which are based on experimental data gained from the materials response to the probing signal. Calibration of Material Models for Human Cervical Spine Ligament Behaviour Using a Genetic... 757 Search space ν and a function g [39] are defined on domain R: :g R (12) The problem is described by the following expression [21, 39]: i iarg min or arg maxa g a g      = = (13) In this particular case, the goal is set in acquiring the minimal value of objective functions given in Eqs. (10) and (11). Search space (or population domain) as well as the definition of genetic operators is also determined in an endeavour to build an appropriate numerical model for the given problem. 7.1. Population Initialization The first step in the definition of genetic algorithm should be the creation of the initial population. It is suggested that one or more individuals should be selected heuristically, instead of doing it all stochastically [23]. The domains for each of the parameters defined through the mechanical principles are chosen as follows: 1 3 4 5 0;100 ; 0; 2 ; 0;8 ; 0;100 ; 2;3c c c c        (14) It should be noted that these are the initial domains for the initial population. The values of the final solution can be outside of these domains because of the influence of the genetic operators. Another variable is the size of the population (n). Larger population will have more solutions included, but would slow down the process of computing, so the user should find the optimal value. Madjidi et al. used both n=100 and n=200 individuals, with both populations producing good results [23]. In [7], only the population of n=100 was used. Martínez-Martínez et al. used populations of n=50, n=100, n=150 and n=200 [26]. 7.2. Ranking, Selection and Elitism In order to proceed with the selection process, individuals are ranked based on a fitness value. Those with the better rank should have better chances of further selection. The process of selection in genetic algorithm can be done by various schemes [40]. K- tournament selection is here seen as an extension of the “normal” tournament selection from two to K contestants (in other words, “normal” tournament selection can be seen as 2-tournament selection) [41]. This selection randomly picks K number of individuals, all of whom are chosen with the same probability, and the best one is selected for the mating pool (the next generation). The tournaments last until there are as many individuals in the next generation as there were in the previous one. For the presented case of genetic algorithm in material parameter identification problem, K value of four (4) was chosen, based on the previous research in similar problems [37]. Although not a genetic operator in a classical meaning of the word, there is also concept of elitism (also called a population overlap). The concept secures that the best individuals survive into the next generation. There is one variable, elite population ratio that is set by the user [39]. The number of individuals that will be replaced is: 758 M. FRANULOVIĆ, K. MARKOVIĆ, A. TRAJKOVSKI R N N E N= −  (15) where NR is the number of individuals to be replaced, N is the total number of individuals (population size) and E is the elite population ratio. 7.3. Crossover There are also various crossover mechanisms. During the implementation of genetic algorithm for the given problem, intermediate crossover has been used because it is characterized by the children created with the influence of the parents’ weighted average. The application of this kind of crossover is widely used in such cases, but with various crossover ratios. Louchet et al. suggested a crossover ratio of 0.3, which is a very low value considering other sources [14]. In this case, several values were tested, but all of them closer to the value of 0.8, as suggested in [10] and [24] or 0.9, suggested in [19]. Only Martínez-Martínez et al. tested many values (from 0 to 1, with steps of 0.1 and finally set on 0.9), with the best one depending on the each case observed and varying greatly [26]. The new generation’s values are calculated here as follows [37]: 1 2 1 0.9 _ _ (1.1 _ 0.9 _ )Child Parent K Ratio K Parent K Parent K=  +   − (16) where Child represents the new generation values, Parent_K1 and Parent_K2 are the first and second parent’s values. Ratio_K is crossover ratio. One of the possible scenarios is that both parents appear as equal, considering they are randomly selected from the same pool. Since this scenario reduces the variety of genetic material, mutation can be additionally implemented on these individuals and then combined through the crossover procedure. In this particular case, the mutation ratio is set at 25% (0.25). Mutation procedure here is the following, based on [37], but adapted for this particular case while the domains are given in expression (14): 1 1 1 _ _ _ _ _ 0.25 _ 0; ( _ ) Parent K Parent K Ratio M Change M Ratio M Domain Change M Random length Parent K = +  =  = (17) where Ratio_M is mutation ratio and Change_M is the size of change. 7.4. Mutation The last genetic operator to be defined is a mutation. Its importance is significant in the cases where the convergence to unwanted local minima is probable [7]. The most influential variable in the mutation process is mutation ratio. Madjidi et al. suggested a mutation ratio of 0.01 [23], Louchet et al. suggested 0.2 [14] and Wright suggested between 0.005 and 0.3 [10], depending on the algorithm. Bianchi and Solenthaler also suggested adaptive mutation, but with a different mutation ratio function [24]. In order to achieve faster convergence to global minima, the adaptive mutation is applied to the given problem [37]. It means that the mutation possibility becomes lower with each new generation: Calibration of Material Models for Human Cervical Spine Ligament Behaviour Using a Genetic... 759 _ _ _ _ _ _ 1 _ _ 0; ( _ ) Change M Parent M Ratio M Change M Current generation Ratio M Last generation Change M Random length Parent M = +  = − = (18) where Parent_M is the parent whose values would be mutated while Current_generation and Last_generation are the ordinal number of current generation and the last generation (or the number of total generations). 8. STOPPING CRITERIA AND ADAPTATION OF GENETIC ALGORITHM As the computing environment for the development of genetic algorithm procedure for the material parameters identification is MATLAB (version R2017b), default stopping criteria has been used here, including number of generations, time limit, fitness limit, stall generations and stall time limit. In practice, most of the times, the stopping criteria is the number of generations, as the genetic algorithm is mostly effective and able to achieve relatively good results in an early stage of reproduction process. Slightly less often, fitness limit was reached just before the last generations, which results in the best solution. Sometimes, though, the run was terminated by “stall generations” criteria, i.e. the results would not converge and the final result is mostly not good in that case, which can be overcome by the fine tuning of genetic operators and genetic algorithm properties described above for the particular case. Some procedures developed in the published researches use only number of iterations (generations) and threshold value [21, 38]. There are more complex stopping criteria developed, such as those described in [42] or [43], but those were not implemented in this case and present a possible subject for further research. Some of the previously mentioned elements of genetic algorithm have variables (crossover ratio, elite population ratio, population size, domains of the parameters and number of generations) that can be set by the user based on previous experiences. During the adaptation of this genetic algorithm, various combinations of those variables were used and fine tuning of their values was deployed in order to acquire an optimal procedure which will result in convergence to the accurate outcomes. As can be expected, some of the developed genetic operators and genetic algorithm properties influence the results more than the others. The only really significant improvement here was accomplished with the increase of the initial population size and even then only to some upper threshold. Increase above that threshold did not increase the quality of the final results by much. Elite ratio also showed slight influence in the way that ratios of 0.1 and 0.15 mostly showed better results than the ratio of 0.05. Representative results of fitness function values are given in Tables 1 and 2, where Table 1 represents the increase in population size and Table 2 the influence of the elite ratio. Other operators and variables showed even less influence than the elite ratio and are not shown in tables. However, higher values of crossover ratios (0.7-0.9) showed more consistent results than the lower ones. Number of generations did not have big influence, as the problems solved here are relatively simple in terms of optimization. 760 M. FRANULOVIĆ, K. MARKOVIĆ, A. TRAJKOVSKI Table 1 Comparison of the results obtained using different population sizes for material model (3) and objective function (10) Population 100 300 500 1000 1500 2500 Results (Fitness values) 1355.88 1516.84 1405.21 1366.44 1358.36 1319.39 1583.09 1611.93 1430.92 1353.01 1319.80 1323.84 2612.97 1368.23 1322.01 1326.68 1346.62 1329.84 2886.96 1349.44 1328.99 1333.65 1353.89 1318.32 1452.07 1531.86 1422.72 1358.97 1348.87 1334.68 2875.70 1443.55 1467.38 1337.50 1443.37 1319.53 3324.66 1404.70 1341.01 1341.34 1325.27 1323.69 1623.73 2595.61 1324.94 1764.23 1462.29 1322.86 2324.78 1333.25 1322.74 1384.54 1350.09 1320.29 1400.56 1383.87 1350.86 1349.65 1318.26 1319.08 Average 2144.04 1553.93 1371.68 1391.60 1362.68 1323.15 Median 1974.26 1424.13 1345.94 1351.33 1349.48 1321.58 The best result 1355.88 1333.25 1322.01 1326.68 1318.26 1318.32 The worst result 3324.66 2595.61 1467.38 1764.23 1462.29 1334.68 Table 2 Comparison of the results obtained using different population sizes for material model (4) and objective function (11) Population 500 1000 Elite ratio 0.05 0.1 0.15 0.05 0.1 0.15 Results (Fitness values) 1385.29 1326.55 1405.21 1346.61 1374.62 1366.44 1696.33 1352.92 1430.92 1319.72 1352.06 1353.01 1440.30 1403.69 1322.01 1730.98 1322.93 1326.68 1397.48 1997.37 1328.99 1336.99 1331.02 1333.65 1357.30 1368.02 1422.72 1335.89 1318.18 1358.97 1651.64 1360.86 1467.38 1370.84 1335.30 1337.50 1327.00 1333.33 1341.01 1382.76 1336.31 1341.37 1359.75 1541.09 1324.94 1352.01 1344.73 1764.23 1675.67 1324.30 1322.74 1364.64 1388.52 1384.54 1630.75 1894.56 1350.86 1523.31 1349.95 1349.65 Average 1492.15 1490.27 1371.68 1406.38 1345.36 1391.60 Median 1418.89 1364.44 1345.94 1358.33 1340.52 1351.33 The best result 1327.00 1324.30 1322.01 1319.72 1318.18 1326.68 The worst result 1696.33 1997.37 1467.38 1730.98 1388.52 1764.23 Population 1500 2500 Elite ratio 0.05 0.1 0.15 0.05 0.1 0.15 Results (Fitness Values) 1327.36 1339.31 1358.36 1358.25 1332.85 1319.39 1320.30 1320.11 1319.80 1318.42 1323.44 1323.84 1399.55 1381.25 1346.62 1331.24 1319.46 1329.84 1320.11 1328.91 1353.89 1346.60 1325.42 1318.32 1475.95 1326.52 1348.87 1320.28 1355.59 1334.68 1321.64 1322.58 1443.37 1328.33 1366.07 1319.53 1424.25 1325.45 1325.27 1320.50 1321.98 1323.69 1366.10 1345.80 1462.29 1333.57 1333.51 1322.86 1370.43 1411.75 1350.09 1358.28 1321.64 1320.29 1330.50 1327.87 1318.26 1362.46 1319.14 1319.08 Average 1365.62 1342.96 1362.68 1337.79 1331.91 1323.15 Median 1348.30 1328.39 1349.48 1332.41 1324.43 1321.58 The best result 1320.11 1320.11 1318.26 1318.42 1319.14 1318.32 The worst result 1475.95 1411.75 1462.29 1362.46 1366.07 1334.68 Calibration of Material Models for Human Cervical Spine Ligament Behaviour Using a Genetic... 761 9. DISCUSSION AND PARAMETERS IDENTIFICATION After the set-up of the numerical procedure is done, and the adaptation of genetic algorithm finished, the comparison of results for the two previously described material models and two given objective functions follows. It should be noted that it is impossible to numerically compare the results obtained using different objective functions and graphical comparison is used instead. Two material models for the same objective function are easily compared using both numerical value and graphical representation. In order to provide comparison of the chosen material models and selected objective functions, to validate the model for the observed material behavior and thus provide the meaningful discussion on the matter, the experimental data set, as mentioned in Chapter 2, has been used. As stated before, for this calibration, only „toe“ and „linear“ regions were taken into account (stretch values from 1 to 4 as shown in Fig. 1). The first comparison is the one using the objective function from expression (10). Graphical representation is shown in Fig. 2a. It can be seen that both material models offer a good representation of experimental values for low stretch values and higher stretch values. However, it is also visible that the material model (3) significantly outperforms material model (4) for stretch values between cca 2 and 3. The difference in the two material models in this section is also clear when observing the numerical value of fitness. The best fitness value obtained using the material model (3) is 1318.18, while it was impossible to achieve the value better than 3900.21 with the one from (4). In other words, the first material model has almost three times better fitness value than the second, as the smaller value is better. A very similar conclusion can be made by observing the graphical comparison (Fig. 2b) between the two models using the objective function (11). In this case, the difference is even more pronounced than in the previous one, as the material model (4) shows a large deviation from the experimental results even for the low stretch values. This time the ratio in fitness values is even for four times in favor of the material model (3). The best value achieved using it was 15.4658 and the best value achieved using the other model was 60.2334. Fig. 2(a) Stress-stretch chart for objective function (10); (b) and objective function (11) 762 M. FRANULOVIĆ, K. MARKOVIĆ, A. TRAJKOVSKI Finally, the comparison between the two objective functions can be made. As stated before, it can be done only by visually comparing the results obtained using both objective functions (Fig. 2a and 2b). It was already stated that the results for objective function (11) were less accurate than those for the function (10). When comparing the two functions for material model (3), which was proved to be better for both functions, it can be seen that the results are pretty close to the experimental values for low and mid-range stretches. However, function (10) significantly outmatches the other one for higher values of stretches. The main reason for this is the fact that the function (11) divides the difference of squares with the experimental value, thus favoring the lower stress (or stretch) values. The best results achieved for parameters are given in the Table 3. Table 3 Genetic algorithm and material model parameters for the best case Genetic algorithm data Material parameters Population 1000 c3 1.1592 Elite ratio 0.1 c4 1.6962 Crossover ratio (χ) 0.8 c5 46.5241 Material model (3) c6 -107.034 Objective function (10) λ* 2.8647 Fitness value 1318.18 10. FURTHER RESEARCH AND IMPROVEMENTS While the genetic algorithm was proven to be an effective tool for optimizing identification of material parameters in highly nonlinear systems, there are various improvements suggested by other authors that could be a subject for further research and implementation and which can additionally improve process of material model calibration. The follow-up in further research could be according to Harb et al. who suggested the use of the combination of a genetic algorithm and exact analytical optimization in order to make it shorter and more effective [38]. In [24], it is predicted that further evaluation of adaptive mutation is a key to faster convergence. Some works suggested other methods could be used instead of genetic algorithms. Han et al. mentioned simulated annealing as a possible alternative [22]. In [26], three methods were compared: an iterative local optimization, scatter search and genetic algorithm. Genetic algorithm showed the best result for 50% of cases, scatter search for 30% and iterative local optimization for 20%. In [13], simulated annealing and genetic algorithm were utilized and the genetic algorithm was proven to be more versatile for that use. The Marquardt- Levenberg nonlinear least squares method was used by Pandit et al. and the results they have obtained show that there is no significant difference between that method and the genetic algorithm [27]. In this paper, adaptive mutation rate was used. However, certain works, [8] and [23], also suggest the use of adaptive crossover rate. According to Eiben et al, there are some cases in which multi-parent recombination should be considered, although they admitted it was a subject for further investigation [44]. Almost 20 years later, Akbari and Ziarati proposed MLEO (Multilevel Evolutionary Optimization algorithm) that features recombination on individual and group level and, although it showed positive results, they also admit future research is needed [45]. Also, Calibration of Material Models for Human Cervical Spine Ligament Behaviour Using a Genetic... 763 there are many different selection schemes mentioned in [40]. It is reasonable to consider that some of them might outperform current 4-tournament selection. 11. CONCLUSION Evolutionary procedure, in this research genetic algorithm, has been used for material parameter identification in modeling material behavior of soft tissues. In order to provide inverse modeling and thus apply presented method for the material models calibration and consequently validation of the chosen parameters set, detailed theoretical basis has been given in the paper. According to the numerous research studies worldwide, which follow a global trend in similar problems, the application of genetic algorithms or similar evolutionary procedures is proven to be generally accepted because of their advantageous properties in acquiring fast and reliable solutions in a very short time. Therefore, it has been applied to calibrate two well-known non-linear material models used for modeling of soft tissues behavior. Since genetic algorithm efficacy depends on its operators and also on its properties, their detailed representation has been given here in order to describe the fine-tuning for the given problem and to make possible its repeatability for possible additional material models or for the simulation of differently structured soft tissues behavior. For both material models tested, the developed algorithm produced good results with a relatively small calculation time while running on an average one-year-old personal computer. Additionally, while comparing the results obtained using both material models, they showed that the application of the proposed method can improve the decision-making in the process of selection of the material model which should be applied in material behavior modeling. As far as the algorithm itself is concerned, the genetic operators and properties were developed and initially tested. The biggest influence of them on the quality of the results has been shown to be a population size, the property which has been chosen for the detailed testing and fine-tuning of algorithm efficacy. The analysis showed that only the increase up to some population size (around 1000) showed significant influence, which turned out as a threshold for the quality increase. Although other genetic operators and properties have been taken into account here, their detailed influence has not been investigated, but rather built based on the previous experience. Their influence on the results are also expected, but not in the same intensity as the influence on the population size shown here. Some of the possible further improvements were also presented, some of them based on the detailed analyses of the influence of the remaining genetic operators and properties, but all of them dealing with the improvement on the genetic algorithm. They take into account operators (use of adaptive crossover, multi-parent recombination and use of selection other than the K-tournament selection). All considered, most of the improvements suggested in the literature would be further researched and, if at all possible, tested and implemented in future genetic algorithms. Also, more complicated material models would be tested once the effectiveness of the algorithm has been proven. Acknowledgements: This work has been supported by the Croatian Science Foundation under the project number IP-2019-04-3607 and by the University of Rijeka under projects number uniri- tehnic-18-34 and uniri-pr-tehnic-19-21. 764 M. FRANULOVIĆ, K. MARKOVIĆ, A. TRAJKOVSKI REFERENCES 1. Mooney, M., 1940, A theory of large elastic deformation, Journal of Applied Physics, 11(9), pp. 582-592. 2. Rivlin, R.S., 1948, Large elastic deformations of isotropic materials. IV. further developments of the general theory, Philosophical Transactions of the Royal Society and Engineering Sciences, 241(835), pp. 379-397. 3. Chagnon, G., Rebouah, M., Favier, D., 2015, Hyperelastic energy densities for soft biological tissues: a review, Journal of Elasticity, 120, pp. 129-160. 4. Marckmann, G., Verron, E. 2006, Comparison of hyperelastic models for rubber-like materials, Rubber Chemistry and Technology, American Chemical Society, 79(5), pp. 835-858. 5. Ali, A., Hosseini,M., Sahari, B.B., 2010, A review of constitutive models for rubber-like materials, American Journal of Engineering and Applied Sciences, 3(1), pp. 232–239. 6. Seibert, D.J., Schöche, N., 2000, Direct comparison of some recent rubber elasticity models, Rubber Chemistry and Technology, 73(2), pp. 366–384. 7. Kohandel, M., Sivaloganathan,S., Tenti, G., 2008, Estimation of the quasi-linear viscoelastic parameters using a genetic algorithm, Mathematical and Computer Modelling, 47(3-4), pp. 266-270. 8. Sinha, N.K., Gupta, M.M., Zadeh, L.A., 1999, Soft computing and intelligent systems : theory and applications, Academic Press, 639 p. 9. Holland, J.H., 1992, Adaptation in Natural and Artificial Systems, The MIT press, Cambridge MA, 232 p. 10. Wright, A.H., 1991, Genetic algorithms for real parameter optimization, Foundations of Genetic Algorithms, 1, pp. 205-218. 11. Forrest, S., 1993, Genetic Algorithms: Principles of Natural Selection Applied to Computation, Science, 261(5123), pp. 872-878. 12. Schraudolph, N.N., Belew, R.K., 1992, Dynamic Parameter Encoding for Genetic Algorithms, Machine Learning, 9(1), pp. 9–21. 13. Louchet, J., 1994, An evolutionary algorithm for physical motion analysis, Proc. Conference on British Machine Vision, BMVC 94, York, 2, pp. 701-710. 14. Louchet, J., Provot, X., Crochemore, D., 1995, Evolutionary identification of cloth animation models, In Computer Animation and Simulation ’95: Proc. Eurographics Workshop in Maastricht, Springer, Vienna, pp. 44-54. 15. Chatterjee, S., Laudato, M., Lynch, L.A., 1996, Genetic algorithms and their statistical applications: an introduction, Computational Statistics and Data Analysis, 22(6), pp. 633-654. 16. Joukhadar, A., Garat, F., Laugier, C., 1997, Parameter identification for dynamic simulation, Proc. IEEE International Conference on Robotics and Automation, CIRA 1997, Monterey, 3, pp: 1928-1933. 17. Sid, B., Domaszewski, M., Peyraut, F., 2005, Topology optimization using an adaptive genetic algorithm and a new geometric representation, WIT Transactions on The Built Environment, 80, pp. 127-135. 18. Karr, C.L., Yakushin, I., Nicolosi, K., 2000, Solving inverse initial-value, boundary-value problems via genetic algorithm, Engineering Applications of Artificial Intelligence, 13(6), pp. 625-633. 19. Vigdergauz, S., 2001, The effective properties of a perforated elastic plate: Numerical optimization by genetic algorithm, International Journal of Solids and Structures, 38, pp. 8593-8616. 20. Albu, A., Precup, R.E., Teban, T.A., 2019, Results and challenges of artificial neural networks used for decision-making and control in medical applications, Facta Universitatis-Series Mechanical Engineering, 17(3), pp. 285-308. 21. Lago, M.A., Ruperez, M.J., Martinez-Martinez, F., Monserrat, C., 2014, Genetic algorithms for estimating the biomechanical behavior of breast tissues, Proc. International Conference on Biomedical and Health Informatics, BHI 2014, Valencia, pp. 760-763. 22. Han, L., Hipwell, J.H., Tanner, C., Taylor, Z., Mertzanidou, T., Cardoso, J., Ourselin, S., Hawkes, D.J., 2012, Development of patient-specific biomechanical models for predicting large breast deformation, Physics in Medicine and Biology, 57(2), pp. 455-472. 23. Madjidi, Y., Shirinzadeh, B., Banirazi, R., Tian, Y., Smith, J., Zhong, Y., 2009, An improved approach to estimate soft tissue parameters using genetic algorithm for minimally invasive measurement, 2nd Proc. International Conference on Biomedical Engineering and Informatics, BMEI 2009, Tianjin, pp. 1-6. 24. Bianchi, G., Solenthaler, B., Székely, G., Harders, M., 2004, Simultaneous topology and stiffness identification for mass-spring models based on fem reference deformations, Medical Image Computing and Computer-Assisted Intervention, 3217, pp. 293-301. 25. Nair, A.U., Taggart, D.G., Vetter, F.J., 2007, Optimizing cardiac material parameters with a genetic algorithm, Journal of Biomechanics, 40(7), pp. 1646-1650. 26. Martínez-Martínez, F., Rupérez, M.J., Martín-Guerrero, J.D., Monserrat, C., Lago, M.A., Pareja, E., Brugger, S., López-Andújar, R., 2013, Estimation of the elastic parameters of human liver biomechanical Calibration of Material Models for Human Cervical Spine Ligament Behaviour Using a Genetic... 765 models by means of medical images and evolutionary computation, Computer Methods and Programs in Biomedicine, 111(3), pp. 537-549. 27. Pandit, A., Lu, X., Wang, C., Kassab, G.S., 2005, Biaxial elastic material properties of porcine coronary media and adventitia, American Journal of Physiology, Heart and Circulatory Physiology, 288, pp. H2581–H2587. 28. Sverdlik, A., Lanir, Y., 2002, Time-dependent mechanical behavior of sheep digital tendons, including the effects of preconditioning, Journal of Biomechanical Engineering, 124(1), pp. 78-84. 29. Hu, J., Klinich, K.D., Miller, C.S., Nazmi, G., Pearlman, M.D., Schneider, L.W., Rupp, J.D., 2009, Quantifying dynamic mechanical properties of human placenta tissue using optimization techniques with specimen-specific finite-element models, Journal of Biomechanics, 42(15), pp. 2528-2534. 30. Schendel, M., Wood, K., Buttermann, G., Lewis, J., Ogilvie, J., 1993, Experimental measurement of ligament force, facet force, and segment motion in the human lumbar spine, Journal of Biomechanics, 26(4-5), pp. 427-438. 31. Haines, D.W., Wilson, W.D., 1979, Strain-energy density function for rubberlike materials, Journal of the Mechanics and Physics of Solids, 27(4), pp. 345-360. 32. De Pascalis, R., 2010, The Semi-Inverse Method in solid mechanics: Theoretical underpinnings and novel applications, PhD Thesis, Universite Pierre et Marie Curie and Universita del Salento, 140 p. 33. Misra, S., Ramesh, K.T., Okamura, A.M., 2010, Modeling of Nonlinear Elastic Tissues for Surgical Simulation, Computer Methods in Biomechanics and Biomedical Engineering, 13(6), pp. 811-818. 34. Hackett, R.M., 2016, Hyperelasticity primer, Springer, Switzerland, 170 p. 35. Guo, Z., Caner, F., Peng, X., Moran, B., 2008, On constitutive modelling of porous neo-Hookean composites, Journal of the Mechanics and Physics of Solids, 56(6), pp. 2338–2357. 36. Quapp, K.M., Weiss, J.A., 1998, Material characterization of human medial collateral ligament, Journal of Biomechanical Engineering, 120(6), pp. 757–763. 37. Franulović, M., Basan, R., Prebil, I., 2009, Genetic algorithm in material model parameters’ identification for low-cycle fatigue, Computational Materials Science, 45(2), pp. 505-510. 38. Harb, N., Labed, N., Domaszewski, M., Peyraut, F., 2011, A new parameter identification method of soft biological tissue combining genetic algorithm with analytical optimization, Computer Methods in Applied Mechanics and Engineering, 200(4), pp. 208-215. 39. Reeves, C.R., Rowe, J.E., 2002, Genetic Algorithms—Principles and Perspectives, Springer US, New York, 332 p. 40. Blickle, T., Thiele, L., 1996, A comparison of selection schemes used in evolutionary algorithms, Evolutionary Computation, 4(4), pp. 361–394. 41. Xie, H., Zhang, M., 2013, Tuning Selection Pressure in Tournament Selection, IEEE Transactions on Evolutionary Computation, 17(1), pp. 1-19. 42. Reynes, C., Sabatier, R., 2012, A New Stopping Criterion for Genetic Algorithms, Proceedings of the 4th International Joint Conference on Computational Intelligence, IJCCI 2012, Barcelona, Spain, pp. 202-207. 43. Bhandari, D., Murthy, C.A., Pal, S.K., 2012, Variance as a stopping criterion for genetic algorithms with elitist model, Fundamenta Informaticae, 120, pp. 145-164. 44. Eiben, A.E., Raué, P.E., Ruttkay, Z., 1994, Genetic algorithms with multi-parent recombination, Proc. Parallel Problem Solving from Nature, PPSN 1994, Jerusalem, pp. 78-87. 45. Akbari, R., Ziarati, K., 2011, A multilevel evolutionary algorithm for optimizing numerical functions, International Journal of Industrial Engineering Computations, 2(2), pp. 419-430.