A Perrier

V Luboz

M Bucki

F Cannard

N Vuillerme
email: nicolas.vuillerme@agim.eu

Y Payan

Chapter 18: Biomechanical modelling of the foot

Keywords: foot model, biomechanical analysis, heterogeneous soft tissues, clinical applications

The foot is a key structure during gait as it adapts to the ground geometry to insure balance and it stores energy to ease the next step. Because this functionality can be affected, by morphological issues or by aging or disease like diabetes, modelling the foot could help in understanding its correct behavior and the way to avoid or solve pathological behaviors. In this chapter, after presenting the foot anatomy and functionality, we propose a summary of the main foot models published in the literature. We also introduce our biomechanical model, using a Finite element mesh representing several soft tissue layers surrounding rigid bones and cables simulating the ligaments. Two clinical applications, namely ankle arthrodesis and foot ulcer prevention, are also presented.

Introduction: clinical context

Understanding the foot anatomy is essential to analyze its functionalities, its pathologies, and the way they can be prevented and treated with orthosis or surgery for instance. These complementary fields of research have yielded to several significant breakthroughs helping understanding the foot, using models adapted to each clinical applications. The first field concerns the physiology of the foot and more especially the way the foot behaves in conjunction with other structures (for example the leg or the ground) and their modifications to monitor the adaptation of the foot to various conditions. Furthermore, among the different clinical fields, we can mention three other fields studying the foot on the pathological side: diabetes, orthosis, and surgery. Diabetes is a pathology affecting the foot through a lack of sensitivity due to the neuropathy it implies [START_REF] Mueller | evaluation, and treatment of the neuropathic foot[END_REF]. As a consequence, diabetic patients have trouble feeling pain in their foot and can consequently undergo pressure ulcer due to high tissue deformations (because of ischemia or cell destruction). Furthermore, diabetes can create joint lesions because of the neuropathy. The third field of research to understand the foot is devoted to orthoses. It concerns studies investigating the influence of a device on the pressure load around the foot, the geometrical changes applied on the bones, or the behavior of the foot after an impact [START_REF] Lemmon | The effect of insoles in therapeutic footwear-a FE approach[END_REF]. The orthosis influence is analyzed based on its shape and/or the materials it uses. Finally, the fourth field of research aims at improving foot surgeries. It concerns mainly the understanding of a surgical gesture, and the influence of an action on the foot structure and functions. It can also investigate the coupling between a foot joint and prosthesis.

Anatomical description of the foot

