Abbass Ballit
email: abbass-khodor.ballit@utc.fr

Imag Mougharbel

Hassan Ghaziri
email: hassan.ghaziri@gmail.com

Ph.D Tien-Tuan Tuan Dao
email: tien-tuan.dao@utc.fr

Fast Soft Tissue Deformation and Stump-Socket Interaction toward a Computer-Aided Design System for Lower Limb Prostheses

Keywords: Lower Limb Prosthesis, Soft Tissue Deformation, Stump-Socket interaction, Mass-Spring Model, CAD System

published or not. The documents may come

I. INTRODUCTION

Lower-limb prostheses are generally developed with three main parts: the socket, the pylon (leg section), and the foot. In particular, the socket is considered as an important element in the design process of the prosthesis. Note that each socket is a tailor-made product designed to fit with the unique geometry of the patient's residual limb [START_REF] Herbert | A Preliminary Investigation into the Development of 3-D Printing of Prosthetic Sockets[END_REF]. A poorly fitting socket can cause significant trauma. Current practice in the design and fabrication of a prosthetic socket usually starts by creating a wrap cast of the residual limb, which is filled by plaster to form the positive mold of the stump (Figure 1). Using his experience, the prosthetist applies some modifications to the positive mold before a socket is formed over it. These modifications aim to produce the expected pressure distribution on the stump/socket interface during usage of the prosthesis. During the fitting stage, the socket is usually modified several times until achieving a well-fit socket. The decision is based on the feedbacks of the patient and the personal judgement of the involved specialist. This process is time-consuming and expensive requiring incremental refinements for the final prototype of the socket, and especially its accuracy depends strongly on the expertise of a high cost specialist [START_REF] Lorenzelli | Socketmaster: Integrated Sensors System for the Optimised Design of Prosthetic Socket for above Knee Amputees[END_REF]. All these elements lead to the unaffordable high prices of the final products: a simple lower limb prosthesis that allows its user to simply stand and walk in level ground costs between 5000$ and 7000$ [START_REF] Mcgimpsey | Limb Prosthetics Services and Devices[END_REF], and even the most expensive ones are built to withstand only three to five years of wear and tear, meaning they will need to be replaced over the course of a lifetime. Therefore, an alternative effective low-cost socket design and fabrication process is needed in order to reduce the incredibly high price of these products. Many researches have been conducted for improving the socket design and fabrication process using Computer-Aided-Design approach. In 2007, Bonacini et al. [START_REF] Bonacini | 3-D Digital Model Reconstruction: Residual Limb Analysis to Improve Prosthesis Design[END_REF] proposed a socket design process where all phases are computer-aided, but they focused on the reconstruction of the stump model. They analyzed 3 different technologies to obtain the geometry model of the stump, which are based on laser scanning to acquire external geometry of the stump, and computer tomography (CT) and magnetic resonance imaging (MRI) for the inner parts. This study did not consider other steps of the process. Another computer-aided process was proposed by Hsu et al. in 2010 [START_REF] Hsu | The development of a rapid prototyping prosthetic socket coated with a resin layer for transtibial amputees[END_REF] to design and develop the prosthetic socket. This process follows the same steps of the conventional method but using manual 3D modeling tools. However, this study does not help to overcome the need of repetitive modifications, the subjectivity and high production time and cost. Frillici et al., in 2008 [6], introduced the role of finite element (FE) simulation to innovate the prosthesis design process. However, only feasibility study with preliminary results was reported. In 2013, a virtual environment to design and test lower limb prosthesis has been proposed by Colombo et al [START_REF] Colombo | A digital patient for computer-aided prosthesis design[END_REF]. The VSL (Virtual Socket Laboratory) consists of a "Socket Modeling Assistance" system where the prosthetist can model the socket using manual 3D modeling tools, starting from the 3D geometry of the stump. This modeling method is assisted by FE simulation of interaction between the designed socket and the stump. The residual limb model is obtained from MRI scan and physical parameters derived from previous works are applied to the soft tissues. This FE model has been evaluated in another article published in 2016 [START_REF] Colombo | A Method to Improve Prosthesis Leg Design Based on Pressure Analysis at the Socket-Residual Limb Interface[END_REF], where the simulation results were compared to experimental ones. Note that FE simulation of the stump-socket system has been commonly performed in other studies [START_REF] Waldenfels | Computer Assisted Optimization of Prosthetic Socket Design for the Lower Limb Amputees Using 3-D Scan[END_REF] [10] [START_REF] Cagle | A finite element model to assess transtibial prosthetic sockets with elastomeric liners[END_REF] to predict the pressure distribution on the stump-socket interface and make quantitative evaluation of the socket design. Overall, these previous works have several common limitations. The first drawback relates to the very high computational cost of the used finite element method, making it difficult to perform fast dynamic simulations. Moreover, the contact formulation in finite element theory is complex and inappropriate for a fast feedback to correct the design of the socket. In fact, previous works did not mention any idea about the integration of their models into a real-time CAD-CAM platform, as efficient alternative for the conventional design and fabrication process. Consequently, in this present study, we developed a new approach to simulate the stump soft tissue deformation and stump-socket interaction using Mass-Spring System (MSS) approach and a point-to-surface contact formulation.

