V Luboz
email: vluboz@imag.fr

A Perrier

I Stavness
email: stavness@gmail.com

J E Lloyd
email: lloyd@cs.ubc.ca

M Bucki

F Cannard

B Diot
email: b.diot@ids-assistance.com

N Vuillerme
email: nicolas.vuillerme@agim.eu

Y Payan
email: ypayan@imag.fr

Foot Ulcer Prevention Using Biomechanical Modeling

Keywords:

Foot ulcers are a common complication of diabetes and are the consequence of trauma to the feet and a reduced ability to perceive pain in persons with diabetes. Ulcers appear internally when pressures applied on the foot create high internal strains below bony structures. It is therefore important to monitor tissue strains in persons with diabetes. We propose to use a biomechanical model of the foot coupled with a pressure sensor to estimate the strains within the foot and to determine if they can cause ulcer formation. Our biomechanical foot model is composed of a Finite Element mesh representing the soft tissues, separated into four Neo Hookean materials with different elasticity: plantar skin, non-plantar skin, fat and muscles. Rigid body models of the bones are integrated within the mesh to rigidify the foot. Thirty-three joints connect those bones around cylindrical or spherical pivots. Cables are included to represent the main ligaments in order to stabilize the foot. This model simulates a realistic behavior when the sole is subjected to pressures measured with a sensor during bipedal standing. Surface strains around 5 % are measured below the heel and metatarsal heads while internal strains are close to 70 %. This strain estimation, when coupled to a pressure sensor, could consequently be used in a patient alert system to prevent ulcer formation.

Introduction

Diabetes affects the lower limb through long-term chronic complications, including peripheral neuropathy and angiopathy [START_REF] Shaw | The pathogenesis of diabetic foot problems. An overview[END_REF]. When the diabetic foot is subject to a trauma, such as repetitive stress from high pressure, illfitting footwear rubbing on the skin or an object inside the shoe, an ulcer can form. This condition if often followed by amputation of a toe or foot: a limb is lost every 30 seconds on Earth because of diabetes. In 2007, in the US, the cost of ulcers was 48 billions USD [START_REF] Driver | The costs of diabetic foot: the economic case for the limb salvage team[END_REF] and life expectancy after foot amputation was 50% after 5 years [START_REF] Reiber | Epidemiology of Foot Ulcerations and Amputations in Diabetes[END_REF].

The two main factors that promote ulcer development are the excessive intensity and the repetition of pressures on the foot. This is worsened by diabetic neuropathy which reduces or even suppresses sensation in the feet of persons with diabetes. Three mechanisms are linked to the pressures leading to ulcer [START_REF] Mueller | Etiology, evaluation, and treatment of the neuropathic foot[END_REF]: [START_REF] Shaw | The pathogenesis of diabetic foot problems. An overview[END_REF] ischemia caused by increased pressure duration, (2) large tissue strains created by increased pressure magnitude, or (3) tissue fatigue caused by increased number of pressure loads.

Today, prevention of foot ulcers is mainly based on the patient's daily vigilance and on frequent medical doctor monitoring. This approach is consequently dependent on the patient's involvement which tends to decrease when no external signs are visible. Unfortunately, when the first ulcers appear, serious complications have already begun for the diabetic patient because of diabetic angiopathy, which limits tissue healing and increases recurrence.

Reducing the pressure load at the interface between the diabetic foot and the trauma source prevents further ulceration and facilitates wound healing [START_REF] Armstrong | Evidence-based options for off-loading diabetic wounds[END_REF]. Pressure loads can be reduced with devices such as casts, orthotics, insoles or foam bandages.

To aid patient monitoring, devices measuring pressures at the foot/insole interface were introduced few years ago [START_REF] Martinez-Nova | BioFoot® in-shoe system: Normal values and assessment of the reliability and repeatability[END_REF][7] [START_REF] Descatoire | Baropodometric information return device for foot unloading[END_REF].