The foot anatomy being complex, it is essential to comprehend it before moving to one of the four fields mentioned above. The human foot is composed of 28 bones and 33 joints (including five main joint complexes: the ankle, the subtalar joint, the tarsal transverse joint, the tarso-metatarsal joint, and the metatarsophalangeal joint). Several bone regionalization schemes have been proposed, one of the most popular is described in figure 18.1 and divides the foot in four regions: the back foot (talus and calcaneus), the medio foot (cuboid, navicular, and the three cuneiform bones), the metatarsus (from M1 to M5), and the toes (from hallux to quintus). The foot has three extrinsic muscle regions located on the tibia and fibula: the anterior region (with the tibial anterior, the hallux long extensor, and the toes' long extensor muscles), the lateral region (with the long fibula and the short fibula muscles), and the posterior region (with the tibial posterior, the hallux long flexor, the toe long flexor muscles, and the sural triceps). Four muscle regions are intrinsic: the dorsal region (with the short extensor muscles of the hallux and of the toes), the medial region (with the abductor and adductor of the hallux, and the short flexor of the hallux), the lateral region (with the abductor of the quintus, the short flexor of the quintus, and the opposing muscle of the quintus), and the median region (with the short flexor of the toes, the plantar square, the interosseous dorsal muscles, the interosseous plantar muscles, and the lumbrical muscles). These muscles allow a wide range of foot motions: dorsiflexion, plantar flexion, pronation, supination, abduction, adduction, inversion and eversion.

FIGURE 18.1 HERE

The foot is a poly-articular system, surrounded by soft tissues, with several joint chains in parallel and serialized in the four previously described regions (figure 18.1). It is a key structure during gait as it adapts to the ground geometry to insure balance and it stores energy to ease the next step in a non-pathological configuration. The foot objective is to transmit the best force vector to the rest of the lower leg whatever the contact surface type (texture, form, rigidity…) or the aim of the gesture (gait, turn, jump, get up, run…). To achieve these goals, the foot must be able to adapt to the shape of the ground, to rigidify itself (good transmission of the vector fields from the lower leg or the external forces), to be thrifty for an optimized gesture (long walk, swim, running….). The serialized organization of the foot allows force transmission from rear foot to front foot during the first contact of the foot on the ground, and then from front foot to rear foot during the digitigrade phase. The parallel chains are active in the rear to front axis and also in the medio lateral axis. The plantar aponeurosis (or plantar fascia) is one of these soft tissues. It helps preventing the stretching of the foot structures. It can be stretched about 10 percent, and has extra length under each metatarsal head in connection with the first phalange of each toe. It has a joint action with the trans-metatarsal ligament to prevent lateral foot stretching, in conjunction with the interosseous muscles. Other ligaments link each bone with its neighbors to create the joints using a complex mechanism mixing bone contacts and cable like connections between ligaments. The poly-articular chain system allows the foot to modulate its length, width and height during loading. If for some reasons, the chains are rigid, the strains are propagated to the foot extremities instead of the center of the foot, which usually concentrates the constraints. The foot is also composed of fat tissue, mainly to serve as a cushion for the other structures. In particular, below the calcaneus, the fat pad is composed of several nodules of high density fat to compensate the impact of the whole body below the heel at the first contact with the ground. Finally, the foot soft tissues are embedded within the protective skin layer. Most of the gait theory is based on conclusions considering only the bones (position and contacts). To the best of our knowledge, soft tissues have not been integrated in early works. However, to be useful and usable in clinical practice, a foot model therefore needs to integrate this component.

Foot models in the literature

Several reviews have already summarized a list of foot models over the past few years. Among them, that recently published by [START_REF] Telfer | What Has Finite Element Analysis Taught Us about Diabetic Foot Disease and Its Management ? A Systematic Review[END_REF] provides a comprehensive review on diabetic foot modelling and investigates several aspects such as the bone geometry, the internal constraints, the soft tissue material properties, the skin surface pressures and their spreading, and the interactions between the foot and the device around the foot. Another review, from [START_REF] Wang | Computational Models of the Foot and Ankle for Pathomechanics and Clinical Applications: A Review[END_REF], focuses on the shoe and the insole and aims at designing the most adapted device for a given foot. Recently, Perrier (2016) further draws a complete survey of most of the foot models presented in the literature, table 18.1. We propose to summarize this review below with few key articles.

TABLE 18.1 HERE

The first biomechanical model was a 2D model introduced by Lemmon and Cavanagh [START_REF] Lemmon | The effect of insoles in therapeutic footwear-a FE approach[END_REF]. Using baropodometry analyzing the plantar pressures, these authors validated their model results against clinical measurements. This work stemmed research studies in the field of diabetes and orthopedic surgery and is still a reference in terms of validation. It showed that insole thickness has an influence on pressure peaks under metatarsal heads. It used a linear constitutive law to model the bone material with a Young's modulus of 10 GPa and a Poisson's ratio of 0.34. A Mooney-Rivlin constitutive law was used to model the soft tissues as a hyperelastic material. In vivo load-deformation data were collected on five subjects using an ultrasound device in conjunction with a load cell. It lead to a simplified Finite Element model with coefficient values of C 10 = 85,550 N.m -2 , C 01 = -58,400 N.m -2 , C 20 = 38,920 N.m -2 , C 11 = -23,100 N.m -2 , C 02 = 8,484 N.m -2 , D 1 = 0.4370E05 N.m -2 , D 2 = 0.6811E-06 N.m -2 . Another 2D model was introduced by (Gefen, 2000) and validated on six phases of gait, using plantar pressure measurements, in a quasi-static set up. This article stemmed foot functional analysis based on biomechanical models. It was followed in [START_REF] Gefen | Plantar soft tissue loading under the medial metatarsals in the standing diabetic foot[END_REF] by a more complete model demonstrating that pressures measured at skin surface lead to a strain increase in deep tissues. From this conclusion, it was hypothesized that pressure ulcers may originate deep in the tissues and not only at the foot surface. In this article, the bones and cartillages were modelled as linear, elastic and isotropic materials, while ligaments, fascia and the soft tissue fat pad were considered as being non-linear materials using an Ogden constitutive law. Young's modulus and Poisson ratio for bone were taken as being 7,300 MPa and 0.3, respectively, and for cartilage as 10 MPa and 0.4, respectively. The Poisson ratios were taken as 0.4 for the ligaments and fascia and 0.49 for the soft tissue pad. The elastic properties of these soft tissues were derived from stress-strain data published by [START_REF] Nakamura | An Analysis of Soft Tissue Loading in the foot -a Preliminary Report[END_REF] but were not given in the article. One of the first 3D biomechanical models was introduced by [START_REF] Chen | Effects of total contact insoles on the plantar stress redistribution: a finite element analysis[END_REF] and showed that an orthosis with an optimal contact with the foot decreases the skin surface pressures and spread them more evenly below the foot. From this conclusion, it seems that a pressure ulcer may appear on a foot depending on the pressure field applied on its surface. Using a 3D Finite Element model therefore appears to be appropriate to estimate the consequences of various orthoses on the foot. In this paper, the foot bone and soft tissue were assumed to be linear, elastic solids. Their Young's moduli and Poisson ration were respectively 10 GPa and 0.34 for the bones, 1.15 MPa and 0.49 for the soft tissues, and 10 MPa and 0.4 for the cartilage. Five cables were used to model the plantar aponeurosis and three cables for the short plantar ligaments. Their stiffness was 11.5 MPa. An improvement of this 3D model was presented in (Chen et al., 2010), with the integration of active muscles, almost all foot ligaments, and soft tissues modelled with the Mooney Rivlin constitutive law defined in [START_REF] Lemmon | The effect of insoles in therapeutic footwear-a FE approach[END_REF] with the same parameters. It showed that internal strains and pressures below the metatarsal heads are directly influenced by the force given by the triceps. Another 3D model has been proposed by (Cheung et al., 2005). It showed that (1) the shape of an orthosis is more important than its materials when it comes to strain distribution and (2) a higher soft tissue stiffness implies higher skin surface pressures. Bones and the soft tissues were also represented as Finite Element models. Except for the soft tissue, all other tissues were idealized as homogeneous, isotropic and linearly elastic materials. The Young's modulus and Poisson's ratio for the bony structures, the cartilage, ligaments, and plantar fascia were selected from the literature (Gefen, 2000). The soft tissue properties of the whole foot were also inspired by the heel pad properties given in [START_REF] Lemmon | The effect of insoles in therapeutic footwear-a FE approach[END_REF] as a nonlinearly elastic material represented by a Mooney-Rivlin constitutive law. This model was used by [START_REF] Wong | Functional Restoration and Risk of Non-union of the First Metatarsocuneiform Arthrodesis for Hallux Valgus: A Finite Element Approach[END_REF] for clinical applications in orthopedic surgery and more specifically in the arthrodesis of the first tarso-metatarsal joint and its consequence on gait. Although the models presented above tend to integrate complex structures and allow studying various clinical fields, they still miss some of the anatomy of the foot (especially no model represents the foot material heterogeneity as the foot is composed of different layers of soft tissues: skin, fat, muscles), they are not patient specific (as they are all modelled from one patient's data, or from an idealized morphology), and they cannot be used clinically as they require long computations (more than one hour). The main limitation concerns the realism of the model: it must be accurate if useful results are needed and therefore the model should not simplify the morphology of the patient and the way the mechanical properties are estimated. A nonlinear heterogeneous model, i.e. with several layers of soft tissues to represent the skin, fat and muscles, consequently seems to be important to simulate the foot accurately. Nevertheless, the key aspect of a good biomechanical model is that it should be adapted to the function (gait, balance, pressure distribution…) it wants to simulate, using an appropriate mechanical hypothesis and an appropriate level of details for the morphology.

Biomechanical foot model for arthrodesis simulation and foot pressure ulcer prevention

In an attempt to improve the previous foot models in a more realistic manner and with an impact on several applications, we have proposed to generate a detailed biomechanical model to study arthrodesis and foot pressure ulcer prevention. A series of medical image acquisitions has been performed on a healthy volunteer (no trauma, surgery or specific pathology of the foot) to obtain morphological data and create his biomechanical foot model. Two imaging modalities were used for this individual (31 years old, 173cm, 84 kg): a CT scan with slices every 0.6 mm, and an MRI scan with slices every 5mm. The scans were segmented to extract the 3D surfaces of the various foot structures: 28 separated foot bones, tibia and fibula, and four layers of soft tissues, namely the skin, the heel fat pad, the fat for the rest of the foot, and the muscles. To allow fast enough computation time, we chose to represent the bones as rigid bodies to concentrate on the cinematic of those structures, especially before and after arthrodesis. On the other hand, since we want to measure the strains inside the soft tissues, in order to study foot pressure ulcer, the soft tissues have been modelled as a Finite Element mesh. Rigid bodies and Finite Element mesh have been coupled semi implicitly in the ArtiSynth open source framework [START_REF] Lloyd | ArtiSynth: a fast interactive biomechanical modeling toolkit combining multibody and finite element simulation[END_REF]. To be realistic in physiological and dynamic points of view, the mass of the different structures have been estimated from (De [START_REF] Leva | Adjustments to Zatsiorsky-Seluyanov's segment inertia parameters[END_REF], with a foot mass of 1.37 % of the subject mass. A density of 1 850 kg.m -3 was set for the bones, leading to a whole bony mass of 0.443 Kg. The density of the soft tissues was set to 1 200 kg.m -3 to reach a soft tissue mass of 1.13 Kg. The muscle and ligament insertions were defined from the medical imaging data and from the literature, Figure 18.2. Cables were used to model the 210 ligaments of the foot included in our model. These cables constrain the elongation of the structures with a linear elongation stiffness of 395 MPa while the structures are free to move towards each other since the compression stiffness is set to 0 MPa. Furthermore, to simulate the coaptation of the bones, the rest length of the cables was set to 90 % of the initial length of the model since the foot was acquired while the subject was lying down with an initial length that is not representative of the position of the bones when standing. Overall, the joints are modelled as the contacts between the bones and are constrained by the ligaments. Plantar aponeurosis is also simulated in our model. It is composed of five cables connecting the calcaneus to each toe first phalange. Transversal cable connections constrain the aponeurosis in the three directions.

FIGURE 18.2 HERE

Fifteen extrinsic and intrinsic muscles are integrated in our foot model. They are simulated as cables using a Hill constitutive law. Their path has been determined on the MRI scan. Each muscle can be activated independently, i.e. its length can be shortened simulating muscle contraction. Since the foot is relaxed during the MRI acquisition with the patient lying down, its muscles are not contracted. Muscle length at equilibrium during stance had therefore to be defined as a ratio of these muscles' MRI lengths. This was defined manually for each muscle and led to length ratios between 0.73 and 0.91. The maximum force developed by a muscle was set by computing the displacement of half the body weight when the foot is in charge. For example, to raise the heel by contracting the two bi articular muscles of the triceps, a force of 1,800 N is needed. This value corroborates the 1,620 N given in (Chen et al., 2010) in a quasi-static setup. The soft tissue Finite Element mesh was generated by a software from Texisense company and based on [START_REF] Lobos | Techniques for the generation of 3D Finite Element Meshes of human organs[END_REF]. It uses the surfaces segmented from the CT scan to define the different domains in the foot and to generate the continuum. This meshing algorithm tends to maximize the number of hexahedrons in the generated mesh to avoid the locking mechanism linked to tetrahedrons in quasi incompressible materials. The final mesh is composed of 142,060 elements (18,784 hexahedrons, 38,774 wedges, 56,655 tetrahedrons, and 27,847 prisms) with a discretization of 3 mm. It is composed of four domains: the skin, the heel fat pad, the rest of the fat, and the muscles, Figure 18.3. The bones are modelled as rigid bodies and no Finite Elements are included in them.

FIGURE 18.3 HERE

To achieve a nonlinear hyper-elastic behavior, a Neo Hookean constitutive law is used in the Finite Element model, defined by a strain energy density function W such as: W=C 10 (I 1 -3)+(J-1) 2 /D Where I 1 is the first invariant of the left Cauchy-Green deformation tensor, C 10 is a material parameter, J is the determinant of the deformation gradient F and D is a material incompressibility parameter. The Neo Hookean model was chosen because of its simplicity to apprehend its parameters. It indeed uses only two material parameters that can be compared to known material parameters, with D linked to the Poisson ratio ν [D=(1-2ν)/ C 10 )] and C 10 that can provide an estimation of the equivalent Young's modulus E at small strains [E≈6* C 10 ]. Furthermore, Neo Hookean models are good representations of the reality in compression, or at least similar to Mooney Rivlin materials, with fewer parameters and more stability. In addition, since the foot is mainly submitted to compression, this range of deformations is well modelled using a Neo Hookean constitutive law. Each soft tissue layer (or domain) has its own material parameters to simulate different behaviors and derives from earlier work [START_REF] Luboz | Influence of the Calcaneus Shape on the Risk of Posterior Heel Ulcer Using 3D Patient-Specific Biomechanical Modeling[END_REF]. It leads to the following four material properties:

-For the skin, a Young's modulus of 200 kPa and a Poisson ration of 0.485, -For the heel fat pad, a Young's modulus of 100 kPa and a Poisson ration of 0.499, -For the rest of the fat;, a Young's modulus of 30 kPa and a Poisson ration of 0.49, -For the muscles, a Young's modulus of 60 kPa and a Poisson ration of 0.495.

This model was used for two different applications: ankle arthrodesis and pressure ulcer prevention. The first one consisted in investigating ankle mobility during muscle contraction with and without arthrodesis. Dorsal and plantar flexions of the ankle were simulated by activating the foot muscles of our model while assessing the mobility of the talus bone and the foot. Dorsal flexion activation is performed in three movements:

(1) activation of the toe and hallux extensor muscles alone, (2) contraction of the tibialis anterior muscle alone, and (3) combined contraction of the extensors and the tibialis anterior muscles. Plantar flexion activation is divided into four movements: (1) activation of toes and hallux flexor muscles alone, (2) contraction of the tibialis posterior muscle alone, (3) contraction of the triceps, and (4) combined contraction of all those ankle flexor muscles. During those movements, the 3D angle formed by the tibia center, the talus center, and the second metatarsal head center gives an estimation of the foot opening angle. The simulated arthrodesis set the angle between the tibia and talus to 117°; it therefore limits the foot opening angle. This angle range of motion is accepted to be for healthy subjects between 60° and 80°. The goal of this study was to estimate the arthrodesis influence by comparing the foot opening angle with and without arthrodesis.

FIGURE 18.4 HERE

The second simulation consisted in investigating pressure ulcer risk in the foot and improving prevention. Foot pressure ulcers are a common complication of diabetes because of patient's lack of sensitivity due to neuropathy. They mainly result from excessive skin surface pressure intensity (leading to internal strains above 50 % for about 10 minutes) or prolonged compression (leading to internal strains above 20 % for about two hours) [START_REF] Loerakker | The effects of deformation, ischaemia and reperfusion on the development of muscle damage during prolonged loading[END_REF]. Monitoring the internal strains using a biomechanical model combined with skin surface pressure measurements could consequently improve the estimation of pressure ulcer risk and help patient daily prevention. We therefore used our biomechanical foot model to simulate a pressure load applied below the foot during static unipodal stance. The plantar pressure pattern was measured using a commercially available pressure sensor (Zebris platform, http://www.zebris.de/) under the right foot of the subject used to build our model. The maximal pressure appears below the calcaneus with a peak value of 14.5 N.cm -2 . During the simulations, the tibia was fixed while the fibula and the foot were free. The simulation is divided in two phases:

(1) reaching equilibrium of the foot while tendons, ligaments and muscles of the model tend to recover their optimal length (defined by the ratio of their initial length) and thus generate pre-stresses in the FE continuum, and (2) gradual application of the pressure pattern below the foot. Once the full pressure is applied and equilibrium is reached, peak internal Von Mises strains have been recorded as well as the cluster volumes of soft tissues above the 20 % and 50 % strain thresholds and defining the risk of pressure ulcers. A cluster is defined as the contiguous elements having a strain peak above one of those two thresholds [START_REF] Bucki | Clinical workflow for personalized foot pressure ulcer prevention[END_REF]. The aim of this study was to estimate the amount of tissues at risk in this configuration and to determine if our model could help in pressure ulcer prevention.

Results and discussions

Arthrodesis simulation

The simulation of the ankle arthrodesis allows to explore the relationships between adjacent joints and to quantify their respective amplitudes within a foot motion. The combined action of the anterior muscles induces a maximal dorsal flexion of 42° for the foot opening angle without arthrodesis while this angle is limited to 9° with arthrodesis. The combined action of the posterior muscles induces a maximal plantar flexion of 33° for the foot opening angle without arthrodesis while this angle is limited to 14° with arthrodesis. Moreover, it is possible to analyse the amplitude of motion for other joints and to estimate their influence during arthrodesis. For example, the angle between the navicular and the talus projected on the sagittal plane during dorsal flexion has a value of 3°. This rotation is produced by the navicular alone since the talus is fixed. Furthermore, the dorsal flexion being 9°, the remaining 6° are probably provided through the Cuneo-metatarsal and Cuneo-navicular joints which are mobilized by the synergistic contraction of the anterior muscles. A similar analysis can be performed for the plantar flexion, leading to an estimation of the muscles' influence.

Based on this example, a biomechanical foot model, driven by muscle activations reproducing the swing motion of the foot with and without arthrodesis, could be used for arthrodesis planning to find the optimal setting of the foot opening angle in order to optimize swing phase during gait. Naturally, the results above need to be validated on real cases before moving to surgical planning. To this aim, a patient specific process to generate a numerical clone of the patient foot would be needed. Using the method presented in [START_REF] Bucki | Clinical workflow for personalized foot pressure ulcer prevention[END_REF] will help to reach this goal. This method is based on the registration of each patient's skin and bone surfaces (extracted from a medical image dataset) with the surfaces of an atlas model. The registration is divided into three steps: rigid registration, affine registration, and local elastic registration. It avoids the tedious and time consuming generation of a new patient mesh and takes few minutes to create a numerical clone for each patient. It consequently gives the possibility to create new patient models easily and to help in surgery planning.

Pressure ulcer prevention

The simulation of static unipodal stance allows estimating the impact of a pressure pattern below the foot and inside the soft tissues. For example, for our model, the highest cluster volume is located in the heel region with 50.8 cm 3 with Von Mises strains above 20 %. No cluster has a volume higher than 0.2 cm 3 when a threshold of 50 % is considered. This almost negligible value seems normal since the pressure applied below the foot is not supposed to cause short term pressure ulcer as it is a unipodal stance. Analyzing only the peak Von Mises strains, figure 18.5, is not really representative as the highest value is in the toe region (161 %), but the cluster volume around this peak strains is extremely small (less than 0.01 cm 3 ). This peak pressure is probably due to a numerical singularity because of a highly deformed element and not really representing the internal strain in this region. Using the cluster volume is consequently more representative as it bypasses those numerical singularities. Overall, the estimation of the cluster volume in our biomechanical model allows determining the level of deformation in the foot when a pressure load is applied below it. It could therefore be used to determine the risk of pressure ulcer and to complete efficiently the monitoring of the skin surface pressure.

FIGURE 18.5 HERE Nevertheless, to allow daily prevention, a model should be fast enough to compute a simulation interactively and it should be coupled to a sensor usable during daily tasks. The computation time of a simulation is difficult to reduce since a high level of details is needed to achieve realistic simulations. The principal area of research here is in the pre computation of most of the possible scenarios to use them as a library on which a daily analysis is based to determine if a specific case, interpolated from the library, produces internal strains high enough to lead to pressure ulcer. Model reduction [START_REF] Cueto | Real time simulation for computational surgery: a review[END_REF] could also be implemented to speed up the computation.

Although fairly inaccurate when the given pressure scenario is outside of the library range, it is accurate and fast when the library covers a wide enough range of scenarios. In term of pressure sensors usable on a daily basis, several are commercially available (such as the ones proposed by Novel, Tekscan, Vista Medical or Texisense).

The key here is to be wearable in a shoe, as an insole or a sock, and to send pressure pattern information regularly to feed simulations with pressure patterns. Finally, for a patient specific prevention, a numerical clone of the patient is needed as the morphology of the patient has a clear impact on the risk of pressure ulcer as demonstrated by [START_REF] Luboz | Influence of the Calcaneus Shape on the Risk of Posterior Heel Ulcer Using 3D Patient-Specific Biomechanical Modeling[END_REF]. Again, using the method presented in [START_REF] Bucki | Clinical workflow for personalized foot pressure ulcer prevention[END_REF] could help to reach this goal.

Conclusion

The foot is a complex structure, involving 28 bones, 30 joints, 210 ligaments, and several layers of soft tissues (skin, fat, muscles). Its purposes are to adapt to the ground during gait or stance, and to store and use the energy while walking. Understanding its behavior during non-pathological function is a whole field of research. Furthermore, it can be afflicted by various pathologies (such as diabetes) or abnormalities (such as hallux valgus), leading to the use of specific devices (orthosis, prosthesis), or even surgery. Musculoskeletal and biomechanical models have been introduced to help understanding these various fields. These models use more or less complex representations, from purely geometrical to biomechanical models using Finite Element modelling to translate constitutive laws that can be linear elastic or more complex nonlinear hyper-elastic and they can model the foot as homogeneous or represent its heterogeneity and the different layer of soft tissues. Given the complexity of the foot, and the variety of behaviors of its components (especially the soft tissues), an accurate model requires a detailed representation of each constituents. This leads to at least two disadvantages: 1) creating a precise foot model is a tedious and long task, and 2) the computation time of a simulation with such model can be very long (quantified in hours). Most of the models published in the literature are either presenting a generic model (with an average ideal morphology) or a clone of a specific patient, therefore creating only one model and avoiding the repetitiveness of the task for other patients. A solution to build rapidly new numerical clones would be to use a registration technique aiming to adapt an atlas model to the morphology of each new patient, in the same fashion as in [START_REF] Bucki | Clinical workflow for personalized foot pressure ulcer prevention[END_REF]. As for the computation time, simulations are most of the time performed for a specific application, and not repeated through a daily clinical set up. It is consequently not a problem to run several simulations over a long period and have the results to be analyzed later. If these models are needed in daily clinical investigations, with a need of interactions, a solution would be to precompute a set of simulations representing a wide range of scenarios and use them to interpolate the results from this abacus.

