Alexandr Klimchik
email: a.klimchik@innopolis.ru

Anatol Pashkevich
email: anatol.pashkevich@imt-atlantique.fr

Stéphane Caro
email: stephane.caro@ls2n.fr

Benoît Furet
email: benoit.furet@univ-nantes.fr

Calibration of Industrial Robots with Pneumatic Gravity Compensators

Keywords: Industrial robot, stiffness modeling, elastostatic calibration, pneumatic gravity compensator, design of calibration experiments

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

I. INTRODUCTION

Advancements in shipbuilding and aeronautic industries demand high-precision and high-speed machining of huge hulls and fuselage components. For these tasks, industrial robots are more attractive comparing to conventional CNCmachines because of large and easily extendable workspace, capability to process complex-shape parts and high-speed motion capability. However, processing of modern and contemporary materials, which are widely used in these industries, requires high processing forces affecting robot positioning accuracy [START_REF] Zhu | An off-line programming system for robotic drilling in aerospace manufacturing[END_REF][START_REF] Guillo | Impact & improvement of tool deviation in friction stir welding: Weld quality & real-time compensation on an industrial robot[END_REF][START_REF] Denkena | Enabling an Industrial Robot for Metal Cutting Operations[END_REF]. To reduce the external force impact on the positioning accuracy, robotic experts usually apply technique that is based on the compliance error estimation via the manipulator stiffness modelling [START_REF] Dumas | Joint stiffness identification of industrial serial robots[END_REF][START_REF] Nubiola | Absolute calibration of an ABB IRB 1600 robot using a laser tracker[END_REF][START_REF] Alici | Enhanced stiffness modeling, identification and characterization for robot manipulators[END_REF][START_REF] Klimchik | Stiffness modeling for perfect and non-perfect parallel manipulators under internal and external loadings[END_REF] and relevant error compensation in the online or offline mode [START_REF] Schneider | Stiffness modeling of industrial robots for deformation compensation in machining[END_REF][START_REF] Klimchik | Compliance Error Compensation in Robotic-Based Milling[END_REF][START_REF] Klimchik | Compliance error compensation technique for parallel robots composed of non-perfect serial chains[END_REF][START_REF] Slavkovic | A method for off-line compensation of cutting force-induced errors in robotic machining by tool path modification[END_REF]. This approach is very efficient if the stiffness and geometric model parameters of the manipulator as well as the external forces are known. To estimate them, additional experimental studies are usually carried out [START_REF] Dumas | Joint stiffness identification of industrial serial robots[END_REF][START_REF] Klimchik | Identification of the manipulator stiffness model parameters in industrial environment[END_REF][START_REF] Wu | Geometric calibration of industrial robots using enhanced partial pose measurements and design of experiments[END_REF][START_REF] Hollerbach | Model Identification[END_REF], which allow user to obtain an extended geometric model suitable for compliance error compensation.

Another practical solution is based on enhancement the robot stiffness by means of increasing the link cross-sections. However, it leads to increasing of the robot components masses causing additional end-effector deflections, which are usually reduced by means of different types of gravity com-pensators. However, integration of mechanical compensators into manipulator kinematics essentially complicates the stiffness modelling, because conventional serial architecture is transformed into the quasi-serial one that contains a kinematic closed loop.

The stiffness modelling of the industrial manipulators with mechanical gravity compensators is quite a new problem. There are rather limited number of works dealing with the impact of gravity compensators on the manipulator forcedeflection relation [START_REF] Wu | Geometric calibration of industrial robots using enhanced partial pose measurements and design of experiments[END_REF][START_REF] Klimchik | Identification of geometrical and elastostatic parameters of heavy industrial robots[END_REF][START_REF] Klimchik | Stiffness Modeling of Robotic Manipulator with Gravity Compensator[END_REF], while there are some work devoted to the compensator design [START_REF] Arakelian | Gravity compensation in robotics[END_REF][START_REF] Cho | Design of a Static Balancing Mechanism for a Serial Manipulator With an Unconstrained Joint Space Using One-DOF Gravity Compensators[END_REF] and softwarebased balancing solutions [START_REF] De Luca | A PD-type regulator with exact gravity cancellation for robots with flexible joints[END_REF][START_REF] De Luca | PD control with on-line gravity compensation for robots with elastic joints: Theory and experiments[END_REF]. In contrast, for conventional strictly serial manipulators [START_REF] Alici | Enhanced stiffness modeling, identification and characterization for robot manipulators[END_REF][START_REF] Chen | Conservative congruence transformation for joint and Cartesian stiffness matrices of robotic hands and fingers[END_REF][START_REF] Salisbury | Active stiffness control of a manipulator in Cartesian coordinates[END_REF] and strictly parallel mechanisms [START_REF] Klimchik | Stiffness modeling for perfect and non-perfect parallel manipulators under internal and external loadings[END_REF][START_REF] Yan | Stiffness analysis of parallelogram-type parallel manipulators using a strain energy method[END_REF][START_REF] Li | Stiffness analysis for a 3-PUU parallel kinematic machine[END_REF][START_REF] Deblaise | A systematic analytical method for PKM stiffness matrix calculation[END_REF][START_REF] Gosselin | Stiffness analysis of parallel mechanisms using a lumped model[END_REF][START_REF] Merlet | Parallel Mechanisms and Robots[END_REF] there were developed a number of methods for the stiffness analysis. At the same time, only limited number of works deals with stiffness modelling of socalled quasi-serial architectures incorporating internal closed-loops [START_REF] Klimchik | Identification of geometrical and elastostatic parameters of heavy industrial robots[END_REF][START_REF] Klimchik | Serial vs. quasi-serial manipulators: Comparison analysis of elasto-static behaviors[END_REF][START_REF] Subrin | New Redundant Architectures in Machining: Serial and Parallel Robots[END_REF][START_REF] Guo | A multilevel calibration technique for an industrial robot with parallelogram mechanism[END_REF][START_REF] Vemula | Stiffness Based Global Indices for Structural Evaluation of Anthropomorphic Manipulators[END_REF]. To our knowledge, the simplest and efficient way to take into account the influence of gravity compensator is utilization of non-linear virtual springs in the frame of the conventional VJM technique [START_REF] Klimchik | Stiffness modeling for perfect and non-perfect parallel manipulators under internal and external loadings[END_REF][START_REF] Pashkevich | Enhanced stiffness modeling of manipulators with passive joints[END_REF][START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF]. This approach was originally proposed in our previous works [START_REF] Klimchik | Identification of geometrical and elastostatic parameters of heavy industrial robots[END_REF][START_REF] Klimchik | Stiffness Modeling of Robotic Manipulator with Gravity Compensator[END_REF] and successfully applied to the manipulators the springbased gravity compensators. However, some additional efforts are required to adapt it to the case of robots with pneumatic compensators, which progressively replace their counterparts in new models of heavy industrial robots available on the market.

This paper proposes a new modification of the VJMbased stiffness modelling technique for the quasi-serial industrial manipulators with a pneumatic gravity compensator that creates a kinematic closed-loop violating the stiffness modelling principles, which are used for pure serial robot architecture. The main attention is paid to the identification of the model parameters and calibration experiment planning. The developed approach is confirmed by the experimental results that deal with compliance error compensation for robotic cell employed in manufacturing of large dimensional components.

To address these problems the remainder of the paper is organized as follows. Section 2 presents the stiffness modelling for pneumatic gravity compensator. Section 3 is devoted to the elastostatic calibration and design of calibration experiments. Section 4 deals with experimental study. Section 5 summarizes the main contributions of the paper. The mechanical structure and principal components of pneumatic gravity compensator considered here is presented in Fig. 1a; its equivalent model is shown in Fig. 1b. The mechanical compensator is a passive mechanism incorporating a constant cross-section cylinder and a constant volume gas reservoir. The volume occupied by the gas linearly depends on the piston position that defines the internal pressure of the cylinder. It is clear that this mechanism can be treated as a non-linear virtual spring influencing on the manipulator stiffness behavior. It is worth mentioning that in general case the gas temperature has impact on the pressure inside the tank, which defines the compensating force. Nevertheless, one can assume that in the case of continuous or periodical manipulator movements the gas temperature remains almost constant, i.e. the process of the gas compression-decompression can be assumed to be the isothermal one. In the frame of the manipulator model, the compensator is attached to the first and second links that creates a closedloop acting on the second actuated joint. This particularity allows us to adapt the conventional stiffness model of the serial manipulator (with constant joint stiffness matrix θ K ) by introducing the configuration dependent joint stiffness matrix θ () Kq that takes into account the compensator impact and depends on the vector of actuated coordinates q . In this case, the Cartesian stiffness matrix C K of the robotic ma- nipulator can be presented in the following form

Calibration of industrial robots with pneumatic gravity compensators

   1 1T C θ θ θ ()      K J K q J  
where θ J is the Jacobian with respect to the virtual joint coordinates θ (in the case of industrial robots it is usually equivalent to the kinematic Jacobian computed with respect to actuated coordinates q ). Thus, to obtain the stiffness model of the industrial robot with the pneumatic gravity compensator it is required to determine the non-linear joint stiffness matrix θ () Kq describing elasticity of both actuators and the gravity compensation mechanism. It should be mentioned that in the majority of works devoted to the stiffness analysis of the serial manipulators the matrix θ K is assumed to be a constant and strictly diagonal one [START_REF] Salisbury | Active stiffness control of a manipulator in Cartesian coordinates[END_REF][START_REF] Gosselin | Stiffness analysis of parallel mechanisms using a lumped model[END_REF][START_REF] Klimchik | Stiffness matrix of manipulators with passive joints: computational aspects[END_REF].

To find the desired matrix θ () PP varies with the robot motions and non-linearly depends on the angle 2 q . Below, this distances are denoted as follows:

Kq
12 ,