Novel (http://www.novel.de), Tekscan (http://www.tekscan.com), Vista Medical (http://www.pressuremapping.com) and Orpyx (http://orpyx.com/) now propose commercial devices that are composed of few pressure sensors located under the heel and the metatarsal heads. Unfortunately, their price and their need to be constantly connected to monitor the plantar pressures in real time make them difficult to use in a daily long-term prevention routine. Moreover, these devices are limited to measuring pressures at the skin surface and it is well known that deep pressure ulcers stem from internal stresses and strains. Such ulcers usually start in deep tissues and progress outward rapidly, causing substantial subcutaneous damage underneath intact skin. Deep pressure ulcers are therefore particularly dangerous since they may be difficult to detect visually. Devices that measure surface pressures are mainly used to alert persons suffering from diabetes to abnormal pressures that may cause skin damage [START_REF] Pipkin | Effect of model design, cushion construction, and interface pressure mats on interface pressure and immersion[END_REF]. These measurements, however, cannot predict ulcer formation due to internal tissue loading [START_REF] Linder-Ganz | Strains and stresses in sub-dermal tissues of the buttocks are greater in paraplegics than in healthy during sitting[END_REF]. For example, a similar pressure distribution could be observed under the heel of a thin person with blunt calcaneus bone and a heavy diabetic person with sharp calcaneus bone; however, the likelihood of a pressure ulcer forming depends on the calcaneus bone curvature as well as the thickness of the soft tissues as proved in [START_REF] Sopher | Exposure to internal muscle tissue loads under the ischial tuberosities during sitting is elevated at abnormally high or low body mass indices[END_REF] for the ischial tuberosity.

In order to account for these anatomical differences and to quantitatively estimate the internal stresses and strains from the measured external pressures, several studies have proposed [START_REF] Shaw | The pathogenesis of diabetic foot problems. An overview[END_REF] to build a patient-specific biomechanical model of the foot including soft tissues and bony prominence, and (2) to use this numerical model to compute the internal strains and stresses [START_REF] Elsner | Is obesity a risk factor for deep tissue injury in patients with spinal cord injury?[END_REF]. Several Finite Element (FE) models have been introduced to simulate the foot deformations. Several studies such as [START_REF] Yarnitzky | Real-time subjectspecific monitoring of internal deformations and stresses in the soft tissues of the foot: a new approach in gait analysis[END_REF] proposed 2D biomechanical models of the foot with fairly complex material behavior and anatomy, providing very interesting results in terms of stresses and strains within the foot. Nevertheless, the foot behavior is highly dependent on its whole 3D shape and the influence of each of its structures (bones or tendons for example). Also, it seems that clinicians prefer a 3D analysis of the foot to determine the development of a pressure ulcer. In [START_REF] Ledoux | A quasilinear, viscoelastic, structural model of the plantar soft tissue with frequency sensitive damping properties[END_REF], the foot soft tissues (skin, fat and muscles) are modeled as a 3D FE mesh with a homogeneous linear elastic material, the bones are modeled as another set of rigid FE meshes, while contacts between the bones simulate the joints, and cables connect the bones for the ligaments located in the mid foot. The 3D biomechanical model presented in [START_REF] Chen | Effects of internal stress concentrations in plantar soft-tissuepreliminary three-dimensional finite element analysis[END_REF] seems closer to reality: it models the soft tissues in a large deformations framework with a Mooney Rivlin constitutive law and it adds almost all the ligaments of the foot. However, the main drawback of these models is their long computation times: it can take hours to compute simple plantar loadings. For this reason, it would be difficult to integrate previously proposed foot models into a daily foot ulcer prevention process. In this article, we propose and discuss the challenges of defining a simple, yet realistic, biomechanical model that could estimate internal foot stresses and strains in interactive time.

Materials and methods

Our foot biomechanical model has been developed using the 3D simulation platform, ArtiSynth [START_REF] Lloyd | ArtiSynth: a fast interactive biomechanical modeling toolkit combining multibody and finite element simulation[END_REF] (www.artisynth.org). The model is composed of soft tissues, bones, ligaments and joints.

Soft tissue FE mesh

The soft tissues are modeled as a FE mesh and are divided into four layers: the plantar skin (below the heel, the metatarsus, and the toes), the non-plantar skin, muscles and fat. The outer surface of this mesh was based on the skin surface from the Zygote database (www.zygote.com). Using an automatic FE mesh generator [START_REF] Lobos | Techniques for the generation of 3D Finite Element Meshes of human organs[END_REF], the model skin surface was filled with finite elements. This generator aims at creating a mesh with a maximum of hexahedrons (to limit the locking effect observed for tetrahedral elements in quasi-incompressible assumptions) while keeping a controlled number of elements to enable fast FE computation. The accuracy of the surface fitting is made possible by optimally subdividing hexahedrons intersecting the organ surface into prisms, pyramids and/or tetrahedron depending on the geometry of the surface locally. Because we assume the bones to be rigid (see section 2.2), there is no need to represent them as FE meshes in our simulation. We therefore set the mesh generator to leave holes in the soft tissue mesh in order to represent the bones' geometry (also given by the Zygote database).

The FE mesh representing the foot soft tissues is shown in Figure 1 and is composed of 36,895 elements (11,045 hexahedrons, 10,206 pyramids, 9,991 tetrahedrons, and 5,653 wedges) and 22,774 nodes. Neo Hookean materials were chosen for the foot tissues in order to represent large deformations. The mesh has four layers, each with distinct material properties: Young moduli were set to 6 MPa for the plantar skin, 200 kPa for the rest of the skin, 50 kPa for the muscles and 4 kPa for the fat (Figure 2). Assuming these tissues are quasi-incompressible, we set their Poisson ratio to 0.495. With the exception of the non-plantar skin, the values of the material properties were taken from measurements reported by Sopher [START_REF] Sopher | The influence of foot posture, support stiffness, heel pad loading and tissue mechanical properties on biomechanical factors associated with a risk of heel ulceration[END_REF] for an Ogden model and converted here to fit our Neo Hookean model (taking the equivalent Young initial modulus). The most superficial layer of elements simulates the foot skin. Two different materials were chosen for the skin because of the high stiffness of the plantar skin compared to the rest of the skin. The elasticity modulus chosen for the non-plantar skin was measured using our own device: LASTIC (for Light Aspiration device for in vivo Soft TIssue Characterization). The LASTIC device is based on the pipette aspiration principle [START_REF] Stavness | Coupled Hard-Soft Tissue Simulation with Contact and Constraints Applied to Jaw-Tongue-Hyoid Dynamics[END_REF] and aims at characterizing in vivo the elastic modulus of soft tissues. It creates a negative pressure against the tissue surface on which it is placed and consequently deforms the surface of these tissues. The height of the deformation is measured on the images captured during the acquisition with several negative pressures. The corresponding height/negative pressure curve allows estimating the elastic modulus of the tissues through inverse method. LASTIC was used to estimate the elasticity of the foot arch skin of a healthy subject (Figure 3) and gave a value of 200 kPa. This value was used for the whole skin except for the plantar surface which was too stiff to be evaluated with LASTIC (as mentioned above, the 6 MPa proposed by [START_REF] Sopher | The influence of foot posture, support stiffness, heel pad loading and tissue mechanical properties on biomechanical factors associated with a risk of heel ulceration[END_REF] was chosen for this plantar skin). To define the two skin layers, all the elements of the surface of the FE mesh that were below a certain height (in our case, about 2 mm) in the model were considered as plantar skin (see Figure 2). The other surface elements were assigned the material with a Young modulus of 200 kPa for the rest of the skin. The muscle layer was defined from the Zygote database and is limited for now to the muscles of the foot arch. This area was mapped in the FE mesh to find the elements considered as muscles (Figure 2). All other elements were assumed to be part of the fat layer.

Bones and joints

The ArtiSynth framework provides a very efficient computational formulation for the coupled simulation of rigid and deformable structures, with a constraint-based mechanism for attaching bones and soft tissues and the use of semi-implicit time integration [START_REF] Luboz | Validation of a Light Aspiration device for in vivo Soft TIssue Characterization (LASTIC)[END_REF]. The bones are therefore represented as rigid body surfaces in our model. The 26 foot bones and a section of the tibia and fibula are integrated in the model from the Zygote database (Figure 4). By not modeling the bones as FE meshes, we are able to decrease the FE matrix size and speed up the simulation. We worked with an anatomist to define all the boundary conditions between bones and soft tissues. Each rigid body is fixed to the nearby finite element nodes of the soft tissue mesh to naturally rigidify the foot. Each bone is automatically connected to its neighbor by a joint. In our model, those joints are simulated by 32 spherical pivots and one cylindrical pivot. The cylindrical pivot is between the talus and calcaneus bones and is placed below the talus (Figure 5). It constrains the motion of the calcaneus to a single rotation axis following the long axis of the foot. Spherical pivots have three rotation axes. The angles permitted by these pivots vary depending on their location to simulate the possible motions of actual foot joints, (Figure 6). For the phalanges, a maximum rotation angle of 45 degrees is allowed while the maximum allowable rotation of the metatarsal pivots is 30 degrees. The joints for the bones in the mid and back foot were limited to rotation angles between 0 and 5 degrees to simulate proper rigidity of the foot. The main motion of all the joints is around a rotation axis given by the intersection of the frontal and horizontal planes at the location of the corresponding joint. The foot is constrained by large ligaments to support body weight and to ensure stability. Our model integrates the six main ligaments of the foot. The plantar fascia connects the calcaneus to the proximal part of the first phalanx of each toe. In the model, it is separated in two sets: (1) the external fascia for the outer ligament layer (Figure 7a) and (2) the internal fascia for the inner ligaments close to the mid foot (Figure 7b). The transversal metatarsal head ligament connects the proximal part of the first phalanx of each toe to constrain lateral spreading of the fore foot (Figure 7c). The Achilles tendon connects the calcaneus to the calf and to the bottom part of the knee (Figures 8 and9). In our model, the Achilles tendon is also separated into two sets: (1) the external tendon for the part closer to the surface and going towards the distal part of the calf, and

Ligaments

(2) the internal tendon for the deeper part going towards the knee. Three smaller ligaments are also integrated in the model: the triangular ligament (Figure 8) between the navicular, calcaneus and cuboid bones, and two internal ligaments (Figure 9) between the calcaneus and navicular bones, and between the talus and navicular bones. Ligaments are modeled as cables connecting the previously mentioned bony structures. Cables were chosen over anisotropic finite elements to reduce computation time, since this kind of elements would have required a much finer mesh to represent separate ligaments. Ligament attachments and paths have been manually created to guarantee precision of the bony insertion site for the ligaments. These insertion sites were defined by an experienced podiatrist. Based on the work of Gefen [START_REF] Gefen | Plantar soft tissue loading under the medial metatarsals in the standing diabetic foot[END_REF], we assigned different stiffness to the ligaments whether they are elongated or compressed. It has indeed been observed that they have a high stiffness when they are elongated, but no stiffness when they are compressed. Consequently, we assigned a stiffness value to the cables of 200 MPa in extension and of 0 kPa in compression. As a result, the ligaments resist elongation, but not compression.

Boundary conditions

As discussed in the introduction, internal tissues strains are likely the most significant factor for inducing deep pressure ulcers [START_REF] Loerakker | The relative contributions of muscle deformation and ischaemia to pressure ulcer development[END_REF]. Therefore, it seems important to monitor these values, which is possible with our biomechanical model. To study the influence of a given foot position on internal strains, a set of pressures simulating the patient's foot in bipedal standing were applied to our biomechanical model. These plantar pressures were measured using a commercially available pressure sensor (Zebris platform, http://www.zebris.de/) under the right foot of a volunteer while standing. The left foot was also on the ground, but not on the pressure sensor, so that about half the body weight was applied onto each foot. The pressures ranged from 0 to 10.5 N.cm -2 and are spread mainly below the heel and the metatarsal heads, Figure 10. The standing posture was simulated in two steps. First, the tibia and fibula bones were fixed and the rest of the foot was let loose under the influence of gravity for 0.2 s to reach a resting position. Secondly, from t = 0.2 s to 3 s, the previously measured set of pressures was applied to the nodes of the foot sole following a slope (0% at 0.2 s and 100 % at 3 s) to model normal bipedal standing.

Results

Figures 11 and12 show the Von Mises strains corresponding to the simulation at t = 3 s at the time of maximum pressures applied to the sole of the foot. Table 1 summarizes the Von Mises strains at six key points of the foot: below the metatarsal head of each toe and below the heel, both at the skin surface and near the bones. The strains at the foot surface were markedly lower than the internal tissue strains. Maximal strains over the foot ranged from 2.7 % to 8.0 % for the skin surface and from 43.0 % to 96.8 % for the soft tissues near the bony structures. This strain distribution (with at least a 10 fold ratio between internal strains and foot surface strains) is consistent with the realistic behavior that is responsible for the creation of foot ulcers: they appear internally near the bones (because of very high strains) before being visible at the skin surface. Furthermore, the range of strains computed within the foot is consistent with the range of maximal strains observed in [START_REF] Yarnitzky | Real-time subjectspecific monitoring of internal deformations and stresses in the soft tissues of the foot: a new approach in gait analysis[END_REF] for four subjects: between 52 % and 107 %. 11. They result from the fact that the tibia and fibula are fixed while the rest of the bones are moving during the simulation, hence creating fairly large displacements near the ankle joints. During bipedal standing, the tibia and fibula would probably move slightly to adjust the position of the hips and the rest of the body, and consequently these high strains would probably be reduced. The sole of the foot is therefore the main region of interest in this simulation. The simulation takes 22 minutes to compute with ArtiSynth on a standard PC equipped with an Intel i5 2.80 GHz processor and 3.42 Go of RAM.

Discussion and conclusion

In this paper, we have proposed a biomechanical foot model to estimate the strains within the foot and determine whether ulcers may appear for given patient postures and plantar pressures. The goal of this work was to propose an anatomically accurate 3D foot model that can simulate at interactive rates.

Our foot model includes several structures. A FE mesh represents the foot soft tissues with four different Neo Hookean materials for plantar skin, non-plantar skin, the arch muscles and the fat in between those layers. Rigid bodies model the 26 foot bones that are connected with 33 pivots representing the joints. Three large ligaments (plantar fascia, transversal metatarsal head ligament, and Achilles tendon) and three smaller ligaments are also integrated in the model to connect certain bones and to stabilize the foot.

Simulations of foot loading with the model showed a realistic behavior in standing stance in terms of surface and internal strains with pressures measured with a commercial sensor. Average surface strains of 5 % were evaluated below the heel and metatarsal heads while average strains of 70 % were measured internally, next to the bones. These initial results are promising, and we are currently performing a thorough validation study.

A key interest of our 3D model is that different levels of strains were computed at different locations within the foot as shown in Table 1. This clearly demonstrates the advantage of a 3D model over a 2D model: the capability to study the variations of the internal strains in terms of the anatomy of the patient and the pressure pattern applied below the foot. Our analysis of the tissue strains resulting from a prescribed load at the sole of the foot reveals a localization of higher pressures inside the foot than on its surface. Consequently, the model could be coupled to the output of a pressure sensor to determine whether or not ulcers may appear as a prevention tool for diabetic patients.

Several limitations of our model need to be addressed before reaching this prevention goal. First, to be used by a given patient, the model needs to take into account the patient's morphology. It is obvious that the location of the bony structures and their shape is a key point in the process of creation of an ulcer. Therefore, using a patientspecific model would be the first step to achieve in order to claim ulcer prevention. Using the Mesh-Match-and-Repair (MMRep) algorithm [START_REF] Bucki | A Fast and Robust Patient Specific Finite Element Mesh Registration Technique: Application to 60 Clinical Cases[END_REF][24], it is possible to generate patient-specific biomechanical models from an image dataset. From a complex dataset (CT scan or MRI), it gives extremely accurate results in terms of modeling. Nevertheless, these imaging modalities are rarely used during ulcer prevention. Using cheap and easily accessible image modalities, such as ultrasounds, X-rays, or even a Kinect scan of the skin [START_REF] Bucki | Multi-modal framework for subject-specific finite element model generation aimed at pressure ulcer prevention[END_REF], the result of the MMRep algorithm remains accurate to a precision of around 1.4 mm [START_REF] Bucki | Multi-modal framework for subject-specific finite element model generation aimed at pressure ulcer prevention[END_REF]. The algorithm computes a non-rigid transformation between the structures segmented in the medical dataset and the biomechanical model. This transformation is then applied to the original model to create a new model corresponding to the anatomy of the patient. Based on this process, it would be possible to build biomechanical model of the foot for each patient and to use it to quantify the internal pressures that could create ulcer for this specific case as could be used in a foot ulcer prevention. The LASTIC device could be used to characterize the patient's skin and help to make the foot model more patient specific. It takes about two minutes to measure the elastic modulus of the skin below the foot arch. Using that device for different locations on the foot could improve the elasticity values given to the model and therefore be closer to the mechanics of the patient's foot.

The second step towards reaching the goal of ulcer prevention is the daily monitoring of the patient's foot. Using a pressure sensor as heavy and as expensive as the Zebris platform, for instance, is only possible in a laboratory or a clinical setting. A lighter and less expensive pressure sensor would permit equipping every patient and monitoring their foot pressures on a daily basis. Using a similar technology than that recently employed for the conception of the TexiCare device dedicated to the prevention of seated buttock pressure ulcers for people with spinal cord injury [START_REF] Chenu | TexiCare: An innovative embedded device for pressure ulcer prevention. Preliminary results with a paraplegic volunteer[END_REF], the "Smart Sock" sensor [START_REF] Bucki | The TexiSense « Smart Sock » -Textile Pressure Sensor and 3D Real-time Finite Element Model of the Diabetic Foot for a Daily Prevention of Pressure Ulcers[END_REF], also developed by Texisense (http://www.texisense.com/), combines both qualities since it is made of a 100 % textile pressure sensing fabric wirelessly connected to a controller which can record and monitor the pressures all around the foot (not only under the sole). It can be used continuously during daily living activities. When combined to a biomechanical model such as ours, its measured pressures could be used to estimate the internal strains. If such strains overcome a threshold in duration or in intensity [START_REF] Loerakker | The relative contributions of muscle deformation and ischaemia to pressure ulcer development[END_REF], an alert would thus have to be sent to the patient to prevent ulcer formation.

Finally, a more accurate and timely monitoring of the internal tissue pressures will require a model that is capable of dynamic estimation of these pressures. The "Smart Sock" sensor would help in this direction since it would give continuous measurements of the external pressures. Nevertheless, moving from simulating static standing to dynamic motion is the third area that needs to be addressed for complete ulcer prevention. The main challenge will reside in adapting our model to reduce its computation time. Although our model is faster than most realistic foot FE models from the literature, in its current form, it is indeed too long to spend 22 minutes to compute internal stresses from a set of external pressures. For this reason, it is not currently possible to integrate this model in a real-time foot ulcer prevention device. Three possible solutions will be soon investigated to achieve real-time foot simulations: (1) decreasing the time of computation by reducing the number of elements in the areas less involved in ulcer formation (above the ankle for example), (2) precomputing accurate simulations with our current model and using them to evaluate the internal pressures from the continuously measured external pressures, or (3) limiting the FE modeling to a specific part of the foot such as the heel and/or the metatarsal heads. In any case, the foot model introduced in this paper will be useful since it now represents a reference in terms of modeling to which any "reduced" model will be compared.

Figure 1 .

 1 Figure 1. The surface of the FE mesh representing the foot soft tissues.

Figure 2 .

 2 Figure 2. Cross section of the FE mesh representing the foot soft tissues. The plantar skin layer, the muscle layer and some gaps defining the bones' locations can be seen. The rest of the surface elements are part of the softer skin layer, while the other internal elements form the fat layer.

Figure 3 .

 3 Figure 3. Acquisition of the elastic modulus of the skin below the foot arch using LASTIC.

Figure 4 .

 4 Figure 4. The bony structures within the foot model. Some tendons and ligaments are visible as yellow cables.

Figure 5 .

 5 Figure 5. Cylindrical pivot (scaled up for visibility) representing the joint between the talus and calcaneus bones.

Figure 6 .

 6 Figure 6. Pivots simulate the joints connecting the different bones. They allow different maximum angles depending on the foot regions: 45 degrees for the phalanges, 30 degrees for the metatarsi, and 0 to 5 degrees for the rest of the foot.

Figure 7 .

 7 Figure 7. The main foot ligaments are simulated as cables. Posterior view of the models shows, a) the external plantar fascia, b) the internal plantar fascia, and c) the transversal metatarsal head ligament.

