Improving the Process of Adjusting the Parameters of Finite Element Models of Healthy Human Intervertebral Discs by the Multi-Response Surface Method

The kinematic behavior of models that are based on the finite element method (FEM) for modeling the human body depends greatly on an accurate estimate of the parameters that define such models. This task is complex, and any small difference between the actual biomaterial model and the simulation model based on FEM can be amplified enormously in the presence of nonlinearities. The current paper attempts to demonstrate how a combination of the FEM and the MRS methods with desirability functions can be used to obtain the material parameters that are most appropriate for use in defining the behavior of Finite Element (FE) models of the healthy human lumbar intervertebral disc (IVD). The FE model parameters were adjusted on the basis of experimental data from selected standard tests (compression, flexion, extension, shear, lateral bending, and torsion) and were developed as follows: First, three-dimensional parameterized FE models were generated on the basis of the mentioned standard tests. Then, 11 parameters were selected to define the proposed parameterized FE models. For each of the standard tests, regression models were generated using MRS to model the six stiffness and nine bulges of the healthy IVD models that were created by changing the parameters of the FE models. The optimal combination of the 11 parameters was based on three different adjustment criteria. The latter, in turn, were based on the combination of stiffness and bulges that were obtained from the standard test FE simulations. The first adjustment criteria considered stiffness and bulges to be equally important in the adjustment of FE model parameters. The second adjustment criteria considered stiffness as most important, whereas the third considered the bulges to be most important. The proposed adjustment methods were applied to a medium-sized human IVD that corresponded to the L3–L4 lumbar level with standard dimensions of width = 50 mm, depth = 35 mm, and height = 10 mm. Agreement between the kinematic behavior that was obtained with the optimized parameters and that obtained from the literature demonstrated that the proposed method is a powerful tool with which to adjust healthy IVD FE models when there are many parameters, stiffnesses, and bulges to which the models must adjust.


Introduction
The human intervertebral disc (IVD) is a fibrocartilage structure that is located between the vertebrae of the spine and absorbs the shock and pressure that are associated with daily movement. The healthy intervertebral disc provides mobility and spine flexibility during body movement and prevents excessive wear of the facet joints during daily movement of the spine. The behavior of the IVD is analogous to a ball full of water. The IVD is usually considered to be incompressible, similar to elastomers, due to its soft tissue and high water content. It permits limited motion while transmitting loads from one vertebra to another. Its complex structure permits significant mobility between adjacent vertebrae while transmitting considerable compressive loads. The IVD has three distinct regions [1]. These are the nucleus pulposus, the annulus fibrosus and the cartilage endplates. The endplate is composed of hyaline cartilage with a thin structure that covers the entire nucleus pulposus and about one third of the annulus fibrosus. The composition of the endplate is similar to a particular cartilage, but with a lower water content. The firm outer region, which is called annulus fibrosus, maintains the shape of the intervertebral disc. The annulus fibrosus is a complex ring structure that surrounds and protects the nucleus pulposus. It can be found between each pair of adjacent vertebrae. It is composed of concentric layers of fibrous tissue that have an approximate thickness of 1 mm. The fibers are an aggregate of collagen fibers and are alternatively oriented to the transverse plane at an angle of approximately ±30 • [2]. The annulus fibrosus serves to resist the nucleus pressure in the radial and tangential directions [3]. The fibers are surrounded by a hydrated proteoglycam gel, or annulus ground substance [4]. The inner portion of the annulus is fibrocartilage, which gradually blends into the nucleus pulposus. The inner portion, which is called the nucleus pulposus, is the soft spongy tissue that enables the disc. It functions as a shock absorber, absorbing the impact of the body's daily movements. The disc may be subjected to a complex combination of loads, but the nucleus pulposus helps to distribute the pressure evenly across the disc. This prevents the development of stress concentrations that could cause damage to the underlying vertebrae or their endplates. Many researchers have experimentally studied the behavior of the intervertebral disc of human cadavers by standardized tests (e.g., compression, flexion, extension, shear, torsion, and lateral bending) based on stiffness and disc bulges [5][6][7][8][9][10][11]. Today, due to medical ethics, the study of human cadavers is declining. Instead, many researchers carry out their studies with animals whose behavior is similar to that of humans [12]. Other researchers conduct their studies by modeling the behavior of the human body's biomechanics by numerical methods (e.g., Finite Element Method, Finite Volume method, Finite-Difference method, etc.). FEM is known to be a powerful tool to design and optimize mechanical devices, as well as different human body parts or even human prostheses for foot, knee or hip. This is especially true when their behavior is nonlinear [13][14][15][16][17][18][19][20]. Over the years, several researchers have used the FEM to study the mechanical behavior of the human IVD. However, the parameters that define its complex structures differ entirely in each of these studies. Also, the process of achieving the best parameters to correctly define the behavior of the IVD FE models according to the standard test by trial-error method is difficult. This may result in an unacceptable cost. A large number of FE studies have proposed methodologies to adjust the deformations [21], the reaction forces [22] and the intradiscal pressure [23] of the proposed FE models were developed essentially by the trial-error method. However, an optimization algorithm was developed to determine the relationship among the components of the collagen fibers, annulus fibrosus and ground substance of an FE model using the biomaterial properties that were extracted from the literature [24]. In this case, the stiffness of the fibers was varied to approximate Young's modulus of the ground substance in order to fulfil the required range of motions possible with the properties found in the literature. A method that was based on differential evolution was proposed by Ezquerro [25] to calibrate the FE model of a functional spinal unit that consisted of ligaments, nucleus pulposus, vertebral arch, and facet joints. In this case, the standard loads that were used were: flexion, extension, lateral bending and torsion. More recently, some researchers have used FEM and regression models based on support vector machine (SVM) combined with genetic algorithms (GA) to adjust the IVD FE model parameters [26]. Eleven material parameters were considered in this paper to define the behavior of the human IVD FE model. They were adjusted with nine different variations of stiffness and bulges from some of the previously mentioned standardized tests. These were two parameters for the nucleus pulposus (C10 and C0); two parameters for the cartilage endplate (E, µ); five parameters for the annulus fibers (Fiber12, Fiber34, Fiber56, Fiber78, and Fiber910), and two for the annulus ground substance (E, µ). The lowest value of error (around 6%) between the stiffness and bulges that were obtained from the proposed FE models and those obtained from the selected standardized test indicated that the proposed methodology was very effective. To date, no further references have been found in the literature for adjustment of parameters that define the non-linear behavior of human IVD FE models. However, some researchers have used the MRS method combined with FEM in order to study and optimize mechanical problems [27][28][29]. The current paper attempts to show how the combination of the FEM and the MRS methods with desirability functions can be used to obtain the most appropriate material parameters to use in defining the behavior of FE models of the human intervertebral lumbar disc. The proposed adjustment methods were applied to a medium-sized human IVD corresponding to the L3-L4 lumbar level with standard dimensions of width = 50 mm, depth = 35 mm, and height = 10 mm and developed as follows: First, a three-dimensional parameterized FE model that consisted of nucleus pulposus, cartilage endplate, annulus fibers, and annulus ground substance was developed and simulated with different standard tests. As in Somovilla Gómez et al. [26], and based on standard tests (e.g., compression, torsion, shear, flexion, extension, and lateral bending), three-dimensional parameterized FE models were generated that consisted of nucleus pulposus, cartilage endplate, annulus fibers, and annulus ground substance. Eleven parameters were used to define the parameterized FE models. Then, for each of the standard tests, regression models were generated to model six stiffness and nine bulges of the healthy IVD models as the parameters of the FE models were changed. Finally, an optimal combination of the eleven parameters was achieved by applying MRS based on desirability functions according to three different adjustment criteria. The first considered stiffness and bulges to be equally important for adjustment of FE model parameters. The second considered stiffness as most important, whereas the third considered bulges to be most important. An agreement between the kinematic behavior that was obtained with the optimized parameters and those obtained experimentally from the literature demonstrated that the proposed method is a powerful tool to use in adjusting healthy IVD-FE models when the latter have a high number of parameters, stiffness values, and bulges to which the models must adjust.