L P P  , 12 ,
a P P  , 01 , s PP  .
In addition, let us introduce parameters  ,  ,

x a and y a defining relevant locations of points P0, P1, P2 (see Fig. 1b). This allows us to compute the compensator length s using the following expression

 2 2 2 2 • • •co 2 s( ) aL s a L q       
which defines the non-linear function 2 () s q . For this geometry, the impact of the gravity compensator can be taken into account by replacing the considered quasiserial architecture by the serial one, where the second joint stiffness coefficient is modified in order to include elasticity of both the actuator and compensator. To find relevant nonlinear expression for this coefficient, let us present the static torque in the second joint 2 M as a geometric sum of two components. The first of them is caused by the deflection 2 q  in the mechanical transmission of the second actuated joint and can be expressed in a usual way as

 2 2 2 2 sin( •• ) qS L a M K q F q s      
where both the force S F and the compensator length s de- pend on the joint variable 0 2 2 2 qq q   . To find the compensating force S F , let us use the iso- thermal process assumption that yields the relation

P V const 

, where P is the tank pressure, 00 () V A s s V     is the corresponding internal volume, A is the piston area, 0 s is compensator link length correspond- ing to zero compensating force and 0 V is the tank volume corresponding to the atmospheric pressure

() •• ( sin( ) ) q A s s L A s s V a M K q P q s         
and present in a more compact form

 2 2 0 2 2 0 s n( ) • i • q V q M ss a L q s KP s s           where a constant 00 / V s V s A  
is the equivalent distance. Further, after computing the partial derivative

( ) ( ) ) () V V q V P L a s s a L q ss s s s ss q a L q K q K s s                             
 which is obviously highly non-linear with respect to manipulator configuration (here, s is also a non-linear function of 2 q ). Nevertheless, it allows us to compute an relevant stiff- ness coefficient 2 K for the equivalent serial chain and direct- ly apply eq. ( 1) to evaluate stiffness of the quasi-serial manipulator with pneumatic gravity compensator.

It should be mentioned that in practice the compensator parameters 0 ,

V s s and actuator stiffness coefficients

1 2 3 6
, , ,..., K K K K are usually not given in the robot datasheets, so they should be identified via dedicated experimental study. For this reason, the following Section focuses on the identification of this extended set of the manipulator elastostatic parameters.

III. ELASTOSTATIC PARAMETERS IDENTIFICATION

A. Methodology

In the frame of the VJM-based modelling approach developed for serial kinematic chains [START_REF] Salisbury | Active stiffness control of a manipulator in Cartesian coordinates[END_REF][START_REF] Pashkevich | Enhanced stiffness modeling of manipulators with passive joints[END_REF] and adapted here for the case of quasi-serial manipulators with pneumatic gravity compensators, the desired stiffness model parameters describe elasticity of the virtual springs located in the actuated joints of the manipulator, and also compensator parameters 0 ,

V s s of defining preloading of the compensator spring and the equivalent distance for the tank volume , , ,..., K K K K used in previous sec- tion) and the compensator elastic parameters as 0 ,

V s s . To find the desired set of elastic parameters, the robotic manipulator sequentially passes through several measurement configurations where the external loading is applied to the specially designed end-effector presented in Fig. 2 (it allows us to generate both forces and torques applied to the manipulator). Using the absolute measurement system (the laser tracker Leica AT900, the Cartesian coordinates of the reference points are measured twice, before and after loading. To increase identification accuracy, it is reasonable to have several markers on the end effector (reference points) and to apply the loading of the maximum allowed magnitude. It should be mentioned that to avoid singularities caused by numerical routines, the external force/torque directions should not be the same for all calibration experiments (while from the practical view point the mass-based gravity loading is the most attractive). Thus, the calibration experiments yield the dataset that includes values of the manipulator joint coordinates   i q , applied forces/torques  

B. Identification algorithm

To take into account the compensator influence while using classical approach developed for strictly serial manipulators without compensators [START_REF] Alici | Enhanced stiffness modeling, identification and characterization for robot manipulators[END_REF], it was proposed below to use in the second joint an equivalent virtual spring with nonlinear stiffness, which depends on the joint coordinate 2 q (see eq. ( 1)). Using this idea, it is convenient to consider several aggregated compliances 2 i k  corresponding to each different value of angle 2 q . This idea allows us to linearize the identification equations with respect to extended set of model parameters and that can be easily solved using standard least-square technique.

Let us denote this extended set of desired parameters as

TT i i i i ni ni i im    F A F J J J J  
that is usually used in stiffness analysis of serial manipulators. Here, ni J denotes the manipulator Jacobian column, superscript '(p)' stands for the Cartesian coordinates (position without orientation). Transformation from i A to () p i B is ra-ther trivial and is based on the extraction from i A the first three lines and inserting in it several zero columns.

In this case, the elastostatic parameters identification can be reduced to the following optimization problem

 0 ,, ( ) ( ) 1 ( ) ( ) min jc m p T p i i i i k k i F          B k p B k p  
which yields to the following solution

 1 ( ) ( ) ( ) 11 • TT mm p p p i i i i ii                 k B B B p  
where the parameters

1 3 6
, ,..., k k k describe the compliance of the virtual joints #1,#3,...#6, while the rest of them 21 22 , ... kk present an auxiliary dataset allowing to separate the compliance of the joint #2 and the compensator parameters 0 , V s s . Using eq. ( 7), the desired optimization problem can be written as

  2 0 2 2 00 2 2 22 1 2 , 0 2 , () ) ( ) m sin ( cos( in ) sin qV q m i iV i i V i ii s V q K s i P L a s s a L q ss s s s ss q aL q s K s                                    
where q m is the number of different angles 2 q in the exper- imental data. It is obvious that eq. ( 12) is highly non-linear and can be solved numerically only.

Thus, the proposed modification of the previously developed calibration technique allows us to find the manipulator and compensator parameters. An open question, however, is how to find the set of measurement configurations that ensure the lowest impact of the measurement noise.

C. Design of calibration experiments

The goal of the calibration experiment design is to select a set of robot configurations/external loadings   , ii qF that ensure the best identification accuracy. The key issue is to rang plans of experiments in accordance with selected performance measure. This problem is well known in the classical regression analysis; however, the results are not suitable for non-linear case of the elastostatic calibration and require additional efforts. Here, an industry oriented performance measure is used, which evaluates the calibration plan quality [START_REF] Wu | Geometric calibration of industrial robots using enhanced partial pose measurements and design of experiments[END_REF][START_REF] Klimchik | Identification of geometrical and elastostatic parameters of heavy industrial robots[END_REF]. Its physical meaning is the robot positioning accuracy (under the loading), which is achieved after compliance error compensation based on the identified elastostatic parameters.

Assuming that experiments include measurement errors i ε , covariance matrix for the parameters k can be expressed as

      1 ( ) ( ) 1 1 ( ) ( ) ( ) ( ) 11 cov( ) E T TT m pp ii i mm p T p p p i i i i i i ii         k B B B ε ε B B B  
Following independent identically distributed assumption with zero expectation and standard deviation 2  for the measurement errors, expression (13) can be simplified to

   1 2 ( ) ( ) 1 cov( ) T m pp ii i      k B B  
Hence, the impact of the measurement errors on the accuracy of the identified parameters k is defined by the matrix

( ) ( ) 1 T m pp ii i
 BB (in regression analysis it is known as the in- formation matrix).

It is evident that in industrial practice the most important issue is not the parameters identification accuracy, but their impact on the robot positioning accuracy. Considering that the end-effector accuracy varies throughout the workspace and highly depends on the manipulator configuration, it is proposed to evaluate the calibration accuracy in a typical manipulator configuration ("test-pose") provided by the user. For the most of applications of heavy industrial robots, the test pose is usually related to the typical machining configuration 0 q and corresponding external loading 0 F related to the corresponding technological process. For the so-called test-pose the mean square value of the positioning error will be denoted as 2 0  and the matrix () p i A corresponding to it as

() 0 p A .
It should be noted that that the proposed approach operates with a specific structure of the parameters included in the vector k , where the second joint is presented by several components 21 22 , ... kk while the other joints are described by a single parameter , ... k k k . This motivates further re- arrangement of the vector k and replacing it by several vectors ) cov( )

