A Multi-Component Conical Spring Model of Soft Tissue in Virtual Surgery

With the rapid development of optical system and medical robotic technology, the virtual surgery approach has achieved good progress but also faces many challenges. Physical and mechanical models of soft tissue characteristics have been a hot research topic in this field. In the present investigation, a physical model with a new multi-component conical spring is proposed, which replaces the cylindrical spring of the traditional mass-spring model. In the new model, a multi-component spring structure has been incorporated that consists of a conical spring and a linear spring in series. Based on the analysis of the nonlinear characteristics of soft tissue, the mechanical model was established to correspond better with the physical behavior. To determine the parameters for the mechanical model, the Euler method was deployed. The determination of the model parameters was achieved using physical tests. Result simulation experiments confirmed that the behavior of the new multi-component conical spring model proposed in the present paper gives results that are closer to those obtained from real-world physical tests than were data given by either the mass-conical-spring model or the mass-spring-model. Average difference, maximum difference and mean square error values were used for the evaluation. The accuracy of simulations for soft tissue virtual surgery can be improved by application of the new model developed in the present study.


I. INTRODUCTION
With the development of optical systems and computer vision technology, virtual surgery has revolutionized medical treatment, enhancing the accuracy and efficiency of procedures in surgical operations and surgery training. Precision medicine is a current direction of medical development and is a hot spot in medical development research. Virtual surgery is an important part of precision medicine and is an important means by which to improve clinical surgical skills. In consequence, it is an important research topic for many scholars [1]. The approach has been used mainly for simulation and training in high-difficulty and complex surgery procedures, using virtual reality technology, and is composed of soft tissue physics The associate editor coordinating the review of this manuscript and approving it for publication was Abdullah Iliyasu . modeling, collision detection, force feedback calculation, visual rendering and so on.
In order that better results can be obtained using virtual surgery, i.e. that visual and tactile immersion in the operation is more realistic, it is helpful if modeling performance can be improved. Finite element models have high precision, but their calculation is cumbersome and time-consuming, and this affects real-time interaction during virtual surgery. Therefore, it was decided that during the present investigation, analysis and research would focus on the basis of the mass spring model, which has good real-time performance, though more moderate precision. A more refined model could be obtained by improving the model spring module to achieve a more realistic experience during virtual surgery.
In terms of physical modeling, a large number of related studies have been conducted and some progress has been VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ made. Zou et al proposed a deformation model based on point interpolation with the meshless radial basis function and greatly improved the topology simulation of complex changes by using the Kelvin model [2]. Ma et al proposed a new method for designing and manufacturing pneumatic systems with ideal deformation responses [3].
Moerman et al discussed the Odgen superelastic approach from the perspective of tension-compression asymmetry. The proposed method can help researchers control independently the degree of tension-compression asymmetry from the degree of nonlinearity, and in the case of anisotropic materials, it can help when studying the modeling of biological tissues [4]. Qian et al proposed a powered soft tissue anatomical model that was simulated using an edge-based anatomical algorithm [5]. Kim et al proposed a layered simulation model, combining volumetric statistical model driving and finite element method driving to simulate the response of human soft tissue [6]. Latorre et al proposed a continuous model based on an isotropic finite element framework to predict tensile and compressive asymmetry stresses [7]. Hou et al proposed a new constraint-based finite element model for brain tissue that was implemented by introducing implicit integration and Graphics Processing Unit calculation [8]. Pereira da Silva et al proposed a new method based on data-driven parameterization of the mass-spring model to achieve object deformation [9], [10]. Tang and Wan proposed a new strain-based constrained finite element method for simulating effectively nonlinear uniform soft structures [11].
Omar et al introduced the stiffness variation of a conical spring in a conventional mass-spring model, which used only a single formula to simulate different types of deformation modes of soft tissue, thereby simplifying and improving the efficiency of soft tissue deformation modeling [12]. Merker et al simulated the hybrid biped walking system using a mass-spring model and solved the complexity of the system while calculating the corresponding boundary problem [13].
Li et al proposed a method that combined the modal equivalent method with the mass-spring model to construct a two-dimensional equivalent multiple (soil seismic response analysis) model [14]. Sanmiguel-Rojas et al analyzed oscillators consisting of two nonlinear springs connected in parallel and in series using an equivalent stiffness method [15]. Gzaiel et al studied the stress response state of a material with a sharp blade for puncture cutting by using finite element analysis [16]. Singh et al solved the motion planning problem of robots avoiding obstacles in static and dynamic environments by using fuzzy logic, genetic algorithms and three-path concepts [17]. Fuentes-Aguilar et al constructed a soft object deformation model based on the finite element method and differential neural network method [18]. Horvat et al implemented a three-dimensional constrained mixed growth and remodeling simulation using finite element analysis [19]. Rakshit et al proposed using the concept of matrix vectorization and the Kronecker product of two matrices to solve the updating problem in finite element models [20]. Zhang [29]. Based on an extensive research of the literature, it was identified that conical springs have a more compact structure than the classic linear cylindrical springs and have a more stable force when a load is applied. At the same time, studies of the similarity between the characteristics of conical springs and soft tissue, and it has been proposed on the basis of the present research that a multi-component-conical-spring model should replace the cylindrical spring component in the traditional mass spring model. The model could be improved further, compared to earlier research reported by the investigating group on nonlinear physical modeling of soft tissue and its parameter determination.
The main contributions of the present paper are as follows: (1) A new mass spring model is proposed. The new multi-component conical spring model formed by a series of conical springs and linear springs has been derived mainly from the analysis of the mechanical properties of soft and conical springs. (2) The least-squares method and experimental data from actual physical measurements were used to calculate the unknown parameters in the model. The accuracy of the new model in simulating force feedback was evaluated to be higher than that of earlier models.