Figure 8 .

 8 Figure 8. Lateral view of the main ligaments of our model: on the left, the Achilles tendon, and below the arch, horizontally, the two sets of plantar fascia. The triangular ligament linking the navicular, calcaneus and cuboid bones is visible near the ankle joint.

Figure 9 .

 9 Figure 9. Medial view of the main ligaments of our model: on the right, the Achilles tendon, and below the arch, horizontally, the two sets of plantar fascia. The two internal ligaments linking the calcaneus and navicular bones, and the talus and navicular bones are visible near the ankle joint.

Figure 10 .

 10 Figure 10. Distribution of the pressures applied under the right foot. Highest pressures, around 10.5 N.cm -2 , are in red (heel and metatarsal heads) and lower are in green.

Figure 11 .

 11 Figure 11. Left: cross section at the second toe showing the internal Von Mises strains resulting from the pressures applied to the foot sole and simulating a standing position. Right: Von Mises strains below the foot. The color map goes from no strain in blue to strains of 100% in red.

Figure 12 .

 12 Figure 12. Left: cross section at the heel showing the internal Von Mises strains resulting from the pressures applied to the foot sole and simulating a standing position (the "staircase" shape of the foot surface is simply a visual artifact due to the clipping plane at the element boundary). Right: Von Mises strains at the metatarsal heads, front view. The color map goes from no strain in blue to strains of 100% in red.

Table 1 . Von Mises (VM) strains measured under the metatarsal (MT) head of each toe and under the heel, internally (just below the bone structure) and externally (at the skin surface) Location Foot surface VM strain Internal VM strain

 1 Note that high strains are visible near the ankle in Figure

	5 th toe MT	2.7 %	63.3 %
	4 th toe MT	5.2 %	96.8 %
	3 rd toe MT	8.0 %	63.0 %
	2 nd toe MT	4.1 %	84.2 %
	1 st toe MT	5.1 %	43.0 %
	Heel	5.0 %	69.8 %

Acknowledgments

This work is partly funded by the French national project ANR/TecSan 2010 IDS, the French national Investments for the Future program project VHP@interactive and by the CAMI Labex.


Khaled Hammami

Jihene Elloumi
email: elloumijihene@yahoo.fr

Sami Aifa
email: sami.aifa@cbs.rnrt.tn

Christophe Drouet
email: christophe.drouet@ensiacet.fr

Hafed El Feki
email: hafed.elfeki@yahoo.fr

Synthesis and Characterization of Hydroxyapatite Ceramics Organofunctionalized with ATP (Adenosine Triphosphate)

Keywords: Hydroxyapatite, Adenosine triphosphate, Grafting, Biomaterials. Academic Discipline and Sub-Disciplines Inorganic Chemistry

Adenosine triphosphate (ATP) is the main source of energy for cells. In this article we investigated, with a possible future finality in bone tissue engineering, the original association of ATP molecules and hydroxyapatite (HAp) fine powders. We have in particular inspected, in a preliminary way, the effects of the presence of ATP molecules on the physico-chemistry of apatite compounds obtained via co-precipitation. Hybrids were produced by association of increasing amounts of adenosine 5'-triphosphate (ATP) disodium hydrate salt with hydroxyapatite, leading to so-called "organoapatites". The structure and composition of the polycrystalline hybrid compounds obtained were studied by XRD, FTIR, SEM and magic angle spinning (MAS) NMR techniques. The presence of ATP was found to alter the apatite crystallization process, leading to a departure from stoichiometry. MAS-NMR data attest to the presence of strong interactions between functional groups of ATP and apatite particles unveiling the obtainment of class II hybrids with potential future applications in cellloaded bone scaffolds.

INTRODUCTION

Hydroxyapatite (Ca10(PO4)6(OH)2, denoted here "HAp") was used extensively for different medical applications such as bone repair biomaterials due to its excellent biocompatibility with human tissues [START_REF] Oguchi | Long-term histological evaluation of hydroxyapatite cerammics in humans[END_REF]. HAp crystallizes in the hexagonal system (space group P63/m), typically with cell parameters a= b = 9.432 Å and c = 6.881 Å [START_REF] Kay | Crystal Structure of Hydroxyapatite[END_REF]. The arrangement of phosphate ions in the structure provides two types of channels (types I and II) containing calcium ions in two different environments, namely Ca (I) and Ca (II) [START_REF] Feki | structure refinements by the Rietveld method of partially substituted Hydroxyapatite[END_REF]. Ca (I) crystallographic sites then form ionic columns parallel to the c-axis, with each Ca (I) ion coordinated to nine oxygen atoms, i.e. three of each from O (1), O [START_REF] Kay | Crystal Structure of Hydroxyapatite[END_REF] and O (3) atoms. Ca (II) ions, in contrast, form equilateral Ca3-triangles, with each Ca (II) ion coordinated to 7 oxygen atoms, one O (1), one O (2), four O (3) and one oxygen from a hydroxide ion; where hydroxide ions are present on the axis of type II-channels. It may be noted that Ca (II) sites correspond, in this structure, to larger volumes than Ca (I) sites.

The properties of synthetic HAp are determined by its chemical composition (nonstoichiometry being rather common), structure, particle morphology and crystallite size which are strongly dependent on preparation methods. Nowadays, a number of protocols have been used for HAp synthesis, such as co-precipitation [4], hydrothermal reaction [START_REF] Liu | Hydroxyapatite synthesized by a simplified hydrothermal method[END_REF], sol-gel synthesis [START_REF] Liu Dean-Mo | Water-based sol-gel synthesis of hydroxyapatite : proces development[END_REF], microemulsion synthesis [START_REF] Lim | Processing of hydroxyapatite via microémulsion and émulsion routes[END_REF], mechanochemical synthesis [START_REF] Yeong | Mechanochemical synthesis of nanocrystalline hydroxyapatite from CaO and CaHPO4[END_REF], combustion synthesis [START_REF] Wag | Synthesis and sintering of nanocrystalline hydroxyapatite powder by citric acid sol -gel combustion method[END_REF] and sonochemical synthesis [START_REF] Wantae | Sonochemical synthesis of hydroxyapatite from H3PO4 solution with Ca(OH)2[END_REF].

Several works were also carried out on the preparation of organic-mineral hybrid apatite-based systems to show the influence of adsorbed organic molecules or surface functionalizations on structural and thermal properties, showing that organic species are among the most effective crystal growth inhibitors for calcium phosphates [START_REF] Van Der Houwen | The effects of organic ligands on the crystallinity of calcium phosphate[END_REF].

Recently, in order to inspect the type of interaction existing between HAp crystals and living cells, investigations have been undergone to study the surface reactivity of HAp-based nanomaterials exposing different organic functions [START_REF] Long | Structure and dynamics of hydrated statherin on hydroxyapatite as determined by solid-state NMR[END_REF][START_REF] Nagadome | Adsorption of human lysozyme onto hydroxyapatite[END_REF].

On the other hand, molecules exposing phosphonate groups have important implications/applications both in biomineralization processes and for the biomaterials industry [START_REF] Perez | Influence of Organic Phosphonates on Hydroxyapatite Crystal Growth Kinetics[END_REF]. The high reactivity of apatites with organic phosphoruscontaining molecules resides in the presence of charges (calcium and phosphate surface ions) and potentially also of amphoteric sites on their surface [START_REF] Willis | Synthesis and Characterization of Nanoapatites Organofunctionalized with Aminotriphosphonate Agents[END_REF]. These surface ions/sites may then enable interactions with organic additives to form hybrid apatite-based materials [START_REF] Portier | Inoganic-organic-hybrids as precursors to functional materials[END_REF], while the organic moieties can change the surface chemical properties of the solid (hydrophilicity, reactivity), as well as structural and physical properties. The presence of grafted (bio)molecules can open many perspectives in the development of applications, especially in the field of bone tissue engineering (a review is accessible on ref. [START_REF] Gomez-Morales | in surface tailoring of inorganic materials for biomedical application[END_REF]). Adenosine 5'-triphosphate (ATP) is the most common energy carrier for cells [START_REF] Chao | Microwave-assisted hydrothermal rapid synthesis of hydroxyapatite nanowires using adenosine 5'-triphosphate disodium salt as phosphorus source[END_REF]. It is a nucleotide constituted of adenosine (association of adenine and ribose) and of three phosphate groups linked together in a linear arrangement (the PO4 groups are denoted "α", "β" and "γ", alpha being the closest to the ribose). How extracellular ATP may interact with cells and living tissues is still a matter of concern and of questions, whether for cancerous or non-cancerous sites.

Recently, nucleic acid and nucleotide-mediated synthesis of inorganic nanoparticles (not of apatitic nature) was reported [START_REF] Berti | Nucleic acid and nucleotide-mediated synthesis of inorganic nanoparticles[END_REF]. For apatitic systems, association with ATP may find future applications in the field of tissue engineering (and/or allow one to better understand cellular activity around bone cancer sites), with the goal to potentially modify advantageously the local cellular activity by providing "additional" ATP in the vicinity of cells. The objective of the present contribution, dedicated to preliminary physico-chemical investigations on this topic, is to demonstrate the possibility to prepared hybrids associating ATP molecules (in the form of their sodium salt Na2ATP) and hydroxyapatite, via a rapid and easy precipitation method; and to inspect the effects of different Na2ATP contents on physico-chemical characteristics of the hybrids.