II. MATERIALS AND METHODS

II.1. SOFT TISSUE DEFORMATION

A. Mass-Spring System Configuration

Our approach is based on the use of the Mass-Spring System (MSS) applied on hexahedral mesh. Using the MSS model, an elastic object is represented by a set of punctual masses, called nodes. Each node is defined by its mass and its position vector, connected together by elastic connections as shown in Figure 2. A stiffness k, a damping coefficient c, an initial length L0, and the two point masses are also defined. This model corresponds to the Kelvin-Voigt model for the viscoelastic material. The distribution of the nodes in the volume of the modeled object is managed by the type of mesh used. In our present study, we used the hexahedral mesh. This choice has been justified by the fact that the hexahedral mesh gives better response to shear stresses than tetrahedral mesh [START_REF] Baudet | Integrating Tensile Parameters in Hexahedral Mass-Spring System for Simulation[END_REF]. However, this model, illustrated in Figure 3(a), has the potential to react correctly to normal compression forces but not to shearing forces, since the springs do not have shear properties. To overcome this limitation, the hexahedral element allows us to benefit from its 4 internal diagonals as well as the 12 diagonals of these faces. The use of these diagonals as elastic connections, in a well definite way, shows their ability to keep several important properties: the shear, and the non-compressibility properties of the soft tissues. Figure 3(b) shows a regular hexahedral element with the 28 possible elastic connections between its 8 nodes. Note that the hexahedral element is considered as linear elastic model defined by its Young's modulus and density.

B. Parameters Identification

After having the shape subdivision (meshing) and MSS configuration established, the parameters (the mass of each vertex, the stiffness and damping coefficient of each connection) of each item in the model are identified. To find the mass mi of a vertex i, we used the following formula presented by Deussen et al. [START_REF] Oliver Deussen | Using Simulated Annealing to Obtain Good Nodal Approximations of Deformable Bodies[END_REF] as follows:

= ∑ (1)
where H is the group of hexahedral elements to which the vertex i belongs, ρj and Vj are the density and volume of the hexahedral element j that belongs to H respectively.

The stiffness of each spring was calculated using the formulation presented in [START_REF] Golec | Hybrid 3D Mass Spring System for Soft Tissue Simulation. Modeling and Simulation[END_REF]. Reference to the hexahedral element shown in the figure 2, there are 3 types of connections (edges, internal diagonals, face diagonals) with stiffness constants ke, kc and kf respectively. The values of these parameters were identified using the following formulations: (

) 2
where E is the young modulus assigned to the hexahedral element, and le is the length of its edge considering that this is a regular cubical element. A is a parameter to be manually tuned. Golec [START_REF] Golec | Hybrid 3D Mass Spring System for Soft Tissue Simulation. Modeling and Simulation[END_REF] proved that the model shows the same results for any given value to A between 0.1 and 100. If a connection is shared by multiple hexahedral elements, its stiffness is equal to the sum of the corresponding quantities calculated for each hexahedral element according to the equation 2.

