Study on Material Parameters Identification of Brain Tissue Considering Uncertainty of Friction Coefficient

Accurate material parameters are critical to construct the high biofidelity finite element (FE) models. However, it is hard to obtain the brain tissue parameters accurately because of the effects of irregular geometry and uncertain boundary conditions. Considering the complexity of material test and the uncertainty of friction coefficient, a computational inverse method for viscoelastic material parameters identification of brain tissue is presented based on the interval analysis method. Firstly, the intervals are used to quantify the friction coefficient in the boundary condition. And then the inverse problem of material parameters identification under uncertain friction coefficient is transformed into two types of deterministic inverse problem. Finally the intelligent optimization algorithm is used to solve the two types of deterministic inverse problems quickly and accurately, and the range of material parameters can be easily acquired with no need of a variety of samples. The efficiency and convergence of this method are demonstrated by the material parameters identification of thalamus. The proposed method provides a potential effective tool for building high biofidelity human finite element model in the study of traffic accident injury.


Introduction
The finite element (FE) method has shown its unique advantages in simulating irregular geometries, inhomogeneous media, nonlinear material properties, complex load and boundary conditions and so on. And FE simulation can capture the response that is difficult to be measured during the experiments. Therefore, many scholars study the intracranial response based on the FE model of the head, so as to better explore the damage mechanism and tissue level damage threshold. However, obtaining the accurate and reliable parameters of the brain tissue under various strain rates is helpful to improve the biofidelity of the FE model.
Based on the optimization strategy, the inverse method has increasingly become an important medium of obtaining biological tissue material parameters [1][2][3]. In the process of traditional biological tissue parameters identification, the positive problem model is usually built based on deterministic boundary conditions and solved by classical deterministic optimization method. However, due to the complexity of the contact surface in the actual biological tissue material experiment device, there are inevitably some uncertain factors in the modeling process of the biological tissue material parameters identification. Therefore, deterministic inverse method especially for the material parameters identification of biological soft tissue has some limitations. Traditionally, the material mechanical properties of soft tissue can be obtained via unconstrained compression experiments. In the experiment, the soft tissue is placed in the middle of two smooth plates, assuming that the sample maintains a uniform stress / strain state during the compression process. However, if there is friction between the sample and the plate, the sample can not maintain a uniform stress / strain state when compressed. Spilker et al. [4] have shown that the friction between the specimen and the plate had a significant effect on its mechanical response, but the effect of friction on material parameters of soft tissue was not further been quantified. To reduce the impact of friction on the test results, many researchers took various measures to reduce the friction between the sample and the platen [5,6]. However, in the actual unconstrained compression experiments, the friction between the sample and the platen cannot be completely eliminated regardless of the measures taken. Even with a well-lubricated contact, such as articular cartilage and stainless steel, the friction coefficient may exceed 0.2 [7]. Therefore, some data about compression characteristics of the soft tissue may contain errors caused by inaccurate friction.
In order to correctly estimate and evaluate the effect of friction on brain tissue material parameter identification results, this study takes the thalamus as the research object and aims at quantifying the influence of coefficient of friction on the material parameter identification based on the interval analysis method [8][9][10][11][12]. The possible range of material parameters can be quickly obtained via the optimization methods and interval algorithm. This will be useful to provide more accurate material parameters for the brain tissue FE model.

Experiment and Simulation Model Set Up
The test data came from the Wayne State University Bioengineering Centre [13]. Brain tissue samples were obtained from a 45-year-old male body (1 week after death). In the area of the thalamus (gray matter), two brain tissue samples with a size of about 15 mm by 15 mm by 8 mm were separated by surgical scalpel in the direction of the fibre. The smooth polyethylene plate was connected to the Instron Tester (Instron 1321 frame, 8500 controller) to load the sample. The unconstrained compression experiments conducted with low strain rate (0.5s -1 ) and high strain rate (35s -1 ) were both chosen to study the effect of friction on material parameter identification. The amount of compression of the low strain rate condition was 2.55mm, and that of high strain rate was 3.35mm.
Three-dimensional geometries of each sample were obtained using a laser scanner, and the hexahedral FE meshes were generated by ANSYS ICEM CFD. The mesh size is 0.45mm. The sensitivity analysis of the meshing shows that the convergence of the results can be sufficiently ensured by this partition density. The FE meshes of thalamus samples are shown in figure 1.  There are some typical and commonly constitutive models of viscoelastic materials, such as complex constant modulus model, standard rheological model, fractional derivative model, fractional exponential model and micro-vibrator model et al [14][15][16]. In this paper, the characteristics of high strain rate of brain tissue samples were characterized by using the material constitutive model (fractional exponential model) which is linear and viscoelastic. The shear relaxation behavior for this Where G 0 represents the short term modulus and G i represents the long term shear modulus, t is time and β indicates the deacy constant that can be determinted by performaing relaxation tests. For the low strain loading tests, the viscosity of the brain tissue can be neglected and FE simulation can be performed with a higher speed than the actual experimental speed to shorten the simulation time. Thus the behavior of the sample can be estimated by a linear elastic constitutive relationship formulated through the Young's modulus E: where K denotes bulk modulus. Since K»Gi, so E ≈3G i and the Poisson's ratio is defined: The loading plate was moved according to the motion of impactor recorded during the test. The contac t force between plate and sample was also recorded. The initial value of the friction coefficient betwee n the sample and the test plate was defined as 0.2.