EXPERIMENTAL

Material preparation

HAp was prepared by a wet chemical method previously reported [START_REF] Hammari | Chemical modification of porous calcium hydroxyapatite surfaces by grafting phenyl phosphonic and phenylphosphite acids[END_REF]. Briefly, 1.67 mol of Ca(NO3)2,4H2O and 1 mol of (NH4)2HPO4, corresponding to the HAp stoichiometric molar ratio of Ca/P = 1.67, were independently dissolved in distilled water at room temperature. The phosphate-containing solution was quickly added to the calcium solution using 10 ml of concentrated ammonia solution at pH=8 and the mixture was stirred for 1 h at 25 °C. The final gel was aged for 24 h at 80°C and left to mature for a total of 16 hours. It was then filtered and thoroughly washed with distilled water. The recovered precipitate was dried overnight at 80 °C. So-called "organoapatites", associating HAp and adenosine 5'-triphosphate disodium hydrate salt (Na2ATP) and denoted "Hap -xNa2ATP", were prepared by a modified protocol, using a mixture of (1-x/3) mol of (NH4)2HPO4 and x/3 mol of Na2ATP with x= 2.5%, 5%, 10% was first dissolved in a 100 ml volume of an acetone-water mixture (50:50%, v/v) at 25°C,during this process, the pH value was maintained at 7.4 by slow addition of 1M NH4OH. The mixture was stirred in water for 30 min at room temperature. The "phosphate" solution (containing both inorganic phosphate ions "Pi" and phosphate groups from the ATP molecules) was added rapidly to the calcium solution, then stirred for 1 h and the suspension was homogenized by ultrasonic and mechanical agitation for 24 h at 37°C. The gel was decanted, matured during 16 h at 80 ° C and dried at 80° C in a conventional oven under air atmosphere . Calcination at 900 °C in air was also carried out when appropriate (see text).

Physico-chemical characterization

The powder samples were subjected to XRD (after calcination at 900°C in air), and to multinuclear solid state MAS-NMR spectroscopy, FTIR measurements and scanning electron microscopy (SEM) observations:

XRD patterns were recorded using a diffractometer Philips PW 3710 in the 2Ɵ range 5-60° using CuKα1 (λ=1.54051 Å) radiation.

Morphology and particle dimensions/agglomeration state of HAp modified by different amounts of Na2ATP were followed using a scanning electron microscope (JSM-6301, JEOL Inc.) at a voltage of 25 kV and a current of 100 mA. Solid-state 1 H and 31 P NMR spectra were recorded using MAS mode.

An NMR BRUKER 300MHz (7.1T) spectrometer was used, with cross polarization (CP) and magic-angle spinning (MAS) at 63.633 MHz (proton and phosphorus Larmor frequencies of 399.5 and 161.7 MHz, respectively). The MAS rotation rate corresponded to 7 kHz. The chemical shift for phosphorus is reported with respect to a 85% solution of phosphoric acid (H3PO4) as external standard; and tetramethylsilane (TMS) was used as standard for 1 H NMR. High-power proton decoupling technique was applied during the recording of 31 P NMR spectra.

FTIR (JASCO FT-IR 420) was used to study the vibrational features of the compounds. 1 mg of dried sample powder was compacted with 0.1 -0.15 g potassium bromide using a hydraulic pressure. For each spectrum, 32 scans between 400 and 4000 cm1 were recorded.

RESULTS AND DISCUSSION

XRD diffraction

Typical XRD profiles of HAp ceramics (calcined at 900°C in air) synthesized here and of crystalline apatite hybrids with ATP (HAp-x%Na2ATP) are shown in Fig. 1. The diffraction pattern relating to "pure" HAp confirms the apatitic nature of this compound. This phase crystallizes in the hexagonal system (space group: P63/m) with the following unit cell parameters: a = 9.426(3) Å; c = 6.883(5) Å. In the presence of ATP in the precipitation medium however, the calcination at 900°C allowed the partial conversion of the apatite phase contained in hybrids HAp-5%Na2ATP and HAp-10% Na2ATP particles into β-tricalcium phosphate Ca3(PO4)2 (β-TCP). In particular, the diffraction peaks of β-tricalcium phosphate appear at 2theta values of 27.78, 30.61 and 33.76°. Interestingly, the extent of this conversion is found to increase with the associated Na2ATP content. Such Hap/β-TCP mixtures are generally obtained when calcining nonstoichiometric apatites, i.e. corresponding to Ca/P ratios lower than 1.67. These findings thus suggest that, for increasing Na2ATP amounts in the precipitating mixture, the apatite phase obtained increasingly departs from stoichiometry. These findings can most probably be related to an inhibitory effect of ATP molecules on the apatite crystal growth process.

31 P-MAS NMR spectroscopy investigation

In order to obtain further information on the bonding of ATP molecules onto hydroxyapatite, 31 P and 1 H MAS-NMR measurements were performed. In a first stage, 31 P MAS-NMR analyses (on samples dried at 80°C) were carried out (Fig. 2). As expected for a hydroxyapatite phase, the 31 P NMR signal for pure HAp gave rise (see Figs. 2a,2b) to a single phosphorus peak around 2.8 ppm relative to 85% H3PO4 (experimental value: 2.83 ppm), similar to the phosphorus peak found in other synthetic apatites as well as in bone and teeth [START_REF] Wu | Phosphate ions in bone:Identification of a calcium-organic phosphate complex by 31P solid-state NMR spectroscopy at early stages of mineralization[END_REF][START_REF] Wu | Nuclear Magnetic Resonance Spin-Spin Relaxation of the Crystals of Bone, Dental Enamel, and Synthetic Hydroxyapatites[END_REF]. It relates to the position of apatitic phosphorus signal as reported in previous works [START_REF] Fujiwara | adsorption of myoglobin onto various synthetic hydroxyapatite particule[END_REF][START_REF] Bouhaouss | Mixed ionic conductivities in sodium fluoroapatites[END_REF]. It may be noted that, in this study, no protonated phosphate (e.g. HPO4 2ions giving rise to an additional contribution typically around 1.4 [START_REF] Rothwell | High-resolution variable-temperature 31P NMR of solid calcium phosphate[END_REF]) has been clearly detected in any sample. In contrast, in addition to this main peak, the 31 P NMR data obtained for the HAp-x%Na2ATP hybrids pointed out the presence of secondary signals in the form of broad resonances clearly detectable at -7.16 ppm and -21.3 ppm (see Fig. 2c). These positions remind those observed for 31 P in pure Na2ATP (see Fig. 2d), respectively at -7.76 (phosphate "γ") and -19.63 ppm (phosphate "β"). The slight but detectable variation in chemical shifts of these two peaks is also interesting to remark, as this can be interpreted as the formation of HAp-Na2ATP complexes (indicative of strong interactions between HAp crystals and Na2ATP molecules) rather than a simple mechanical mixture between non-interacting entities [START_REF] Zhaoyang | First-principles study of the adsorption of lysine on hydroxyapatite (1 0 0) surface[END_REF]. These findings are also corroborated by the detection of a slight shift of the apatitic 31 P NMR signal from 2.83 down to about 2.76 ppm (Fig. 2b and Table 2) for Na2ATP contents between 0 and 10%. P nuclei, which may be related to the discussion of Vasavada et al [START_REF] Vasavada | 31P NMR lineshapes of beta-P (ATP) in the presence of Mg2+ and Ca2+: estimate of exchange rates[END_REF] on the interaction of calcium ions with nucleotides. An interaction is indeed likely to take place between surface calcium ions from apatite particles and phosphate groups of the ATP molecules. Phosphate-coordinated metal ions (e.g. Ca2+ in our case) may also bind to the N-7 nitrogen in the purine residue (innersphere-type of complex, Fig. 3a), or interact with N-7 via a water molecule forming an outersphere-type of complex (Fig. 3b). Such calcium ions-ATP interactions probably play a role in the formation of Ca2+-ATP soluble complexes during the precipitation process, and possibly also directly until the grafting of ATP molecules on the apatite surface. It is worth noting here also that no trace of secondary phase has been detected in this study on such apatite-ATP hybrids, therefore suggesting that the ATP molecules are interacting directly with the apatite surface (rather than a scenario where some insoluble calcium ATP slat would form and co-precipitate beside apatite).

1 H-MAS NMR spectroscopy investigation

In order to probe hydrogen environments in such organoapatite hybrids, 1 H MAS NMR data were recorded. The results have been added in Table 2. Fig. 4 reports the 1 H NMR spectra obtained, pointing out two intense broad peaks centered on -0.05 and 6.81 ppm, respectively.

The NMR peak near -0.05 ppm can be assigned to apatitic OH -ions from the Hap phase [START_REF] Yesinowski | High-Resolution Structural Insights into Bone: A Solid-State NMR Relaxation Study Utilizing Paramagnetic Doping[END_REF]; and the peak observed around 6.81 ppm may presumably be assigned to surface-adsorbed water on apatite crystals, as was previously found for hydroxyapatite ceramics with colloidal-like dimensions [START_REF] Miquel | Solid State NMR to Study Calcium Phosphate Ceramics[END_REF].

It is clear from Table 2 that the chemical shifts and line widths (FWHM) of the proton NMR spectra vary as we compare the profiles for different amounts of Na2ATP. The line widths (FWHM) for apatitic OH-ions increase from 0.77 ppm (for pure HAp) to 0.96, 0.93 and 0.92 ppm for the 2.5%, 5% and 10% Na2ATP-Hap hybrids respectively. As for phosphorus, such spectral modifications upon ATP association suggest the existence of strong interactions between the organic molecules and the apatite surface like in "class II" hybrids [START_REF] Sanchez | Hybrid organic-inorganic materials: a land of multidisciplinarity[END_REF] (which are defined as hybrids involving strong binding between an organic entity and an inorganic substrate, as opposed to class I hybrids which exhibit only a weak interaction). In the present case, this situation may be related to the binding of some functional groups from the ATP molecules and some accessible calcium sites from the apatite surface. As mentioned previously, Ca 2+ ions are likely to interact with phosphate groups from ATP and may present also some affinity for the aromatic site of nucleobase residues. In any case, the interaction of phosphate groups and calcium ions from the surface of apatite particles probably determines a large part of the stability of such hybrids [START_REF] Szab´ | Multinuclear NMR studies of the interaction of metal ions with adenine-nucleotides[END_REF].

The intensity of the 1 H band related to apatitic OH -is also found to decrease when the ATP content increases (Fig. 3). This is in agreement with our above findings (see Fig. 1 and related text) pointing out a greater departure from stoichiometry for apatite phases precipitated in the presence of ATP molecules. In this context, the configuration of Na2ATP molecules and the orientation of the nitrogen may play a role in mineral binding affinity [START_REF] Russell | Mechanisms of action of bisphosphonates: similarities and differences and their potential influence on clinical efficacy[END_REF]. These complex processes appear to slow down the restructuration of immature apatite crystals towards more stable states, possibly by blocking "active" growth sites of the apatite phase.