Figure captions:

Figure 18.1: Transverse regionalization of the foot, with four regions: the back foot (talus and calcaneus), the medio foot (cuboid, navicular, and the three cuneiform bones), the metatarsus (from M1 to M5), and the toes (from hallux to quintus). Figure 18.2: The musculoskeletal model with 30 bones, 15 activated muscles (in red) and 210 ligaments (in green). The insertions of those structures have been determined on the MRI scan and from the literature. Figure 18.3: The complete biomechanical foot model with the bony structures as rigid bodies, and the soft tissues as a Finite Element mesh. The soft tissues are divided in four layers: the skin (in pink), the heel fat pad (in green), the rest of the fat (in yellow), and the muscles (in red). Figure 18.4: Dorsiflexion under three different muscle activations: on the left, activation of the toe and hallux extensor muscles alone; on the center, contraction of the tibialis anterior muscle alone; and on the right, combined contraction of the extensors and the tibialis anterior muscles. Figure 18.5: Peak Von Mises strains in the foot during unipodal stance pressure load. The blue area corresponds to a strain value close to 0 % while the red area is for the peak strain at 161 %. Left, with 20 % of the total pressure, the strains remain below the 20 % strain threshold. Right, with 100 % of the pressure load, the peak strain reaches 161 % (area in red in the toes). It can be seen that the cluster volume (contiguous elements with a color different than blue) is larger in the heel region.