Use of FEM for Modeling the Lumbar Intervertebral Disc
A majority of the authors consider the nucleus pulposus to be an incompressible material. Its behavior has been modeled with FEM. It is considered to be a non-linear incompressible solid [30] as well as an incompressible fluid [11]. This non-linear behavior is due to its soft tissue condition and high water content [31,32]. In contrast, other authors have modeled the nucleus pulposus considering it to be a linear isotropic incompressible solid, with a Young's modulus of 0.1 MPa [2,[33][34][35], 1 MPa [36][37][38][39][40] and in a range of 0.5 to 1 MPa [41,42] or in a range of 1 to 4 MPa [43], whereas its Poisson's module was 0.4999 or 0.5. Other authors have considered the incompressible elastomer formulation based on the Mooney-Rivlin model for modeling nucleus pulposus behavior. In this case, C10 and C0 are the empirical constants parameters that are most frequently used. For example, Smit et al. [41] considered the nucleus pulposus as incompressible and used the Mooney-Rivlin model with the parameters C10 = 0.12 and C0 = 0.09 and µ = 0.4999 to model it. Some authors also have considered the Mooney-Rivlin model with values of C10 = 0.12 and C0 = 0.03 [3,30,44], whereas others have used values of C10 = 0.10 and C0 = 0.09 [45,46]. Similarly, Ibarz et al. [47] considered parameters of C10 = 0.0343 MPa and C0 = 0.1369 MPa (E = 1.0 MPa and µ = 0.49).
On the other hand, Lu et al. [48], Martínez et al. [49], and González Gutiérrez et al. [42] considered the endplates' cartilage behavior. They assumed an isotropic formulation with an elastic modulus E and a Poisson ratio µ of 20 MPa and 0.3, respectively. In a similar fashion, Kim et al. [36] and Dicko et al. [37], used the values for the endplates cartilage of E = 24 MPa and µ = 0.4, whereas Rohlmann et al. [46] and Schmidt et al. [44] used values of E = 23.8 MPa and µ = 0.8. However, Denozière [2] considered values of E = 12 MPa and µ = 0.3. The literature provides a Young's modulus for the annulus ground substance of 2 MPa to 10 MPa. However, the most commonly used value is 4.2 MPa. Denozière [50], Dietrich et al. [51], and Baroud et al. [30] selected values of 4.2 and 10 MPa, respectively, for Young's modulus. For Poisson's ratio, they used values of 0.45 to 0.35. Other researchers considered the annulus ground substance as an incompressible, and used the Mooney-Rivlin model with the parameters of C10 = 0.348, C0 = 0.3, and µ = 0.45 [45,46,52]. Schmidt et al. [44] used values of C10 = 0.10, C0 = 0.05 and µ = 0.45. Dicko et al. [37] considered C10 = 0.18, C0 = 0.04, and µ = 0.45. Ayturk et al. [38] considered the annulus ground substance as a Yeoh hookean hyperelastic material that is a model for the deformation of nearly incompressible nonlinear realistic materials such a rubber. For the collagen fiber layers, most researchers have chosen linear properties in their studies, which provide an average Young's modulus of approximately 360 to 550 MPa and a Poisson' ratio of 0.3 to 0.45 [2,30,33,35,36,38]. However, Grauer et al. [39] considered values of 175 MPa for all fiber layers. Other authors assumed a nonlinear stress-strain behavior of the annulus collagen fiber layers based on a nonlinear function [37,[44][45][46]52] that was obtained from previous reports that were based on the stress-strain curve of Shirazi-Adl et al. [11]. Ayturk et al. [38] considered the annulus fibers as Yeoh hyperelastic material with C10 = 0.0146, C20 = −0.0189, C30 = 0.04; a3 = 0.03, and b3 = 120. Table 1 shows the large number of parameters that some researchers have used during the last 30 years to develop intervertebral disc FE models. The latter adjusted the intervertebral disc FE models according to standard tests by trial and error, although this is difficult and very costly. To date, no author has made a methodical adjustment of these parameters. Only Somovilla Gómez et al. [26] have used FEM and Data Mining Techniques to adjust intervertebral disc FE models. In this paper, the authors propose the adjustment of the parameters that best define the behavior of the intervertebral disc. They propose the use of FEM and MRS with desirability functions and consider three different adjustment criteria. The first considers stiffness and bulges to be equally important, the second considers stiffness as most important, and the third considers bulges as most important. - [11] Incompress. Fluid ------4.2 0.45 σ = 23,000 × ε 1.9 Yeoh material. Material coefficients: C10 = 0.0146, C20 = −0.0189, C30 = 0.041; a3 = 0.03, b3 = 120.0 (b3 is unitless). C10 = 0.0343 MPa; C0 = 0.1369 MPa. An elastic analysis with a Young modulus of 1.0 MPa and Poisson ratio of 0.49 was conducted with similar results and a volume change of less than 0. * is Stress-strain curve by Shirazi et al. [11]: in this case, ε is the value for the deformation and σ is the stress.

Biomechanics of the Intervertebral Disc
The intervertebral disc (IVD) is one of the most important constitutive elements of a functional spine unit (FSU). An FSU consists of two adjacent vertebrae, the intervertebral disc, and all adjoining ligaments between them. The FSL excludes other connecting tissues such as muscles. Its complex structure allows significant mobility between two adjacent vertebrae while transmitting considerable compressive loads from one vertebra to another. Over time, the normal aging process causes the intervertebral discs to degenerate. This reduces the water contained in the intervertebral disc and the ability to absorb the impacts that are associated with spinal movements. Excessive pressure, strain, or injury to a rigid disc can cause the disc to tear or bulge. The degeneration of the disc's size and the loss of its functionality will bring the adjacent vertebrae closer to one another. This causes impingement and compression of a spinal nerve root. Nerve impingement can result in intermittent low back pain or leg pain, depending on the level of impingement [3]. This problem can be relieved by a surgical procedure called artificial disc replacement. In this procedure, the damaged disc is replaced by an artificial prosthesis (arthroplasty). Figure 1a,b show, respectively, an FSU with a healthy disc and an FSU with a herniated disc.