Infrared spectroscopy investigation and Scanning Electron Microscopy (SEM)

Fig. 5. shows the FTIR spectra of the various hybrids (heated at 80°C). The infrared band positions and their assignments are summarized in Table 3. These data confirm in particular the presence of ATP molecules associated to the apatite phase: broad bands appear in particular at 717 and 912 cm-1, characteristic of ATP. These two bands are assignable respectively to asymmetric and symmetric stretching vibration of the P-O-P bridge. The band at 912 cm-1 may for instance be attributed to the Pβ-O-Pγ vibration in ATP. For the ATP-apatite hybrids, additional peaks related to P-O-P from adenine-nucleotides species are also detectable at 1244 cm_1, which is characteristically generated by the ν asPO 2vibration in ATP [33-34-35].

Characteristic frequencies assigned to PO4 groups in HAp occur at 1096 cm-1, 1040 cm-1, 963 cm-1 and at 603 cm-1, 568 cm-1 as discussed in previous publications [START_REF] Hammari | Crystallinity and fluorine substitution effects on the proton conductivity of porous hydroxyapatites[END_REF][START_REF] Bouhaouss | Cationic conductivity and structural studies in the Pb8 K2-x Nax (PO4 )6 system Solid State Ion[END_REF]. A broad band attributed to water O-H stretching vibrations appears at 3500 cm-1 for all samples. It is also associated to the deformation band of water molecules that appears around 1640 cm-1. These observations point out the high hydration state of such polycrystalline apatites. Weaker bands attributable to apatitic OH-ions are also visible around 3570 and 632 cm-1. The decrease in intensity of these two contributions as the amount of ATP increases is in perfect agreement with our previous findings indicating the departure from stoichiometry of apatite in these conditions. SEM analyses (Fig. 6) indicated that the materials consist of aggregated particles. Interestingly, with 2.5% and 5% of Na2ATP, one can observe a rather open structure with individual particles easy to distinguish, whereas the density of aggregation appears to increase at 10% Na2ATP. This fact points out some effect on apatite particle aggregation played by adsorbed ATP molecules, thus likely to modify electrostatic interparticle interactions.

Concluding remarks

Apatite precipitation carried out in the presence of ATP molecules (in the form of Na2ATP hydrate salt) was investigated in this work. The obtainment of class II hybrids, closely associating the organic and inorganic subparts, was evidenced by way of various complementary analyses, especially solid state NMR. The effects of the presence of ATP-grafted molecules on the physico-chemical characteristics of the apatite phase were also investigated. In particular, a departure from stoichiometry was evidenced as the ATP content increased in the precipitating medium, which suggests a role of crystal growth inhibition played by ATP.

Fig 1 :

 1 Fig 1: X-ray diffraction pattern for HAp-x%Na2ATP synthesized by co-precipitation (after calcination at 900°C).

Fig 2 : 31 P

 231 Fig 2: 31 P NMR spectra of the organoapatite (Hap-x%Na2ATP) (dried at 80 °C): a) full spectra, b) main 31 P peak, c) zoom on basis of main peak, d) spectrum for pure Na2ATP.

Fig 3 :

 3 Fig 3: Structures for chelated ATP complexes formed by inner-(a) and outer-sphere (b) coordination of the purine residue (the central ion "•" denotes a calcium ion.

Fig 4 : 1 H

 41 Fig 4:1 H NMR spectra of the organoapatite (HAp-x%Na2ATP) (dried at 80 °C).

Fig 5 :

 5 Fig 5: Infrared spectra of HAp -x%Na2ATP with different Na2ATP contents.

Fig 6 :

 6 Fig 6: SEM images of HAp (samples dried at 80°C), (scale bar = 20 μm) and Hap-2.5% Na2ATP, Hap-5% Na2ATP, Hap-10%Na2ATP (scale bar = 20 μm) (dried at 80 °C).

. Starting synthesis conditions for the preparation of ATP-HAp organoapatites.

  Table .1 summarizes the starting synthesis conditions.

	Table1HAp-x%Na2ATP	Ca(NO3)2,4H2O	(NH4)2HPO4 (=inorganic phosphate, Pi)	Na2ATP.3H2O	Ca/P
		Mole	g	Mole	g	Mole	g
	Hap			0,100	13.2	0	0
	HAp-2.5%Na2ATP	0.1667	39.35	0.0975	12.87	2.5*1/3.10-3 0.50	1.667
	HAp-5%Na2ATP			0,095	12.54	5*1/3.10-3	1
	HAp-10%Na2ATP			0,090	11.88	10*1/3.10-3	2.01

Table 2 . 1 H and 31 P chemical shifts and MAS line widths of HAp functionalized with Na2ATP at various percent.

 2

			1 H MAS NMR					31 P MAS NMR
	Sample (dried at	Assignment	Chemical	line width	Assignment	Chemical	line width
	80°C)		shift (ppm)	(ppm)				shift (ppm)	(ppm)
		OH-	-0.05	0.77
	HAp	Surface H2O	6.81	2.18	PO4	3-	2.83	0.90
	HAp-2.5%Na2ATP	OH-	-0.24	0.96
		Surface H2O	6.34	3.33	PO4	3-	2.77	0.91
	HAp-5%Na2ATP	OH-	-0.20	0.93
		Surface H2O	6.65	2	PO4	3-	2.77	0.91
	HAp-10%Na2ATP	OH-	-0.20	0.92
		Surface H2O	6.65	1.3	PO4	3-	2.76	0.91

Table3. Assignments of the observed vibrational frequencies of HAp-x%Na2ATP hybrids (dried at 80°C).

	Assignments			Observed vibrational frequencies (cm-1)
		HAp		2.5%Na2ATP-HAp	5%Na2ATP-HAp	10%Na2ATP-HAp
		position intensity	position intensity	position intensity	position intensity
	PO4 bend ν2	475	low	475	low	474	low	473	low
	PO4 bend ν4	568	strong	568	strong	568	strong	568	strong
	PO4 bend ν4	605	strong	602	strong	602	strong	602	strong
	Structural OH	635	low	632	low	631	low	632	very low
	P-Ô-P bend νs		_	717	very low 717	low	717	average
	P-Ô-P bend νas		_	912	low	912	low	912	average
	PO4 stretch ν1	963	average 962	average 962	average 961	average
	PO4 bend ν3	1040	strong	1038	strong	1038	strong	1037	strong
	PO4 bend ν3	1096	average 1094	average 1093	average 1092	average
	H2O adsorbed (ν2) 1647	average 1642	average 1642	average 1641	average
	H2O adsorbed	3427	average 3425	average 3425	average 3439	average
	apatitic OH-	3571	low	3571	low	3570	low	3569	very low

This study adds to the existing literature on calcium phosphate apatite-based hybrids and shows the possibility to associate closely ATP molecules and apatite surfaces which may find, in the future, some biomedical applications such as in tissue engineering, by directly exploiting the presence of ATP which represents a source of energy for cells.


Vincent Luboz

Antoine Perrier

Marek Bucki

Bruno Diot

Francis Cannard

Nicolas Vuillerme

Yohan Payan

Influence of the calcaneus morphology on the risk of posterior heel ulcer creation

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Influence of the calcaneus morphology on the risk of posterior heel ulcer creation Vincent Luboz, Antoine Perrier, Marek Bucki, Bruno Diot, Francis Cannard, Nicolas Vuillerme, Yohan Payan Abstract: Forty percent of the reanimation or geriatric patients develop a pressure ulcer, of which 40 % are posterior heel ulcers. The main suspected causes are the excessive pressure intensity (leading to internal strains above 50 % for about 10 minutes) and 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]. Prevention based on daily clinical examination lacks efficiency because deep tissue injuries result from internal strains and when visual symptoms appear it is often too late for treatment. It is consequently important to monitor the internal strains on a patientspecific basis as the shape of the calcaneus in the lying patient might affect the outcome. Based on our previous foot model [START_REF] Luboz | Foot Ulcer Prevention Using Biomechanical Modeling[END_REF], we study the influence of the calcaneus shape on internal strains. The Finite Element (FE) model used here represents the soft tissues of the lower leg and is composed of four different sub-domains each modelled using a Neo Hookean material with Young moduli and Poisson ratios of 200 kPa and 0.495 for the skin, 30 kPa and 0.49 for the fat, 1 GPa and 0.495 for the tendon, and 60 kPa and 0.495 for the muscles [START_REF] Sopher | Exposure to internal muscle tissue loads under the ischial tuberosities during sitting is elevated at abnormally high or low body mass indices[END_REF], see figure. Bones are modelled as rigid solids. This leg model rests on a FE model of a cushion with three compartments of varying stiffness: under the calf, the Achilles tendon, and the heel back. Simulations carried out on 18 models each implementing a different calcaneus demonstrate the influence of this bone's morphology on the magnitude of internal strains with a standard deviation (STD) of the Von Mises strains of 19.3 percentage points. Additionally, the model was used to identify the supporting cushions' stiffness combination that minimizes the ulceration risk for each individual. The pressure magnitudes below the different sections of the lower leg can be related to the speed of the pressure ulcer creation: higher cushion stiffness below the heel lead to a mean internal strain of 56.9 % with a STD of 13.2 % and could result in ulcers in less than 10 minutes, while mild cushion stiffness maintain the internal strains between 20 and 50 % (with a mean internal strain of 27.7 % with a STD of 2.0 %) and the risk of ulcer creation around two hours [START_REF] Loerakker | The effects of deformation, ischaemia and reperfusion on the development of muscle damage during prolonged loading[END_REF]. Additionally, the maximal strains are located at the interface between the calcaneus and the fat layer when the cushion is more inflated under the heel, therefore increasing the risk of pressure ulcer when maintained for a long time.

Figure

  Figure


Antoine Perrier

Vincent Luboz

Marek Bucki

Francis Cannard

Nicolas Vuillerme

Yohan Payan

Evaluation of a musculoskeletal finite element model of the foot

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Evaluation of a musculoskeletal finite element model of the foot

Antoine Perrier, Vincent Luboz, Marek Bucki, Francis Cannard, Nicolas Vuillerme, Yohan Payan

This work aims at developing a patient-specific Finite Element (FE) model of the foot in the context of pressure ulcer prevention. Starting from our previous model {1} we propose here to improve the realism of the anatomy by introducing all the detailed bony structures, ligaments and muscles segmented from the CT and MRI exams of the same patient. The model has been developed using the 3D biomechanical simulation platform Artisynth (artisynth.org). It includes 30 bones modeled as articulated rigid-bodies connected with cables that simulate the 210 segmented ligaments in their actual positions and therefore define the articulations. The aponeurosis is modeled with fives parallel multipoint ligaments. A set of 15 extrinsic and intrinsic muscles based on linear axial muscle activation law is anatomically positioned in order to allow a natural movement of the foot. A FE mesh of the soft tissue was created by applying an automatic FE mesh generator {2} to the surfaces resulting from MRI segmentation. The FE mesh has 95883 elements and 56806 nodes .Three soft tissue layers with Neo Hookean materials (Young moduli, Poisson Ratio) were created to represent a 1mm skin layer (200, 0.485), the fat (30, 0.49) and muscle (60, 0.495) tissues. The foot model was evaluated by comparing its simulations with the pressure mat collected with the patient standing onto a Zebris FDM platform (www.zebris.de). For this, the patient weight was applied onto the foot model put in contact with anhorizontal plate (figure). Another evaluation of the model aimed at comparing the kinematics of the model resulting from simulated muscle activations with an actual 3D motion analysis of the patient foot {3} coupled with an EMG monitoring of the tibialis posterior, tibialis anterior, soleus, medial gastrocnemius, lateral gastrocnemius and peroneus muscles. Globally, this new version of the foot model provides realistic simulations in both static and dynamic frameworks. In addition to its use for pressure ulcer prevention, it could become relevant for the simulation of neuro-orthopedic surgical interventions, orthotic devices analysis or educational purposes.