Table captions:

Table 18.1: Survey of the foot models presented in the literature with: the first author of the reference, the team leader, the year of publication, the application, the constitutive law, and the material parameters. Lig = ligaments, SoftTis = Soft Tissues, Li = Linear material (all rows in blue in the table), NoLi = Non Linear material (all rows in yellow in the table) with a non-specified constitutive law, Og = Ogden material (all rows in orange in the table), MoRi = Mooney Rivlin material (all rows in green in the table), Neo Hook heterogeneous material (rows in red), El = Elastic material, HyEl = Hyper-elastic material, Ho = Homogeneous material, Is = Isotropic material.

Acknowledgements

Competing interests:

Some authors are involved with the TexiSense Company (http://www.texisense.com/home_en). Funding: This work was partly funded by the 2010 ANR TecSan IDS project, by the CAMI Labex (ANR-11-LABX-0004) and by the Institut Universitaire de France. Ethical approval: the study was approved by the French institutional review board of Grenoble (IRB 5891 (CECIC) for Rhône-Alpes-Auvergne).


Seyed Ali Elahi

Nathanael Connesson
email: nathanael.connesson@gmail.com

Yohan Payan

A volume-based aspiration method to estimate in-vivo soft tissues stiffness: evaluation of the device with silicone samples

Computer Assisted Medical Interventions (CAMI) involving soft tissues require devices that model such tissues in order to estimate the way they are going to be deformed and/or resected by the surgical gesture. The corresponding biomechanical patient-specific models require an estimation of the constitutive behavior of the soft tissues. Since the mechanical behavior of living tissues varies between in-vivo and ex-vivo conditions, it is important for the CAMI devices to offer an in-vivo estimation with non-traumatic measurements that should undergo sterile conditions. For this purpose, among all the methods proposed in the literature, aspiration/suction is the most widely used technique due to its simplicity and robustness. For such a technique, a device with a hole is put in contact with the soft tissue while a negative pressure aspires part of this tissue. Knowing the relationship between the negative pressure and the aspired tissue height, an inverse problem is then solved to identify the material mechanical properties. In the literature, the apex height is usually measured with a camera and a mirror or a prism, which induces design difficulties, in particular in regards on the required sterilization process for in-vivo measurements. This paper introduces a new method that replaces the optical apex height measurement with a measurement of the aspired tissue volume. The method is referred to as "rate-based method". The camera, mirror and all electronic parts are not required, which makes our device the simplest, lightest and cheapest one could achieve. In particular, this simplification enables the system to meet the severe sterilizations constraints that can be present for some surgeries. Indeed, the proposed device is only composed of a cylindrical plastic chamber (with a hole at its basis), some connection tubes, a syringe pump, a manometer and an aspiration chamber. The idea is to aspirate the tissue inside the chamber using the syringe pump, while measuring the negative pressure and the corresponding removed volume. Such a volume is due to both the aspired tissue volume inside the chamber and the volume changes in the device (air expansion and elasticity of the connections, tubes, syringe, etc.). An off-line calibration process is defined to differentiate the volume changes due to tissue aspiration from the volume changes due to the compressibility of the device (tubes, connections, syringe).

Our device was evaluated with the estimation of the constitutive behavior of two bulk cylindrical silicone samples, one being slightly softer than the other one. For each silicone, thin

flat samples were extracted to characterize their constitutive behavior using a classical traction experiment (MTS Criterion machine, Model 41) and an equibiaxial extension (bulge test). Such experiments respectively provide Young modulus of 16.82kPa and 24.87kPa for silicones #1 and #2. The aspiration tests on each silicone sample were repeated nine times to check the reproducibility of the results. The inverse method to estimate the constitutive behavior assumed a hyper-elastic Gent model with material constants computed using an updated FE simulation. The results validate the ability of our method to discriminate the two silicones with equivalent Young moduli estimated to 18.35kPa and 26.52kPa for silicones #1 and #2.

Keywords: Aspiration Method, In, vivo Measurement, Soft Tissues Characterization, Experimental Mechanics, Inverse Characterization


Anna Mîra

Yohan Payan

Ann-Katherine Carton

Pablo Milioni De Carvalho

Zhijin Li

Viviane Devauges

Serge Muller

Simulation of breast compression using a new biomechanical model

Keywords: Mammography, breast compression, biomechanical model, finite element analysis

 . Therefore, the characterization of the pain perceived during breast compression is of potential interest to compare different compression approaches.

The aim of this work is to develop simulation tools enabling the characterization of existing breast compression techniques in terms of patient comfort, dose delivered to the patient and resulting image quality. A 3D biomechanical model of the breast was developed providing physics-based predictions of tissue motion and internal stress and strain intensity. The internal stress and strain intensity are assumed to be directly correlated with the patient discomfort. The resulting compressed breast model is integrated in an image simulation framework to assess both image quality and average glandular dose. We present the results of compression simulations on two breast geometries, under different compression paddles (flex or rigid).

INTRODUCTION

Today, mammography is the primary imaging modality for breast cancer screening and plays an important role in cancer diagnostics. Screening with mammography has been shown to reduce breast cancer mortality in women over 40-50 years in age. However, some women are still refusing to undergo regular mammography exams. Fleming et al [START_REF] Fleming | Intermittent attendance at breast cancer screening[END_REF] have shown, in a study of 2500 women, that 15% of those who skipped the second appointment, cited an unpleasant or painful first mammogram. The scope of this work is to compare the pain perceived by the women when their breasts are compressed using different designs of compression paddles, and therefore to identify the best compression paddle that may help improving the patient comfort and therefore the adherence to breast cancer screening.

The scope of this work is to assess the differences in terms of patient comfort, image quality and patient dose exposure when using a rigid or a flex paddle. Several studies showed that, the pain experienced by women during the mammographic exam depends on psychologic factor [START_REF] Aro | Pain and discomfort during mammography[END_REF] (technician behavior, patient anxiety), sociologic factors [START_REF] Dullum | Rates and correlates of discomfort associated with mammography[END_REF] (ethnicity, education level) as well as physiologic [START_REF] Poulos | Breast compression in mammography: how much is enough?[END_REF] factors (compression level, breast size). Here, the psychologic and sociologic factors are neglected. The study focuses on physiological factors as the compression force or structural specifications of the compression paddle to characterize the patient comfort.

In this purpose, MR images of two subjects are used to create patient specific finite elements breast models. The mechanical behavior of soft tissues under compression is computed for both subjects and for both paddle designs. The perceived pain for a given paddle design is quantitatively characterized by contact pressure, internal stress and strain distributions. After compression, three sets of macrocalcifications are inserted into breast volumes. The latter are then subject to a Monte-Carlo based simulation (CatSim 8 ) enabling to simulate the image acquisition of the compressed breast with a mammography system. Then, the diagnosis quality is assessed by measuring the signal-difference-to-noise-ratio (SDNR), signal-to-noise-ratio (SNR) and the average glandular dose (AGD).

MATERIALS AND METHODS

Subject specific data

To assess the entire breast volume as well as the surrounding soft tissues, the patient specific data is acquired using a specific MRI coil for breast imaging. The two volunteers taking part to this study agreed to participate in an experiment part of a pilot study approved by an ethics committee (MammoBio MAP-VS pilot study). The volunteers are 59 and 58 years old and have A-cup (volunteer 1) and F-cup (volunteer 2) breast size respectively. The images were acquired with a Siemens 3T MRI scanner using T2 weighted image sequences (Figure 1). The in-plane image resolution was 0.5x0.5 mm, and the slice thickness was 0.6 mm. The volunteers were also asked to provide the compression force and breast thickness as measured on their most recent mammograms. Such data are summarized in Table 1.

Volunteer

Biomechanical breast model

Breast deformation under compression is computed using the improved version of our previously developed breast biomechanical model [START_REF] Milioni De Carvalho | Low-dose 3D quantitative vascular X-ray imaging of the breast[END_REF] . The model shows a good fidelity to breast anatomy considering breast heterogeneity, anisotropy, personalized hyper-elastic properties of breast tissues and sliding boundary conditions between breast and pectoral muscle. The 3D geometrical breast shape is generated from segmented MR images of the volunteer using ITK-snap 6 software. The 3D breast geometry is discretized in a tetrahedral Finite Element (FE) mesh and is subject to a hyper-elastic quasistatic deformation using the ANSYS software framework [START_REF] Poulos | Breast compression in mammography: how much is enough?[END_REF] . Five types of tissues are considered: muscle, adipose and glandular tissues (here referred as breast tissues), skin, fascia and ligaments (Figure 2).

The muscle, breast tissues, skin and fascia are modeled as homogeneous, hyper-elastic quasi-incompressible materials (Poisson ratio υ = 0.499). Mechanical response of soft tissues to external stimuli was modeled using the Gent form of strain-energy potential [START_REF] Mîra | Breast Biomechanical Modeling for Compression Optimization in Digital Breast Tomosynthesis[END_REF] . Considering small strain behavior, the Gent expression reduces to the Neo-Hookean one [START_REF] Yushkevich | User-guided 3D active contour segmentation of anatomical structures: significantly improved efficiency and reliability[END_REF] . Therefore the Young's moduli (E) of the tissues were chosen in conformity with the Neo-Hookean constitutive models described in our previous work [START_REF] Milioni De Carvalho | Low-dose 3D quantitative vascular X-ray imaging of the breast[END_REF] and are summarized in Table 2. The limiting value J m of the first invariant of the left Cauchy-Green deformation tensor was fitted to comply with the compression data from Table 2

Compression paddles

During mammography, a qualified radiologic technologist positions the breast of the patient between a stationary image receptor and a movable paddle (Figure 3). The technologist gradually compresses the breast in order to even out the breast thickness and to spread out the soft tissues. Nowadays, two types of compression paddles are widely available: rigid compression paddles (RCP) and flex compression paddles (FCP).

The RCP is fixed to its frame and is constrained to move in the up-down direction. This paddle has some flexibility because of material mechanical properties and can slightly bend when compressing the breast, while remaining globally flat and parallel to the image receptor. On the other hand, the FCP is attached to its frame by rotational joints and therefore, presents an additional rotational degree of freedom enabling the paddle to tilt with respect to the image receptor plane (Figure 3.c).

With a rigid paddle, the breast under compression presents a nearly uniform thickness all over the contact surface. Contrariwise, with a flex paddle, the compressed breast thickness decreases quasi linearly from the chest wall to the nipple. Flex paddles are used to better conform the breast contours and thereby to improve compression. However, Broeders and collegues [START_REF] Ansys® | Academic Research Mechanical, Help System, Modeling and Meshing Guide[END_REF] have shown that such compression paddle may decrease the diagnostic quality of mammograms as the breast tissues may be pushed out to the chest wall resulting in less retro-glandular tissue visible on the image.

Breast compression

The prone configuration of the breast is used to simulate the effect of the breast pulling gesture performed by the technologist during breast positioning. Neglecting the gravity effects, the breast geometry is then tilted to the body upright position. To simulate the cranio-caudal incidence, the image receptor is positioned at the inframammary ligament level while the paddle compresses the breast by a downward movement. The minimal thickness and maximal compression force were consistent with the subject real acquisition data (Table 1).

A finite element model was defined for each type of compression paddles based on the geometry specifications from a Senographe Pristina TM mammography unit (GE-Healthcare, Buc, France). On both cases, the paddle flexibility due to the material properties was neglected. The RCP was therefore defined as a fully rigid body with only one translational degree of freedom in the downwards direction. For the FCP, a rotation around the longitudinal axis is added. An additional degree of freedom was modeled using a rotational-only joint type of element. The joint stiffness was computed by fitting the force-deflection curve of a flex paddle from a Senographe Pristina TM unit.

The compression was then computed for both right breast volumes. To compare the two types of compression, a nominal breast thickness was defined. The nominal thickness for the RCP was defined as the distance between the paddle and the image receptor. The nominal thickness for the FCP was defined as the distance between the compression paddle and the image receptor at mid-length of the compressed breast. Then, we compared the contact pressure and the internal strain/stress distributions. The nominal breast thickness was computed and compared at different compression levels.

Image quality and average glandular dose

To assess the impact of breast compression on image quality, we inserted a set of microcalcifications into each compressed breast volume. The smallest breast volume contains 21 microcalcifications arranged in a matrix of 7 rows and 3 columns (Figure 4a). The largest breast volume contains 56 microcalcifications arranged in a matrix of 7 rows and 8 columns. The matrix of calcifications is parallel with the entrance surface of the image receptor and positioned at the breast mid thickness (Figure 4b). The distance between two consecutive columns or rows is 10mm. We assumed a uniform breast-equivalent material composed of glandular/adipose tissue with a 20/80 ratio. A mammogram was simulated using typical clinical acquisition parameters obtained with the standard automatic optimization of parameters (AOP) mode. Two simulations were performed with microcalcifications of 0.2 mm and 0.3mm in diameter. The signaldifference-tonoise ratio (SDNR) per pixel of these microcalcifications was measured. Additionally, the signal-to-noise ratio (SNR) was computed on the same pixels excluding the microcalcifications. The average glandular dose (AGD) was derived using the approach proposed by Dance et al12 regardless the paddle type. In practice, it is very difficult to accurately measure the exact breast thickness. Thus, the nominal breast thickness was used to compute conversions factors which relate measurements of incident air kerma to the delivered mean glandular dose.

The force versus breast thickness curves are plotted in Figure 5 for both volunteers with the rigid and flex paddles. We observed a nominal compression thickness roughly equals for both paddles. The resulting internal stress and strain distributions, as well as contact pressure maps were derived at compressive forces of 22 N for the first volunteer (Figure 6) and 95 N for the second one (Figure 7). As concerns the small breast volume (Figure 6), there is no significant difference between FCP and RCP in pressure distribution over the skin surface or in internal stress/strain intensity distributions. For both compression paddles, high pressure at the skin surface is concentrated in the juxtathoracic region with a maximum pressure of 77.7 kPa. Several clinical studies [START_REF] Ansys® | Academic Research Mechanical, Help System, Modeling and Meshing Guide[END_REF][START_REF] Chagnon | A comparison of the Hart-Smith model with Arruda-Boyce and Gent formulations for rubber elasticity[END_REF] sustained this result of no significant difference in experienced pain when using FCP or RCP. In addition, the FE simulations confirm that in small breasts the paddle tilt is too small to impact the tissues compression in the middle part of the breast.

FCP applied on large breast volumes (Figure 7) results in significantly lower intensities of pressure at the skin surface in contact with the compression paddle, with a maximal pressure of 37 kPa, compared to 56 kPa when using RCP. No significant difference in the measured maximal intensities of strain and stress was observed, however strain and stress distribution patterns are different. When the breast is compressed with a rigid paddle, maximal strain and stress are concentrated in the retromammary space and decrease considerably toward the nipple. When a flex paddle is used, stress and strain are more uniformly distributed over the breast volume with the highest values in the middle third of the breast.

The areal pressure distribution patterns has already been demonstrated in the work by Dustler and collegues [START_REF] Chagnon | A comparison of the Hart-Smith model with Arruda-Boyce and Gent formulations for rubber elasticity[END_REF] . The authors have studied the pressure distribution patterns of 103 women undergoing breast compression with a rigid paddle at different compression levels. Four groups have been differentiated: a) skin pressure widespread over the breast (29%); b) skin pressure concentrated on the central part of the breast (8%); c) skin pressure concentrated on the juxtathoracic region (16%); d) skin pressure concentrated along a narrow zone at the juxtathoracic region (26%). The pressure distribution patterns observed for our first and second volunteers correspond to the group d and a respectively. The nominal breast thickness may vary by about 2mm between rigid and flex paddle for both volunteers (Table 3). Accordingly, no significant difference was found between the estimated AGD, while a dose reduction of 2% for the smaller breast and 4% for larger breast was observed. The SNR and SDNR have been estimated and compared between flex and rigid paddles. When using a flex paddle instead of a rigid paddle on the largest breast (volunteer 2), we observe a statistically significantly higher SNR. The same trend is observed on SDNR for both 200 and 300 µm microcalcifications, while not statistically significant. We did not observe any statistically significant difference in SNR or SDNR for microcalcification of any size when considering the compression of the smallest breast by a rigid or a flex paddle. Therefore, despite a breast thickness varying linearly from chest wall to nipple when the flex compression paddle is used, the image quality is preserved or improves compared to the image quality obtained with the rigid compression paddle.