Biomechanics of the Intervertebral Disc
The intervertebral disc (IVD) is one of the most important constitutive elements of a functional spine unit (FSU). An FSU consists of two adjacent vertebrae, the intervertebral disc, and all adjoining ligaments between them. The FSL excludes other connecting tissues such as muscles. Its complex structure allows significant mobility between two adjacent vertebrae while transmitting considerable compressive loads from one vertebra to another. Over time, the normal aging process causes the intervertebral discs to degenerate. This reduces the water contained in the intervertebral disc and the ability to absorb the impacts that are associated with spinal movements. Excessive pressure, strain, or injury to a rigid disc can cause the disc to tear or bulge. The degeneration of the disc's size and the loss of its functionality will bring the adjacent vertebrae closer to one another. This causes impingement and compression of a spinal nerve root. Nerve impingement can result in intermittent low back pain or leg pain, depending on the level of impingement [3]. This problem can be relieved by a surgical procedure called artificial disc replacement. In this procedure, the damaged disc is replaced by an artificial prosthesis (arthroplasty). Figure 1a,b show, respectively, an FSU with a healthy disc and an FSU with a herniated disc.

Spine Movement
Spinal movement during daily activities is complex because it involves a combination of plane movements with several motions in axial, coronal, and sagittal planes. In the coronal plane, the spinal movement takes place when an individual bends forward or backward. This typically occurs when one is exercising or performing heavy tasks. Forward bending of the spine is defined as flexion, whereas backward bending is called extension. Movement in the sagittal plane occurs when bending laterally to the right or to the left. Bending the spine to the right is called right bending and to the left is called left bending. Movement of the spine in the axial plane occurs after applying

Spine Movement
Spinal movement during daily activities is complex because it involves a combination of plane movements with several motions in axial, coronal, and sagittal planes. In the coronal plane, the spinal movement takes place when an individual bends forward or backward. This typically occurs when one is exercising or performing heavy tasks. Forward bending of the spine is defined as flexion, whereas backward bending is called extension. Movement in the sagittal plane occurs when bending laterally to the right or to the left. Bending the spine to the right is called right bending and to the left is called left bending. Movement of the spine in the axial plane occurs after applying torsion clockwise and counterclockwise. Torsion of the spine is defined as torsion. Figure

Bulge
In this work, the radial bulge of the lower disc was examined in the healthy model. The bulge was calculated at only three characteristic locations of the disc (anterior, lateral, and posterior) as shown in Figure 3, when axial compression, flexion, extension and lateral bending were simulated. The bulges for the torsion and shear tests were not considered. The disc bulged anteriorly during flexion, posteriorly during extension, and toward the concavity of the spinal curve during lateral bending. One of the mechanisms of nerve root irritation and herniated disc is root impingement by disc bulge. Clinical interest has led to several in vitro studies that have carefully measured the disc bulge. In the early studies, only the compression load was applied [54], but in later studies, the effects of various other loads also were considered in several directions, such as flexion, extension, and lateral bending [55][56][57][58][59].

Bulge
In this work, the radial bulge of the lower disc was examined in the healthy model. The bulge was calculated at only three characteristic locations of the disc (anterior, lateral, and posterior) as shown in Figure 3, when axial compression, flexion, extension and lateral bending were simulated. The bulges for the torsion and shear tests were not considered. The disc bulged anteriorly during flexion, posteriorly during extension, and toward the concavity of the spinal curve during lateral bending. One of the mechanisms of nerve root irritation and herniated disc is root impingement by disc bulge. Clinical interest has led to several in vitro studies that have carefully measured the disc bulge. In the early studies, only the compression load was applied [54], but in later studies, the effects of various other loads also were considered in several directions, such as flexion, extension, and lateral bending [55][56][57][58][59].

Experimental Behavior of the Intervertebral Disc
For decades, most authors have studied the behavior of the human intervertebral disc by standard tests using human cadavers based on stiffness and disc bulges. The standard tests examine compression, flexion, extension, shear, lateral bending, and torsion. One of the first authors to study the behavior of the intervertebral disc was Virgin [60]. He investigated the elastic properties of the intervertebral disc in order to ascertain their contribution to the function of the vertebral column.