Figure

  Figure


Christian Baudet
email: christian.baudet@univ-lyon1.fr

Ulisses Dias
email: udias@ic.unicamp.br

Zanoni Dias
email: zanoni@ic.unicamp.br

Length and Symmetry on the Sorting by Weighted Inversions Problem

. Despite the fact that symmetry and length of the reversed segments seem very important, they have not been considered together in any problem in the genome rearrangement field. Here, we define the problem of sorting genomes (or permutations) using inversions whose costs are assigned based on their lengths and asymmetries. We present five procedures and we assess these procedure performances on small sized permutations. The ideas presented in this paper provide insights to solve the problem and set the stage for a proper theoretical analysis.

Introduction

Among various large-scale rearrangement events that have been proposed to date, inversions were established as the main explanation for the genomic divergence in many organisms [START_REF] Darling | Dynamics of genome rearrangement in bacterial populations[END_REF][START_REF] Dias | A simulation tool for the study of symmetric inversions in bacterial genomes[END_REF][START_REF] Hannenhalli | Transforming cabbage into turnip: polynomial algorithm for sorting signed permutations by reversals[END_REF]. An inversion occurs when a chromosome breaks at two locations, and the DNA between those locations is reversed.

In some families of bacteria, an 'X'-pattern is observed when two circular chromosomes are aligned [START_REF] Dias | A simulation tool for the study of symmetric inversions in bacterial genomes[END_REF][START_REF] Eisen | Evidence for symmetric chromosomal inversions around the replication origin in bacteria[END_REF]. Inversions symmetric to the origin of replication (meaning that the breakpoints are equally distant from the origin of replication) have been proposed as the primary mechanism that explains the pattern [START_REF] Eisen | Evidence for symmetric chromosomal inversions around the replication origin in bacteria[END_REF]. The justification relies on the fact that one single highly asymmetric inversion affecting a large area of the genome could destroy the 'X'-pattern, although short inversions may still preserve it.

Darlink, Miklós and Ragan [START_REF] Darling | Dynamics of genome rearrangement in bacterial populations[END_REF] studied eight Yersinia genomes and added evidence that symmetric inversions are "over-represented" with respect to other types of inversions. They also found that inversions are shorter than expected under a neutral model. In many cases, short inversions affect only a single gene, as observed by Lefebvre et al. [START_REF] Lefebvre | Detection and validation of single gene inversions[END_REF] and Sankoff et al. [START_REF] Sankoff | The distribution of inversion lengths in bacteria[END_REF], which contrasts with the null hypothesis that the two endpoints of an inversion occur by random and independently.

Despite the importance of symmetry and length of the reversed segment, both have been somewhat overlooked in the genome rearrangement field. Indeed, the most important result regarding inversions is a polynomial time algoritm presented by Hannenhalli and Pevzner [START_REF] Hannenhalli | Transforming cabbage into turnip: polynomial algorithm for sorting signed permutations by reversals[END_REF] that considers an unit cost for each inversion no matter its length or symmetry. When gene orientation is not taken into account, finding the minimum number of inversions that transform one genome into the other is a NP-Hard problem [START_REF] Caprara | Sorting permutations by reversals and Eulerian cycle decompositions[END_REF].

Some results have considered at least one of the concepts. There is a research line that considers the total sum of the inversion lengths as the objective function of a minimization problem. Several results have been presented both when gene orientation is considered [START_REF] Arruda | Heuristics for the sorting by length-weighted inversions problem on signed permutations[END_REF][START_REF] Swidan | Sorting by lengthweighted reversals: Dealing with signs and circularity[END_REF] and when it is not [START_REF] Arruda | Heuristics for the sorting by length-weighted inversion problem[END_REF][START_REF] Bender | Improved bounds on sorting by length-weighted reversals[END_REF][START_REF] Pinter | Genomic sorting with length-weighted reversals[END_REF]. Recently, Arruda et al. [START_REF] Arruda | Heuristics for the sorting by length-weighted inversions problem on signed permutations[END_REF] developed a randomized approach that starts with a scenario with the minimum number of inversions, but allows the number of inversions to increase if the outcome is a reduction in the total sum of the inversion lengths.

Regarding symmetry, the first results were presented by Ohlebusch et al [START_REF] Ohlebusch | The median problem for the reversal distance in circular bacterial genomes[END_REF]. Their algorithm uses symmetric inversions in a restricted setting to compute an ancestral genome and, therefore, is not a generic algorithm to compute the rearrangement distance using only symmetric inversions. In 2012, Dias et al. presented an algorithm that considers only symmetric and almost-simmetric inversions [START_REF] Dias | Sorting genomes using almostsymmetric inversions[END_REF]. They later included unitary inversions to the problem and provided a randomized heuristic to compute scenarios between two genomes that uses solely these operations [START_REF] Dias | Greedy randomized search procedure to sort genomes using symmetric, almost-symmetric and unitary inversions[END_REF].

Here we propose a new genome rearrangement problem that combines the concepts of symmetry and length of the reversed segments. Whereas previous works restricted the set of allowed operations by considering only inversions that satisfy constrains like symmetry or almost-symmetry [START_REF] Dias | Greedy randomized search procedure to sort genomes using symmetric, almost-symmetric and unitary inversions[END_REF][START_REF] Dias | Sorting genomes using almostsymmetric inversions[END_REF], here we allow all possible inversions. The problem we are proposing aims at finding low-cost scenarios between genomes when gene orientation is not taken into account. The results obtained are the first steps in exploring this interesting new problem.

Definitions

Formally, a chromosome is represented as a n-tuple whose elements represent genes. If we assume no gene duplication, then this n-tuple is a permutation π = (π 1 π 2 . . . π n ), 1 ≤ π i ≤ n and π i = π j ↔ i = j. Because we focus on bacterial chromosomes, we assume permutations to be circular, and π 1 is the first gene after the origin of replication.

The inverse of a permutation π is denoted by π -1 , for which π -1 πi = i for all 1 ≤ i ≤ n. The composition between two permutations π and σ is similar to function composition in such way that π • σ = (π σ1 π σ2 . . . π σn ).

An inversion ρ(i, j),

1 ≤ i ≤ j ≤ n, is a rearrangement that transforms π into π • ρ(i, j) = (π 1 . . . π i-1 π j π j-1 . . . π i+1 π i π j+1 . . . π n ).
Given two permutations α and σ, the inversion distance d(α, σ) is the size t of the minimum sequence of operations ρ 1 , ρ 2 , . . . , ρ t such that α

•ρ 1 •ρ 2 •. . .•ρ t = σ.
Let ι = (1 2 . . . n) be the identity permutation, sorting a permutation π = (π 1 π 2 . . . π n ) is the process of transforming π into ι and the inversion distance between them is denoted as

d(π, ι) = d(π). Note that σ • σ -1 = σ -1 • σ = ι n . Therefore, the inversion distance d(α, σ) is equivalent to transform a permutation π into a permutation ι if we take π = σ -1 •α, because d(α, σ) = d(σ -1 •α, σ -1 •σ) = d(π, ι) = d(π)
. Therefore, we hereafter consider the sorting by inversions problem which aims at finding the sorting distance d(π) for an arbitrary permutation π.

The following functions can be applied to identify any element i in the permutation π. Position:

p(π, i) = k ⇔ |π[k]| = i, p(π, i) ∈ {1, 2, . . . , n}. Slice: slice(π, i) = min{p(π, i), n -p(π, i) + 1}, slice(π, i) ∈ {1, 2, . . . , ⌈ n
2 ⌉}. We will now define the cost function for each inversion ρ(i, j) and then we explain two cases that arise. Our cost function is given by cost(ρ(i, j)) = |slice(ι, i)slice(ι, j)| + 1. Figure 1 illustrates both cases.

Case 1: i, j ≤ ⌈ n 2 ⌉ or i, j ≥ ⌈ n 2 ⌉.
In this case, the cost function can be simplified to cost(ρ(i, j)) = abs(i-j)+1, which means that it is proportional to the number of elements in the reversed segment. This cost is what one would expect from a length-weighted inversion distance in such a way that larger inversions cost more than short inversions. Case 2: i > ⌈ n 2 ⌉ and j < ⌈ n 2 ⌉, or j > ⌈ n 2 ⌉ and i < ⌈ n 2 ⌉. In this case, the cost function is penalizing the asymmetry instead of the number of elements in the reversed segment. In effect, if the inversion ρ(i, j) is perfectly symmetric (meaning that i and j are equally distant from the origin of replication), then the cost is given by cost(ρ(i, j)) = 1.

(a) i, j ≤ ⌈ n 2 ⌉ or i, j ≥ ⌈ n 2 ⌉ and (b) i > ⌈ n 2 ⌉ and j < ⌈ n 2 ⌉, or j > ⌈ n 2 ⌉ and i < ⌈ n 2 ⌉.

Algorithms

This section presents several algorithms that take advantage from the characteristics of the cost functions.

Greedy Function

This algorithm uses a greedy function in order to estimate how good an inversion might be. The greedy function uses two concepts: breakpoints and slicemisplaced pairs. Definition 1. Breakpoints: let us extend the permutation π to include the elements π 0 = 0 and π n+1 = n + 1. The pair π i and π i+1 , for

0 ≤ i ≤ n, is a breakpoint if |π i+1 -π i | = 1.
We use b(π) to represent the number of breakpoints in a permutation and

∆ b (π, ρ) = b(π • ρ) -b(π)
to represent the variation in the number of breakpoints caused by an inversion ρ.

Definition 2. Slice-Misplaced pairs: let π i , π j be two elements in π such that slice(π, π i ) > slice(π, π j ). We say that π i , π j corresponds to a slice-misplaced pair if slice(ι, π i ) < slice(ι, π j ). We use m(π) to represent the number of slicemisplaced pairs in π and ∆ m (π, ρ) = m(π • ρ)m(π) to represent the variation in the number of misplaced-pairs caused by an inversion ρ.

The identity permutation ι is the only permutation with no breakpoints. Therefore, an inversion that decreases the number of breakpoints indirectly leads up to the identity permutation. The identity permutation also has no misplaced pairs, so we will create a function that combines both concepts in order to estimate how good an inversion is.

Since an inversion can affect only two breakpoints, we know that ∆ b (π, ρ) ∈ {-2, -1, 0, 1, 2}. The variation in the number of slice-misplaced pairs ∆ m is not well-behaved like ∆ b . Thus, we decided to favour breakpoints reduction in our greedy function:

h(π, ρ) = ∆ b (π, ρ) + ∆m(π,ρ) n 2
. Therefore, the benefit of an inversion is given by δ(π, ρ) = h(π,ρ) cost(ρ) . We construct a sequence of inversions that sorts π by iteratively adding an inversion with the best benefit among all possible inversions.

Left or Right Heuristic

We first divide the elements in π in two groups. The first group refers to the elements that are in slices classified as sorted and the second group comprises those elements that are in unsorted slices. A slice s is in the sorted group if p(π, s) = s, p(π, ns + 1) = ns + 1, and the slices {1, 2, . . . , s -1} are also in the sorted group. Otherwise, s is in the unsorted group.

First, the Left or Right heuristic selects the least slice in the unsorted group. Then, we determine the element that should be moved first to that slice: the left or the right. The left (right) side is composed of the elements which are in positions that have indices bigger (lower) than the middle position.

To make this choice, we use an auxiliary function f 1 which determines the total weight to put a given element in its right place. The function considers only inversions ρ such that cost(ρ) ≤ 2. The minimum number of such inversions that is necessary to move the element onto its right position is counted and the weight is computed.

If the slice has only one element that does not belong to it, we just find the element that should be in that slice and perform an inversion to place it closer to its right position. We use inversions ρ whose cost(ρ) ≤ 2 following the principle of always applying non-expensive movement.

If the slice has two elements that are misplaced, we compute the total weight to place the left element and the total weight to place the right element. Then, we choose the element that has the smallest total weight. In case of tie, we move the right element.

Lock Heuristic

This heuristic adds a new step that will be called before using the Left or Right heuristic. We first search for a lock inversion, which puts an element directly into its final position in the unsorted slice with the least value.

We also consider the possibility of indirect lock inversions that put the element in its final position in two steps: first, it moves the element to its final slice, but in the opposite side. Then, it moves the element to its right place with a perfectly symmetric inversion. In the case where the opposite side is already locked (opposite side has already the right element), we consider moving the element one slice above and then applying an almost-symmetric inversion.

An auxiliary function f 2 is used to evaluate the inversions. This function takes the elements π i and π j that are in the endpoints of the inversion ρ(i, j) and uses the function f 1 to compute the total weight to put them in their final positions. We compute the total weight before and after applying the inversions and compute the gain. If the gain is positive (i.e., the total weight decreased) and this gain is equal to or bigger than the inversion weight, ρ(i, j) is considered a valid inversion.

The heuristic evaluates the inversions for placing the right and left elements, and then applies the inversion with least cost. If a lock inversion does not exist, then the Left or Right heuristic is used.

Reduce Heuristic

This heuristic computes all possible inversions that remove at least one breakpoint and is considered valid by the function f 2 as described in Section 3.3. We filter the inversions having minimum weight according to f 2 , and if more than one inversion remains, we select one that removes two breakpoints. If no inversion reduces the number of breakpoints, we apply the Left or Right heuristic (Section 3.2). This new heuristic is called Reduce heuristic.

It is also possible to combine the Reduce heuristic with the lock heuristic to generate other approach called Lock or Reduce. This new approach searches for a valid lock inversion ρ 1 and for a valid inversion ρ 2 which reduces the number of breakpoints. If just one of ρ 1 or ρ 2 exists, the heuristic simply applies it. If both exist, the heuristic selects the one that has the smallest weight (in case of tie, apply ρ 1 ). Finally, if none exists, the Left or Right heuristic is used.

We generated a dataset with all possible instance π, such that |π| ≤ 12. Since the permutations in the dataset are small, we were able to compute a minimum cost solution for comparison purposes. The minimum cost solution was calculated using a graph structure G n , for n ∈ {1, 2, . . . , 12}. We define G n as follows. A permutation π is a vertex in G n iff π has n elements. Let π and σ be two vertices in G n , we build an edge from (π, σ) iff there is an inversion ρ that transforms π into σ. The weight assigned to this edge is cost(ρ). Finally, we calculate the shortest path from ι to each vertex in G n using a variant of Dijkstra's algorithm for the single-source shortest-paths problem. This variant gives us the minimum cost to sort permutations in G n , as well as an optimum scenario of inversions.

Let heuristic cost(π) be the cost for a sorting sequence of π and cost(π) be the optimum cost, we can compute the approximation ratio as heuristic cost(π) cost(π)

. The first graph in Figure 2 shows how often each heuristic returns the optimum cost. The second and third graphs exhibit, respectively, the average and maximum ratios observed among permutations of same size. The graphs in figures 2(a), 2(b) and 2(c) are maximum or average values. Thus, they may not answer the question: for a single instance π, is there any algorithm that is likely to provide the best answer? The fourth graph discuss this question by assessing the number of times each algorithm provides the least costly sequence. We observe that Greedy leads to the best results. In fact, this heuristic is consistently better than the other heuristics in every aspects we plot in Figure 2. The heuristics Reduce and Lock or Reduce are very similar when we consider how often each one returns the optimum cost and the average ratio. However, the Lock or Reduce heuristic seems more appropriate for those permutations that are hard to optimize, since the maximum ratio curve for Lock or Reduce (shown in the third graph) is significantly better than that for Reduce.

We added new curves to see if one could take advantage of using more than one heuristic other than Greedy. The curve labeled as A+E selects for each instance the less costly result between those produced by the Greedy heuristic and the Lock or Reduce heuristic. As we can see, using both heuristics is consistently better than using solely the Greedy heuristic in every aspect we are studying.

A final test checks if any profit is gained from running all possible heuristic, which is what we consider in the curve All. The increase in performance produced by All is consistent when compared to A+E. In every graph, the All curves are far from the others, which leads to the conclusion that using a combination of more than one heuristic accomplishes satisfactory results.

Conclusions

In this work, we have defined a new genome rearrangement problem based on the concepts of symmetry and length of the reversed segments in order to assign a non-unit cost for each inversion. The problem we are proposing aims at finding low-cost scenarios between genomes when gene orientation is not taken into account. We have provided the first steps in exploring this problem.

We presented five heuristics and we assessed their performances on small sized permutations. The ideas used in order to develop these heuristics together with the experimental results set the stage for a proper theoretical analysis.

As in other inversion sorting problems, we would like to know the complexity of determining the distance between any two genomes using only the operations we defined. That seems to be a difficult problem that we intend to keep studying. We plan to design approximation algorithms and more effective heuristics.

  cost is proportional to the number of elements in the re-asymmetry between i and j. A symmetric inversion costs 1 unit

Fig. 1 .

 1 Fig. 1. Effect of the cost function when (a) i, j ≤ ⌈ n 2 ⌉ or i, j ≥ ⌈ n 2 ⌉ and (b) i > ⌈ n 2 ⌉ and j < ⌈ n 2 ⌉, or j > ⌈ n 2 ⌉ and i < ⌈ n 2 ⌉.

Fig. 2 .

 2 Fig.2. In (a) we show how often each heuristic returns a minimum cost solution. In (b) and (c) we show the average and the maximum ratio, respectively. In (d) we show how often each heuristic succeeds in providing the best answer among all the heuristics.

Acknowledgments

This work was supported by a Postdoctoral Fellowship from FAPESP to UD (number 2012/01584-3), by project fundings from CNPq to ZD (numbers 306730/ 2012-0, 477692/2012-5 and 483370/2013-4), and by French Project ANR MIRI BLAN08-1335497 and the ERC Advanced Grant SISYPHE to CB. The authors also thank the Center for Computational Engineering and Sciences at Unicamp for financial support through the FAPESP/CEPID Grant 2013/08293-7.


Paul Alain Nana

Pierre Ngassam

Geneviève Bricheux

Philippe Bouchard

Bernard Vigues

Télesphore Sime-Ngando

Zéphyrin Fokam

Description of two new species of Sicuophoridae and Nyctotheridae (Heterotrichina), endocommensal in the rectal ampulla of Bufo regularis (Amphibia: Anura) from the Northwest of Cameroon

Keywords: Bambili, Bambui, Bufo regularis, Nyctotheroides bambuiensis n. sp., Prosicuophora cyclostomatus n. sp

scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The amphibian digestive tract harbors in its posterior part a diverse fauna of ciliated Protozoa, among others including a great number of heterotrichs that have been examined by many investigators [START_REF] Mackinnon | An introduction to the study of Protozoa[END_REF][START_REF] Albaret | Etude systématique et cytologique sur les ciliés hétérotriches endocommensaux. Mé moires du Musé um National[END_REF]de Puytorac and Öktem, 1967;[START_REF] Shorr | Populational changes of the enteric protozoans Opalina spp. and Nyctotherus cordiformis during the ontogeny of anuran tadpoles[END_REF][START_REF] Vojtkova | Parasites of the frogs and toads in Europe. Part I: Protozoa[END_REF].