CONCLUSION

Breast compression with flex and rigid paddle have been simulated using the finite elements theory applied to segmented MRI images acquired on 2 volunteers under different geometries. Appling the Gent form of strain-energy potential, instead of the Neo-Hookean form, allowed to obtain compression force magnitudes comparable with the real subject data.

After simulating the breast compression, the SDNR of microcalcifications and the AGD, delivered during the acquisition of the corresponding simulated mammography, have been computed. The simulations have been repeated for two different breast volumes (cup sizes A and F) with a rigid and a flex paddle. The four configurations have been analyzed to compare patient perceived pain (measured as strain and stress) and image quality (measured as SNR, SDNR and AGD). The results of our simulations indicate that, for the smallest breast, there is no significant difference for the patient perceived pain when using the rigid or the flex paddle. The shape of the breast under compression does not present significant changes between the two paddle designs. We did not observe any statistically significant difference in SNR or SDNR for microcalcification of any size when considering the compression of the smallest breast by a rigid or a flex paddle. Therefore, our results suggest that using a flex paddle should not significantly impact image quality and delivered dose in small breasts, and should not reduce significantly the perceived pain.

For the largest breast, our simulations indicate that using a flex paddle may reduce the maximal pressure intensity on the skin surface by about 30% compared to the rigid paddle. The tissues deformation is more uniformly distributed inside the breast volume, and the highest deformation is occurring in the middle breast region corresponding to the supposed location of dense tissues. Moreover, our simulations have shown that flex paddle have no significant impact on the average glandular dose and improves image quality compared to the rigid paddle.