Finally, to ensure the best behavior consistency for different and combined resolutions [START_REF] Golec | Hybrid 3D Mass Spring System for Soft Tissue Simulation. Modeling and Simulation[END_REF] [21] [START_REF] Lloyd | Identification of spring parameters for deformable object simulation[END_REF] [23], the damping coefficient ci of a spring of index i and stiffness ki connecting two punctual masses mi1 and mi2 with initial length li0, are expressed in the following equation:

' = "( ) (* )+ * ), ) )- (3)

C. Volume Conservation

Soft tissues are incompressible materials, i.e. their volume remains constant under external and internal loadings. For this kind of material, the Poisson's Modulus tends to 0.5, while the system defined in equation 2 deals only with a Poisson's modulus equal to 0.25 [17] [20]. To solve this problem, a novel approach, called MSS-CS (Mass-Spring System Corrective Springs), was developed in this study. In fact, additional springs, named "Corrective Springs" were integrated to the MSS model. Thus, the proposed configuration consists of adding, in parallel with the internal diagonal springs with stiffness kc and initial length lc (resp. the face diagonal springs with stiffness kf and initial length lf), new springs having stiffness k'c and initial length l'c (resp. stiffness k'f and initial length l'f). Constitutive equations of our MSS-CS are based on the following two conditions, which are described below.

Condition 1: The initial lengths of these additional springs are greater than the initial geometric lengths as follows:

. = . = / 012ℎ / > 1 (4)
This condition allows added springs to be always extended to conserve the volume (Figure 4).

Condition 2: The stiffness of each of these springs is calculated in terms of the volume variation as follows:

67 = 8 98 - 8 - (5)
where V0 and Vc are the initial and actual volume of the hexahedral element respectively.

The corrective spring's stiffness is defined as follows:

′ ! = ′ = ; <= % (67) (6)
where C is a constant to be identified.

Supposing that the hexahedral element is regular, i.e. a cube, then:

= % = √ > (7)
Therefore:

′ ! = ′ = ; < ( > ( 8 98 - 8 - ) (8)
Finally, our proposed model is then characterized by the equations ( 2), ( 4) and ( 8), with three constants = 0.5 ; ; = 25 ; / = √2.

D. Simulation Algorithm

Let N be the total number of vertices in the MSS-CS model. Two types of forces are applied on each vertex of index i: internal forces, noted F i int, applied by the related elastic connections, and external forces, noted F i ext. Let mi and ai be the mass and acceleration of the vertex i respectively. Then, according to the second law of Newton, the equilibrium equation is defined as follows:

C DE + C %FE = G (10)
Knowing that the internal forces are spring and damping forces calculated in terms of the position Xi and its first derivative X'i respectively, and that ai is equal to the second derivative Xi", we need to apply an integration on the equation ( 10) to find the new position of the vertex i in each simulation step ∆t. We chose the implicit Euler integration method that gives a more stable simulation than the explicit one. A first order linearization is performed since the forces are non-linear according to the position variable. Applying these steps on all N vertices of the system, we get a system of 3N linear equations with 3N unknowns. The matrix of such system is always symmetric positive definite. Therefore, we can solve the system using Conjugate Gradient method.

II.2. STUMP-SOCKET INTERACTION

A. Socket modeling

To simplify the implementation and reduce the computational cost, the socket is considered as a rigid body. Usually, sockets are made of perlon stockinette hardened by resin, which makes it firmly stiff and non-deformable. Beside its triangular 3D surface (Figure 5), the rigid socket is characterized by the following physical parameters: mass Ms, moment of interia Is, center of mass CM, and center of inertia Cs. The behavior of the socket in the 3D space is described by the following second law of Newton:

C H IIIJ = K H G H IIIJ GLM N H IIIJ = O H P H IIIIJ ( 11
)
where Fs is the sum of all forces applied to the socket, τs is the sum of all torques applied to the socket, as and ωs are the acceleration and angular speed of the socket.

B. Contact modeling

Modeling contact between two 3D objects includes two numerical processing components: 1) a collision detection algorithm and 2) a formulation of physical response.

1) Collision detection algorithm