A. THE CONICAL SPRING
Conical springs, also known as conical coil springs, have varying pitches. Their characteristic lines are non-linear, and their construction is compact, making the system more stable when loads are applied [30]- [32]. The conical spring can be considered to exhibit two phase relationships, which are the linear phase and the nonlinear phase, when loaded during the pressure process [33], [34], as shown in Figure 1. Such a force phase relationship is similar to that of soft tissue in virtual surgery.

B. THE MECHANICAL MODEL OF A CONICAL SPRING
In the ideal state of the spring system, its load can be considered to be as close as possible to the axial load. When the coil of the conical spring reaches the maximum physical deformation under the action of an external force, deformation will no longer occur, and it is no longer the active coil of the model. Through analysis of the conical spring, the T point is defined as the ''Transition Point'', which is the endpoint of the initial linear phase of loading. After that the deformation and force of the coil exhibit nonlinear behavior. When the spring coil is deformed by force to make the adjacent coils tight, the tight cone spring coil is called a passive coil and no longer participates in subsequent deformation. Therefore, as the force increases, the number of coils that can be changed in the system model will be fewer and fewer, and the index of the spring will become smaller and smaller, which makes the spring stiffer. In the limit, the number of active coils is zero and deformation is stopped. As such, the point C is defined as the ''Maximum Point''. Therefore, the characteristic curve of a conical spring is a mixture of linear and nonlinear. The main structural parameters of one conical spring are shown in Figure 2 [35].
In Figure 2, n is the number of active coils, d is the diameter of the spring wire, L 0 is the total length of the spring, and D 1 and D 2 are the diameters of the smallest and largest coils, respectively.

1) LINEAR PHASE
During the initial linear phase of the conical spring curve (as shown in Figure 1), the maximum coil can be considered as the active coil with the largest load, the characteristic relationship between load and length is linear, and the spring stiffness k is constant, and therefore: where G is the shear modulus of the material. For the constant k, the expression of the spring system formula can be obtained in the linear phase as follows: where: After transformation, this becomes: In the formulae (2) and (3), L is the length after spring deformation in the linear phase, and L 0 is the initial length of the spring. F T and L T are the load force and the displacement of the transition point of the spring from the linear phase to the nonlinear phase, respectively.