Material Parameters Deterministic Identification of Brain Tissue Based on Specimen-Specific Fe Model
Combined with thalamic sample unconstrained compression experiments, the objective function of the viscoelastic material parameters identification of thalamic can be defined as following: Where the subscript m represents the measurement data in the compression experiment, and c represents the predicted data from the simulation model, and the superscript l represents the low strain rate test, and h represents the high strain rate test, and n represents the number of points on the impact force and deformation curve obtained from the experiment and simulation. It can be seen that the fitting degree of the experimental measurement data and the simulation data can be evaluated by the objective function. The function is closer to 0, the result of material parameters identification of biological tissue is better.
According to the literature [17][18][19][20] and E≈3G i , the ranges and initial values of linear viscoelastic material parameters (G 0 , G i , β) and elastic material parameter E are listed in table 1. Based on the FE method and sequence response surface method (FEM-SRSM) [21], the measurement data of compression experiments with two different strain rate were considered and the viscoelastic material parameters of the thalamus were investigated. The elastic material parameter E is used as a non-independent variable in the inverse process, and E is input into the three-dimensional FE model of low-strain rate compression experiment according to E≈3G i . The result is taken as the initial value of the inverse considering uncertainty.

Material Parameters Identification of Brain Tissue Considering Uncertain Parameters
Aiming at the uncertainty of the friction coefficient, a mathematical model of the material parameter identification of the biological tissue considering the uncertain boundary condition is established. The inverse problem is transformed into an optimization problem, and the following objective function is constructed: Where F represents the error between the response curve of the experimental measurement and the response curve of the simulated prediction; M is the parameter in the constitutive equation of the material corresponding to the specific biological tissue, such as the instantaneous shear modulus G 0 , the long-term shear modulus G i , and the reduction coefficient β in constitutive equation of viscoelastic material, etc.; U is the uncertainty parameter, such as the friction coefficient f in the boundary condition.
The objective function is nonlinear and continuous. Since it is a continuous function with respect to M and U, and U is an interval vector. Therefore, for any determined experimental measurement response curve, the obtained material parameters may take a range of values, rather than a constant values. Therefore, the above problems can not be solved by the traditional deterministic optimization method.

The Interval Model of Material Parameters Identification.
Based on the interval mathematical theory [22], the interval number is defined as an ordered pair of real numbers: In equation (6), the superscript I denotes the interval, the superscript L and R are the lower bound and the upper bound of the interval, respectively. When A L =A R , the interval degrades to a real number. For the material parameters inverse problem with uncertainty of the boundary conditions, the interval vector can be used: , , [ , ], 1, 2, , The biological tissue material parameters identified under this boundary condition are no longer a single real number but a bounded interval corresponding to all possible uncertainty parameters of boundary condition: ( ) , In equation (8)

Interval Analysis Method for Uncertain Boundary Conditions.
The interval number f I can be written as equation (10) based on the interval mathematics theory [22]: , 1,1 (10) In equation (10), f c and f w represent the midpoint and radius of the interval f I . The uncertainty level γ (f I ) of the interval f I is defined as [23]: (11) So the uncertainty vector f can be written as: ∈ − = f f (13) Assuming that the uncertainty level of each parameters in f is small, the first-order Taylor expansion of the uncertain material parameter M(f) can be performed at f c ( 14) The range of structural displacement response could be calculated by performing a natural interval extension on equation (14) 0.3 under well controlled test conditions. Therefore, in the boundary condition, the friction coefficient is set as an uncertainty variable, and the median value is 0.2. The viscoelastic material parameter of the tissue sample is inversed at 10% uncertainty level (i.e. f I = [0.18, 0.22]). The inverse problem is transformed into an optimization problem, and the following objective function is constructed: 100≤β≤900 The solution flow of viscoelastic material parameters of brain tissue considering the uncertainty of friction coefficient in the material experiment is shown as figure 2. To solve the positive problem model, the FEM -SRSM is adopted, and the genetic algorithm is used as the optimization solver.

Deterministic Inverse Results
The material parameters were obtained by the sequence linear response surface method (shown in table 2). The FE simulation was conducted using the above material parameters. And the comparison of the force and deformation curve obtained from experiment and simulation respectively is shown in figure 3. It can be seen that the two curves fit well, which proves the feasibility of the inverse method based on the FEM-SRSM and the correctness of the material parameters obtained by this method.

Uncertain Inverse Results
Considering the uncertainty of the boundary condition (friction coefficient), the inverse results of the parameters of the thalamus based on the interval analysis method are listed in table 3. It is obvious that the uncertainty level of the parameter β is the largest, namely the friction coefficient has the greatest influence on the material parameter β. Meanwhile, the inverse result can provide an optional parameter interval for the establishment of the FE model. The interval analysis method is based on the approximation of first-order Taylor expansion, which means the method can be only applied to the case where the uncertainty level of the parameter is small. If the uncertainty level is large, the accurate agent model can be constructed in the material parameter space, and then the upper and lower bounds of material parameters can be acquired using the uncertainty optimization method based on the agent model.

Conclusion
Considering the uncertainty of friction coefficient, an identification method for parameters estimation of brain tissue based on interval analysis theory is proposed in this paper. Firstly, the intervals are used to quantify the friction coefficient in the experimental boundary condition. The interval analysis method based on first-order Taylor expansion is adopted to transform the material parameter identification problem under the uncertain boundary condition into two types of deterministic inverse problems. The intelligent optimization algorithm is used to solve the two types of deterministic inverse problems quickly and accurately, and the upper and lower bounds of the material parameters are obtained by interval operation. This method can be used to evaluate the effect of uncertainty parameters on the inverse results of biological tissue parameters. This provides an effective method for the development of human FE models with high biological fidelity and high applicability. At the same time, this method avoids the time-consuming nesting optimization problem, which greatly improves the computational efficiency.

5.Acknowledgments
This research was funded by the National Natural Science Foundation of China (11402296).