Compression
In this work, Virgin [60] obtained a stiffness value of 2500 N/mm in a compression tests when he applied a load of 4500 N. Other authors used loads of 5500 N to 1000 N for the compression test and obtained stiffness values of 3000 to 700 N/mm [3,54,[60][61][62][63][64]. In other works, for lower compression values with loads between 850 and 500 N, stiffness values between 2420 and 510 N/mm were obtained [5,7,[65][66][67]. Others authors applied loads with values of less than 400 N and obtained stiffness values of 800 and 250 N/mm [8,[68][69][70]. In regard to the bulges in the standard compression test, many authors analyzed the anterior, posterior, and lateral bulges of the intervertebral disc. In this regard, Shirazi-Adl et al. [11] used loads of 500, 720, and 1000 N and obtained values of 0.5-0.7-0.8 mm for the anterior bulge, 0.75-1.0-1.5 mm for the posterior bulge, and 0.35-0.4-0.6 mm for the lateral bulge. However, Denozière [2] considered a compression load of 2500 N and obtained results for the anterior, posterior and lateral bulges of 0.5, 0.7, and 0.4 mm. Table 2 summarizes the loads that were used in the compression test and the different stiffness and bulge values that were obtained by different authors.

Experimental Behavior of the Intervertebral Disc
For decades, most authors have studied the behavior of the human intervertebral disc by standard tests using human cadavers based on stiffness and disc bulges. The standard tests examine compression, flexion, extension, shear, lateral bending, and torsion. One of the first authors to study the behavior of the intervertebral disc was Virgin [60]. He investigated the elastic properties of the intervertebral disc in order to ascertain their contribution to the function of the vertebral column.

Flexion and Extension
One of the first studies of the behavior of Flexion and Extension of the human intervertebral disc was [61]. In this work, a load of 20 Nm was considered for both flexion and extension. Stiffness values of 4.5 Nm/ • were obtained in both tests. Similarly, Brown et al. [69] considered the same load (20 Nm) but obtained a stiffness value of 2 Nm/ • in both tests. Another study, Miller et al. [71], used higher loads (70 Nm), and obtained stiffness values of 5.51 and 7.60 Nm/ • for the Flexion and Extension tests respectively. Finally, Adams et al. [72] used a load of 80 Nm and obtained stiffness values of 7.3 Nm/ • for both tests. However, most researchers have used loads of 5 Nm to 10 Nm for their studies. For example, the work by Schultz et al. [9] considered a load of 10.6 Nm for a study of the behavior of the intervertebral disc, and obtained stiffness values of 1.92 Nm/ • and 3.53 Nm/ • for the Flexion and Extension test respectively. Adams et al. [73] considered a load of 10.7 Nm, and the stiffness obtained was 1.34 Nm/ • . Other authors used a load of 10 Nm for the Flexion and Extension test. They obtained stiffness values that ranged between 0.8 and 2.04 Nm/ • for Flexion, and between 2 and 3.53 Nm/ • for Extension [7,53,67]. Guan et al. [74], Busscher et al. [75] and [76], González Gutiérrez [3], and Patwardhan et al. [77] considered a load range of 4 to 8 Nm for both the Flexion and Extension tests. The stiffness values that they obtained from these tests varied from 0.8 to 1.18 Nm/ • for the Flexion test and from 0.8 to 1.53 Nm/ • for the Extension test. In regard to the study of anterior, posterior and lateral bulges by the Flexion and Extension tests, Reuber et al. [58] considered loads of 3.9 Nm to 7.9 Nm and obtained values of 0.73 and 1.11 mm for the posterior bulge, respectively. For the lateral bulge, they obtained values of 0.07 and 0.21 mm respectively. Finally, Denozière [2] used a load of 10 Nm and obtained values of 1.3, 1.9, and 2.6 mm for the anterior, posterior, and lateral bulges, respectively. Table 3 summarizes the stiffness and bulge values for the Flexion and Extension standard tests by the various authors.  [3,7,9,53,[74][75][76]78]. In addition, one of the few authors who were found in the extensive literature to have studied the anterior, posterior and lateral IVD bulges under a Lateral Bending load was Reuber, M. et al. [58]. In this work, two different loads of 9.8 Nm and 3.9 Nm were used. They resulted in bulge values of 0.49 and 1.13 mm for the posterior bulge respectively, whereas the lateral bulge values were 0.83 and 2.11 mm respectively (See Table 4).

Shear and Torsion
Only a few studies in the literature have studied the stiffness of the IVD in a shear test. One of these was conducted by Schultz et al. in [9,61], who used loads of 1000 N and 980 N and obtained, respectively, 685 N/m and 1000 N/m for stiffness values. Similarly, Weisse et al. [79] obtained a stiffness value of 830 N/m when they applied a load of 950 N. Another research team to conduct this type of test was Liu et al. [6]. In this work, the load applied was 450 N and the stiffness value that was obtained was 300 N/m. Other authors, such as Miller et al. [71] and Markolf [62], applied loads of 150 N for the same shear test but obtained different stiffness values (115 N/m and 260 N/m). Finally, Moroney et al. [68] used a much lower load than those used by the previously mentioned authors (20 N) and obtained a stiffness value of 60 N/m.
When the disc is subjected to a torsion load, the stress distribution in the IVD depends on the disc's degree of degeneration and whether the posterior elements are intact. Miller et al. [71] was one of the groups of investigators in the literature, who considered the highest torsion load. It was 70 Nm and produced a stiffness value of 10.9 Nm/ • . Farfan et al. [ [42,75,76]. Table 5 summarizes the stiffness values and the loads applied in shear and torsion tests.  [72] 10. 9 70 As described in the previous paragraphs, the published studies of the behavior of IVD are extensive. In this regard, the current work selected only fifteen different stiffness values and bulges, as well as the corresponding loads, for a search for the optimal combination of the disc's FE model parameters. Table 6 summarizes the 15 different stiffness and bulge values that were selected (six for stiffness and nine for bulges) and the load that each author, whose work was selected, applied in the corresponding standard test. These values were used to adjust the parameters that define the behavior of the FE model.

FE Model for Modeling the Intervertebral Disc Proposed
Like the authors' disc that was modeled by FEM and described in Section 2 (See Table 1), this paper considered that intervertebral discs have four main parts. They are the annulus fibrosus, the nucleus pulposus, and the two endplates that link it to the vertebrae. Figure 4 shows a 3D view of the proposed human healthy IVD FE model. In the lumbar spine, the nucleus is located between the middle and the posterior third of the sagittal diameter. Its volume is 30-50% of the total disc volume [53]. In some finite element analyses of the lumbar spine, the authors concluded that the nucleus was 43% of the total disc volume, which is the volume that the current paper assumes.

FE Model for Modeling the Intervertebral Disc Proposed
Like the authors' disc that was modeled by FEM and described in Section 2 (See Table 1), this paper considered that intervertebral discs have four main parts. They are the annulus fibrosus, the nucleus pulposus, and the two endplates that link it to the vertebrae. Figure 4 shows a 3D view of the proposed human healthy IVD FE model. In the lumbar spine, the nucleus is located between the middle and the posterior third of the sagittal diameter. Its volume is 30-50% of the total disc volume [53]. In some finite element analyses of the lumbar spine, the authors concluded that the nucleus was 43% of the total disc volume, which is the volume that the current paper assumes. Also, the nucleus pulposus was considered to be a non-linear incompressible and hyper-elastic Mooney Rivlin solid material. It was formulated according to empirical constants C10 and C0 [30,41]. It was assumed that the nucleus pulposus is surrounded by the annulus fibrosus. In order to prevent volumetric blockage, as the nucleus was considered to be practically incompressible, a full integration with Herrmann formulation was used for the eight-node solid FE elements that constituted the nucleus pulposus [83,84]. The annulus fibrosus was assumed to be a homogenous ground substance of hydrated proteoglycan gel that is reinforced by collagen fibers. The current FE model includes five layers of fibers that were oriented at an angle of ±30° relative to the transverse plane. The fibers were simulated by three-dimensional unidirectional line FE elements. The homogenous ground Also, the nucleus pulposus was considered to be a non-linear incompressible and hyper-elastic Mooney Rivlin solid material. It was formulated according to empirical constants C10 and C0 [30,41]. It was assumed that the nucleus pulposus is surrounded by the annulus fibrosus. In order to prevent volumetric blockage, as the nucleus was considered to be practically incompressible, a full integration with Herrmann formulation was used for the eight-node solid FE elements that constituted the nucleus pulposus [83,84]. The annulus fibrosus was assumed to be a homogenous ground substance of hydrated proteoglycan gel that is reinforced by collagen fibers. The current FE model includes five layers of fibers that were oriented at an angle of ±30 • relative to the transverse plane. The fibers were simulated by three-dimensional unidirectional line FE elements. The homogenous ground substance was simulated as eight-node isoparametric solid FE elements. For the five layers (Fiber12, Fiber34, Fiber56, Fiber78, and Fiber910), a range of E of 550 to 360 MPa, respectively, was assumed for their elastic modulus values. The endplates, which are composed of hyaline cartilage, link the disc to the vertebrae. This work considered cartilage endplates with an isotropic formulation (an elastic modulus E and a Poisson ratio of µ) for the FE model). A total of 11 material parameters were selected to define the behavior of the human disc FE model. They were adjusted after considering the above mentioned standard test. Table 7 summarizes these eleven different proposed material parameters.

Configuration of the Finite Element Analysis and Mesh Size
Due to the large displacements and deformations that the intervertebral disc suffered in this study, as well as the hyper-elastic behavior of the nucleus pulposus, a nonlinear analysis that involved large strain procedure with an updated Lagrangian formulation was used. Also, the mesh size that was established for each of the four parts of the intervertebral disc was based on the extensive literature available [85,86]. In this case, the largest mesh size for hexahedral elements of the endplate and annulus fibrosus was 2.2 mm, whereas the lowest was 0.5 mm, which corresponds to the thickness of the endplate itself (See Figure 5a). Also, the largest mesh size for the nucleus pulposus was 1.6 mm, and the smallest size was 0.3 mm (See Figure 5b). The openings in all of the mesh sizes that were established for the proposed FE models were smaller than each of the four parts of the intervertebral disc than those sizes of the FE models proposed in the literature. In other words, the discretization that is chosen can obtain satisfactory results without depending on further mesh refinement. substance was simulated as eight-node isoparametric solid FE elements. For the five layers (Fiber12, Fiber34, Fiber56, Fiber78, and Fiber910), a range of E of 550 to 360 MPa, respectively, was assumed for their elastic modulus values. The endplates, which are composed of hyaline cartilage, link the disc to the vertebrae. This work considered cartilage endplates with an isotropic formulation (an elastic modulus E and a Poisson ratio of μ) for the FE model). A total of 11 material parameters were selected to define the behavior of the human disc FE model. They were adjusted after considering the above mentioned standard test. Table 7 summarizes these eleven different proposed material parameters.

Configuration of the Finite Element Analysis and Mesh Size
Due to the large displacements and deformations that the intervertebral disc suffered in this study, as well as the hyper-elastic behavior of the nucleus pulposus, a nonlinear analysis that involved large strain procedure with an updated Lagrangian formulation was used. Also, the mesh size that was established for each of the four parts of the intervertebral disc was based on the extensive literature available [85,86]. In this case, the largest mesh size for hexahedral elements of the endplate and annulus fibrosus was 2.2 mm, whereas the lowest was 0.5 mm, which corresponds to the thickness of the endplate itself (See Figure 5a). Also, the largest mesh size for the nucleus pulposus was 1.6 mm, and the smallest size was 0.3 mm (See Figure 5b). The openings in all of the mesh sizes that were established for the proposed FE models were smaller than each of the four parts of the intervertebral disc than those sizes of the FE models proposed in the literature. In other words, the discretization that is chosen can obtain satisfactory results without depending on further mesh refinement.

Intervertebral Disc Dimensions
For decades, many authors have studied the kinematic behavior of human intervertebral discs of patients of different ages, sex, and stature, as well as different states of degeneration of the intervertebral disc itself. Many of these studies have been based exclusively on experimental studies using intervertebral discs from cadavers. Other works have been based mainly on FE models of discs for the study of their kinematics. In the studies that were based on FE models, the proposed models have been experimentally validated using data obtained by the authors themselves, or have been experimentally validated using data obtained from the literature. For a very broad range of patients that were studied, all lumbar intervertebral discs for the L1-L5 segment analyzed have had very similar widths, depths, and heights (See Figure 6a). For example, Zhou et al. [87] studied 126 patients with differing degrees of disc degeneration according to the Friberg and Hirsch scale [88]. Fifty-five male patients with a mean age of 50 ± 13.6 and 71 female patients with a mean age of 49 ± 12.04 and an age range of 22-80 years were studied. The approximate dimensions of the intervertebral discs analyzed in this case were, according to Figure 6a, a width of 46.1 ± 3.1 mm, a depth of 34.1 ± 2.6 mm, and a height of 12.4 ± 1.7 mm. Rostedt et al. [5] experimentally studied six lumbar FSUs from four subjects with different degrees of disc degeneration. The less degenerate discs (with a degree of degeneration equal to 1) were those of individuals of 45 years of age and a disc height of 12 mm. Schultz et al. [9] studied 16 UVF from 10 females and 26 UVF from 13 males with an age range of 21 to 60 and a mean age of 43 years with different degrees of degeneration. In this case, the less damaged discs were those of men of 35, 40, and 53 years of age with a disc surface of 15.9, 16.8, and 15 cm 2 , respectively. Other authors, such as Panjabi [89], studied UVF with healthy intervertebral discs from a sample of 60 patients of age 19-59 years (a mean of 46.3) and disc dimensions of width of 48.1 mm and a depth of 34.7 mm. Eijkelkamp [90] also studied a group of 60 patients of age 18 to 65 years, and observed that the mean intervertebral disc height was 13.5 mm. With a group of 157 patients and ages between 20 and 38 years, Nissan and Gilad [91] determined an average depth and height of the disc of 34.6 mm and 10.8 mm, respectively. With a smaller group (11 individuals) of age between 25 and 36 years, Tibrewal and Pearcy [92] found a depth and height of the disc of 33 mm and 9.8 mm, respectively. However, Wolf et al. [93] determined a mean width and depth of 44.1 mm and 31.7 mm respectively for a disc in a group of 55 individuals of 20 to 90 years of age. A more extensive study was that of Amonoo-Kuofi [94]. In the latter case, a group of 305 males and 310 females of 10 to 64 years of age and 10 to 61 years respectively were studied. The mean dimensions of the intervertebral discs were a depth of 42.8 mm and a height of 13.5 mm for males and a depth of 39.9 mm and a height of 13 mm for females. Other authors proposed FE models and validated them later with their own experimental data. For example, Schmidt et al. [24] developed a method to determine the individual contribution of the fibers and the ground substance to bending moments with four different magnitudes. They proposed a 3D, non-linear FE-model of an intervertebral disc from the L4-L5 FSU. In order to determine the appropriate geometry of the FE-model, six specimens were tested. The median anteroposterior dimension of the inferior endplate of L4 was 37.4 mm, whereas the median left to right dimension was 58.7 mm. The geometry used in the FE-model was determined from one of the tested specimens. The latter had dimensions nearest to the median with an anteroposterior distance of 38.1 mm and a left to right distance of 59.0 mm. Kim et al. [36] studied with the FEM the influence of the disc degeneration. A 3D nonlinear FE model of the lumbar spine was developed that consisted of four lumbar vertebrae, three intervertebral discs, and the associated spinal ligaments. Geometrical details of the human lumbar spine (L2-L5) were obtained from the high-resolution computed tomography (CT) images of a 46-year-old male patient who had no spinal deformities. In this case, the average area of the intervertebral disc was 1119 mm 2 . González et al. [42] studied the biomechanics of the intervertebral disc with different degrees of degeneration in compression using a FE model that was validated by experimental data. The study was performed on 5 patients aged 65-75 years for L2-L3 and L4-L5 lumbar levels. The height of the healthy intervertebral disc that was studied was 9.9 mm with an area of 1845 mm 2 for the L2-L3 level. For the L4-L5 level, it was 10 mm with an area of 1951 mm 2 . Shirazi-Adl et al. [11] proposed an FE model for analyzing the stress of the UVF para el nivel lumbar L2-L3 in compression for a 29 years old female patient. The geometry of the FE model that was analyzed was based in in vitro measurements. It had a width of 49.2 mm, a depth of 34 mm, a heigth of 11 mm and an area of 1371 mm 2 . Smit et al. [41] proposed a complete UVF FE model and studied the stresses on the cortical bone of the vertebrae. In that case, the intervertebral disc that was proposed in the FE model had a width of 42 mm and a depth of 35 mm. Finally, other authors proposed FE models and later validated them with experimental data from other authors. For example, Ibarz et al. [47] proposed a FE model for the Lumbar Spine to study the degeneration of the intervertebral discs. To verify the FE model, radiological images (X-rays) were taken of a group of 25 healthy male individuals who had an average age of 27.4 and average weight of 78.6 Kg. The results obtained from the simulation and from radiological measurements were compared to the results obtained from White and Panjabi [95]. Good agreement was obtained for the different movements. However, Ayturk [38] proposed an FE model to study the changes in lumbar spine mechanics due to degenerative disc disease. The experimental study was conducted on a 49-year-old female patient. The results obtained from the FE model were also compared to the studies that were conducted by Heuer et al. [96]. Other authors such as Weisse et al. [79] proposed an FE model for the determination of the translational and rotational stiffness of an L4-L5 FSU using the FEM. In this case, the geometry of the disc was obtained by computed tomography, and the results were compared subsequently to the studies by several of the aforementioned authors [87,[89][90][91][92][93][94]. Once the proposed FE model was simulated, the results were compared experimentally to results obtained from the literature [96][97][98]. Denozière and Ku [50] proposed a 3D FE model of a two-level ligamentous lumbar segment. The proposed FE model was validated with the clinical data obtained from the literature [98][99][100][101]. These authors developed an FE model of the intervertebral disc according to Figure 6a. In this case, the intervertebral disc had a width of 50 mm, a depth of 35 mm, a height of 10 mm and an area of 1440 mm 2 . Table 8 summarizes the anatomical dimensions of the L1-L5 lumbar segment and the differents intervertebral discs that were studied. As has been observed from the references in Table 8, the dimensions of the intervertebral discs are very similar for the range of patients studied. Also, it can be seen in this table that many of the mentioned authors validate their FE models with the experimental data obtained from the studies by other authors. For this reason, the FE model of the intervertebral disc that is studied in the current paper has been elaborated on the basis of the study by Denozière and Ku [50], which had the following dimensions according to where Fmax is the value of the maximum torque, which is 450 N in this case (see Table 6), and n is the number of nodes to which the load Fs is applied. The value of the stiffness of shear was obtained from the Y displacement of these nodes. In addition, for all previously proposed simulations, a boundary condition of embedding all nodes was imposed on the lower support of steel.

Design of the Experiments and Design Matrix
The response surface method (RSM) needs to establish a design of experiments (DoE) in order to obtain accurate models with the least amount of data to support the initial hypotheses [102]. Several methods have been proposed to develop DoE. However, they generally require the construction of a design matrix (inputs) to measure the outputs or responses of the experiments. One of the most widely used methods is the full factorial design method [103]. This method is based on the adoption all possible combinations of each of the values (or levels) and each of the factors that are considered in the DoE. Another DoE that is widely used is the 2k factorial design. This DoE is a

Boundary Conditions
Numerous authors have studied the behavior of intervertebral discs with the FEM. In doing do, they applied the contour conditions to these models and thus reproduced, as rigorously as possible, the corresponding experimental analysis [2,3,50]. In the current paper, the boundary conditions applied to the proposed FE models were applied in a similar way to the work that was developed by Denozière [2] and Denozière and Ku [50]. The proposed intervertebral disc FE model consisted of the intervertebral disc itself (nucleus pulposus, endplate, annulus fibrosus, and the fiber layers), as well as two steel supports to which the disc was glued (See Figure 6a). The boundary conditions applied to the proposed FE model to obtain the compression stiffness were applied as follows (see Figure 6b). On the upper steel support the pressure P c was applied. It was calculated by the following equation: where F max is the maximum compression test load, whose value was 500 N (See Table 6). S is the upper surface of the steel support, which has an area of 1440 mm 2 . The compression stiffness value was obtained from the "Z" displacement. The boundary conditions to obtain the stiffness value for the bending test were applied as follows (See Figure 6c): the corresponding maximum torque of the bending test (B max ) was applied on the anterior area of the steel upper support, which, in this case, is 5 Nm (See Table 6). The pressure P b was calculated as follows: where S a is the anterior area of the intervertebral disc, namely 745 mm 2 , and CG is the distance from the geometric center of the disc surface or instantaneous axis of rotation (IAR) to the center of gravity of S a . The bending stiffness value was calculated once the angle α had been obtained. The boundary conditions of the FE model that were required to obtain the lateral bending stiffness were applied in similar fashion to the process to obtain the bending stiffness (see Figure 6d). Over half area of the upper steel support was applied the pressure P lb , which was calculated according to the following equation: where LB max is the maximum torque of the bending test, which in this case was to 10.6 Nm (see Table 6). S/2 is the area of the intervertebral disc that is considered for the lateral bending, namely 720 mm 2 and CG is the distance from the geometric center of the disc to the center of gravity of S/2 area. The lateral bending stiffness was obtained from angle β. The boundary conditions necessary to obtain the torsional stiffness (see Figure 6e) are: on a group of nodes that belong to the periphery of the steel upper plate and a force F t (in this case in the Y direction) that is applied in a way that twists the disc. In order that the load F t is always at an angle of 90 • to a given R, a "following force" analysis was considered. The value of this force F t was calculated with the following equation: where T max is the maximum torque value. In this case, it is 10 Nm (see Table 6), Also, n is the number of nodes at which the load F t was applied and R is one the half of the width of the intervertebral disc (in this case is 25 mm). The value of the stiffness of lateral flexion is obtained from the angle η. Finally, the boundary conditions that are necessary to obtain shear stiffness (see Figure 6f) are as follows: on all of the nodes on the upper surface of the steel support, an F s load was applied in the Y direction. The value of this force F s was calculated with the following equation: where F max is the value of the maximum torque, which is 450 N in this case (see Table 6), and n is the number of nodes to which the load F s is applied. The value of the stiffness of shear was obtained from the Y displacement of these nodes. In addition, for all previously proposed simulations, a boundary condition of embedding all nodes was imposed on the lower support of steel.

Design of the Experiments and Design Matrix
The response surface method (RSM) needs to establish a design of experiments (DoE) in order to obtain accurate models with the least amount of data to support the initial hypotheses [102]. Several methods have been proposed to develop DoE. However, they generally require the construction of a design matrix (inputs) to measure the outputs or responses of the experiments. One of the most widely used methods is the full factorial design method [103]. This method is based on the adoption all possible combinations of each of the values (or levels) and each of the factors that are considered in the DoE. Another DoE that is widely used is the 2k factorial design. This DoE is a full factorial design that has two levels and generates 2k experiments in which k is the number of factors. In using this method to adjust the parameters of the FE model, the number of factors is k = 11 (C0, C10, Fiber12, Fiber34, Fiber56, Fiber78, Fiber910, Annulus_E, Annulus_µ, Cartil_E, and Cartil_µ) and the number of experiments or FE simulations needed is 2048. Similarly, a 3k factorial design generates 3k experiments or 177147 FE simulations. For both methods, this amount of data may be sufficient to cover completely the entire range of possibilities. However, it has the disadvantage of generating a large number of experiments or FE simulations. Another of the methods that are used to develop a DoE is Central Composite Design (CCD). This method is considered to be a fractional three-level design that is useful in obtaining regression models, thereby reducing the number of experiments in comparison to a factorial design 3k. A CCD generates 2070 experiments. This is very similar to the number that a 2k factorial generates. However, a Box-Behnken design (BBD) generates 176 experiments. This is an advantage because it significantly reduces the number of FE simulations [104]. In the current paper, the DoE was employed using a two-level fractional factorial design. In this case, 128 FE simulations were needed. This reduced amount of data may be sufficient to completely cover the entire range of possibilities and subsequently to search for the optimal parameters that define the behavior of the intervertebral disc. Table 9 shows some of the 128 combinations of the material parameters (C0, C10, Fiber 12, AnnulusE, etc.) that are generated with a two-level fractional factorial design when the ranges that appear in Table 7 are considered. The design matrix and their corresponding combination of parameters that define de behavior of the human lumbar Intervertebral disc FE models are generated by using the "R" statistical software (R 2014).

Response Surface Method for Modeling and Optimizing Problems
RSM is a method that attempts to determine the relationships between input variables and one or more response variables. The method was introduced by Box and Wilson in [104] and was used first to obtain a regression model from experimental data and an optimal response. Originally, RSM was developed to model experimental responses. However, it has been used recently to optimize products and industrial processes [105,106]. RSM is a group of statistical techniques that utilize a regression model that is based on a polynomial function (Equation (6)).
where Y is the response variables for the experiments and (x 1 , x 2 , x 3 , . . . , x k ) are the input variables, e is an error, and f is a function that consists of cross-products of the terms that form the polynomial.
To optimize the response Y, it is necessary to find an approximation functional relationship between the inputs and the response surface. A polynomial is the functional relationship used by RSM ((Equation (7)) in this kind of work.
where the first summation is the linear part, the second is the quadratic part and the third is the product of pairs of all inputs of the polynomial. The values of the coefficients b 0 , b i , b ii , and b ij are calculated by use of a regression analysis to determine the relationship between inputs and outputs. Also, the terms that are selected to form the equation are chosen according to their levels of significance [107]. One computes this level by the analysis of variance (ANOVA) and selecting terms according to the p-value obtained. When a problem has more than one output, as does this work, the optimization problem can be solved by using MRS [108]. In this regard, Harrington developed the desirability functions to obtain a compromise between the different outputs (Equations (8) and (9)) and the overall desirability. The latter is defined as the geometric mean of the desirability of each output (Equation (10)), [109].
In these equations, parameters A and B correspond to the limits of the input ranges, S is an exponent that determines the importance of reaching the target value, X corresponds to the input vector, and f r is the polynomial function that is used to predict the response. A second polynomial degree should be used to optimize one or several responses [108]. The desirability approach involves transforming each estimated response into a unitless utility that is bounded by 0 < d r < 1, where a higher value of d r indicates that the response value is more desirable. The optimization, which was developed with R statistical package, searches for a combination of importance factors that simultaneously satisfy the optimization criteria of each response and input [110].

Combining FEM and MRS to Optimize Mechanical Problems
Over the years, many researchers have used the FEM as an alternative to reduce costs during the design and optimization phases of mechanical problems. The FEM is widely used in industry to analyze engineering problems that are too complicated to solve by classical analytical methods. The object of its application is to reduce the expense incurred by experimental tests. FEM has the disadvantage of requiring a high computational cost especially when the process of fitting the proposed FE models is based only on the experience of the designer through simulations and trial and error. This adjustment process is greatly amplified when the FE model tries to solve non-linear problems, such as mechanical contacts, large deformations and hyper-elastic material behavior. However, the computational cost is reduced with the RSM, which provides an approximation of the same problem that is modeled with FEM [111]. Building regression models based on RSM is a good strategy to reduce the number of simulations that are necessary to model and optimize FE models. These models learn from the most characteristic samples that are obtained from FEM simulations and use their outputs. In this regard, the combination of RSM and FEM has been used widely to optimize many industrial processes. For example, Lin et al. [112] proposed a functionally graded material (FGM) as a potential upgrade to some conventional implant materials like titanium in prosthetic dentistry. In that work, computational bone remodeling and design optimization are used. Based on the results of remodeling, RSM has been adopted to develop a multi-objective optimal design for FGM implantation. Similarly, Sadollah and Bahreininejad [113] investigated the development of an optimal design of an FGM dental implant to promote long-term success. Also, Rungsiyakull et al. [114] combined RSM and FEM to establish a relationship between the surface morphology induced micromechanics and bone remodeling responses to a solid bead coated porous implant. In this case, the RSM was used to relate the major implant coating parameters to the bone responses. More recently, Bahraminasab et al. [115] studied the use of a functionally graded material (FGM) for the femoral component of knee implants. This is attractive because the properties can be designed to vary in a certain pattern to meet the desired requirements at different regions of the knee joint system, thereby decreasing the loosening problem. Therefore, a multi-objective design optimization of a FGM femoral component is conducted in this study by use of the FEM and RSM. The results of using an optimized FGM are then compared to those of using a standard Co-Cr alloy in a femoral component knee implant to demonstrate relative performance. In the current paper, the RSM used the results of FE analysis to search for the optimal parameter that defines correctly the behavior of healthy human lumbar intervertebral disc (IVD) FE models. Eleven parameters with which to define the parameterized FE models were selected. They were C0, C10, Fiber12, Fiber34, Fiber56, Fiber78, Fiber910, Annulus_E, Annulus_µ, Cartil_E, and Cartil_µ. For each of the standard tests, regression models were generated to model the six stiffness and nine bulges of the healthy IVD models when the parameters of the FE models were changed. The optimal combination of the eleven parameters was achieved by applying MRS based on desirability functions according to three different adjustment criteria.

FE Models' Results
After the parameterized healthy IVD FE models were created, an automatic procedure conducted the 128 simulations according to the design matrix that appears in Table 9. Table 10 shows the stiffness and bulges that were obtained from the FE simulations. They are named as follows and grouped according to the corresponding standard test. For the compression test, the anterior, posterior and lateral bulges values (Comp_bulgeA, Comp_bulgeL, Comp_bulgeP), as well as the stiffness (Comp_stiff), were obtained. The same applies to Shear, Extension, Lateral Bending, Flexion and Torsion tests: Shear_stiff, Exte_bulgeL, Exte_bulgeP, Exte_stiff, LBend_bulgeL, LBend_bulgeP, LBend_stiff, Flex_bulgeL, Flex_bulgeP, Flex_stiff, and Tors_stiff. These 128 values formed the training dataset and were used to generate the regression models. The following subsections show the process to generate and optimize these models. Table 10. Results of the simulation of FE models when a combination of 128 material parameters are considered in Table 9.

Run
Comp BulgeA

Analysis of Variance
Equation (6) was fitted with the data that appear in Tables 9 and 10 to obtain the regression equations for all responses by the use of the RMS "R" package [116]. Each of the responses or outputs from Equations (6) through (20) is shown. They were obtained by a combination of polynomials that are formed by input variables. The responses are: Comp_bulgeA, Comp_bulgeL, Comp_bulgeP, Comp_stiff, Shear_stiff, Exte_bulgeL, Exte_bulgeP, Exte_stiff, LBend_bulgeL, LBend_bulgeP, LBend_stiff, Flex_bulgeL, Flex_bulgeP, Flex_stiff, and Tors_stiff.
Comp_bulge L = 0.169454 − 0.084948 × C10 − 0.087169 × C0 − 5 × 10 −0.5 ×Fiber12 − 0.006264 × Annulus_E + 0.088348 × Annulus _µ −0.001015 × Cartil_E + 0.031074 × Cartil _µ (12) In order determine whether the variables that were used in the linear regression models are statistically significant, an ANOVA test was developed. Tables 11-25 show the ANOVA results. The tables show that most of the variables have a p-value of less than 0.01. This means that the models are statistically significant. In addition, the multiple correlation coefficient (R 2 ) was calculated as a measure of the variation around the mean that the regression model produced. The results show that all values of R 2 are close to one. This indicates that these models possess good predictive capacity.               Also, MAE and RMSE are calculated to determine the generalization capacity of the regression models that were obtained by using the results of the design matrix (inputs) from Table 9 and their corresponding outputs from Table 10 according to Equations (26) and (27).   Figure 7 shows the relationship between some of the actual values that were obtained experimentally with FEM (Table 10)  Additionally, 30 new FE models were created to test the proposed regression models with parameters that had not been used previously to generate regression models. Table 27 shows the errors incurred during the testing stage, when the maximum error corresponds to Tor_stiff (MAE equal to 10.06% and RMSE equal to 19.50%) and the minimum error corresponds to Flex_bulgeP (MAE equal to 2.86% and RMSE equal to 3.19%). The errors show that the adjustment between the regression models and the results obtained from the FE models is almost accurate. This also demonstrates its good capacity for generalization.     correlations are high, whereas the residuals that were obtained are small, which indicate that these regression models are adequate for the prediction of IVD behavior.

Multiple Response Optimization
Tables 28-30 show a combination of the eleven material parameter (inputs) that were studied in searching for the optimal behavior of human intervertebral lumbar disc FE models when using the desirability functions with the "R" package [110]. Three different adjustment criteria were considered.
The first column of the tables shows the material parameters (inputs) and stiffness values and bulges (outputs) that were studied. The second column shows the goals that were established in the goal setting process for both inputs and outputs when different adjustment criteria are considered. The third column shows the optimal values to define the behavior of the FE model and the last column shows the desirability values. Table 28 shows the results when all material parameters, as well as stiffness and bulges, were considered to have the same level of importance (equal to "inRange"). In this case, the value of the overall desirability was 0.625. Also, the table shows that some of the values that were obtained are very close to the targets that were proposed. For example, the target proposed for the Exte_bulgeL was 0.1 and the optimal value that was obtained was 0.100 with a desirability value of 1. In contrast, the proposed target for the Comp_bulgeL was 0.35 and the optimal value that was obtained was 0.0970 with a desirability value of 0.244.  Table 29 shows the results when the stiffness was considered to have a higher level of importance than bulges. This table shows that the material parameters that were obtained for all different goals required in the first criteria are very similar, although the value of overall desirability was 0.817, which is higher than the first criteria. Also, the target for the parameters of Shear_stiff and Tors_stiff were 300 and 2.1, and the values obtained were 300.000 (desirability = 1) and 3.456 (desirability = 0.428), respectively. Table 30 shows the results that were obtained with the third criterion when the bulges were considered to be the most important. Analogously to the other adjustment criteria that were previously studied, the material parameters that have been obtained are very similar for all different goals. In this case, the value of overall desirability was 0.554, which is the lowest result of the three criteria that are studied in this paper. Also, it is observed that the targets for the parameters of Exte_bulgeL and Comp_bulgeL were 0.1 and 0.35, and the values obtained were 0.103 (desirability = 0.947) and 0.101 (desirability = 0.257), respectively. Finally, three new FE models were simulated with the eleven different optimal material parameters that considered the three different adjustment criteria. These FE models were simulated again by the same standard test (Compression, Flexion, etc.) in order to compare the methodology proposed for adjust the optimal parameters. Table 31 shows a comparison of the results with the FE models using the optimized parameters, the optimal results from the regression models using the RMS with desirability functions and the experimental standard test. In order to compare the different errors that were obtained using three different adjustment criteria that were studied in this case, different MAE were obtained from the normalized data. The data is commonly normalized in statistical processes to transform all variables to the same scale (from 0-1). The transformation was achieved in this case by subtracting the minimum value from each original value and dividing the result by the range of each variable according to Equation (28).
where Y k,norm were the normalized outputs that were obtained from the results of the FE models with the optimized parameters and the outputs that were obtained experimentally from the standard test. The first column of the table shows the stiffness and bulges (outputs) that were studied. The second, third, and fourth columns show, respectively, the results obtained from the FE models with the optimized parameters for each of the different criteria that were considered. The fifth column shows the values of the standard test, and the last column shows the normalized MAE that was obtained. It can be seen in the table that the normalized MAE of the three criteria considered are very similar (Criteria 1 = 0.2782, Criteria 2 = 0.2795 and Criteria 3 = 0.2788). In contrast, the normalized MAE obtained for each of the outputs is lower when predicting the Shear_stiff (MAE = 0.01) and greater when predicting Comp_bulgeL (MAE = 0.73). The reason for this difference may be that the proposed FE model is generally less accurate in predicting bulges than stiffness. In addition, all of the values of MAE that were obtained for each of the different outputs or stiffness and bulges are in acceptable agreement.

Conclusions
This paper sets out a fully automated method that combines the finite element method (FEM) and multi-response surface method (MRS) based on desirability functions. This work looks for the optimal parameters that will correctly define the behavior of a medium-sized healthy human lumbar intervertebral disc (IVD) models based on FEM. First, based on standard tests (compression, flexion, extension, shear, lateral bending, and torsion), three-dimensional parameterized finite element (FE) models were generated. Then, 11 parameters were selected to define the parameterized intervertebral lumbar disc FE models. For each of the standard tests, regression models were generated for modeling the six stiffness and nine bulges of the healthy IVD models when the parameters of the FE models were varied. The optimal combination of the 11 parameters was achieved by applying MRS based on desirability functions according to three different adjustment criteria. The first criterion considered all of the material parameters (inputs), as well as the stiffness and bulges (outputs), with the same level of importance (equal to "inRange"). The second criterion considered, with a higher level of importance the stiffness, then bulges (goal of stiffness = target), whereas the third criteria considered, the bulges with a higher level of importance, then stiffness (goal of bulges = target). The best fit of the FE model parameters was achieved with the proposed second criterion. This produced a value for the normalized MAE of 0.2795. However, the results were very similar for the first and the third criteria that were considered. These were, respectively, MAE = 0.2782 and MAE = 0.2788. In contrast, the normalized MAE obtained for each of the outputs was lower when predicting the Shear_stiff (MAE = 0.01) and greater when predicting Comp_bulgeL (MAE = 0.73). This reason for this difference may be that the proposed FE model is generally less accurate in predicting bulges than stiffness. The MAE obtained from each criterion that was studied demonstrated that the proposed method is a powerful tool for adjusting healthy IVD FE models when there are many parameter and the stiffness and number of bulges to which the models must be adjusted, are high.

Conflicts of Interest:
The authors declare no conflict of interest.