2) NONLINEAR PHASE
The spring is in a nonlinear phase when the load applied to the spring system is greater than the transition load F T . The nonlinear phase has a corresponding nonlinear load-displacement function, and the nonlinear load-displacement of each spring is different [36]. Based on Wahl's basic cylindrical assumption, the conical spring considered in the present study was analyzed [37]. The coil is discretely composed of several arc-measuring parts. Under a given load, the total deformation of the conical spring is determined by the additive deformation of each part, and the individual deformation is constrained by its maximum geometric value, which is illustrated in Figure 3 [38], [39].
The relationship between the number of free coils n f and the number of initial active coils n is as follows (where Lc is the incompressible length of the spring): The total deformation of the spring system is approximately the sum of the deformations of the free coils f and the solid coils c , as shown in the following formula: where δ f is the axial deformation component of the free coil and δ c is the axial deformation component of the solid coil corresponding to formulae (8), (9) and (10): In formula (8), D(n D ) is the variable diameter of the conical spring, and n D is the number of variable coils: The deformation of a spring model is reflected mainly in the spring length variation. In the present study, the length was chosen as a function of the load in the nonlinear phase. Formula (11) is obtained by integrating the total deformation of the springs of formulae (6), (7), (8), (9), and (10). Therefore, when the parameters in the spring system are given, the relationship between the deformation length of the conical spring and the load-bearing force F can be expressed as: After substituting formula (5) into formula (11), the stiffness coefficient of the nonlinear phase of the conical spring is obtained by the derivative transformation of the formula, as shown in formula (12):

C. MULTI-COMPONENT CONICAL SPRING MODEL
The multi-component spring model in [40] is the main structure for constructing the 'multi-component conical spring model' proposed in the present paper. This model is composed of a conical spring S e and a pre-variable spring S p in series, where i and j represent two adjacent nodes. The prevariable spring S p has a stiffness of K p and a pre-tightening force of F T . 2 indicates the total displacement of the spring model and 1 indicates the displacement of the conical spring, as shown in Figure 4. For the new model proposed in the present paper, the load applied to the spring model will delay the deformation of the pre-variable spring S p in the spring model. The conical spring and the pre-variable spring will participate in the deformation together until the applied load is greater than or equal to the pre-tightening force F T . In view of the real-time requirements of virtual surgery, it was specified during the present study that the conical spring S e would work alone when it is in its linear phase. When the two spring types participate together, the deformation is produced by the nonlinear phase of the conical spring S e together with the pre-variable spring S p . The working process of the spring model is mainly as follows: When the applied load F is smaller than F T , only the conical spring participates in the deformation in the spring model due to the preload force, and the displacement and load relationship is as follows: where k can be obtained from formula (1). When the applied force F is greater than the preload force F T , the conical spring and the pre-variable spring in the spring model will participate in the deformation together, i.e., the displacement will continue to change on the basis of 1 . In such cases, the relationship between displacement and load in the spring model is as follows: In formula (14), 2 is the total displacement of the entire spring, G P is the shear elastic modulus of the pre-variable spring, n p is the effective number of turns of the pre-variable spring, D P is the linear spring mean diameter, and d P is the linear spring wire diameter. By simulating the mechanical behavior of the spring model, the combined curve is smooth, as shown in Figure 4(b) above, and is not segmented. It is also consistent with the nonlinear characteristics of soft tissue.