A node-to-triangle discrete collision detection algorithm was established. In this method, a collision occurs when the point penetrates into the triangle, i.e. it passes through the triangle from the positive side of the triangle's plane to the negative one (Figure 6). Thus, in each simulation time-step, the involved point is checked to verify if this point is behind the triangle's plane and whether its orthogonal projection to this plane is inside the triangle (Figure 7). In other words, we trace a ray from the point in the 3D space and verify whether it intersects with the triangle and whether the point is behind this triangle. Noting that, in our method, the ray is to be parallel to the triangle's normal vector. Practically, a collision occurs between a point Pi of the stump and a triangle Tj of the socket if these following equations are true:

Q ,S T U IIIIIIIIIJ . V S IIIJ W 0 0 W X W 1 0 W 7 W 1 0 W X + 7 W 1 (12)
where V0,j , V1,j , V2,j are the three points of the triangle Tj. u and v are the barycentric coordinates of the projection Prij of the point Pi with respect to the triangle Tj. These coordinates define the position of the projection point with respect to the three triangle's vertices as follows:

TY = + X Z - \ + 7 ( " - ) (13)
Figure 7: collision detected between the point Pi and the triangle Tj

In the equation ( 12), the first inequality means that the point is behind the triangle. The remaining three inequalities confirm that the projection point to the triangle's plane is inside the triangle. If all these inequalities are true, the point Pi and the triangle Tj are considered as "contact pair". In case if these conditions are satisfied between Pi and multiple triangles, the nearest triangle to the point is selected.

Basically, the test is applied between each point of the stump surface and each triangle of the socket. Having N points and M triangles, the complexity in this case will be O(NxM), which increases exponentially when the resolution of meshes increases. The final step in the collision detection phase is placing the point of each contact pair on a well selected collision-free position. In our method, we choose the projection of this point to the triangle as the corrected position.

2) Formulation of physical response

Physical response or reaction of the collision consists of all the forces applied to the socket and the stump due to the collision. Basically, this mechanical reaction consists of two types of forces: the normal contact force, and the friction force. According to the third law of Newton, the forces applied to the elements of the contact pair are equal in magnitude and opposite in direction.

Thus, having selected Nc contact pairs in the first phase, for each contact pair Ck (k ∈ {1, 2, 3, …., Nc}), the point Pk and the triangle Tk are each subject to a contact force of magnitude C ]DE and a friction of magnitude C !^ . Thus, the force applied by the point of the stump surface to the contact surface of the socket is equal to the sum of the inner forces C J DE,_ ` applied to the point due to the material mechanical properties along with the external forces C J %FE,_ `. This total force applied to the point is noted C J _ ` and expressed in the following equation:

C J _ `= C J DE,_ `+ C J %FE,_ ` (14)
Talking about a normal contact force, we will consider the forces perpendicular to the triangle and disregard the forces aligned with it. Therefore, the magnitude of the contact force is defined as follows:

C ]DE = -C J _ `. V a IIIIIJ (15)
The minus sign in the equation ( 15) is added to obtain a positive value, since the force C J _ ` applied to the point should be pushing the triangle, having opposite direction of the triangle normal vector V a IIIIIJ (Figure 8(a)).

c C J ]DE,_ `= C ]DE . V a IIIIIJ C J ]DE,b `= -C ]DE . V a IIIIIJ (16)
Once a contact force exists between the contact pair elements, the friction force is computed. This force basically depends on the relative velocity between the two elements. Let IJ be the relative velocity vector of the point Pk moving with a velocity IJ _ ` with respect to the triangle Tk moving with a velocity IJ b `:

IJ = IJ _ `- IJ b ` (17)
If this relative velocity is null, the friction is ignored, otherwise, its magnitude is defined as follows:

C !^ = d C ]DE ( 18
)
where µ is the friction coefficient assigned to the stump-socket interface. The direction of the friction vector C J !^ ,_ ` applied to the point is opposite to its relative velocity vector. The friction vector C J !^ ,b ` applied to the triangle is in the other direction as follows:

Q C J !^ ,_ `= -C !^ . 8 I IJ è8 I IJ `e C J !^ ,b `= C ]DE . 8 I IJ è8 I IJ `e (19)
The forces generated by the collision between the elements of each contact pair are thus calculated. To introduce them into the simulation, we apply the following steps:

Step 1: For each point Pk of the contact pair Ck and existing on the surface of the stump model, the contact force C J ]DE,_ ` and the friction C J !^ ,_ ` are added to the external forces in equation [START_REF] Surapureddy | Predicting Pressure Distribution Between Transfemoral Prosthetic Socket and Residual Limb Using Finite Element Analysis[END_REF].

Step 2: Regarding the socket, two vectors are necessary to simulate its behavior in the 3D space: the total force C H IIIJ and the total torque N H IIIJ, as shown in the equation [START_REF] Cagle | A finite element model to assess transtibial prosthetic sockets with elastomeric liners[END_REF]. The total force is simply found by summation of the forces applied to the triangles Tk, k ∈ {1, 2, 3, …., Nc}, along with any other external force C J f,%FE applied to the socket as follows:

C H IIIJ = C J f,%FE + ∑ (C J ]DE,b `+ C J !^ ,b `) g h i (20)
For the torque, let first consider for each contact pair Ck a torque τk that is generated by the interaction of Pk and Tk and defined as follows:

N IIIJ = (C J ]DE,b `+ C J !^ ,b `). T ; j IIIIIIIIJ (21)
where CI is the center of inertia of the socket. The total torque is then the summation of all these particular torques with any other external torque τs,ext that may exist:

N H IIIJ = N H,%FE IIIIIIIIIJ + ∑ N IIIJ g h i (22)

II.3. SIMULATION ENVIRONMENT DEVELOPMENT

A simulation environment was developed by using C++ programming language. Graphical user interface was designed and developed with "Visual Studio / Windows Form Application". OpenGL library [START_REF] Architecture | OpenGL Programming Guide: The Official Guide to Learning OpenGL[END_REF] was used for the 3D visualization. A virtual transtibial residual limb was developed using the available CAD model library (free3D.com) (Figure 9(a)). A hexahedral mesh of the model (Figure 9(b)) was generated using the open-source tool IA-FEMesh [START_REF] Grosland | IA-FEMesh: An open-source, interactive, multiblock approach to musculoskeletal finite element model development[END_REF], from which the code automatically generates the related MSS-CS model (Figure 9(c)) and calculate its parameters from the predefined physical parameters of the soft tissues. A Young's modulus of 50KPa [START_REF] Colombo | A Method to Improve Prosthesis Leg Design Based on Pressure Analysis at the Socket-Residual Limb Interface[END_REF], [START_REF] Kot | Elastic Modulus of Muscle and Tendon with Shear Wave Ultrasound Elastography: Variations with Different Technical Settings[END_REF] was used. The density for soft tissue was set up as 1060Kg/m3 [START_REF] Urbancheka | Specific Force Deficit in Skeletal Muscles of Old Rats Is Partially Explained by the Existence of Denervated Muscle Fibers[END_REF]. A geometrical process was developed to create the socket model (Figure 5). Stump geometry was used to extrapolate the socket geometry with a thickness of 10 mm. Finally, the soft tissue deformation and stump-socket interaction algorithms were integrated into the developed environment.

II.4. ACCURACY ANALYSIS