T j j j   k k k , the performance measure 2 0  can be presented as    1 ( ) ( ) 2 00 2 ) ( 0 1 1 () trace T q T m j p j p m pp ii i j            A A A A  
Based on this performance measure, the calibration experiment design can be reduced to the following optimization problem

   1 ( ) ( ) 00 1 { , } ( ) ( ) 1 trace min T q T ii m j p j p i pp j i m i          qF A A A A   subject to max , 1.. i F i m  F
whose solution gives a set of the desired manipulator configurations and corresponding external loadings. It is evident that its analytical solution can hardly be obtained and a numerical approach is the only reasonable one.

IV. EXPERIMENTAL STUDY

The developed technique was applied to the elastostatic calibration of robot Kuka KR-120. The parameters to be identified were the compliances j k of the actuated joints and the gravity compensator parameters 0 , V s s . To generate de- flections in the actuated joints, the gravity forces 120 kg were applied to the robot end-effector (see Fig 3). The Cartesian coordinates of three markers located on the tool (see Fig, 2) have been measured before and after the loading. To find optimal measurement configurations for calibration, the design of experiments was used for six different angles 2 q that are distributed between the joint limits. For each 2 q from three to seven optimal measurement configurations were found, which satisfy joint limits and physical constraints related to the possibility carry out experiments. In total 31 different measurement configurations and 186 measurements were considered for the identification, from which 7 physical parameters were obtained. The obtained experimental data have been processed using the identification algorithm presented in Section 3. Identified values for the extended set of joint compliances (for 6 different angles q2) and their confidence intervals are presented in Table 1. As follows from this results, wrist compliances were identified with lower accuracy. The reason for it is smaller shoulder from the applied external forces comparing with manipulator joints. Relatedly small accuracy of first joint due to a smaller number of measurements in the experiments in which the deflections were generated in the first joint. Further, obtained compliances k21… k26 were used to estimate pneumonic compensator parameters by solving optimization problem [START_REF] Klimchik | Identification of the manipulator stiffness model parameters in industrial environment[END_REF]. The identified joint compliances can be used to predict robot deformations under the external loading. V. CONCLUSIONS

Gravity compensator

The paper presents a new approach for the modelling and identification of the elastostatic parameters of heavy industrial robots with the pneumatic gravity compensator. It proposes a methodology and data processing algorithms for the identification of elastostatic parameters of gravity compensator and manipulator. To increase the identification accuracy, the design of experiments has been used aimed at proper selection of the measurement configurations. In contrast to other works, it is based on the industry oriented performance measure that is related to the robot accuracy under the loading. The advantages of the developed techniques are illustrated by experimental study of the industrial robot Kuka KR-120, for which the joint compliances and parameters of the gravity compensator have been identified.

Figure 1 .

 1 Figure 1. Pneumatic gravity compensator and its mode

  , let us consider the compensator geometry in detail. As follows from Fig 1b, the compensator geometrical model contains three principal node points P0, P1, P2., where P0, P1 define the passive joint rotation axes and P2 defines the second actuated joint axis. In this model, two distances 12 ,

K 2 M

 2 is the stiffness coefficient. The second component can be presented as• sin force generated by the gravity compensator. It is clear that sin can be computed from the triangle 12 allows us to express the torque in the following form:

  Taking into account that compensating force S F depends on the internal and external pressure difference and is computed as

0V

  . In the frame of this model, let us denote the manipulator joint com-

iFFigure 2 .

 2 Figure 2. End-effector used for the elastostatic calibration experiments and it model

F

  them in the vector k . In this case the linearize force-deflection relation with respect to this vector can be present in the following form () p ii  p B k  where i p is the vector of the end-effector displacements under the external loading i

Figure 3 .

 3 Figure 3. Experimental setup for the identification of the elastostatic parameters.

TABLE I .

 I ELASTO-STATIC PARAMETERS OF ROBOT KUKA KR-120

	Parameter	value	CI
	k1, [rad×μm/N]	1.13	±0.15 (13.3%)
	k21, [rad×μm/N]	0.34	±0.004 (1.1%)
	k22, [rad×μm/N]	0.36	±0.005 (1.4%)
	k23, [rad×μm/N]	0.35	±0.005 (1.4%)
	k24, [rad×μm/N]	0.28	±0.007 (2.6%)
	k25, [rad×μm/N]	0.32	±0.011 (3.6%)
	k26, [rad×μm/N]	0.26	±0.007 (2.8%)
	k3, [rad×μm/N]	0.43	±0.007 (1.8%)
	k4, [rad×μm/N]	0.95	±0.31 (31.8%)
	k5, [rad×μm/N]	3.82	±0.27 (7.0%)
	k6, [rad×μm/N]	4.01	±0.35 (8.7%)

ACKNOWLEDGMENTS

The work presented in this paper was partially funded by Innopolis University and the project Partenariat Hubert Curien SIAM 2016 France-Thailand.

*The work presented in this paper was partially funded by Innopolis University and the project Partenariat Hubert Curien SIAM 2016 France-Thailand


Jiuchun Gao
email: jiuchun.gao@ls2n.fr

Anatol Pashkevich

Sté Phane Caro

Optimal Trajectories Generation in Robotic Fiber Placement Systems

The paper proposes a methodology for optimal trajectories generation in robotic fiber placement systems. A strategy to tune the parameters of the optimization algorithm at hand is also introduced. The presented technique transforms the original continuous problem into a discrete one where the time-optimal motions are generated by using dynamic programming. The developed strategy for the optimization algorithm tuning allows essentially reducing the computing time and obtaining trajectories satisfying industrial constraints. Feasibilities and advantages of the proposed methodology are confirmed by an application example.

Introduction

Robotic fiber placement technology has been increasingly implemented recently in aerospace and automotive industries for fabricating complex composite parts [START_REF] Gay | Composite materials: design and applications[END_REF][START_REF] Gallet-Hamlyn | Multiple-use robots for composite part manufacturing[END_REF]. It is a specific technique that uses robotic workcell to place the heated fiber tows on the workpiece surface [START_REF] Peters | Handbook of composites[END_REF]. Corresponding robotic systems usually include a 6-axis industrial robot and a one-axis positioner (see Figure 1), which are kinematically redundant and provides the user with some freedom in terms of optimization of robot and positioner motions.

To deal with the robotic system redundancy, a common technique based on the pseudo-inverse of kinematic Jacobian is usually applied. However, as follows from relevant studies, this standard approach does not satisfy the real-life industrial requirements of the fiber placement [START_REF] Kazerounian | Redundancy resolution of serial manipulators based on robot dynamics[END_REF][START_REF] Buss | Introduction to inverse kinematics with jacobian transpose, pseudoinverse and damped least squares methods[END_REF]. In literature, there is also an alternative technique (that deals with multi-goal tasks) that is based on conversion of the original continuous problem to the combinatorial one [START_REF] Gueta | Coordinated motion control of a robot arm and a positioning table with arrangement of multiple goals[END_REF][START_REF] Gueta | Hybrid design for multiple-goal task realization of robot arm with rotating table[END_REF], but it only generates trajectories for point-to-point motions, e.g. for spot welding applications. A slightly different method was introduced in [START_REF] Dolgui | Manipulator motion planning for high-speed robotic laser cutting[END_REF][START_REF] Pashkevich | Multiobjective optimization of robot motion for laser cutting applications[END_REF][START_REF] Zhou | Off-line programming system of industrial robot for spraying manufacturing optimization[END_REF], and it was successfully applied to laser-cutting and arc-welding processes where the tool speed was assumed to be constant (which is not valid in the considered problem). Another approach has been proposed in [START_REF] Debout | Tool path smoothing of a redundant machine: Application to Automated Fiber Placement[END_REF], where the authors concentrated on the tool path smoothing in Cartesian space in order to decrease the manufacturing time in fiber placement applications. For the considered process, where the tool speed variations are allowed (in certain degree), a discrete optimization based methodology was proposed in our previous work [START_REF] Gao | Manipulator Motion Planning in Redundant Robotic System for Fiber Placement Process[END_REF]. It allows the user to convert the original continuous problem to the combinatorial one taking into account particularities of the fiber placement technology and to generate time-optimal trajectories for both the robot and the positioner. Nevertheless, there are still a number of open questions related to selection of the optimization algorithm parameters (i.e., its "tuning") that are addressed in this paper and targeted to the improvement of the algorithm efficiency and the reduction of the computing time.

Robotic system model

In practice, the procedure of off-line programming for robotic fiber placement is implemented in the following way. The fiber placement path is firstly generated and discretized in CAM system. Further, the obtained set of task points is transformed into the task graph that describes all probable configurations of the robot and the positioner joints. The motion generator module finds the optimal trajectories that are presented as the "best" path on the graph. Finally, the obtained motions are converted into the robotic system program by the post processor. The core for the programming of this task is a set of optimization routines addressed in this paper.

To describe the fiber placement task, let us present it as a set of discrete task frames

n i F i task ,... 2 , 1 , ) (



, in such a way that the X-axis is directed along the path direction and Z-axis is normal to the workpiece surface pointing outside of it (see Figure 1). Using these notations, the task locations can be described by 4 4 homogenous transformation matrices and the considered task is formalized as follows: where all vectors of positions and orientations are expressed with respect to the workpiece frame (see superscript "w"). To execute the given fiber placement task, the robot tool must visit the frames defined by (1) as fast as possible.

) ( ) ( ) 2
The considered robotic system, shown in Figure 1, is composed of an industrial robot and an actuated positioner. Their spatial configurations can be described by the joint coordinates R q and P q respectively. The task frames can be presented in two ways using the robot and positioner kinematics that are expressed as

) ( R R g q and
) ( P P q g , respectively. To obtain the kinematic model of the whole system that is expressed as a closed loop containing the robot, the workpiece and the positioner, a global frame 0 F is selected. Then, the tool frame , tool F and task frame ) (i task F can be aligned in such a way that: (i) the origins of the two frames coincide; (ii) Z-axes are opposite; (iii) X-axes have the same direction. Due to the foregoing closed-loop, two paths can be followed to express the transformation matrices from the global frame to the task frames, namely,

n i q g g i task w i P P Pbase task Tool i R R Rbase ,... 2 , 1 ; ) ( ) ( ) ( ) ( 0 ) ( 0       T T T q T (2)
Equation ( 2) does not lead to a unique solution for R q and P q as the robotic system, i.e., robot and positioner, is kinematically redundant. Therefore, the optimum robot and positioner configurations can be searched based on specific criteria.

Algorithm for trajectories generation

To take advantage of the kinematic redundancy, it is reasonable to partition the desired motion between the robot and the positioner ensuring that the technology tool executes the given task with smooth motion as fast as possible.

To present the problem in a formal way, let us define the functions

) (t R q and ) (t q P
describing the robot and positioner motion as a function of time

] , 0 [ T t  . Additionally, a sequence of time instants } ,... , { 2 1 n t t t
corresponds to the cases where the tool visits the locations defined by (1), and

T t t n   , 0
1 . As a result, the problem at hand is formulated as an optimization problem aiming at minimizing the robot processing time

) ( ), ( min t q t P R T q  (3)
This problem is subjected to the equality constraints [START_REF] Gallet-Hamlyn | Multiple-use robots for composite part manufacturing[END_REF] and some inequality constraints associated to the capacities of the robot/positioner actuators that are defined by upper bounds of the joint velocities and accelerations. Besides, the collision constraints verifying the intersections between the system components are also taken into account.

) ( 0 0 )) ( ( )) ( ( i task w i P P Pbase task Tool i R R Rbase t q g t g T T T q T      defined in
For this considered problem aiming at finding desired continuous function of

) (t R q and ) (t q P
, there is no standard approach that can be applied to straightforwardly. The main difficulty here is that the equality constraints are written for the unknown time instants

} ,... , { 2 1 n t t t
. Besides, this problem is nonlinear and includes a redundant variable. For these reasons, this paper presents a combinatorial optimization based methodology to generate the desired trajectories.

For the considered robotic system, there is one redundant variable with respect to the given task. It is convenient here to treat P q as the redundant one since it allows us to use the kinematic models of the robot and the positioner independently and to consider the previous equality constraints.

To present the problem in a discrete way, the allowable domain of ] , [ max min P P P q q q  is sampled with the step

P q  as m k k q q q P P k P ,... 1 , 0 ; min ) (     
, where

P P P q q q m    ) ( min max
. Then, applying sequentially the positioner direct kinematics and the robot inverse kinematics, a set of possible configuration states for the robotic system can be obtained as

n i t q g g t tool task i task w i k P P Pbase Rbase R i k R ,... 2 , 1 ); , )) ( ( ( ) ( ) ( ) ( 0 0 1 ) (        T T T T q
, where μ is a configuration index vector corresponding to the robot posture. Therefore, for i th task location, a set of candidate configuration states can be obtained, i.e., T in joint space as above, the original task can be converted into the directed graph shown in Figure 2. It should be noted that some of the configuration cells should be excluded because of violation of the collision constraints or the actuator joint limits. These cases are denoted as "inadmissible" in Figure 2, and are not connected to any neighbor. Here, the allowable connection between the graph nodes are limited to the subsequent configuration states

) , ( ) , ( i k task i k task     L L
, and the edge weights correspond to the minimum robot processing time restricted by the maximum velocities and accelerations of the robot and the positioner.

Using the discrete search space above, the considered problem is transformed to the classic shortest path searching and the desired solution can be represented as the sequence , where

} ...{ } { } { n) , (k ,2) (k ,1) (k
6 ,... 1 , 0 ; ) max( ) , ( max ) ( 1 , ) ( , ) 1 , ( ) , ( 1 1
       j q q q dist j k i j k i j i k task i k task i i i i  L L
. It should be mentioned that the above expression takes into account the velocity constraints automatically and the acceleration constraints should be considered by means of the following formula:

max 1 1 ) ( 1 , ) ( 1 ) ( ) ( 1 , ) ( ) ( ) ( 2 1 1 j i i i i k i j k i k k i j i q t t t t q t q t i i i i                     i j, i j, q q (5) where ) , ( ) 1 ( ) ( 1 1      ,i k task ,i k task i i i dist t L L and ) , ( ) 1 ( ) ( 1     ,i k task ,i k task i i i dist t L L .
By discretizing the search space, the original problem is converted to a combinatorial one, which can be solved by using conventional way, e.g. However, this straightforward approach is extremely time-consuming and can be hardly accepted for industrial applications. For example, it takes over 20 hours to find a desired solution in a relatively simple case (two-axis robot and one-axis positioner), where the search space is built for 100 task points and the discretization step 1° (processor Intel® i5 2.67 GHz) [START_REF] Gao | Manipulator Motion Planning in Redundant Robotic System for Fiber Placement Process[END_REF]. Besides, known methods are not able to take into account the acceleration constraints that are necessary here. For these reasons, a problem-oriented algorithm taking into account the particularities of the graph based search space is proposed in this paper.

The developed algorithm is based on the dynamic programming principle, aiming at finding the shortest path from

} , { 1 ) 1 ( 1 k , k task  L to the current } , { ) ( i ,i k task k i  L
. The length of this shortest path is denoted as i k d

, . Then, the shortest path for the locations corresponding to the next

} , { ) 1 ( k k,i task   L
can be obtained by combining the optimal solutions for the previous column

} , { ) ( k ,i k task    L
and the distances between the task locations with the indices i and i+1,

; } ) , ( { ) ( ) 1 , ( , 1 , ,i k task i k task i k k i k dist d d        L L min (4)
This formula is applied sequentially from the second column of the task graph to the last one, and the desired optimal path can be obtained after selection of the minimum length . This proposed algorithm is rather time-efficient since it takes about 30 seconds [START_REF] Gao | Manipulator Motion Planning in Redundant Robotic System for Fiber Placement Process[END_REF] to find the optimal solution for the above mentioned example.

Tuning of trajectories generation algorithm

For the proposed methodology, the discretization step for the redundant variable is a key parameter, which has a big influence on the algorithm efficiency. An unsuitable discretization step may lead either to a bad solution or high computational time. For this reason, a new strategy for the determination of the discretization step is proposed thereafter to tune the optimization algorithm.

Influence of the discretization step

Let us consider a simple case study that deals with a three-axis planar robotic system executing a straight-line-task (see Figure 3). For this problem, the fiber placement path is uniformly discretized into 40 segments. Relevant optimization results are presented in Table 1. It is clear that here (as well as in other cases) smaller discretization step should provide better results but there exists a reasonable lower bound related to an acceptable computing time.

   2 P q    1 P q    75 . 0 P q    5 . 0 P q    25 . 0 P q    1 . 0 P q
Robot processing time are not acceptable because they lead to a robot processing time 20-50% higher than the optimal one. Moreover, in the case of

   2 P q
, the optimization algorithm generates a bad solution that does not take advantage of the positioner motion capabilities (positioner is locked,

var const q R P   q ,
). The reason for this phenomenon is that the discretization step here is so large that the positioner step-time is always higher than the robot moving time between the subsequent task points.

Another interesting phenomenon can be observed for slightly smaller discretization steps, where the algorithm may produce non-smooth intermittent rotation of the positioner (start-stop motion: (without acceleration constraints), the optimization algorithm generates solution that includes only several steps where the positioner is not locked.

In addition, it is noteworthy that in the case of with acceleration constraints, the discretization step reduction from 2° to 1° leads to even worse solution, where the robot processing time is about 10% higher. This phenomenon can be explained by heuristic integration of the acceleration constraints into the optimization algorithm, which may slightly violate the dynamic programming principle. Nevertheless, further reduction of P q  allows to restore the expected algorithm behavior. Hence, to apply the developed technique in practice, users need some simple "rules of thumb" that allows setting an initial value of P q  . Then, the optimization algorithm can be applied several times (sequentially decreasing P q  ) until the objective function convergence. To reduce computing time in the case of small P q  , some local optimization techniques have been also developed by the authors.

Initial tuning of the optimization algorithm

To find a reasonable initial value of the discretization step, let us investigate in details robot and positioner motions between two sequential task locations. It is clear that for smooth positioner motions, it is required that corresponding increments of the coordinate P q should include at least one discretization step P q  . To find the maximum value of P q  , let us denote   as the increment of P q for the movement between two adjacent task locations (Pi-Pi+1) and s  as the length of the path segment. It is clear that s  can be also treated as the arc length between Pi and Pi+1 around the positioner joint axis. Let us assume that the distance from a path point to the rotational axis is r, and rmax represents the furthest task location with respect to the positioner axis. To avoid undesired intermittent positioner rotations, the following constraint      max r s should be verified, since the positioner velocity is usually smaller than the velocity of the robot. The latter inequality can be rewritten in terms of the robot/positioner motion time as max max max

) ( P R q v r s         
, which can be further

    P q
where number of the positioner steps is no less one. Hence, the initial value of P q  should be at least equal to ) ( max max max max r q v s q P R P       in order to provide acceptable motions of the robot and positioner. For instance, for the previous case study, this expression gives the discretization step about 0.5° that allows to generate trajectories that are very close to the optimal ones, namely, the robot processing time is only 1% higher than the minimum value.

Conclusions

This paper contributes to optimization of robot/positioner motions in redundant robotic systems for the fiber placement process. It proposes a new strategy for the optimization algorithm tuning. The developed technique converts the continuous optimization into a combinatorial one where dynamic programming is applied to find time-optimal motions. The proposed strategy of the optimization algorithm tuning allows essentially decreasing the computing time and generating desired motions satisfying industrial constraints. Feasibilities and advantages of the presented technique are confirmed by a case study. Future research will focus on application of those results in real-life industrial environments.

Figure 1 .

 1 Figure 1. Typical robotic fiber placement system (6-axis robot and one-axis positioner).

Figure 2 .

 2 Figure 2. Graph-based representation of the discrete search space After presenting

  function (robot processing time) can be presented as the sum of the edge weights

  the final column. Therefore, the desired path is described by the recorded indices

Figure 3 .

 3 Figure 3. Three-axis planar robotic system and straight-line-task

  To estimate the reasonable discretization step for the considered fiber placement problem, let us analyze Table1in more details. From Table1, the discretization steps

		T=1.90s	T=1.84s	T=1.54s	T=1.30s		T=1.29s	T=1.29s
	(without acc-constraint)	(38s comp.)	(2m comp.)	(4m comp.)	(9m comp.)	(47m comp.)	(5h comp.)
	Robot processing time	T=1.90s	T=2.11s	T=1.60s	T=1.30s		T=1.29s	T=1.29s
	(with acc-constraint)	(67s comp.)	(4m comp.)	(8m comp.)	(17m comp.)	(1.2h comp.)	(9h comp.)
						 P		,  2 {	1 	,	}  75 . 0
						q

Acknowledgments

This work has been supported by the China Scholarship Council (Grant N° 201404490018). The authors also acknowledge CETIM for the motivation of this research work.


Charles Desmarchelier

Véronique Rosilio
email: veronique.rosilio@u-psud.fr

David Chapron

Ali Makky

Damien Preveraud

Estelle Devillard

Véronique Legrand-Defretin

Patrick Borel

Damien P Prévéraud

Molecular interactions governing the incorporation of cholecalciferol and retinyl-palmitate in mixed taurocholate-lipid micelles

Keywords: bioaccessibility, surface pressure, bile salt, compression isotherm, lipid monolayer, vitamin A, vitamin D, phospholipid

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

Introduction

Retinyl esters and cholecalciferol (D 3 ) (Figure 1) are the two main fat-soluble vitamins found in foods of animal origin. There is a renewed interest in deciphering their absorption mechanisms because vitamin A and D deficiency is a public health concern in numerous countries, and it is thus of relevance to identify factors limiting their absorption to tackle this global issue. The fate of these vitamins in the human upper gastrointestinal tract during digestion is assumed to follow that of dietary lipids [START_REF] Borel | Vitamin D bioavailability: state of the art[END_REF]. This includes emulsification, solubilization in mixed micelles, diffusion across the unstirred water layer and uptake by the enterocyte via passive diffusion or apical membrane proteins [START_REF] Reboul | Proteins involved in uptake, intracellular transport and basolateral secretion of fat-soluble vitamins and carotenoids by mammalian enterocytes[END_REF].

Briefly, following consumption of vitamin-rich food sources, the food matrix starts to undergo degradation in the acidic environment of the stomach, which contains several enzymes, leading to a partial release of these lipophilic molecules and to their transfer to the lipid phase of the meal. Upon reaching the duodenum, the food matrix is further degraded by pancreatic secretions, promoting additional release from the food matrix, and both vitamins then transfer from oil-in-water emulsions to mixed micelles (and possibly other structures, such as vesicles, although not demonstrated yet). As it is assumed that only free retinol can be taken up by enterocytes, retinyl esters are hydrolyzed by pancreatic enzymes, namely pancreatic lipase, pancreatic lipase-related protein 2 and cholesterol ester hydrolase [START_REF] Desmarchelier | The distribution and relative hydrolysis of tocopheryl acetate in the different matrices coexisting in the lumen of the small intestine during digestion could explain its low bioavailability[END_REF].

Bioaccessible vitamins are then taken up by enterocytes via simple passive diffusion or facilitated diffusion mediated by apical membrane proteins (Desmarchelier et al. 2017). The apical membrane protein(s) involved in retinol uptake by enterocytes is(are) yet to be identified but in the case of D 3 , three proteins have been shown to facilitate its uptake: NPC1L1 (NPC1 like intracellular cholesterol transporter 1), SR-BI (scavenger receptor class B member 1) and CD36 (Cluster of differentiation 36) [START_REF] Reboul | Proteins involved in uptake, intracellular transport and basolateral secretion of fat-soluble vitamins and carotenoids by mammalian enterocytes[END_REF]. Both vitamins then transfer across the enterocyte towards the basolateral side. The transfer of vitamin A is mediated, at least partly, by the cellular retinol-binding protein, type II (CRBPII), while that of vitamin D is carried out by unknown mechanisms. Additionally, a fraction of retinol is re-esterified by several enzymes (Borel & Desmarchelier 2017). Vitamin A and D are then incorporated in chylomicrons in the Golgi apparatus before secretion in the lymph.

The solubilization of vitamins A and D in mixed micelles, also called micellarization or micellization, is considered as a key step for their bioavailability because it is assumed that the non-negligible fraction of fat-soluble vitamin that is not micellarized is not absorbed [START_REF] Desmarchelier | The distribution and relative hydrolysis of tocopheryl acetate in the different matrices coexisting in the lumen of the small intestine during digestion could explain its low bioavailability[END_REF]. Mixed micelles are mainly made of a mixture of bile salts, phospholipids and lysophospholipids, cholesterol, fatty acids and monoglycerides [START_REF] Hernell | Physical-chemical behavior of dietary and biliary lipids during intestinal digestion and absorption. 2. Phase analysis and aggregation states of luminal lipids during duodenal fat digestion in healthy adult human beings[END_REF]). These compounds may form various self-assembled structures, e.g., spherical, cylindrical or disk-shaped micelles [START_REF] Walter | Intermediate structures in the cholate-phosphatidylcholine vesicle-micelle transition[END_REF][START_REF] Leng | Kinetics of the micelle-to-vesicle transition ; aquous lecithin-bile salt mixtures[END_REF] or vesicles, depending on their concentration, the bile salt/phospholipid ratio [START_REF] Walter | Intermediate structures in the cholate-phosphatidylcholine vesicle-micelle transition[END_REF], the phospholipid concentration, but also the ionic strength, pH and temperature of the aqueous medium (Madency & Egelhaaf 2010;[START_REF] Salentinig | Self-assembled structures and pKa value of oleic acid in systems of biological relevance[END_REF][START_REF] Cheng | Mixtures of lecithin and bile salt can form highly viscous wormlike micellar solutions in water[END_REF]. Fat-soluble micronutrients display large variations with regards to their solubility in mixed micelles [START_REF] Sy | Effects of physicochemical properties of carotenoids on their bioaccessibility, intestinal cell uptake, and blood and tissue concentrations[END_REF][START_REF] Gleize | Form of phytosterols and food matrix in which they are incorporated modulate their incorporation into mixed micelles and impact cholesterol micellarization[END_REF] and several factors are assumed to account for these differences (Desmarchelier & Borel 2017, for review).

The mixed micelle lipid composition has been shown to significantly affect vitamin absorption. For example, the substitution of lysophospholipids by phospholipids diminished the lymphatic absorption of vitamin E in rats [START_REF] Koo | Phosphatidylcholine inhibits and lysophosphatidylcholine enhances the lymphatic absorption of alpha-tocopherol in adult rats[END_REF]. In rat perfused intestine, the addition of fatty acids of varying chain length and saturation degree, i.e. butyric, octanoic, oleic and linoleic acid, resulted in a decrease in the rate of D 3 absorption [START_REF] Hollander | Vitamin D-3 intestinal absorption in vivo: influence of fatty acids, bile salts, and perfusate pH on absorption[END_REF]. The effect was more pronounced in the ileal part of the small intestine following the addition of oleic and linoleic acid. It was suggested that unlike short-and medium-chain fatty acids, which are not incorporated into micelles, long-chain fatty acids hinder vitamin D absorption by causing enlargement of micelle size, thereby slowing their diffusion towards the enterocyte.

Moreover, the possibility that D 3 could form self-aggregates in water [START_REF] Meredith | The Supramolecular Structure of Vitamin-D3 in Water[END_REF], although not clearly demonstrated, has led to question the need of mixed micelles for its solubilization in the aqueous environment of the intestinal tract lumen [START_REF] Rautureau | Aqueous solubilisation of vitamin D3 in normal man[END_REF][START_REF] Maislos | Bile salt deficiency and the absorption of vitamin D metabolites. In vivo study in the rat[END_REF].

This study was designed to compare the relative solubility of D 3 and RP in the aqueous phase rich in mixed micelles that exists in the upper intestinal lumen during digestion, and to dissect, by surface tension and surface pressure measurements, the molecular interactions existing between these vitamins and the mixed micelle components that explain the different solubility of D 3 and RP in mixed micelles.

Materials and methods

Chemicals

2-oleoyl-1-palmitoyl-sn-glycero-3-phosphocholine (POPC) (phosphatidylcholine, ≥99%; Mw 760.08 g/mol), 1-palmitoyl-sn-glycero-3-phosphocholine (Lyso-PC) (lysophosphatidylcholine, ≥99%; Mw 495.63 g/mol), free cholesterol (≥99%; Mw 386.65 g/mol), oleic acid (reagent grade, ≥99%; Mw 282.46 g/mol), 1-monooleoyl-rac-glycerol (monoolein, C18:1,-cis-9, Mw 356.54 g/mol), taurocholic acid sodium salt hydrate (NaTC) (≥95%; Mw 537.68 g/mol) ), cholecalciferol (>98%; Mw 384.64 g/mol; melting point 84.5°C; solubility in water: 10 -4 -10 -5 mg/mL; logP 7.5) and retinyl palmitate (>93.5%; Mw 524.86 g/mol; melting point 28.5°C; logP 13.6) were purchased from Sigma-Aldrich (Saint-Quentin-Fallavier, France). Chloroform and methanol (99% pure) were analytical grade reagents from Merck (Germany). Ethanol (99.9%), n-hexane, chloroform, acetonitrile, dichloromethane and methanol were HPLC grade reagents from Carlo Erba Reagent (Peypin, France). Ultrapure water was produced by a Milli-Q ® Direct 8 Water Purification System (Millipore, Molsheim, France). Prior to all surface tension, and surface pressure experiments, all glassware was soaked for an hour in a freshly prepared hot TFD4 (Franklab, Guyancourt, France) detergent solution (15% v/v), and then thoroughly rinsed with ultrapure water. Physico-chemical properties of D 3 and RP were retrieved from PubChem (https://pubchem.ncbi.nlm.nih.gov/).

Micelle formation

The micellar mixture contained 0.3 mM monoolein, 0.5 mM oleic acid, 0.04 mM POPC, 0.1 mM cholesterol, 0.16 mM Lyso-PC, and 5 mM NaTC [START_REF] Reboul | Lutein transport by Caco-2 TC-7 cells occurs partly by a facilitated process involving the scavenger receptor class B type I (SR-BI)[END_REF]. Total component concentration was thus 6.1 mM, with NaTC amounting to 82 mol%. Two vitamins were studied: crystalline D 3 and RP.

Mixed micelles were formed according to the protocol described by [START_REF] Desmarchelier | The distribution and relative hydrolysis of tocopheryl acetate in the different matrices coexisting in the lumen of the small intestine during digestion could explain its low bioavailability[END_REF]. Lipid digestion products (LDP) (monoolein, oleic acid, POPC, cholesterol and Lyso-PC, total concentration 1.1 mM) dissolved in chloroform/methanol (2:1, v/v), and D 3 or RP dissolved in ethanol were transferred to a glass tube and the solvent mixture was carefully evaporated under nitrogen. The dried residue was dispersed in Tris buffer (Tris-HCl 1mM, CaCl 2 5mM, NaCl 100 mM, pH 6.0) containing 5 mM taurocholate, and incubated at 37 °C for 30 min. The solution was then vigorously mixed by sonication at 25 W (Branson 250W sonifier; Danbury, CT, U.S.A.) for 2 min, and incubated at 37 °C for 1 hour. To determine the amount of vitamin solubilized in structures allowing their subsequent absorption by enterocytes (bioaccessible fraction), i.e. micelles and possibly small lipid vesicles, whose size is smaller than that of mucus pores [START_REF] Cone | Barrier properties of mucus[END_REF], the solutions were filtered through cellulose ester membranes (0.22 µm) (Millipore), according to [START_REF] Tyssandier | Processing of vegetable-borne carotenoids in the human stomach and duodenum[END_REF]. The resulting optically clear solution was stored at -20 °C until vitamin extraction and HPLC analysis. D 3 and RP concentrations were measured by HPLC before and after filtration. For surface tension measurements and cryoTEM experiments, the mixed micelle systems were not filtered.

Self-micellarization of D 3

Molecular assemblies of D 3 were prepared in Tris buffer using the same protocol as for mixed micelles. D 3 was dissolved into the solvent mixture and after evaporation, the dry film was hydrated for 30 min at 37°C with taurocholate-free buffer. The suspension was then sonicated. All D 3 concentrations reported in the surface tension measurements were obtained from independent micellarization experiments -not from the dilution of one concentrated D 3 solution.

Surface tension measurements

Mixed micelle solutions were prepared as described above, at concentrations ranging from 5.5 nM to 55 mM, with the same proportion of components as previously mentioned. The surface tension of LDP mixtures hydrated with a taurocholate-free buffer, and that of pure taurocholate solutions were also measured at various concentrations. The solutions were poured into glass cuvettes. The aqueous surface was cleaned by suction, and the solutions were left at rest under saturated vapor pressure for 24 hours before measurements. For penetration studies, glass cuvettes with a side arm were used, allowing injection of NaTC beneath a spread LDP or vitamin monolayer. Surface tension measurements were performed by the Wilhelmy plate method, using a thermostated automatic digital tensiometer (K10 Krüss, Germany). The surface tension g was recorded continuously as a function of time until equilibrium was reached. All experiments were performed at 25 ±1°C under saturated vapor pressure to maintain a constant level of liquid. The reported values are mean of three measurements. The experimental uncertainty was estimated to be 0.2 mN/m. Surface pressure (π) values were deduced from the relationship π = γ 0 -γ, with γ 0 the surface tension of the subphase and γ the surface tension in the presence of a film.

Surface pressure measurements

Surface pressure-area π-A isotherms of the LDP and LDP-vitamin mixtures were obtained using a thermostated Langmuir film trough (775.75 cm 2 , Biolin Scientific, Finland) enclosed into a Plexiglas box (Essaid et al. 2016). Solutions of lipids in a chloroform/methanol (9:1, v/v) mixture were spread onto a clean buffer subphase. Monolayers were left at rest for 20 minutes to allow complete evaporation of the solvents. They were then compressed at low speed (6.5 Å 2 .molecule -1 .min -1 ) to minimize the occurrence of metastable phases. The experimental uncertainty was estimated to be 0.1 mN/m. All experiments were run at 25 ±1°C. Mean isotherms were deduced from at least three compression isotherms. The surface compressional moduli K of monolayers were calculated using Eq. 1:

(Eq. 1)

Excess free energies of mixing were calculated according to Eq. 2:

(Eq. 2) with X L and A L the molar fraction and molecular area of lipid molecules, and X VIT and A VIT the molar fraction and molecular area of vitamin molecules, respectively [START_REF] Ambike | Interaction of self-assembled squalenoyl gemcitabine nanoparticles with phospholipid-cholesterol monolayers mimicking a biomembrane[END_REF].

Cryo-TEM analysis

A drop (5 µL) of LDP-NaTC micellar solution (15 mM), LDP-NaTC-D 3 (3:1 molar ratio) or pure D 3 "micellar suspension" (5 mM, theoretical concentration) was deposited onto a perforated carbon-coated, copper grid (TedPella, Inc); the excess of liquid was blotted with a filter paper. The grid was immediately plunged into a liquid ethane bath cooled with liquid nitrogen (180 °C) and then mounted on a cryo holder [START_REF] Da Cunha | Overview of chemical imaging methods to address biological questions[END_REF]. Transmission electron measurements (TEM) measurements were performed just after grid preparation using a JEOL 2200FS (JEOL USA, Inc., Peabody, MA, U.S.A.) working under an acceleration voltage of 200 kV (Institut Curie). Electron micrographs were recorded by a CCD camera (Gatan, Evry, France).

K = -A dπ dA       T ∆ G EXC = X L A L + X VIT A VIT ( ) 0 π ∫ dπ 2.7. Vitamin analysis 2.7.1. Vitamin extraction.
D 3 and RP were extracted from 500 µL aqueous samples using the following method [START_REF] Desmarchelier | The distribution and relative hydrolysis of tocopheryl acetate in the different matrices coexisting in the lumen of the small intestine during digestion could explain its low bioavailability[END_REF]: retinyl acetate was used as an internal standard and was added to the samples in 500 µL ethanol. The mixture was extracted twice with two volumes of hexane.

The hexane phases obtained after centrifugation (1200 × g, 10 min, 10°C) were evaporated to dryness under nitrogen, and the dried extract was dissolved in 200 µL of acetonitrile/dichloromethane/methanol (70:20:10, v/v/v). A volume of 150 µL was used for HPLC analysis. Extraction efficiency was between 75 and 100%. Sample whose extraction efficiency was below 75% were re-extracted or taken out from the analysis.

2.7.2. Vitamin HPLC analysis.

D 3 , and RP and retinyl acetate were separated using a 250 x 4.6-nm RP C18, 5-µm Zorbax

Eclipse XDB column (Agilent Technologies, Les Ulis, France) and a guard column. The mobile phase was a mixture of acetonitrile/dichloromethane/methanol (70:20:10, v/v/v). Flow rate was 1.8 mL/min and the column was kept at a constant temperature (35 °C). The HPLC system comprised a Dionex separation module (P680 HPLC pump and ASI-100 automated sample injector, Dionex, Aix-en-Provence, France). D 3 was detected at 265 nm while retinyl esters were detected at 325 nm and were identified by retention time compared with pure (>95%) standards. Quantification was performed using Chromeleon software (version 6.50, SP4 Build 1000) comparing the peak area with standard reference curves. All solvents used were HPLC grade.

Statistical analysis

Results are expressed as means ± standard deviation. Statistical analyses were performed using Statview software version 5.0 (SAS Institute, Cary, NC, U.S.A.). Means were compared by the non-parametric Kruskal-Wallis test, followed by Mann-Whitney U test as a post hoc test for pairwise comparisons, when the mean difference using the Kruskal-Wallis test was found to be significant (P<0.05). For all tests, the bilateral alpha risk was α = 0.05.

Results

Solubilization of D 3 and RP in aqueous solutions rich in mixed micelles

D 3 and RP at various concentrations were mixed with micelle components (LDP-NaTC).

D 3 and RP concentrations were measured by HPLC before and after filtration of aggregates with a diameter smaller than 0.22 µm (Figure 2). D 3 and RP solubilization in the solution that contained mixed micelle solution followed different curves: D 3 solubilization was linear (R²=0.98, regression slope = 0.71) and significantly higher than that of RP, which reached a plateau with a maximum concentration around 125 µM.

The morphology of the LDP-NaTC and LDP-NaTC-D 3 samples before filtration was analyzed by cryoTEM. In Figure 3, micelles are too small to be distinguished from ice. At high LDP-NaTC concentration (15 mM) small and large unilamellar vesicles (a), nano-fibers (b) and aggregates (c) are observed (Figure 3A). Both nano-fibers and aggregates seem to emerge from the vesicles. In the presence of D 3 at low micelle and D3 concentration (5 mM LDP-NaTC + 1.7 mM D 3 ) (Figures 3B and3C), the morphology of the nano-assemblies is greatly modified.

Vesicles are smaller and deformed, with irregular and more angular shapes (a'). There are also more abundant. A difference in contrast in the bilayers is observed, which would account for leaflets with asymmetric composition. Some of them coalesce into larger structures, extending along the walls of the grid (d). Fragments and sheets are also observed (figure 3B). They exhibit irregular contour and unidentified membrane organization. The bilayer structure is not clearly observable. New organized assemblies appear, such as disk-like nano-assemblies (e) and emulsion-like droplets (f). At higher concentration (15 mM LDP-NaTC + 5 mM D 3 in figure 3D), the emulsion-like droplets and vesicles with unidentified membrane structure (g) are enlarged. They coexist with small deformed vesicules.

Compression properties of LDP components, the LDP mixture and the vitamins

To better understand the mechanism of D 3 and RP interaction with LDP-NaTC micelles, we focused on the interfacial behavior of the various components of the system. We first determined the interfacial behavior of the LDP components and their mixture in proportions similar to those in the micellar solution, by surface pressure measurements. The π-A isotherms are plotted in Figure 4A. Based on the calculated compressibility modulus values, the lipid monolayers can be classified into poorly organized (K < 100 mN/m, for lyso-PC, monoolein, and oleic acid), liquid condensed (100 < K < 250 mN/m, for POPC and the LDP mixture) and highly rigid monolayers (K > 250 mN/m, for cholesterol) [START_REF] Davies | Interfacial phenomena 2nd ed[END_REF].

The interfacial behavior of the two studied vitamins is illustrated in Figure 4B. D 3 shows a similar compression profile to that of the LDP mixture, with comparable surface area and surface pressure at collapse (A c = 35 Å 2 , π c = 38 mN/m) but a much higher rigidity, as inferred from the comparison of their maximal K values (187.4 mN/m and 115.4 mN/m for D 3 and LDP, respectively). RP exhibits much larger surface areas and lower surface pressures than D 3 . The collapse of its monolayer is not clearly identified from the isotherms, and is estimated to occur at π c = 16.2 mN/m (A c = 56.0 Å 2 ), as deduced from the slope change in the π-A plot.

Self-assembling properties of D 3 in an aqueous solution

Since D 3 showed an interfacial behavior similar to that of the lipid mixture, and since it could be solubilized at very high concentrations in an aqueous phase rich in mixed micelles (as shown in Figure 2), its self-assembling properties were more specifically investigated. Dried D 3 films were hydrated with the sodium taurocholate free-buffer. Surface tension measurements at various D 3 concentrations revealed that the vitamin could adsorb at the air/solution interface, and significantly lower the surface tension of the buffer to g cmc = 30.6 mN/m. A critical micellar concentration (cmc = 0.45 µM) could be deduced from the γ-log C relationships and HPLC assays. Concentrated samples D 3 samples were analyzed by cryo-TEM (Figure 3E and3F).

Different D 3 self-assemblies were observed, including circular nano-assemblies (h) coexisting with nano-fibers (i), and large aggregates (j) with unidentified structure. The analysis in depth of the circular nano-assemblies allowed to conclude that they were disk-like nano-assemblies, rather than nanoparticles.

Interaction of LDP with NaTC

To better understand how the two studied vitamins interacted with the mixed micelles, we compared the interfacial behaviors of the pure NaTC solutions, LDP mixtures hydrated by NaTC-free buffer, and LDP mixtures hydrated by the NaTC buffered solutions (full mixed micelle composition). The LDP mixture composition was maintained constant, while its concentration in the aqueous medium was increased. The concentration of NaTC in the aqueous phase was also calculated so that the relative proportion of the various components (LDP and NaTC) remained unchanged in all experiments. From the results plotted in Figure 5A, the critical micellar concentration (cmc) of the LDP-NaTC mixture was 0.122 mM (γ cmc = 29.0 mN/m), a concentration 50.8 times lower than the concentration used for vitamin solubilization.

The cmc values for the LDP mixture and the pure NaTC solutions were 0.025 mM (γ cmc = 24.0 mN/m), and 1.5 mM (γ cmc = 45.3 mN/m), respectively.

Experiments modeling the insertion of NaTC into the LDP film during rehydration by the buffer suggested that only few NaTC molecules could penetrate in the condensed LDP film (initial surface pressure: π i = 28 mN/m) and that the LDP-NaTC mixed film was not stable, as

shown by the decrease in surface pressure over time (Figure 5B).

Interaction of D 3 and RP with NaTC

The surface tension of the mixed NaTC-LDP micelle solutions was only barely affected by the addition of 0.1 or 1 mM D 3 or RP: the surface tension values increased by no more than 2.8 mN/m. Conversely, both vitamins strongly affected the interfacial behavior of the NaTC micellar solution, as inferred from the significant surface tension lowering observed (-7.0 and -8.1 mN/m for RP and D 3 , respectively).

Interaction of D 3 and RP with lipid digestion products

The interaction between the vitamins and LDP molecules following their insertion into LDP micelles was modeled by compression of LDP/D 3 and LDP/RP mixtures at a 7:3 molar ratio. This ratio was chosen arbitrarily, to model a system in which LDP was in excess. The π-A isotherms are presented in Figures 6A and6B. They show that both vitamins modified the isotherm profile of the lipid mixture, however, not in the same way. In the LDP/D 3 mixture, the surface pressure and molecular area at collapse were controlled by LDP. For LDP/RP, despite the high content in LDP, the interfacial behavior was clearly controlled by RP. From the isotherms in Figures 6A and6B, compressibility moduli and excess free energies of mixing were calculated and compared (Figures 6C, and6D). D 3 increased the rigidity of LDP monolayers, whereas RP disorganized them. The negative ∆G EXC values calculated for the LDP-D 3 monolayers at all surface pressures account for the good mixing properties of D 3 and the lipids in all molecular packing conditions. Conversely for RP, the positive and increasing ∆G EXC values with the surface pressure demonstrate that its interaction with the lipids was unfavorable.

Discussion

The objective of this study was to compare the solubility of RP and D 3 in aqueous solutions containing mixed micelles, and to decipher the molecular interactions that explain their different extent of solubilization. Our first experiment revealed that the two vitamins exhibit very different solubilities in an aqueous medium rich in mixed micelles. Furthermore, the solubility of D 3 was so high that we did not observe any limit, even when D 3 was introduced at a concentration > 1mM in the aqueous medium. To our knowledge, this is the first time that such a difference is reported. Cryo-TEM pictures showed that D 3 dramatically altered the organization of the various components of the mixed micelles. The spherical vesicles were deformed with angular shapes. The nano-fibers initiating from the vesicles were no longer observed. Large irregular in shape vesicle and sheets, disk-like nano-assemblies and emulsionlike droplets appeared in LDP-NaTC-D 3 mixtures, only. The existence of so many different assemblies would account for a different interaction of D 3 with the various components of mixed micelles, and for a reorganization of the components. D 3 could insert in the bilayer of vesicles and deform them, but also form emulsion-like droplets with fatty acids and monoglyceride. It is noteworthy that these emulsion-like droplets were not observed in pure D 3 samples, nor mixed micelles. Since previous studies have shown that both bile salts and some mixed micelle lipids, e.g. fatty acids and phospholipids, can modulate the solubility of fatsoluble vitamins in these vehicles [START_REF] Yang | Vitamin E and vitamin E acetate solubilization in mixed micelles: physicochemical basis of bioaccessibility[END_REF], we decided to study the interactions of these two vitamins with either bile salts or micelle lipids to assess the specific role of each component on vitamin solubility in mixed micelles.

The characteristics of pure POPC, Lyso-PC, monoolein, and cholesterol isotherms were in agreement with values published in the literature [START_REF] Pezron | Monoglyceride Surface-Films -Stability and Interlayer Interactions[END_REF]Flasinsky et al. 2014;[START_REF] Huynh | Structural properties of POPC monolayers under lateral compression: computer simulations analysis[END_REF]. For oleic acid, the surface pressure at collapse was higher (π c = 37 mN/m) and the corresponding molecular area (A c = 26 Å 2 ) smaller than those previously published [START_REF] Tomoaia-Cotisel | Insoluble mixed monolayers[END_REF], likely due to the pH of the buffer solution (pH 6) and the presence of calcium.

The interfacial properties of D 3 were close to those deduced from the isotherm published by [START_REF] Meredith | The Supramolecular Structure of Vitamin-D3 in Water[END_REF] for a D 3 monolayer spread from a benzene solution onto a pure water subphase. The molecular areas at collapse are almost identical in the two studies (about 36 Å 2 ), but the surface pressures differ (30 mN/m in Meredith and coworkers' study, and 38 mN/m in ours). Compressibility modulus values show that D 3 molecules form monolayers with higher molecular order than the LDP mixture, which suggests that they might easily insert into LDP domains.

As could be expected from its chemical structure, RP exhibited a completely different interfacial behavior compared to D 3 and the LDP, even to lyso-PC which formed the most expanded monolayers of the series, and displayed the lowest collapse surface pressure. The anomalous isotherm profile of lyso-PC has been attributed to monolayer instability and progressive solubilization of molecules into the aqueous phase [START_REF] Heffner | Thermodynamic and kinetic investigations of the release of oxidized phospholipids from lipid membranes and its effect on vascular integrity[END_REF]. The molecular areas and surface pressures for RP have been compared to those measured by [START_REF] Asai | Formation and stability of the dispersed particles composed of retinyl palmitate and phosphatidylcholine[END_REF] for RP monolayers spread from benzene solutions at 25°C onto a water subphase. Their values are much lower than ours, accounting for even more poorly organized monolayers. The low collapse surface pressure could correspond to molecules partially lying onto the aqueous surface, possibly forming multilayers above 16 mN/m as inferred from the continuous increase in surface pressure above the change in slope of the isotherm. The maximal compressibility modulus confirms the poor monolayer order. The significant differences in RP surface pressure and surface area compared to the LDP mixture might compromise its insertion and stability into LDP domains.

The dogma in nutrition is that fat-soluble vitamins need to be solubilized in bile salt micelles to be transported to the enterocyte and then absorbed. It is also well known that although NaTC primary micelles can be formed at 2-3 mM with a small aggregation number, concentrations as high as 10-12 mM are usually necessary for efficient lipid solubilization in the intestine [START_REF] Baskin | Bile salt-phospholipid aggregation at submicellar concentrations[END_REF]. Due to their chemical structure bile salts have a facial arrangement of polar and non-polar domains (Madency & Egelhaaf 2010). Their selfassembling (dimers, multimers, micelles) is a complex process involving hydrophobic interaction and cooperative hydrogen bonding, highly dependent on the medium conditions, and that is not completely elucidated. The cmc value for sodium taurocholate in the studied buffer was particularly low compared to some of those reported in the literature for NaTC in water or sodium chloride solutions (3-12 mM) [START_REF] Kratohvil | Concentration-dependent aggregation patterns of conjugated bile-salts in aqueous sodiumchloride solutions -a comparison between sodium taurodeoxycholate and sodium taurocholate[END_REF][START_REF] Meyerhoffer | Critical Micelle Concentration Behavior of Sodium Taurocholate in Water[END_REF][START_REF] Madenci | Self-assembly in aqueous bile salt solutions[END_REF]. At concentrations as high as 10-12 mM, NaTC molecules form elongated cylindrical "secondary" micelles [START_REF] Madenci | Self-assembly in aqueous bile salt solutions[END_REF][START_REF] Bottari | Structure and composition of sodium taurocholate micellar aggregates[END_REF]. The cryoTEM analysis did not allow to distinguish micelles from the ice. In our solubilization experiment, the concentration of NaTC did not exceed 5 mM. Nevertheless, the micelles proved to be very efficient with regards to vitamin solubilization.

When bile salts and lipids are simultaneously present in the same environment, they form mixed micelles [START_REF] Hernell | Physical-chemical behavior of dietary and biliary lipids during intestinal digestion and absorption. 2. Phase analysis and aggregation states of luminal lipids during duodenal fat digestion in healthy adult human beings[END_REF]. Bile salts solubilize phospholipid vesicles and transform into cylindrical micelles [START_REF] Cheng | Mixtures of lecithin and bile salt can form highly viscous wormlike micellar solutions in water[END_REF]. [START_REF] Walter | Intermediate structures in the cholate-phosphatidylcholine vesicle-micelle transition[END_REF] suggested that sodium cholate cylindrical micelles evolved from the edge of lecithin bilayer sheets. Most published studies were performed at high phospholipid/bile salt ratio. In our system, the concentration of the phospholipids was very low compared to that of NaTC. We observed however the presence of vesicles, and nano-fiber structures emerging from them. In their cryoTEM analysis, [START_REF] Fatouros | Colloidal structures in media simulating intestinal fed state conditions with and without lypolysis products[END_REF] compared bile salt/phospholipid mixtures to bile salt/phospholipid/fatty acid/monoglyceride ones at concentrations closer to ours. They observed only micelles in bile salt/phospholipid mixtures. However, in the presence of oleic acid and monoolein, vesicles and bilayer sheets were formed. This would account for a reorganization of the lipids and bile salts in the presence of the fatty acid and the monoglyceride.

We therefore decided to study the interactions between bile salts and LDP. The results obtained

show that the surface tension, the effective surface tension lowering concentration, and cmc Solubilization experiments and the analysis of vitamin-NaTC interaction cannot explain why the LDP-NaTC mixed micelles solubilize D 3 better than RP. Therefore, we studied the interfacial behavior of the LDP mixture in the presence of each vitamin, to determine the extent of their interaction with the lipids. The results obtained showed that D 3 penetrated in LDP domains and remained in the lipid monolayer throughout compression. At large molecular areas, the π-A isotherm profile of the mixture followed that of the LDP isotherm with a slight condensation due to the presence of D 3 molecules. Above 10 mN/m, an enlargement of the molecular area at collapse and a change in the slope of the mixed monolayer was observed.

However, the surface pressure at collapse was not modified, and the shape of the isotherm accounted for the insertion of D 3 molecules into LDP domains. This was confirmed by the surface compressional moduli. D 3 interacted with lipid molecules in such manner that it increased monolayer rigidity (K max = 134.8 mN/m), without changing the general organization of the LDP monolayer. The LDP-D 3 mixed monolayer thus appeared more structured than the LDP one. D 3 behavior resembles that of cholesterol in phospholipid monolayers, however without the condensing effect of the sterol [START_REF] Ambike | Interaction of self-assembled squalenoyl gemcitabine nanoparticles with phospholipid-cholesterol monolayers mimicking a biomembrane[END_REF]. The higher rigidity of LDP monolayer in the presence of D 3 could be related to the cryo-TEM pictures showing the deformed, more angular vesicles formed with LDP-NaTC-D 3 . The angular shape would account for vesicles with rigid bilayers [START_REF] Kuntsche | Cryogenic transmission electron microscopy (cryo-TEM) for studying the morphology of collidal drug delivery systems[END_REF].

For RP, the shape of the isotherms show evidence that lipid molecules penetrated in RP domains, rather than the opposite. Indeed, the π-A isotherm profile of the LDP-RP monolayer is similar to that of RP alone. The insertion of lipid molecules into RP domains is also attested by the increase in the collapse surface pressure from 16 to 22 mN/m. Partial collapse is confirmed by the decrease in the compressibility modulus above 22 mN/m. Thus, RP led to a destructuration of the LDP mixed monolayer and when the surface density of the monolayer increased, the vitamin was partially squeezed out from the interface. The calculated ∆G EXC values for both systems suggest that insertion of D 3 into LDP domains was controlled by favorable (attractive) interactions, whereas mixing of RP with LDP was limited due to unfavorable (repulsive) interactions, even at low surface pressures. According to [START_REF] Asai | Formation and stability of the dispersed particles composed of retinyl palmitate and phosphatidylcholine[END_REF], RP can be partially solubilized in the bilayer of phospholipids (up to 3 mol%), and the excess is separated from the phospholipids, and dispersed as emulsion droplets stabilized by a phospholipid monolayer.

On the whole, the information obtained regarding the interactions of the two vitamins with NaTC and LDP explain why D 3 is more soluble than RP in an aqueous medium rich in mixed micelles. Both vitamins can insert into pure NaTC domains, but only D 3 can also insert into the LDP domains in LDP-enriched NaTC micelles.

Furthermore, the results obtained suggest that this is not the only explanation. Indeed, since it has been suggested that D 3 could form cylindrical micelle-like aggregates [START_REF] Meredith | The Supramolecular Structure of Vitamin-D3 in Water[END_REF], we hypothesize that the very high solubility of D 3 in the aqueous medium rich in mixed micelles was partly due to the solubilization of a fraction of D 3 as self-aggregates. Indeed, we observed that D 3 at concentrations higher than 0.45 µM, could self-assemble into various structures including nano-fibers. To our knowledge, no such structures, especially nanofibers, have been reported for D 3 so far. Rod diameter was smaller than 10 nm, much smaller than for the rods formed by lithocholic acid, for example [START_REF] Terech | Self-assembled monodisperse steroid nanotubes in water[END_REF]. They were similar to those observed in highly concentrated LDP-NaTC mixtures, which seemed formed via desorganization of lipid vesicles. Disk-like and aggregates with unidentified structure, also observed in concentrated D 3 samples, could be related to these nano-fibers.

In our solubilization experiments, which were performed at much higher D 3 concentrations, both insertion of D 3 molecules into NaTC and LDP domains, and D 3 self-assembling could occur, depending on the kinetics of insertion of D 3 into the NaTC-DLP mixed micelles.

Conclusion

The solubilization of a hydrophobic compound in bile salt-lipid micelles is dependent upon its chemical structure and its ability to interact with the mixed micelles components. Most hydrophobic compounds are expected to insert into the bile salt-lipid micelles. The extent of the solubilizing effect is, however, much more difficult to predict. As shown by others before us, mixed micelles components form a heterogeneous system with various molecular assemblies differing in shape and composition. The conditions of the medium (pH, ionic strength and temperature) affect the formation of these molecular assemblies, although we did not study this effect on our system. Our results showed that D 3 displayed a higher solubility in mixed micelle solutions than RP. This difference was attributed to the different abilities of the two vitamins to insert in between micelle components, but it was also explained by the propensity of D 3 , contrarily to RP, to self-associate into structures that are readily soluble in the aqueous phase. It is difficult to predict the propensity of a compound to self-association. We propose here a methodology that was efficient to distinguish between two solubilizing behaviors, and could be easily used to predict the solubilization efficiency of other hydrophobic compounds. Whether the D 3 self-assemblies are available for absorption by the intestinal cells needs further studies.

  values were very much influenced by LDP. The almost parallel slopes of Gibbs adsorption isotherms for pure NaTC and mixed NaTC-LDP suggest that LDP molecules inserted into NaTC domains, rather than the opposite. This was confirmed by penetration studies, which showed that NaTC (0.1 mM) could hardly penetrate in a compact LDP film. So, during lipid hydration, LDP molecules could insert into NaTC domains. The presence of LDP molecules improved NaTC micellarization.After having determined the interfacial properties of each micelle component and measured the interactions between NaTC and LDP, we assessed the ability of D 3 and RP to solubilize in either NaTC or NaTC-LDP micelles. Surface tension values clearly show that both vitamins could insert in between NaTC molecules adsorbed at the interface, and affected the surface tension in the same way. The interfacial behavior of the molecules being representative of their behavior in the bulk, it is reasonable to think that both D 3 and RP can be solubilized into pure NaTC micelles. For the mixed NaTC-LDP micelles, the change in surface tension was too limited to allow conclusions, but the solubilization experiments clearly indicated that neither vitamin was solubilized to the same extent.

Figure 1 :

 1 Figure 1: Chemical structures for D 3 and RP.

Figure 2 :

 2 Figure 2: Solubilization of D 3 and RP in aqueous solutions rich in mixed micelles: (l),

Figure 3 :

 3 Figure 3: Cryo-TEM morphology of (A) 15 mM mixed LDP-NaTC micelles, (B) and (C) 5

Figure 4 :

 4 Figure 4: Mean compression isotherms for (A) the pure micelles components and the LDP

Figure 5 :

 5 Figure 5: (A) Adsorption isotherms for LDP hydrated in NaTC-free buffer (○), LDP hydrated

Figure 6 :

 6 Figure 6: π-A isotherms (A,B), compressibility moduli (C) and excess free energies (D) for

Ø

  Figure 1

Acknowledgements:

The authors are grateful to Dr Sylvain Trépout (Institut Curie, Orsay, France) for his contribution to cryoTEM experiments and the fruitful discussions.

Funding: This study was funded by Adisseo France SAS.

Conflicts of interest: DP, ED and VLD are employed by Adisseo. Adisseo markets formulated vitamins for animal nutrition.