In conclusion, our simulations confirm that using the flex paddle used for breast compression may improve the patient comfort without affecting the image quality and the delivered average glandular dose. Moreover, despite a breast thickness varying linearly from chest wall to nipple, when a flex compression paddle is used on large breasts, the image quality seems to be preserved or improved compared to the image quality obtained with a rigid compression paddle.

Figure 1 :

 1 Figure 1 : Patient breast MR images in prone body position: left -volunteer 1; right -volunteer 2.

Figure 2 :

 2 Figure 2: Biomechanical breast model: left panel -breast components; right panel breast suspensory matrix.

Figure 3 :

 3 Figure 3: Breast compression between the paddle (up) and the receiver (down): a) Rigid paddle; b) Flex paddle with flexion angle α.

Figure 4 :

 4 Figure 4: Microcalcification distribution over the smallest breast volume (d = 10mm): a) axial view, b) sagittal view.

Figure 5 :

 5 Figure 5: Resulting breast thickness for a given compression force

Figure 6 :

 6 Figure 6: Stress, strain and contact pressure distribution for the first subject.

Figure 7 :

 7 Figure 7: Stress, strain and contact pressure distribution for the second subject.

Table 1 :

 1 Compression force and breast thickness for both volunteers for a cranio-caudal mammogram.

			1	Volunteer 2
		Right breast Left breast Right breast Left breast
	Force (N)	21.9	40.9	94.8	56.6
	Breast thickness (mm)	47	42	50	49