D. MECHANICAL ANALYSIS BASED ON THE NEW MODEL
For the model proposed in the present paper, the 'finverse' function in MATLAB was chosen to inverse the displacement and load relationship of formulae (13) and (14), to obtain the function of the elastic force F on the displacement of the new multi-component conical spring model. Then, according to Newton's law, the mechanical equation of motion of each particle affected by the load can be expressed as: where m i is the mass of the mass point i, and x i is the coordinate of the mass i, which is a three-dimensional vector.
x i andẍ i are the velocity and acceleration of the mass, F ci is the damping of the mass, and F i is the spring force of the other mass on the mass i. f i is the external force on the mass point i. The specific forms of F i and F ci are as follows: In the above formula (16), |x j − x i | is the relative distance between two mass points, and |x j − x i | 0 is the initial length of the spring before deformation. K is the spring coefficient, C is the damping coefficient, j 1 , j 2 , . . . , j q represents the mass point connected to the particle i.   (15) is established for each mass point of the whole model; that is, the relationship between force and displacement in the structure, and they are combined to form the system of motion formulae. Suppose there are N mass points in the model and the position vector of each mass point is represented by a vector, if a 3N-dimensional displacement vector x is introduced, formula (15) can be rewritten in the matrix form as: In formula (18), M and C are 3N×3N vector matrices, respectively representing the mass diagonal matrix and the damping coefficient matrix. Both F and f are 3N-dimensional vectors, representing the spring force vector of the mass point and all external forces. When the mass point has no external force, according to the principle of force balance, the initial coordinates of each mass point can be calculated. When the mass point is subjected to an external force, the force state of the model is in an unbalanced state. According to the magnitude of the external force, the kinematic formula can be calculate to obtain the new coordinates of the mass point in the equilibrium stress state, such that the simulation of the deformation of the model can be realized. The model solves the deformation process of soft tissue under the condition of a known force. The physical quantity to be solved by the new model is the position and velocity information for a mass point in the soft tissue model. In formula (18), as the position vector in the mechanical equation is not explicitly given, it is necessary to discretize the formula in time to obtain a numerical solution. Formula (18) is converted to the following first-order differential formulae: where v represents the velocity vector of the entire mass model. In the present study, an improved Euler method was used to solve the differential formulae. Due to the low accuracy of the initial value of the model, the velocity v i is VOLUME 8, 2020 calculated by the explicit Euler method as the initial value, and the calculation result for v i is used directly to calculate x i in the implicit Euler method, which reduces the iteration process for the implicit Euler method. This approach not only reduces the amount of calculation but also guarantees the accuracy of the implicit Euler method.
where v t i is the velocity of the mass point i at time t, v t+ t i is the velocity of mass point i at the time t + t, x t i is the position vector of mass point i at time t, and x t+ t i is the position vector of the mass point at the time t + t. By cyclic calculation of the above formulae (21) and (22), the variation with time of velocity, position, and force of other mass points in the force region can be obtained. A series of discrete position vectors x(0), x( t), x(2 t), and so on, can be obtained by numerical solution, where t is the time step of the force deformation process. These discrete position vectors are dynamically and continuously output to achieve simulation of the deformation of the soft tissue (the pseudo-code of the model calculation is shown in Figure 6). This section first introduces the stress characteristics of the conical spring and its linear and nonlinear phases. Through the mechanical analysis of the conical spring, the new multi-component conical spring model proposed in the present paper was derived. Then, based on Newton's mechanical equation, the corresponding mechanical model was established. Finally, the deformation calculation of the model was completed by the modified Euler method. In the next section, we introduce mainly the parameter solution in the mechanical formula of this section, and the comparative physical tests, using the method introduced in this section.

A. TEST METHODOLOGY
In order to verify the proposed new multi-component conical spring model, a micro-controlled electronic universal testing machine (model WDW3300) produced by Changchun Kexin Test Instrument Co., Ltd. was used to measure and analyze the liver tissue of pigs. The liver test samples were cubes with a length, width, and height of 1 cm, to which a force of 1 N, at a test speed of 0.001m/s was applied. The physical test equipment measured the curve of the external force of the sample and the vertical axial displacement.
Because the precondition for modeling was that the soft tissue had body invariance, the physical test result is a weighted average of multiple measurements, and the results were converted into a curve of the pressure and displacement variation ratio. The displacement variation ratio is the ratio of the displacement to the original axial length. The tests included compression stretching tests. However, the discrimination of the stretching data was low, which is not conducive to the determination of model parameters and model comparison. Therefore, only the compression data were selected for subsequent analysis. The measurement process is shown in Figure 7, and the measurement result data were organized as shown in Table 1.

B. PARAMETER DETERMINATION
According to the data shown in Table 1, the parameters of the new multi-component conical spring model proposed in this study were determined by the least-squares method (its pseudo code is shown in Figure 8), and the results are shown in Table 2.

C. THE CONTROLLED MODELS
There are two controlled models, which are the mass-conicalspring model and the mass-spring model.
The formula of the conical-spring model is as follows: (23) and K 1 , K 2 , K 3 are calculated using the following formulae: The formula of the mass-spring model is as follows: where T s is the spring model stress, λ is the elongation ratio, and C 1 , C 2 , C 3 are the constant coefficients of the mass spring model. According to the test data shown in Table 1, the parameters of the mass-conical-spring model [12] also are obtained  using the least-squares method and the data fitting method, as shown in Table 3 below. At the same time, according to the physical test data in Table 1, the mass-spring model [41] was calculated and fitted. The results are also shown in Table 3.
According to the data in Table 1, the mass-spring model, the new model, and the mass-conical-spring model were used by MATLAB to simulate soft tissue deformation with the same external force. The simulation results for the load and displacement are shown in Table 4.
This section first introduces the physical test equipment and the measured real soft-tissue data and then applies the data to determine the relevant parameters of the new model proposed in the present paper. The new model then is compared with the real data, and with the other two models through simulation, and the evaluation of the comparison data is shown in the next section.