The first study of heterotrich cilates in central Africa was conducted by de Puytorac and Öktem (1967) on nyctotherans of the genera Nyctotherus [START_REF] Leidy | New genera of entozoa[END_REF] and Prosicuophora (de [START_REF] Puytorac | Observations cytologiques sur les Nyctothères des Genres Nyctotherus Leidy et Prosicuophora N. Gen., ciliés parasites de batraciens anoures du Gabon[END_REF] collected in anuran amphibians from Gabon. In Cameroon, similar studies have been carried out by [START_REF] Affa | Nyctotheroides nouveaux ou peu connus du Cameroun (première série)[END_REF][START_REF] Affa | Nyctotheroides nouveaux ou peu connus du Cameroun (deuxième série)[END_REF], [START_REF] Affa | Progrès récents dans la connaissance des Nyctothères (Protozoaires, Ciliés Hétérotriches) associés aux anoures[END_REF]. However, works on the nyctotheran ciliates in amphibians in Cameroon are scarce. For this reason, surveys and comparisons of the ciliate fauna of the amphibian digestive system from the hosts collected in different regions should provide new information about the ciliated Protozoa fauna of Cameroon. Here we describe two new species, Prosicuophora cyclostomatus n. sp. and Nyctotheroides bambuiensis n. sp., obtained from the Bufo regularis in Tubah district of the Northwest region of Cameroon.

Material and methods

Our study areas in the north west of Cameroon (Bambili and Bambui in Tubah district) were situated within the latitude 5 0 59' 0''N, and the longitude 10 0 15' 0''E, at an altitude of about 1,350 m (Fig. 1).

Rectal ampulla samples were dissected from Bufo regularis collected in the vicinity of Bambili and Bambui in 2012. Portions of the rectal ampulla were opened either in physiological Ringer's solution or in commercial mineral water (Volvic™ in France or Supermont™ in Cameroon). Ciliates were observed under a stereomicroscope and the actively swimming ones were selected by micropipette and first examined in vivo. Then the ciliates were stained with Fernandez Galiano's pyridinated silver impregnation (1994). All of the cell measurements were made with a calibrated ocular micrometer. To obtain morphometric data, groups of 30 ciliates of each species were examined. Drawings of the cells were made with the aid of a camera lucida attached to a Wild M20 microscope. A digital camera was used for light micrographs. Identification and classification was done according to the key of [START_REF] Albaret | Etude systématique et cytologique sur les ciliés hétérotriches endocommensaux. Mé moires du Musé um National[END_REF].

Results

PROSICUOPHORA CYCLOSTOMATUS SP. N.

This species occurs in the rectal ampulla of Bufo regularis where it lives in cohabitation with various species of the genera Nyctotheroides, Opalina and Protoopalina.

Diagnosis. Size about 101 × 74 µm. Body ovoid. Anterior pole pointed, posterior pole rounded. Peristome starts near the anterior pole as a spiral, extends posteriorly, and joints the infundibulum in the equatorial part of the body. Infundibulum bends inwards forming nearly one spire in the central part of the cell and ends with a forked cytopharynx. Macronucleus elongated, with a pointed posterior end and a wide anterior part, located in the anterior half of the body, close to the infundibulum. Globulous micronucleus surmounts the macronucleus. The cilia pattern shows at the poles two suture lines extending both to the left and right surfaces and forming two secant systems.

Type host. Bufo regularis Reuss, 1833. Type locality. Bambili and Bambui, Tubah District, Bamenda, Cameroon.

Habitat. Rectal ampulla. Type specimens. Slides with paratypes are deposited at the Department of Biology (Higher Teacher Training College, The University of Bamenda, Cameroon).

Description. The ciliate ovoid with a broadened posterior end; anterior end slightly pointed (Fig. The study area (Helvetas Cameroon, 2002;modified, 2013).

Zéphyrin Fokam et al.

2). Body measurements and proportions: length, 101 ± 13.4 µm; width, 74.6 ± 6.3 µm (Table 1). Nuclear apparatus in the anterior part of the body. Macronucleus, 22.3 ± 3.49 × 14.4 ± 4.2 µm, major axis parallel to the axis of the cell. Macronuclear shape, irregular lozenge-like with a tapered posterior end, typically characteristic of the species. Micronucleus spherical, located in the concavity at the anterior side of the macronucleus (Fig. 2).

Buccal apparatus begins externally with a welldeveloped oblique peristome (34.5 ± 0.1 µm) and continues internally with the infundibulum (26.8 ± 1.4 µm). The terminal part of the infundibulum is prolonged by a spiral cytostome and a forked cytopharynx (Fig. 3A). In most of the specimens studied, the peristome and infundibulum are linked by a joint forming an obtuse angle (130.1 ± 13.4)°. However, in a few cases this angle is flattened instead (Fig. 3B). The vestibular opening (oral opening) is close to the equatorial plane of the cell. The buccal kinetosomes form adoral membranelles composed of short striated kineties covering the peristome and the infundibulum.

The 178-206 somatic kineties are arranged in meridian and tight bipolar rows throughout the body. The ciliary topography shows on the right side two suture lines prolonged by secant systems (Fig. 4). The apical suture line is shorter and located in front of the peristome. Posteriorly, the caudal suture line is slightly discrete.

NYCTOTHEROIDES BAMBUIENSIS SP. N.
This ciliate co-inhabits with the previous species the rectum of Bufo regularis.

Diagnosis. Body shape oval. Both poles rounded, but the posterior slightly larger. Peristome begins at a short distance from the apex and extends backward in a continuous circling flat curve leading to the oral opening. Infundibulum, very long, runs laterally in the median plane of the cell, then twists and Notes. AIPe -angle infundibulum-peristome in degree, In -infundibulum, Max -maximum, Min -minimum, Mn -macronucleus, MnPAmacronucleus posterior angle in degree, Pe -peristome, SD -standard deviation; N=30. Description. Body shape oval, the apical part narrower than the bottom. Cell size: 146.50 ± 4.24 µm long and 108.90 ± 24.75 µm in the widest part (Fig. 5, Table 2). Left side flattened with a hollow structure constituting an oral groove; right side convex.

The buccal apparatus starts close to the apex with a curved peristome (Figs 5,6). The peristome is largely opened in the equatorial part of the cell by an oral opening. The oral groove, relatively long, appears as a U-shaped infolding forming the infundibulum (Fig. 6).

Nuclear apparatus located in the anterior half of the body, above the infundibulum. Macronucleus, kidney-shaped, about 58.10 ± 5.66 µm long and 29.12 ± 3.53 µm in the widest part, left side slightly pointed, right side broader. Macronucleus oriented transversally to the cell axis with the convex side facing the top of the body and the flattened side facing the bottom. Micronucleus ellipsoid, located in a depression in the middle of the flattened side (Figs 5,6).

Buccal apparatus divided into two parts: the peristome and infundibulum. Peristome begins close to the apex on the left lateral face and ends with an oral opening in the mid-body. Infundibulum U-shaped, begins at the oral opening, plunges into the cytoplasm, bends and runs backwards parallel to its anterior part, ends with a cytostome prolonged by a short cytopharynx (Fig. 6). Coiled thick structure, the monitorium, supports the internal base of the infundibulum bend.

Ciliary topography: about 126-162 bipolar kineties, unequally distributed throughout the body. The preoral secant system and the caudal secant system are extensions of the anterior and posterior suture lines respectively. In the axial region of the left side, the confrontation of kineties forms posteriorly, a very short caudal secant system (Fig. 7).

Discussion and conclusion

According to its general morphology (shape, buccal and nuclear apparatus), the first ciliate described in this work is typical to the genus Prosicuophora with P. basoglui as a typical species described by de Puytorac and Oktem (1967). The pattern of the kineties is classical, with two secant systems, the preoral and the caudal one. However, the cell described here differs drastically from P. basoglui, as it lacks a karyophore (a sort of a sac-like membrane enclosing the macronuleus). Additionally, the infundibulum of P. basoglui bears numerous short fibers, distinguishing it from our specimen.

The first specimen also shows some similarities with Nyctotheroides spirotomatus examined by [START_REF] Albaret | Etude systématique et cytologique sur les ciliés hétérotriches endocommensaux. Mé moires du Musé um National[END_REF] in terms of morphology (the infundibulum structure). However, our ciliate differs from N. spirotomatus in the structure of the nuclear apparatus, notably the absence of karyophore. Presumably, the differences mentioned above may be due to the difference in hosts, since P. basoglui has been recorded in Bufo superciliaris.

Nyctotheroides bambuiensis possesses typical characteristics of the genus Nyctotheroides (general morphology, pattern of the buccal and somatic kineties) as stated by [START_REF] Affa | Nyctotheroides nouveaux ou peu connus du Cameroun (deuxième série)[END_REF]. This ciliate shares certain features with the species N. heterotomus described by [START_REF] Affa | Nyctotheroides nouveaux ou peu connus du Cameroun (première série)[END_REF], namely, the shape of the cell, the position of the oral opening, the form of the buccal apparatus and the nuclear apparatus. However, the presence of a monitorium, a long and curved infundibulum, a voluminous cytostome and a short cytopharynx set these two species apart. The second studied ciliate also resembles N. cordifomis [START_REF] Albaret | Etude systématique et cytologique sur les ciliés hétérotriches endocommensaux. Mé moires du Musé um National[END_REF] in the body shape, the position of the nuclear apparatus and its shape, the position of the oral opening, the buccal apparatus, and buccal and somatic ciliature, determining two suture lines (apical and caudal). Alongside with these similarities, there are certain distinctions between our specimen of Nyctotheroides and N. cordiformis. The oral or vestibular opening, which is in the median part of the cell in N. bambuiensis, is closer to the posterior end in N. cordiformis. The area occupied by infundibulum in N. bambuiensis is larger than in N. cordiformis. The curve formed by the infundibulum is wider in N. cordiformis; in this species the infundibulum, after passing the middle part of the cell and forming a wide curve, turns to the posterior end. However, in our ciliate, the infundibulum forms a small curve somewhere near the dorsal side of the cell and thereafter turns towards the ventral side of the cell, where it is extended by a voluminous cytopharynx. Another distinct feature to discriminate the two species is the presence of a monitorium supporting the curve of the infundibulum in N. bambuiensis and its absence in N. cordiformis. One of the main distinguishing characters of N. cordiformis is the caryophore, which is absent in N. bambuiensis.

The differences we observed could be explained by some adaptations or resistance to the environmental factors of the hosts, and from evolutionary perspective, such situation can create adaptive radiation leading to diversification in the same genus. This diversification can also be explained by the difference of hosts, since N. heterotomus is

Fig

  Fig 1.The study area(Helvetas Cameroon, 2002; modified, 2013).

Fig. 2 .

 2 Fig. 2. General morphology of Prosicuophora cyclostomatus n. sp. A -light micrograph after silver staining (×400), B-drawing; Ap -apical pole, Cv -contractile vacuoles, Cy-cytopyge, Cypcytopharynx, In-infundibulum, Ma-macronucleus, Mi-micronucleus, Pe-peristome.

Fig. 3 .

 3 Fig. 3. Buccal apparatus (BA) of Prosicuophora cyclostomatus n. sp. A -detail, B -linear disposition (×400); BK -buccal kineties, In -Infundibulum, Pe -peristome.

Fig. 4 .

 4 Fig. 4. Ciliary pattern of Prosicuophora cyclostomatus n. sp. CSL -caudal suture line, PrSL -preoral suture line.

Fig. 5 .

 5 Fig. 5. General morphology of Nyctotheroides bambuiensis n. sp. A -light micrograph after silver staining (×400), B -drawing; Ap -apical pole, Cy -cytopyge, Cyp -cytopharynx, Ma -macronucleus, Mi -micronucleus, Oo -oral opening.

Fig. 6 .

 6 Fig. 6. Detail of the buccal and nuclear apparatus of Nyctotheroides bambuiensis n. sp. Cypcytopharynx, Ma -macronucleus, Mi -micronucleus; Oo -oral opening.

Fig. 7 .

 7 Fig. 7. Ciliary pattern of of Nyctotheroides bambuiensis n. sp. CSL -caudal suture line, Cycytopyge, Inf -infundibulum, Pe -peristome, PrSL -preoral suture line, Oo -oral opening.

Table 1 .

 1 Morphometric characters of Prosicuophora cyclostomatus n. sp.

		Cell length (µm)	Cell width (µm)	Mn length (µm)	Mn width (µm)	Pe length (µm)	In length (µm)	MnPA (°)	AIPe (°)
	Max	140	103	31	22	45	38	48	148
	Mean	101	74.6	22.3	14.4	34.5	26.8	41.7	130.1
	Min	70	52	15	8	22	17	27	104
	SD	13.4	6.3	4.9	4.2	0.1	1.4	0.7	13.4

Table 2 .

 2 Morphometric characters of Nyctotheroides bambuiensis n. sp.

		Cell length	Cell width	Mn length	Mn width
		(µm)	(µm)	(µm)	(µm)
	Max	175	133	68	35
	Mean	146.5	108.9	58.1	29.1
	M i n	6 0	4 5	2 5	1 1
	SD	4.24	24.75	5.66	3.53

Notes. Max -maximum, Min -minimum, SD -standard deviation; N=30.

a commensal ciliate of Hylarana albolaris and the specimen described here resides in Bufo regularis. Besides that, we can hypothesise that morphological differences observed in N. bambuiensis may have developed as a climatic adaptation.

Aknowledgements

The authors are indebted to the 2013 team of young researchers in the unit of protistology (Department of Biology, Higher Teacher Training College, University of Bamenda, Cameroon) for their assistance in collecting the ciliate's host.