Two simulation cases were developed to evaluate the accuracy of our soft tissue deformation and stump-socket interaction algorithms. We started first by a simple cubical elastic object subject to normal compression due to the weight of a rigid box (Figure 10). We repeated the simulation with different weight values, and we recorded, for each one, the displacement and contact pressure on elastic cube's upper surface as well as its volume variation, and we compared them to analytical results. We also implemented the MSS-CF model developed by Golec [START_REF] Golec | Hybrid 3D Mass Spring System for Soft Tissue Simulation. Modeling and Simulation[END_REF] for comparison. The dimensions of the elastic cube are 50mm. Young' modulus of 50KPa and a density of 1060Kg/m 3 were used. In the second case, we performed a dynamic simulation for the socket donning process (Figure 11), and we mapped the contact pressure distribution on the stump-socket interface. The same physical properties (Young' modulus of 50KPa and a density of 1060Kg/m 3 ) were assigned to the stump soft tissues. The socket moves according to a displacement constraint. Finally, a FE model of the stump-socket interaction was developed using Abaqus for the comparison with our MSS-CS model. A C++ code was written to convert the MSS-CS model to an Abaqus model, so that the same meshes are used in both simulations.

A convergence analysis was conducted to determine the optimal mesh configuration (7393 nodes of which 778 are related to the bones and are fixed, and 6270 hexahedral elements (Figure 12)) of the involved models. The same physical properties (Young' modulus of 50KPa and a density of 1060Kg/m 3 ) were assigned to the stump soft tissues. A Poisson ratio of 0.45 was used [START_REF] Chawla | Mechanical properties of soft tissues in the human chest, abdomen and upper extremities[END_REF]. All parameters used in the FE simulation are given in table 1. The simulation was performed on an HP workstation with 16GB RAM and Intel® Core™ i7-7700HQ 2.8GHz CPU. Obtained results show that the MSS-CF is very accurate regarding volume conservation but does not reflect the real linear elastic behavior with error reaching 20% for a displacement of 5%, compared to the analytical results. In other hand, our MSS-CS approach shows more satisfying elastic behavior with errors that don't exceed 5% for up to 50% displacement. Even though the Figure 15 shows that the MSS-CS is not as efficient as the MSS-CF regarding the volume conservation, we consider it acceptable with error equal to 3.1% for 10% displacement. Figure 13 confirms that our contact model has a very good precision level and the contact pressures on the interface between the elastic and rigid bodies are very close to the analytical solutions. These simulations were performed in real time and the average calculation time of each time step is equal to 15.29µs.

Parameters

B. Stump-socket interaction simulation

The results shown in the Figure 16 were captured at the end of the socket donning process simulation by using Abaqus software and our MSS-CS approach. Overall, the two pressure patterns have a satisfactory correspondence in a qualitative manner. To make a quantitative evaluation, two profiles were tracked on the front vertical line as show in figure 17(a), and those existing on the circumference of a circle as shown in figure 17(b). The first points set has been chosen since the line cut a critical point in the stump, the patellar tendon, and the points are represented by their height, i.e. z coordinate. The second point set has been chosen in a way to evaluate the pressure on tabular flares as well as medial, lateral and posterior pressure zones. We represented these points by their angular coordinate with respect to the vertical axis passing through the middle of the stump. The contact pressure distributions for these 2 points set are given in Figures 18 and19. We observed a very good fitting pattern between our MSS-CS approach and Abaqus simulation. However, a relative error ranging from 25% to 50% was found at the peak contact pressure. In addition, the simulation time using Abaqus is around 5 hours 25 minutes and 10 second, while our simulation time is around 69 seconds. Although the simulation time is considerably reduced, the real time simulation is not reached because of a time step reduction for assuring simulation stability. The time step was set to 1ms while it's needed calculation time was equal to 69ms almost divided equally between the stump deformation calculation and the contact calculation (30ms and 33.9ms respectively, and the rest is for rendering and other calculations).

IV. DISCUSSION

Prosthesis design has become an attractive field in the recent years thank to the progress of 3D printing capacity and biomechanical researches. However, traditional manual design and fitting process has been commonly performed in many hospitals and clinics over the world. This time-consuming and inefficient process leads to a high medical cost and requires significantly human resources [START_REF] Lorenzelli | Socketmaster: Integrated Sensors System for the Optimised Design of Prosthetic Socket for above Knee Amputees[END_REF]. Research studies have been conducted to improve this process. Virtual prototyping studies have been proposed. Numerical modeling based on finite element methods has been usually used to provide objective and quantitative indicators for evaluating and optimizing the design process [START_REF] Colombo | A digital patient for computer-aided prosthesis design[END_REF] [8] [START_REF] Waldenfels | Computer Assisted Optimization of Prosthetic Socket Design for the Lower Limb Amputees Using 3-D Scan[END_REF] [10] [START_REF] Cagle | A finite element model to assess transtibial prosthetic sockets with elastomeric liners[END_REF]. Even if FE-based approach leads to a detailed level of feedback, simulation outcomes require significantly computational resources and complex modeling skills. To deal with this problem, in this present study, a novel approach was proposed to estimate the soft tissue deformation and stumpsocket interaction during the design process. A Mass-Spring System coupled with corrective springs (MSS-CS) was developed and evaluated according to a finite element model. Thanks to its simplicity and low computational cost, MSS-based approach attracts the interest of many researchers in different domains to model the elastic objects like in video games and computer graphics [START_REF] Bender | Position-Based Simulation of Continuous Materials[END_REF] [13] [START_REF] Liu | Quasi-newton methods for real-time simulation of hyperelastic materials[END_REF] as well as in surgical simulators [15] [16] where researchers work on developing virtual training environment. Unfortunately, the literature analysis shows that this approach is not yet introduced in the prosthetic domain. This approach was adopted in this study to provide a more efficient computational approach for soft tissue deformation. Regarding the MSS-based studies, many researchers proposed to add a new type of forces to each vertex, called "Corrective Forces", which are functions of Poisson's modulus and volume variation. Inspired by previous works, Golec [START_REF] Golec | Hybrid 3D Mass Spring System for Soft Tissue Simulation. Modeling and Simulation[END_REF] proposed a formulation of these volume correction forces, giving a new model called MSS-CF (Mass-Spring System Corrective Forces). Golec's results showed the validity of these additional forces to perfectly correct the volume but did not show their impact on the elastic behavior. In this study, we already showed that MSS-CF model would negatively impact the elastic response of the model to an applied stress. We noticed that the problem resides in the direction of the resistive forces to an applied stress when the material is compressed. By comparing to this approach, our MSS-CS showed more consistent and accurate results (shown in Figure 14).

In addition, a formulation of stump-socket contact was also proposed in this study. Collision detection (CD) is the computational problem of detecting the geometric intersection of objects in the 3D space. Many CD methods for rigid and deformable bodies have been proposed and used in different applications such as video games [START_REF] Ericson | Real-time collision detection[END_REF], robotics [START_REF] Balan | Real-time 3D collision avoidance method for safe human and robot coexistence[END_REF], surgical simulators [START_REF] Courtecuisse | Gpu-based real-time soft tissue deformation with cutting and haptic feedback[END_REF] [27] [START_REF] Fazioli | Implementation of a soft-rigid collision detection algorithm in an open-source engine for surgical realistic simulation[END_REF]. Generally speaking, there are two basic types of collision detection algorithms: discrete collision detection (DCD) and continuous collision detection (CCD). DCD detects the collision between the two objects in a discrete time, so the exact collision time is missed, and the detection is made after the penetration of an object into the other. The collided items then are moved to corrected collision-free positions. In other hand, CCD searches for the exact collision moment in the interval bounded by two time-steps. It has been proposed for collision between rigid bodies [START_REF] Redon | Fast continuous collision detection between rigid bodies[END_REF] [30] and then used for detection of self-collision of deformable objects [31] [32]. CCD is more robust than DCD but requires very high computational cost, while the latter is much simpler and cheaper. Usually, it is extremely difficult to compare the different collision detection methods and choose one of them since they are very sensitive to specific scenario [START_REF] Trenkel | A benchmarking suite for static collision detection algorithms[END_REF]. In our case, a high simulation speed is required, that's why the DCD approach was used. More precisely, an algorithm derived from ray-casting-based method proposed by Fukuhara et al. in 2014 [START_REF] Fukuhara | Proposition and evaluation of a collision detection method for real time surgery simulation of opening a brain fissure[END_REF] was developed with few modifications (inverse way of collision (deformable object to rigid one) and ignorance of false contact pairs) to make it compatible with our application. By comparing with a detailed finite element contact simulation, our formulation gives a similar contact behavior with a perfect qualitative matching. However, quantitative deviation of contact pressure at peak values remains high. Note that this quantitative difference may be due to the different contact formulations between Abaqus (surface-to-surface contact) and our approach (point-to-surface contact). Thus, further studies will be investigated to improve our contact formulation. Even if our MSS-CS is more efficient than Abaqus simulation in term of computational time, our computation cost remains high. A possible simple solution to reduce the calculation time is to automatically assign to each point, at the beginning of the simulation, the potential contact triangle, so the algorithm will test each point with the assigned triangles only. For instance, we remarked that the points of the posterior aspect of the stump never collide with the triangles in front of the socket. In addition, all triangles of the external surface of the socket can be discarded. Thus, further studies will be investigated to improve this point.

Another limitation of our study relates to the use of linear elastic material to model soft tissue. However, this simplification is considered acceptable in our application because the deformations of stump soft tissues are small and within the linear region.

Previous work [START_REF] Colombo | Physically based modelling and simulation to innovate socket design[END_REF] experimented both linear and hyper-elastic models, and results show similar pressure distribution for both cases but a considerably higher simulation speed in case of linear model.

V. CONCLUSIONS AND PERSPECTIVES

A novel approach was proposed to model soft tissue deformation and stump-socket interaction toward a Computer-Aided Design (CAD) system for lower limb prostheses. A Mass-Spring System Corrective Spring (MSS-CS) model was developed and evaluated. Nonlinear "Corrective" springs were introduced to respect the non-compressibility of the soft tissues. A node-to-surface contact formulation was also integrated and evaluated. Obtained results confirmed the efficiency and accuracy of the proposed method comparing to a complex and detailed finite element model. As perspectives, this present approach will be used to simulate the stump-socket interaction in real case studies. Thus, the evaluation of the proposed approach will be enhanced with more datasets. Finally, this modeling approach will be integrated into a CAD-CAM platform providing a fast feedback to the prosthetist during the socket manufacturing process.

Figure 1 :

 1 Figure 1: Current practice in the fabrication of the prosthetic socket

Figure 2 :

 2 Figure 2: Schematic representation of the connection between 2 punctual masses.

Figure 3 :

 3 Figure 3: (a) Hexahedral MSS model and (b) hexahedral element with the 28 possible connections: 12 edge connections in blue, 12 faces diagonals in red, and 4 internal diagonals in green

Figure 4 :

 4 Figure 4: Cubical elastic object subject to normal compression pressure: (a) Compressible material modeled using MSS without corrective springs. (b)Uncompressible material modeled using MSS with corrective springs performing horizontal expansion to conserve its

Figure 5 :

 5 Figure 5: 3D triangular surface of a transtibial prosthetic socket (a) shaded view, (b) wireframe view

Figure 6 :

 6 Figure 6: Schematic representation of the collision taking place between the point Pi and the triangle Tj

Figure 8 :

 8 Figure 8: (a) Forces applied to the point are toward the triangle, so the contact force exists. (b) The forces applied to the point are not toward the triangle, so no contact force is considered In the opposite case, if this equation gives a negative value, this would mean that the point is "pulled" instead of being "pushed" by the internal and external forces (Figure 8(b)). In this case, we discard it and set the contact force to null. Finally, the vectors of contact forces C J ]DE,_ ` and C J ]DE,b ` applied to the point Pk and the triangle Tk respectively are defined as follows:

Figure 9 :

 9 Figure 9: Stump model: (a) original 3D model, (b) hexahedral mesh, (c) MSS-CS model

Figure 10 :

 10 Figure 10: Compression of the pink elastic cube due to the weight of the blue rigid box: (a) initial state, (b) final state

Figure 11 :

 11 Figure 11: Schematic simulation of the socket donning process: (a) initial state, (b) final state

Figure 12 :

 12 Figure 12: Optimal meshed stump model used in socket donning simulation

Figure 13 :Figure 14 :Figure 15 :

 131415 Figure 13: Contact pressure on the upper surface of the elastic cube in function of the rigid box weight

Figure 16 :

 16 Figure 16: Contact pressure distribution on the stump-socket interface from both MSS-CS and FE simulations at the end of the socket donning process

Figure 17 :Figure 18 :

 1718 Figure 17: the two points sets chosen for quantitative evaluation

Figure 19 :

 19 Figure 19: Contact pressure distribution on the horizontal points set

Table 1 :

 1 Parameters of the FE simulation performed by Abaqus

	Values

ACKNOWLEDGMENT

The authors would like to acknowledge the support and funding of "Beirut Research and Innovation Center" (BRIC), Beirut, Lebanon.