IV. DISCUSSION
To compare the accuracy between different models, it was necessary to determine the simulation model parameters based on unified physical test data firstly, then obtain the displacement results of the models under the external forces; and then to evaluate the performance of the models according to closer evaluation criteria, based on the physical test results. VOLUME 8, 2020  According to the data obtained from Figure 7, the various simulations were carried out in the MATLAB environment, and the results are shown in Figure 9, which is a comparison diagram according to the chosen evaluation criteria based on the above approach. It can be seen from Figure 9, compared with the other two models, the simulation results from the new model match the deformation data of the real soft tissue test measurements better, hence they reflect better the characteristics of the soft tissue.
As shown in Table 5, the simulation data were analyzed quantitatively, and the simulation data of the three models were evaluated for their mean difference, maximum difference and mean square error, as indicators. According to the data in Table 5, when the vertical axis pressure is constant, the average difference between the mass-spring model and the liver soft tissue physical test displacement variation ratio was 0.13, and the maximum difference was 0.18. The average difference between the mass-conical spring-model and the liver soft tissue displacement variation ratio from the physical tests was 0.028, and the maximum difference was 0.066. The average difference between the new model and the liver soft tissue displacement variation ratio in the physical tests was 0.018, and the maximum difference was 0.050. The data show that the new model has a better prediction capability than either of the other two models in terms of mean square error, maximum difference and average difference, which means better accuracy in the stability and model simulations. Therefore, based on the above analysis of physical test data and simulation data, the new multi-component conical spring model was demonstrated to improve upon the deficiencies of the existing models of the soft tissue structure. Meanwhile, the real liver tissue was subjected to interrelated axial compression animation simulations using the finite element analysis software ABAQUS 6.14 on the Personal Computer of Intel (R) Core (TM) i3-8100 GHz CPU and 8G memory which is equipped with AMD Radeon R7M265 Series graphics card. For the physical model, we chose the same size as the liver tissue used in the real measurement experiment, and set the material properties, assembly, loads, and other parts in different models. As shown in Figure 10, the bottom surface of the physical models were fixed during the axial loading process, and the top surface were continuously loaded an axial pressure of 0 to 14 kPa with the interval of 0.01 kPa. Through the comparative analysis of different deformation stages from Figure 10, it is obvious that the new model proposed in this paper has more similarities in mechanical behavior to the real liver tissue. VOLUME 8, 2020 TABLE 5. Comparison of the mean difference, the maximum difference, and the mean square error between the three simulation models and the real experiment of soft tissue.  In order to evaluate the applicability of the proposed method, OpenGL was used for the mainstream graphics application programming interface, which not only has strong portability and high-performance graphics rendering capabilities but also is suitable for a variety of professional computers, such as medical reality devices. Therefore, the OpenGL environment was configured based on VC++6.0, and different loading tests were performed on the proposed new multi-component conical spring model. The model and the deformation effect are shown in Figure 11.
Through the comparison of the above modeling effects, the model was stable under the application of different loads, the deformation effect was clearly evident, and the local variation was balanced. This verified the applicability of the new multi-component conical spring model proposed in the present paper.

V. CONCLUSION
As a multi-disciplinary subject that combines optical engineering and computer technology, virtual surgery plays a vital role in the field of clinical medicine. Compared with HD video technology, which is relatively mature, real haptic feedback is crucial to its successful application. Therefore, accurate soft tissue modeling methods will be key to achieving these desired objectives.
The soft tissue new multi-component conical spring model proposed in the present paper integrates a conical spring into the multi-component spring model after analyzing the force characteristics and the mechanical model of the conical spring. In the subsequent simulation tests, the data fitting curves for the mass spring model, for liver soft tissue derived from physical measurements and for the mass-conical spring model were compared. It was verified that the new multicomponent conical spring model reflected accurately the deformation of compressed soft tissue. This capability will improve the effectiveness and the quality of simulated force feedback in relevant virtual surgery procedures.
In future research, it is envisaged that a more generalized solution could be developed for modeling other soft tissue applications.