Table 2 :

 2 .

		E muscle (kPa) E breast (kPa) E skin (kPa) E fascia (kPa)	J m
	Subject 1	10	0.4	3	120	1
	Subject 2	10	0.5	10	160	2

Constitutive parameters for the used materials: E -Young's modulus; Jmthe limiting values of the first invariant of the left Cauchy-Green deformation tensor.

Table 3 :

 3 Breast nominal thickness (BNT), average glandular dose (AGD), signal-to-noise-ratio (SNR) and signal-difference-to-noise (SDNR) for both volunteers and both compression paddle types.

		Rigid Paddle	Flex Paddle			Rigid Paddle	Flex Paddle
		Mean SNR StdDev SNR Mean SNR StdDev SNR p-Values		BNT (mm)	AGD (mGy) BNT (mm) AGD (mGy)
	Volunteer 1	82,90	43,09	83,70	37,72	0,706		46	1,15	44	1,12
	Volunteer 2	126,89	8,75	137,21	10,73	0,000		48	1,20	46	1,15
				200 um					300 um
		Rigid Paddle	Flex Paddle		Rigid Paddle	Flex Paddle
		Mean SDNR StdDev SDNR Mean SDNR StdDev SDNR p-Values Mean SDNR StdDev SDNR Mean SDNR StdDev SDNR p-Values
	Volunteer 1	0,74	0,68	0,79	0,54	0,689	2,01	1,28	1,85	1,02	0,224
	Volunteer 2	1,14	0,57	1,13	0,53	0,885	2,96	0,76	3,15	0,92	0,093