A Pashkevich
email: anatol.pashkevich@emn.fr

A Klimchik

S Briot

D Chablat
email: damien.chablat@irccyn.ec-nantes.fr

Performance evaluation of parallel manipulators for milling application

Keywords: Performance evaluation, Kinetostatic modeling, Elastic errors, Parallel manipulators, Milling application

This paper focuses on the performance evaluation of the parallel manipulators for milling of composite materials. For this application the most significant performance measurements, which denote the ability of the manipulator for the machining are defined. In this case, optimal synthesis task is solved as a multicriterion optimization problem with respect to the geometric, kinematic, kinetostatic, elastostostatic, dynamic properties. It is shown that stiffness is an important performance factor. Previous models operate with links approximation and calculate stiffness matrix in the neighborhood of initial point. This is a reason why a new way for stiffness matrix calculation is proposed. This method is illustrated in a concrete industrial problem.

INTRODUCTION

Currently, parallel manipulators have become more and more popular for a variety of technological processes, including high-speed precision machining [START_REF] Brogardh | Present and future robot control developmentan industrial perspective[END_REF] [START_REF] Chanal | A study of the impact of machine tool structure on machining processes[END_REF]. This growing attention is inspired by their essential advantages over serial manipulators, which have already reached the dynamic performance limits. In contrast, parallel manipulators are claimed to offer better accuracy, lower mass/inertia properties, and higher structural rigidity (i.e. stiffness-to-mass ratio) [START_REF] Merlet | Parallel Robots[END_REF]. These features are induced by their specific kinematic structure, which resists to the error accumulation in kinematic chains and allows convenient actuators location close to the manipulator base. This makes them attractive for innovative robotic systems, but practical utilization of the potential benefits requires development of efficient stiffness analysis techniques, which satisfy the computational speed and accuracy requirements of relevant design procedures. Generally, the stiffness analysis evaluates the effect of the applied external torques and forces on the compliant displacements of the end-effector. Numerically, this property is defined through the "stiffness matrix", which gives the relation between the translational/rotational displacement and the static forces/torques causing this transition [START_REF] Koseki | Kinematic Analysis of Translational 3-DOF Micro Parallel Mechanism Using Matrix Method[END_REF]. Similar to other manipulator properties (kinematical, for instance), the stiffness essentially depends on the force/torque direction and on the manipulator configuration [START_REF] Alici | Enhanced stiffness modeling, identification and characterization for robot manipulators[END_REF]. Several approaches exist for the computation of the stiffness matrix, such as the Finite Element Analysis (FEA), the matrix structural analysis (MSA), and the virtual joint method (VJM). Each of methods has their own assumptions and in general approximate/linearized the stiffness model with adequate accuracy for the specific problem. Let us consider the recent modification of the VJM proposed by the authors [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF]. It allows extending it to the over-constrained manipulators and applying it at any workspace point, including the singular ones. The method is based on a multidimensional lumpedparameter model that replaces the link flexibility by localized 6-dof virtual springs that describe both the linear/rotational deflections and the coupling between them. The spring stiffness parameters are evaluated using FEA modelling to take into account real shape of the manipulator components. This gives almost the same accuracy as FEA but with essentially lower computational effort because it eliminates re-meshing through the workspace.

In contrast to previous works, the novelty of this paper appears in the computation of the stiffness matrix for the loaded mode in the neighborhood of static equilibrium with true dimensions links. This paper focuses on the performance evaluation of the milling process for the composite materials via enhancing the precise stiffness model for the compliance errors estimation. VJM is used for the stiffness modeling. It incorporates accurate stiffness properties of links via improved FEA-based stiffness calculation algorithm. This leads to adequate estimation of deflections and allows improving accuracy of manufacturing via accurate planning of the technological process. The reminder of the paper is organized as follows. In the Section 2, it is defined the set of problems which are considered in the paper. Section 3 presents performance evaluation factors. Section 4 proposes a new method for the stiffness modeling of parallel manipulator. Section 5 illustrate the efficiency of the propose technique on the Orthoglide manipulator by specializing it for milling of bathroom component. And, finally, section 6 summarizes the main contributions of the paper.

PROBLEM STATEMENT

The main problem to be solved in this article is the accuracy evaluation during the manufacturing step for specific parallel manipulators and technological process. Position accuracy and performances characteristics for the industrial manipulators depend on the technological process, kinetostatic and geometry properties of both manipulator and workpiece. Therefore specializing of the manipulator for the current manufacturing process with respect to the required performances properties allows increasing the position accuracy of the tool. Mechanical stiffness is one of the most important properties, which have influence on the position accuracy of low mass parallel mechanisms. As it is mentioned in [START_REF] Pashkevich | Enhanced stiffness modeling of serial manipulators with passive joints, Intech, Austria[END_REF] Cartesian stiffness is a nonlinear function of the external loading and depends on the configuration of the manipulator and, as result, of the position of end effector in the workspace [START_REF] Pashkevich | Stiffness analysis of multichain parallel robotic systems with loading[END_REF]. The computed Cartesian stiffness both depends on the stiffness of the links and stiffness model accuracy. So the aim of this article is to specialize a parallel manipulator for the milling of products made of composite material with known geometric sizes. The idea is to improve the accuracy for the high-speed milling via mass reduction and improving of the stiffness model.

PERFORMANCE EVALUATION

Let us present here the most essential technologyoriented performance measures that are used in the design process of a machine-tool. They may be classified into two groups:

 the measures based on geometric, kinematic and kinetostatic properties; these measures are based on simple models that evaluate the geometric, kinematic and kinetostatic properties of a mechanism. Their specificity is that they only use the primary geometric parameters of the mechanism, i.e., the size of the base, the length of links, etc.

 the measures based on dynamic, elastostatic and elastodynamic properties; these measures mostly use more complicated models. Their specificity is that they not only use the primary geometric parameters of the mechanism, but also the secondary geometric parameters, as for example, the cross-section of links. Let us now describe them.

Geometric, kinematic and kinetostatic performances

The most commonly used geometric, kinematic and kinetostatic performances are:

 the velocity transmission factors;

For appreciating the speed capability of a manipulator, several kinematic performance indices are defined using the Jacobian matrix J (see [START_REF] Merlet | Jacobian, manipulability, condition number, and accuracy of parallel robots[END_REF]), such as the condition number, the largest/smallest singular value (also called the maximal/minimal transmission factor and denoted as , respectively -Figure 1a), the dexterity, the manipulability, etc. However, as mentioned by Merlet in [START_REF] Merlet | Jacobian, manipulability, condition number, and accuracy of parallel robots[END_REF], the previously cited indices does not take into account the "technological reality" of the mechanism, as they are based on the use of the Euclidian norm of the input velocity vector Φ  (|| Φ  || being considered equal to 1) while it is clear that each actuator may have a velocity

i    [-max i   , max i  
], where i   and max i   are the actual and maximal velocities for the actuator i (i = 1 to n). Thus, it is necessary to redefine the transmission factors.

Two types of transmission factors may be used: (i) the velocity transmission factors along all directions of the workspace; in such a case, the unit square is mapped into a parallelepiped (Figure 1b) and (ii) the velocity transmission along some particular directions of the workspace; in such a case, the unit square is mapped into a hexagon (Figure 1c) [START_REF] Briot | Optimal Technology-Oriented Design of Parallel Robots for High-Speed Machining Applications[END_REF].

 the accuracy transmission factors;

The accuracy of a mechanism may be due to different factors. However, as pointed out by Merlet [START_REF] Merlet | Computing the worst case accuracy of a PKM over a workspace or a trajectory[END_REF], activejoint errors are the most significant source of errors in a properly designed, manufactured, and calibrated parallel robot. The classical approach consists in considering the first order approximation that maps the input error to the output error:

Φ J p δ δ  (1)
where  represents the vector of the active-joint errors, p the vector of end-effector errors. This method will give only an approximation of the end-effector maximum error, but at an optimisation stage, this model may be sufficient. It is possible to show that, due the use of the expression (1) as the accuracy model, the accuracy transmission factors are similar to the velocity transmission factors. More information about the computation of these transmission factors is proposed in [START_REF] Briot | Optimal Technology-Oriented Design of Parallel Robots for High-Speed Machining Applications[END_REF].

 the force transmission factors;

For 3-DOF translational parallel manipulators, the input efforts  are related to the forces f applied on the platform by the following relation:

τ J f T   (2)
Looking at expression (2), it appears that the relationship for the accuracy is similar to the relationship for the velocity, but J -T is used instead to J. Moreover, it is also clear that the effort i of one actuator is comprised between - max and + max ,  max being the maximal effort admissible by the actuated pair. So, the force transmission factors may be computed in the same manner as the velocity transmission factors. The maximal and minimal force and moment transmission factors will be denoted as min respectively. To have more information about the computation of these transmission factors, the reader is proposed to refer to [START_REF] Briot | Optimal Technology-Oriented Design of Parallel Robots for High-Speed Machining Applications[END_REF].

 the size of the workspace;

Using the set of geometric parameters, the workspace W may be generated using the kinematic equations and the (passive and active) joint limits [START_REF] Merlet | Parallel Robots[END_REF]. Since, for the considered application, the desired regular workspace is a parallelepiped W0 of size {a0 × b0 × c0}, the relevant measure may be defined by the largest similar object

W abc = {µa0 × µb0 × µc0} inscribed in W, i.e.   W W T T W T W 0 T 0    ) ( max arg ) , ( ); ( ,      abc (3)
where µ, T are respectively the scalar scaling factor and the coordinate transformation operator in the Cartesian space. An algorithm that is able to compute the size largest parallelepiped inside a given workspace is presented in [START_REF] Briot | Optimal Technology-Oriented Design of Parallel Robots for High-Speed Machining Applications[END_REF] Let us now introduce the dynamic, elastostatic and elastodynamic performances.

Dynamic, elastostatic and elastodynamic performance

The most commonly used elastostatic and elastodynamic performances are:

 the total mass of the robot; This is probably the simplest performance measure. The masses of the different elements of a robot have a direct influence on its dynamic behaviour. For one given structure with fixed value of links length, the robot that will have the smallest mass will be the one with the smallest input efforts, and as a result, the highest acceleration capacity.

 the input efforts;

The input efforts, denoted as , depend on the mass and axial moment of inertia of links, such as friction in joints and position, velocity and acceleration of the robot. There expression is given by:

f J q G q q q C q q M τ T     ) ( ) , ( ) (     (4)
where M is the mass matrix, C the matrix of the Coriolis, centrifugal and viscous friction effects, G is the matrix of the gravity and Coulomb friction effects. q, q  and q  represent the vectors of the positions, velocities and accelerations of the actuators. f is an external force applied on the platform. Generally, it is considered that, given a desired trajectory, the robot that has the lowest input efforts along the trajectory has the best performance.

 the maximal deformations;

For parallel manipulators, elasticity is an essential performance measure since it is directly related to the positioning accuracy and the payload capability. Mathematically, this benchmark is defined by the stiffness matrix K, which describes the relation between the linear/angular displacements t of the end-effector and the external forces/torques f applied on the tool:

t K f δ  (5)
It is obvious that elasticity is highly dependent upon geometry, materials and link shapes that are completely defined with the CAD model. The stiffness matrix may be computed using three methods: the finite element analysis (FEA) [START_REF] El-Khasawneh | Computation of Stiffness and Stiffness Bounds for Parallel Link Manipulators[END_REF][START_REF] Bouzgarrou | Rigidity Analysis of T3R1 Parallel Robot Uncoupled Kinematics[END_REF], the matrix structural analysis (SMA) [START_REF] Deblaise | Systematic Analytical Method for PKM Stiffness Matrix Calculation[END_REF][START_REF] Ghali | Structural Analysis: A Unified Classical and Matrix Approach[END_REF] and the virtual joint method (VJM) [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF][START_REF] Gosselin | Stiffness Analysis of Parallel Mechanisms Using a Lumped Model[END_REF][START_REF] Majou | Parametric Stiffness Analysis of the Orthoglide[END_REF]. Whatever will be the way to compute the deformations of the robot, the mechanism that will have the best elastostatic performances will be the one that will have the smallest deformations under a given force, along a specific trajectory or in the totality (or some portions) of the workspace. External forces are caused by the coupling of the tool and machining piece. But previous stiffness models calculate stiffness matrix in the neighborhood of initial point and operate with links approximations. Such models can"t guarantee good accuracy. Therefore let us consider the stiffness model improvement and deflections evaluation.

STIFFNESS MODEL

Problem of stiffness modeling

To evaluate the manipulator stiffness, let us apply the VJM method that assumes that the traditional rigid model is extended by adding virtual joints, which describe stiffness of the actuator and links. Thus, the end-effector position for each chain of the manipulator can be described as ( , )  t g q θ [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF] where (...) g is the geometry function which depends on the passive q and virtual joint θ coordinates, the vectors 12 ( , , ..., ) T n q q q  q includes all passive joint coordinates, the vector 12 ( , , ...,

) T m  θ    collects all virtual joint coordinates,
n is the number of passive joins, m is the number of virtual joints. This expression includes both traditional geometric variables (passive and active joint coordinates) and stiffness variables (virtual joint coordinates). To evaluate the manipulator ability to respond to external forces and torques, it is necessary to introduce additional equations that define the virtual joint reactions to the corresponding spring deformations. For analytical convenience, corresponding expressions may be collected in a single matrix equation

 θθ τ K θ (7)
where

  ,1 , 2 , , , ..., T m  θ τ      
is the aggregated vector of the virtual joint reactions,

  , , ..., diag  θ θ,1 θ, 2 θ, m K K K K
is the aggregated spring stiffness matrix of the size mm, and ,i θ K is the spring stiffness matrix of the corresponding link. For the compliant link, the matrix K can be computed using the FEA-based techniques, which usually produce rather accurate result. Using the FEA, the stiffness matrix K (or its inverse k) is evaluated from several numerical experiments, each of which produces the vectors of linear and angular deflections (t, ) corresponding to the applied force and torque (F, M). Then, the desired matrix is computed from the linear system

mm mm               F F t t k MM φφ (8)
where []  is pseudoinverse of the rectangle matrix, m is the number of experiments (

6  m
) and the matrix pseudoinverse is replaced by the inverse in the case of

6  m
. It is obvious that this case with special arrangement of the forces and torques is numerically attractive (for more detail see [START_REF] Pashkevich | Accuracy Improvement for Stiffness Modeling of Parallel Manipulators[END_REF]). In order to increase accuracy it is worth to improve the deflection estimation technique. It is proposed to evaluate (t, ) from the displacement field describing transitions of rather large number of nodes located in the neighborhood of reference point (RP).

To formulate this problem strictly, let us denote the displacement field by a set of vector couples

{ , | 1, 2, ..., } ii in  pp
where the first component i p define the node initial location (before applying the force/torque), i p  refers to the node displacement due to the applied force/torque, and n is the number of considered nodes. Then, assuming that all the nodes are close enough to the reference point, this set can be approximated by a "rigid transformation"

t p R p p      i i i ) ( , 1, 2,..., in  (9)
that includes as the parameters the linear displacement t and the orthogonal 33 matrix R that depends on the rotational displacement . Then, the problem of the deflection estimation can be presented as the best fit of the considered vector field by equation ( 4) with respect to six scalar variables incorporated in t, R.

In general case, the desired stiffness model is defined by a non-liner relation [START_REF] Merlet | Jacobian, manipulability, condition number, and accuracy of parallel robots[END_REF] that describes resistance of a mechanism to deformations Δt caused by an external force/torque F [START_REF] Brogardh | Present and future robot control developmentan industrial perspective[END_REF]. It should be noted that the mapping Δt  F is strictly mathematically defined and physically tractable in all cases, including under-constrained kinematics and singular configurations of the manipulator. However, the converse is not true.

()  F ξ Δ t
In engineering practice, function (...) ξ is usually linearized in the neighborhood of the static equilibrium ( , ) q θ corresponding to the end-effector position t and external loading F . For the unloaded mode, i.e. when  F0 and  0 θ0 the stiffness model is expressed by a simple relation

( , )  00 F K q θ Δ t ( 11
)
where K is 66 ""stiffness matrix" and the vector 0 01 02 0 ( , , ..., ) T n q q q  q defines the equilibrium configuration corresponding to the end-effector location 0 t , in accordance with the manipulator geometry. However, for the loaded mode, stiffness model have to be defined in the neighborhood of the static equilibrium that corresponds to another manipulator configuration ( , ) q θ , which is caused by external forces F . In this case, the stiffness model describes the relation between the increments of the force δF and the position δt ( , )  F K q θt  [START_REF] El-Khasawneh | Computation of Stiffness and Stiffness Bounds for Parallel Link Manipulators[END_REF] where  0 qq Δq and  0 θ θ Δ θ denote the new position of the manipulator, Δq and Δθ are the deviations of the passive joint and virtual spring coordinates. Hence, the problem of the stiffness modeling in the loaded mode may be divided into two sequential subtasks: (i) stiffness model identification from the vector field of displacements and (ii) linearization of relevant force/position relations in the neighborhood of the loaded configuration. Let us consider these two sub-problems consequently.

Stiffness model identification

To estimate the desired deflections (t, ), let us apply the least square technique that leads to minimization of the sum of squared residuals

t R t p R p p , n i i i i f min ) ( 1 2          (13)
with respect to the vector t and the orthogonal matrix R representing the rotational deflections . The specificity of this problem (that does not allow direct application of the standard methods) are the orthogonally constraint

I R R  T
and non-trivial relation between elements of the matrix R and the vector . To reduce the computational efforts, let us linearize the rotational matrix R [START_REF] Pashkevich | Accuracy Improvement for Stiffness Modeling of Parallel Manipulators[END_REF]. This allows to rewrite equation of the "rigid transformation" (9) in the form

n i i i , 1 ;      t p p  (14)
that can be further transformed into a linear system of the following form

  n i i i , 1 ;          p t P I  (15)
where i P is a skew-symmetric matrix corresponding to the vector i p . Then, applying the standard least-square technique and shifting the origin of the coordinate system to the point

    n i i c n 1 1 p p leading to expression                                                    n i i T i n i i n i i T i n 1 1 1 1 1 ˆp P p P P 0 0 I t  (16)
that requires inversion of the matrix of size 33. Here, following the adopted notation i P ˆ is a skew-symmetric matrix corresponding to the vector

c i i p p p   ˆ.
By its general principle, the FEA-modeling is an approximate method that produces some errors caused by the discretization. Beside, even for the perfect modeling, the deflections in the neighborhood of the reference point do not exactly obey the equation ( 4). Hence, it is reasonable to assume that the "rigid transformation" (9) incorporates some random errors

n i i i i i , 1 ; ) (        ε t p R p p  (17)
that are supposed to be independent and identically distributed Gaussian random variables with zero-mean and standard deviation .

In the frame of this assumption, the expression for the deflections ( 16) can be rewritten as

                    n i i T i n i i T i n i i n 1 1 1 ο 1 1 o ; ε P P P ε t t   (18)
where the superscript "o" corresponds to the "true" parameter value. This justifies usual properties of the adopted point-type estimator [START_REF] Gosselin | Stiffness Analysis of Parallel Mechanisms Using a Lumped Model[END_REF], which is obviously unbiased and consistent. Furthermore, the variancecovariance matrices for t,  may be expressed as

    1 1 2 2 ĉov ; cov              n i i T i n P P I t    (19)
allowing to evaluate the estimation accuracy using common confidence interval technique. Another practical question is related to detecting zero elements in the compliance matrix or, in other word, evaluating the statistical significance of the computed values compared to zero. Relevant statistical technique [START_REF] Akin | Finite Element Analysis With Error Estimators: An Introduction to the FEM and Adaptive Error Analysis for Engineering Students[END_REF] operates with the p-values that may be easily converted in the form a k  , where k is usually from 3 to 5 and the subscript "

a " refers to a particular component of the vectors t, .

To evaluate the standard deviation  describing the random errors ε , one may use the residual-based estimator obtained from the expression

2 1 2 ) 6 3 ( ) ( E                   n n i i i i t p R p p  . ( 20
)
The latter may be easily derived taking into account that, for each experiment, the deflection filed consist of n three-dimensional vectors that are approximated by the model containing 6 scalar parameters. Moreover, to increase accuracy, it is prudent to aggregate the squared residuals for all FEA-experiments and to make relevant estimation using the coefficient

2 ) 6 3 (  m n 
, where m is the experiments number. In addition, to increase accuracy and robustness, it is reasonable to eliminate outliers in the experimental data. They may appear in the FEA-field due to some anomalous causes, such as insufficient meshing of some elements, violation of the boundary conditions in some areas of the mechanical joints, etc. The simplest and reliable method that is adopted in this paper is based on the "data filtering" with respect to the residuals.

Stiffness model in the loaded mode

To compute the desired stiffness matrix, let us consider the neighborhood of the loaded configuration and assume that the external force and the end-effector location are by some values F , t . Besides, let us assume that a new configuration satisfies the equilibrium conditions [START_REF] Pashkevich | Enhanced stiffness modeling of serial manipulators with passive joints, Intech, Austria[END_REF]. Hence, it is necessary to consider simultaneously two equilibriums corresponding to the manipulator state variables ( , , , ) Fqθt and

( , , , )         F F q q θ θ t t .
Relevant equations of statics may be written as

;0 TT q       F J K θ F J (21)
and

          ; 0 T T                θ θ θ qq F F J J K θθ F F J J (22) where ( , ) q Jq θ and ( , )
 Jq θ are the differentials of the Jacobians due to changes in ( , ) q θ . Besides, in the neighborhood of ( , ) q θ , the kinematic equation may be also presented in the linearized form:

( , ) ( , ) q      δt J q θ δθ J q θ δq , (23)
Hence, after neglecting the high-order small terms and expending the differentials via the Hessians of the function ( , ) T gqθF , equations ( 21), ( 22) may be rewritten as

( ) ( ) ( ) ( ) ( ) ( )
T F F q T F F q qq q                          J q, θ F H q, θ q H q, θ θ K θ J q, θ F H q, θ q H q, θ θ 0 (24)
and the general relation for the stiffness matrix in the loaded mode can be presented as

1 F T F F Fq qq T F F T F F F F q q qq q q                          KK J k J J J k H J H k J H H k H (25)
Hence, the presented technique allows computing the stiffness matrix in the presence of the external load and to generalize previous results both for serial kinematic chains and for parallel manipulators.

APPLICATION EXAMPLE

Industrial problem

Let us demonstrate the efficiency of our design approach on a concrete problem coming from the industrial sector of the region of Nantes (France). One of the most important activity areas of this region is the manufacturing of bathroom components (shower cabin, washbasin, bathtub, etc.). Most of parts used during the assembly process are made of thermosetting materials.

The main operations achieved on these parts is trimming, i.e. the suppression of the edges of the parts in order to obtain a good surface roughness. The tools used for milling are specific mills, which are composed of a large number a diamond glued on its surface. Therefore the number of tooth is supposed to be infinite. However, the use of these specific tools allows the simplification of the model for computing the milling forces.

Let us consider the trimming operation depicted at Figure 2. The milling force may be decomposed into two components, denoted as Ft (the tangential force) and Ff (the radial force) (the vertical component can be neglected for such kind of operation

Architecture of the manipulator

For the milling process let us specialize the Orthoglide manipulator (Figure 4) [START_REF] Chablat | Architecture Optimization of a 3-DOF Parallel Mechanism for Machining Applications, the Orthoglide[END_REF]. This architecture was built in Institut de Recherche en Communications et Cybernetique de Nantes (IRCCyN) and satisfies the following design objectives: cubic Cartesian workspace of size 200200200 mm 3 (while for selection treatment required workspace 200200200 mm 3 ), Cartesian velocity and acceleration in the isotropic point of 1.2 m/s and 14 m/s 2 ; payload of 4 kg; transmission factor range 0.5-2.0. The legs nominal geometry is defined by the following parameters: L = 310 mm, d = 100 mm, r = 31 mm where L, d are the parallelogram length and width, and r is the distance between the kinematic parallelogram and the tool centre point. Stiffness model (elements and whole model) of Orthoglide manipulator is presented in [START_REF] Pashkevich | Accuracy Improvement for Stiffness Modeling of Parallel Manipulators[END_REF].

Performance evaluation

Since the workspace of the Orthoglide manipulator is lower than required for milling, the length and crosssection of the bar element was increased. Using proposed technique for the stiffness model identification, the new stiffness matrix for the bar element was obtained (Table 1). Now let us compute the displacements caused by the forces during milling of the composite material. The three subtasks are considered: estimation of the force size, estimation of the displacement through the whole workspace and analyzing of the force direction. The maximum compliance is 20% higher than for the single force (Fx, Fy) along the principal directions x and y of the frame. For the work point x=y=z=300, a force oriented along the principal axes of the frame causes deformations 20% less than for the worst direction and 50% more than for the strongest direction. It should be noted that in the workpoint (x=-200, y=300, z=0), forces along x and y directions cause different deformations. Along the x direction, the deformations are close to the minimum, while along the y direction, they are close to the maximum. The rate between the maximal and minimal deformations is about 7. Force size analyses. Results for the same four workpoints are presented on Figure 6. The forcedeflection relationship for the force less than 1000 N is linear, but depends on the manipulator configuration. Another conclusion is that in the work point x=y=z=0, the required accuracy can be satisfied for the force up to 600 N, while for all other tested workpoints, it can be satisfied only for forces inferior to 100 N. Moreover, taking into account the force direction, maximum compliance errors for the force 300 N may raise up to 1 mm and more. 7a and z = 300 mm Figure 7c) and for the "zero plane" (z = 0 mm Figure 7b) of workspace are presented. It should be noted that the position accuracy depends on the configuration of the manipulator and vary from 0,1 mm to more than 0.5 mm for z=-200 and z=300. The accuracy of the "zero plane" is satisfied for the operation, while for z=-200 mm, guaranteed accuracy is only 0.4 mm, and for z=300, 0.5 mm. But milling in the strongest direction of the manipulator leads to an accuracy more than 0,2 mm. So, optimizing the milling process for the specialized Orthoglide manipulator may improve the accuracy more than 3 times.

SUMMARY

The accuracy of milling of composite materials depends on the number of factors such as accurate kinematic and stiffness modeling, performance evaluation, force control, planning of milling process and others. This paper contributes to the methodology, which is used for the accurate stiffness modeling for the manipulator and it links for the estimation of the deflections errors. It allows evaluating compliance errors caused by the technological process and detecting strongest and lowest directions for the compliance errors. Using this information for the planning of the milling process allows increasing the accuracy of processing of the pieces made of composite materials. The method is efficiently illustrated on the milling of bathroom components produced in the Nantes region. For the manufacturing, Orthoglide manipulator was revised, and deformations for different forces and work points were presented. The results allow estimating the accuracy of technological process and improving it by the accurate planning of the milling. While analyzing the modeling results, the several directions for prospective research activities were identified. They include accurate modeling of milling, improving stiffness of the manipulator for the working direction and compliance deformation compensation in the milling process.

  (a) mapping of the unit sphere (b) mapping of the unit cube. (c) mapping of the velocity transmission in the xy plane.

Figure 1 :

 1 Figure 1: Mapping, using the Jacobian matrix.

Figure 2 :

 2 Figure 2: Trimming operation for composite material. The machines tools that are used for the trimming of the bathroom components proposed on the Figure 3 must be

Force

  direction analysis. Here let us plot errors for the constant force (300 N) rotating it from -180 to 180 deg for four typical work points: original configuration (x=y=z=0), two opposite corners of workspace x=y=z=-200 and x=y=z=300 and point for nonsymmetrical configuration of the manipulator (x=-200, y=300, z=0). Results are presented on Figure 5. Figures 5a-d show the size of deflections for each direction of the force and represent it in polar coordinates, while Figures 5e-h show ellipses of deflections in the Cartesian coordinates. As we can see the compliance of the manipulator depends on the direction of external force and on the workpoint. Only for the isotropic work point (x=y=z=0) it is constant for each direction. The lowest and the highest compliant directions of the manipulator differ on /2. For the work point x=y=z=-200, the manipulator is 7 times stronger in the direction 135/-45 deg than in the direction 45/-135 deg.

Figure 3 :

 3 Figure 3: Typical examples of bathroom components manufactured in the region of Nantes (France).

Figure 4

 4 Figure 4: CAD model of Orthoglide manipulator

Figure 5 :

 5 Figure 5: Norm of the deformations of the end-effector as a function of the direction  phi of the planar force F= 300 N: (a)-(d) in the polar coordinates (

Figure 6 :Figure 7 :

 67 Figure 6: Force-Deflections relationship in the four test points

Table 1 :

 1

Stiffness of bar element

Workspace analysis. For the accuracy control through the whole workspace, error maps for opposite planes (z = -200 mm Figure

ACKNOWLEDGMENTS

The work presented in this paper was partially funded by the Region "Pays de la Loire", France and by the EU commission (project NEXT).


Pierre-Emmanuel Morant

Quentin Thommen

Benjamin Pfeuty

Constant Vandermoere

Florence Corellou

François-Yves Bouget

Marc Lefranc
email: marc.lefranc@univ-lille1.fr

A robust two-gene oscillator at the core of Ostreococcus tauri circadian clock

for the presence of a robust two-gene oscillator at the core of Ostreococcus clock by showing that a minimal lightindependent model can reproduce very accurately microarray data and luminescent reporter data recorded in different experiments.

The microscopic green alga Ostreococcus tauri is rapidly emerging as a promising model organism in the green lineage. In particular, recent results by Corellou et al. [Plant Cell 21, 3436 (2009)] and Thommen et al. [PLoS Comput. Biol. 6, e1000990 (2010)] strongly suggest that its circadian clock is a simplified version of Arabidopsis thaliana clock, and that it is architectured so as to be robust to natural daylight fluctuations. In this work, we analyze time series data from luminescent reporters for the two central clock genes TOC1 and CCA1 and correlate them with microarray data previously analyzed. Our mathematical analysis strongly supports both the existence of a simple two-gene oscillator at the core of Ostreococcus tauri clock and the fact that its dynamics is not affected by light in normal entrainment conditions, a signature of its robustness.

In order to anticipate periodic environmental changes induced by Earth rotation, many organisms have evolved a circadian clock, a genetic oscillator which generates biochemical rhythms with a period around 24 hours. Exact synchronization with the day/night cycle requires that one or more clock components sense daylight (for example, a protein degrades faster in the light). However, daylight intensity can highly fluctuate from day to day, or during the day, due to environmental factors such as sky cover. The question then arises as to how circadian clocks can keep time without being continously reset by the signal that should entrain them? A mathematical analysis of Ostreococcus tauri clock, whose molecular basis has been identified recently [1], has unveiled a simple and elegant mechanism which exploits the dynamical properties of the core clock oscillator to make it robust to daylight fluctuations [2]. When the clock is on time, coupling to light is activated precisely when the oscillator does not respond to external perturbations, making it blind to light and its fluctuations. If the clock has drifted and needs resetting, however, light affects the oscillator in a different part of its cycle, where it reacts so as to recover the entrainment phase. In this work, we provide strong evidence

I. INTRODUCTION

Biochemical oscillations are widespread biological phenomena involved in many important cellular processes such as signaling, development, motility or metabolism [START_REF] Goldbeter | Biochemical Oscillations and Cellular Rhythms[END_REF][START_REF] Hasty | [END_REF]. Understanding how such a simple dynamical pattern has been implemented repeatedly to support a great diversity of biological functions is appealing for scientists seeking to uncover unifying principles [START_REF] Winfree | The Geometry of Biological Time[END_REF]. By identifying the molecular machinery behind cellular rhythms, molecular biology has first laid the groundwork for the development of strategies toward their comprehensive understanding [START_REF]Molecular Genetics of Biological Rhythms[END_REF]. In a second stage, how a collective dynamics can emerge in networks of interacting molecular actors and how it can be harnessed robustly has been under focus. Besides the design of synthetic genetic circuits with specific oscillatory abilities [START_REF] Elowitz | [END_REF]8], a common strategy to gain insight into this question has been to construct quantitative dynamical models constrained by experiments [START_REF] Goldbeter | Biochemical Oscillations and Cellular Rhythms[END_REF][9][10][11]. Indeed, the nonlinear dynamical behaviors that underlie oscillations can only be fully captured through a mathematical description.

However, the ever-growing amount of experimental data obtained with genetic transformations and real-time monitoring presents extraordinary challenges to modelers. How to design relevant mathematical models that not only reproduce the data but give insight into the architecture of a biological system, and avoid pitfalls such as overfitting which can ascribe biological meaning to ex-perimental artifacts? As we hope to illustrate here, it is important to combine careful data analysis and minimal modeling, and to check the consistence of the different data sets at all stages of model building.

Circadian clocks are systems of choice for quantitative studies of the function and design of biochemical oscillators [12]. They operate in many species to keep track of the most regular environmental constraint, the alternation of daylight and darkness caused by Earth rotation, so as to finely control the cellular physiology accordingly [13]. This basic function is vital in many organisms such as plants, which need to timely coordinate their photosynthesis, and thus their growth and division, to daily changes in light intensity [14][15][16]. However, a precise entrainment of circadian rhythms to the diurnal cycle is potentially challenged by many sources of variability including molecular fluctuations [START_REF] Gonze | Proc. Nat. Acad. Sci. USA[END_REF][START_REF] Barkai | [END_REF] and temperature variations [START_REF] Pittendrigh | Proc. Natl. Acad. Sci[END_REF][START_REF] Rensing | [END_REF], but also fluctuations of daylight intensity from day to day or during the day, due to environmental factors such as sky cover [2,[21][22][23]. To account for the remarkable ability of circadian oscillators to run autonomously and to be precisely and robustly entrained, experimental efforts aimed at unraveling their complex architecture [24][25][26] have motivated studies trying to adjust mathematical models to experimental data [2,[START_REF] Forger | Proc. Nat. Acad. Sci. USA[END_REF][START_REF] Forger | [END_REF][29][30][31][32][33], an approach that is increasingly necessary as models become more complex, featuring generally several feedback loops.

In this paper, we use such a quantitative modeling approach to investigate the dynamics of the circadian rhythms of the smallest free-living eukaryote known to date, Ostreococcus tauri. This microscopic green alga displays a very simple cellular organization and a small and compact genome [34][35][START_REF] Derelle | Proc. Nat. Acad. Sci. USA[END_REF]. Very recently, the molecular basis of its circadian clock has been extensively characterized by Corellou et al., who carried out an extensive work of genetic transformation, leading to transcriptional and translational fusion lines allowing one to monitor transcriptional activity and protein dynamics in living cells [1,[START_REF] Djouani-Tahri | [END_REF]. Their results point to a core archictecture comprising two genes, similar to Arabidopsis central clock genes TOC1 and CCA1 [38]. These two genes display rhythmic expression both under light/dark alternation and in constant light conditions. Thus, the unicellular green alga Ostreococcus tauri has emerged as a promising organism model to study the circadian rhythm in single photosynthetic eukaryotic cell combining experimental and modeling approaches.

The goal of this modeling study is to check carefully whether experimental data obtained through various channels support the hypothesis that the circadian clock of Ostreococcus Tauri contains a simple two-gene transcriptional loop serving as a core oscillator [1,2], and that the dynamics of this oscillator is not affected by light in normal entrainment conditions [2]. An important point of the present work is the simultaneous adjustment of microarray data and luminescent reporter data recorded in two different experiments [1,15] under different conditions. While this may not be sensible for a general system whose parameters will typically change with experimental conditions, it is expected here that the core clock oscillator is robust to environmental changes. As advocated by Thommen et al. [2], it should thus deliver similar signals in different experiments carried out with the same photoperiod. Therefore, being able to reproduce with the same model signals recorded with different techniques in different experiments provides strong evidence that the core Ostreococcus clock oscillator is indeed robust.

Before adjusting a minimal model to experimental data, a careful data analysis of both microarray and luminescence time series has been performed. It allowed us to detect and remove an experimental artifact that would otherwise have prevented an optimal adjustment and thus would have falsely called for a more complex model. The combined use of microarray data and luminescence reporter data provided a unique opportunity to calibrate the latter with respect to the former, given that very little is known on luciferase dynamics in Ostreococcus. This is all the more important as an extensive collection of experimental data has been acquired with these reporters and will form the basis of future quantitative studies of Ostreococcus clock.

We not only found that a simple two-gene transcriptional feedback loop model reproduces perfectly the CCA1 and TOC1 transcript and protein profiles obtained from data analysis, but that it does so with no model parameter depending on light intensity whereas the data have been recorded under light/dark alternation. This confirms our previous observation based only on microarray data of limited time resolution [2]. As we have proposed previously, this counterintuitive finding indicates that the coupling to light does not influence the oscillator in normal entrainment conditions, shielding it from daylight fluctuations [2]. Besides confirming the two-gene loop hypothesis and the oscillator robustness, this work also identifies unambiguously the mechanistic origin of oscillation in the transcriptional negative feedback loop involving CCA1 and TOC1, delayed by the saturated degradation of CCA1 mRNA and TOC1 protein.

II. DATA ANALYSIS

In this section, we first present the experimental data and explain how we determine mRNA and protein concentrations from them. We then discuss the simplest model of the TOC1 -CCA1 negative feedback loop, which consists of four differential equations describing the time evolution of the TOC1 and CCA1 mRNA and protein concentrations. To match the numerical time profiles generated by this model, we need to reconstruct target profiles from the experimental data. RNA profiles will be interpolated from microarray data whereas protein profiles will be extracted from luminescence time series in the form of a Fourier series describing the dynamics at long time scales. Both types of profile will be characterized by their times of passage at 20% and 80% of the oscillation amplitude, which provide a uniform characterization of the two types of data considered here. Model adjustment will be carried out by minimizing discrepancies between the experimental and numerical passage times.

A. Experimental data

The microarray data used here come from the study in Ref. [15]. They were recorded at three-hour time intervals under 12:12 light/dark alternation. The mRNA time profiles that will serve as targets for model adjustment are the same as in our previous study [2], where we showed that they could be accurately approached by solutions of a simple set of differential equations, an indication of the high quality of these data. The fact that microarray data reflect mRNA level without ambiguity was also checked by quantitative RT-PCR.

The luminescent reporter data used here are those presented by Corellou et al. in [1] to provide evidence of the central role of the TOC1 and CCA1 genes in Ostreococcus clock. They had been recorded at one-hour intervals under a 12:12 light/dark cycle from two types of transgenic cell lines.

In transcriptional fusion lines pTOC1:luc and pCCA1:luc, the sequence inserted into the genome consists of the promoter of one of the two clock genes fused to the coding sequence of firefly luciferase [1], so that TOC1 or CCA1 transcriptional activity drives luciferase expression. Luciferase catalyzes the transformation of the substrate luciferin into oxyluciferin, in which one photon is emitted [39]. The luminescence signal provides information about the quantity of luciferase synthesized and thus on the transcriptional activity of the promoters.

In translational fusion lines TOC1:luc and CCA1:luc, both the coding sequences of one of the clock proteins (TOC1 or CCA1) and of luciferase are fused to the promoter [1], so that a fusion protein combining the clock protein with luciferase is synthesized by the additional gene with a trancriptional activity similar to the original clock gene. In this case, the luminescence signal also provides information about the protein dynamics, and in particular on its degradation kinetics. It is usually assumed that luciferase remains complexed and inactive after the reaction and does not recover its activity before being degraded. Thus two limiting cases can be considered depending on whether the luciferase reaction time (defined as the average time from synthesis to photon emission) is much shorter or much longer that the protein lifetime.

In the former case, luciferase reacts immediately after being synthesized. Since only freshly made proteins then contribute to the luminescence signal, the latter is obviously proportional to the protein synthesis rate, hence to cytoplasmic fusion mRNA concentration. If we as-sume for simplicity that the kinetic constants (synthesis and degradation rates) of the fusion proteins TOC:luc and CCA1:luc and of their mRNA are similar to that of their native counterparts, then the concentrations of the former are proportional to the concentrations of the latter. In this case, the luminescence signal tracks mRNA concentration.

In the opposite case where the probability of photon emission by luciferase is very low, and with the same assumption of identical constants for fusion and native molecular actors, it is easily seen that the luminescence signal is proportional to protein concentration. In this limit, indeed, the photon emission probability is constant throughout protein lifetime, thus luminescence intensity is proportional to the number of TOC1:luc or CCA1:luc fusion proteins, which is in turn proportional to the number of native clock proteins. In this case the luminescence signal can be used as an indicator of protein concentration. In the general case, it should be intermediate between the RNA and protein time courses.

A time lag of at least two hours between the maxima of RNA concentration as indicated by microarray data and the maxima of translational fusion lines can be observed, which indicates that photon emission probability is low and that luciferase lifetime is large compared to other protein lifetimes. This is consistent with recent experiments in which translation was blocked in Ostreococcus cell cultures with emetine dihydrochloride [START_REF] Djouani-Tahri | [END_REF] and time evolution of TOC:luc luminescence was monitored. It was observed that the signal would slowly decay over more than 12 hours, implying that the luciferase reaction time is at least of this order. Thus, we will use in the following translational fusion line signals as indicators of protein concentrations, which allows us to keep our mathematical model as simple as possible.

B. Model of the TOC1 -CCA1 oscillator used for data adjustment

The minimal model for the TOC1 -CCA1 transcriptional feedback loop, where TOC1 activates CCA1 and CCA1 represses TOC1, consists of the following four ordinary differential equations:

ṀT = µ T + λ T 1 + (P C /P C0 ) nC -δ MT K MT M T K MT + M T (1a) ṖT = β T M T -δ PT K PT P T K PT + P T (1b) ṀC = µ C + λ C (P T /P T 0 ) nT 1 + (P T /P T 0 ) nT -δ MC K MC M C K MC + M C (1c) ṖC = β C M C -δ PC K PC P C K PC + P C (1d)
Eqs. (1) describe the time evolution of mRNA concentrations M C and M T and protein concentrations P C and P T for the CCA1 and TOC1 genes, respectively, as they result from mRNA synthesis regulated by the other protein, translation and enzymatic degradation. TOC1 transcription rate varies between µ T at infinite CCA1 concentration and µ T + λ T at zero CCA1 concentration according to the usual gene regulation function with threshold P C0 and cooperativity n C . Similarly, CCA1 transcription rate is µ C (resp., µ C + λ C ) at zero (resp., infinite) TOC1 concentration, with threshold P T 0 and cooperativity n T . TOC1 and CCA1 translation rates are β T and β C , respectively. For each species X, the Michaelis-Menten degradation term is written so that δ X is the low-concentration degradation rate and K X is the saturation threshold.

A more detailed model would take into account compartmentalization, with each actor having separate nuclear and cytoplasmic concentrations, as well as the fact that the luminescence signal is linked to a third gene artificially inserted in the genome, with the kinetic constants of its mRNA and protein possibly different from those of its native counterpart. With this model, it is the profile predicted for the total fusion protein concentration that would be adjusted to the experimental time series rather than the native protein profile. However, the important point is that such a biochemically detailed model taking into account all native and inserted molecular actors reduces to Eqs. (1) in the limiting case where fusion proteins and mRNA have the same kinetic constants as the native molecules (e.g., TOC1:luc and TOC1 degradation rates are equal), luciferase reaction time is large and nucleocytoplasmic transport is fast. When the minimal model already adjusts the data with excellent accuracy, as will be the case here, there is no point in using the sophisticated model, which can only fit the data better. In fact, its higher flexibility could lead to overfit the data and ascribe incorrectly biological meaning to experimental artifacts.

As in [2], the free-running period (FRP) of the uncoupled oscillator is chosen equal to 24 hours, which was the mean value observed in experiments [1]. In this case, this oscillator is adjusted to experimental data without modulation since our goal in this work is to show that the average oscillation does not carry any signature of coupling to light. As a control, we will also consider adjustment to oscillators of FRP 23.8 and 25 hours under day conditions. In these two cases, a small modulation of some parameters is required to achieve frequency locking. To keep the FRP fixed during adjustment, we rescale kinetic parameters after each optimization step so that the period of the free-running limit cycle matches the desired value. This subsequently applies to the parameter set which the adjustment converges to.

C. Reconstruction of target profiles and determination of passage times

Experimental RNA target profiles were estimated from the same microarray data points as used in [2]. Due to the high quality of data, noise reduction was not required. However, the relatively low temporal resolution (one sample every three hours) made it difficult to estimate accurately passage times and positions of expression peaks. Thus the profiles were obtained by interpolating between data points. Time series in the logarithm of RNA concentration showed a very smooth behavior and were well approximated by simple cubic splines (Fig. 1). These interpolating curves were thus considered to provide an optimal approximation to the variations of the logarithm of mRNA concentration, and mRNA linear target curves were generated by exponentiation of the interpolating curves.

That the largest data point is less than half of the maximum of the reconstructed TOC1 mRNA profile should not be surprising. In fact, this profile is the most probable one given the data points if we assume sufficient smoothness, or equivalently that the same dynamics acts over the 24-hour cycle. Any other curve would imply fast, transient, processes not linked to the TOC1 -CCA1 loop. Since RNA concentration rises from almost zero at ZT6 (Zeitgeber Time 6, i.e., 6 hours after dawn) to almost the maximum level measured at ZT9 and decays from the same level at ZT12 to almost zero at ZT15, it is indeed obvious that the true maximum located between ZT9 and ZT12 must be much higher than the two largest data points. Incorrectly assuming that the largest data point is the profile maximum would in fact bias the analysis.

We checked that interpolating either the CCA1 mRNA concentrations or their logarithms lead to two almost superimposing profile curves, showing that both sequences of values provide the same information. However the profile interpolated directly from TOC1 mRNA concentrations rather than from their logarithms was highly nonphysical, with interpolated concentrations becoming clearly negative on each side of the expression peak. This is a direct consequence of the fact that only two data points are well above the zero level. Implicit in the interpolation procedure is the reconstruction of successive time derivatives of a function from its values at different points. The long sequence of almost zero data points makes the correspondence between function values and time derivatives nearly singular, so that the reconstruction is highly unstable. During the time interval near zero, all time derivatives up to high order appear to be zero, then jump to high values at the beginning and the end of the expression peak. Since the interpolation error is proportional to the maximum value of the fourth time derivate in the interval, the resulting approximation is poor.

The series of logarithms of TOC1 mRNA concentrations does not have this problem and allows optimal re- construction of the TOC1 mRNA profile. The fact that the logarithm of a time series with long intervals near zero provides more information about the dynamics than the original time series has also been noted in the analysis of chaotic time series [40,41]. Furthermore, it should be noted that many statistical analyses of microarray data use the logarithms of mRNA concentrations because they are more evenly distributed and provide more information.

Compared to a more usual method based on leastsquare-fitting the data points (as used in [2]), this interpolation procedure can only make adjustment more difficult, since it constrains tightly the shape of the profile. However, it is fully consistent with the good agreement found with a low-dimensional ODE set in our previous study [2]. The passage times were then determined as the times at which the interpolated profile would cross the appropriate level (Fig. 1).

Protein profiles were reconstructed from translational fusion luminescence signals used as indicators of protein concentrations, as discussed above. The reconstruction was more involved than for mRNA because these time series display significant amplitude variations from peak to peak across the experiment, which may be for instance linked to variations in the number of cells contributing to light emission or other experimental factors.

The first 48 recording hours, considered as transient, were discarded and traces were then least-squares fitted by the product p(t)F (t) of a Fourier series F (t) = N k=0 a k cos(2πkt/T + φ k ), with T = 24 hours, by a loworder polynomial p(t) accounting for slowly varying experimental conditions (Fig. 2). The order of the polynomial p(t) was generally chosen to be 3 or 4, roughly equal to the number of periods in the fitted segment to avoid over-fitting. The Fourier series F (t) represents the average periodic biochemical oscillation associated with circadian oscillations while p(t) models a slowly varying gain (due for example to a variable number of cells contributing to light emission).

Because we are interested in long lasting mechanisms sustaining the oscillation, the number N of harmonics of the Fourier series was fixed to 5 (an odd number, since the spectrum of the square wave forcing by light only contains odd harmonics). We also tried slightly higher numbers of harmonics, with little difference in the results. Fourier fitting naturally smoothes out acute responses at day/night and night/day transitions that can be seen in the raw signals and which may reflect rapid mechanisms involved in clock resetting. This separation of dynamical processes occurring on a 24-hour time scale on the one hand, and confined to short time intervals on the other hand, must be understood with reference to our previous finding that the average dynamics of the TOC1 -CCA1 oscillator is identical to that of a free-running oscillator, with coupling to light and resetting occuring transiently in specific time intervals [2]. Confirming this key observation with a careful analysis of the luminescence time series is one of the main goals of this paper. It will be reached if we find that the average protein profile does not carry any signature of a coupling to light.

The Fourier series obtained from the fits of all bioluminescence traces (possibly corresponding to different transgenic lines) recorded in the experiment were then averaged, and the average was used as surrogate for the protein temporal profile after suitable normalization. The resulting profile is shown in Fig. 2(b), along with individual luminescence traces from different wells, renormalized using the slowly varying polynomial function p(t). The variability observed is partially due to the fact that a few different transgenic lines, corresponding to different insertions in the genome, were used. The times of passage at 20% and 80% of the amplitude were then determined, as well as the minimal level reached, for subsequent model adjustment.

We also tried fitting each luminescence trace to p(t)F (t) + b where b represents a possible constant bias. Surprisingly, the fit was consistently better when b was equal to the minimum luminescence level. In other words, the Fourier fitting procedure tended to remove the floor level. Because this could indicate the existence of a bias in the luminescence level, we examined more closely the raw data, and found that indeed the zero luminescence Large dots indicate passages through levels corresponding to 20% and 80% of maximum amplitude and serve as target points for adjustment.

level does not correspond to the zero protein level. In Fig. 3, we show two individual traces monitoring luminescence intensities over time in two wells of the luminometer plate. These two wells contain genetically identical cell cultures monitored simultaneously in the same experiment. Besides an excellent overall reproducibility, the comparison of the two traces reveals the existence of a significant bias. Assuming that identical clocks run in different cells, and that there is a well-defined average number of fusion proteins per cell, the luminescence signals should be proportional to numbers of cells in each well. Therefore, the two signals should be approximately proportional to each other.

Contrary to this, we observe that both for TOC1:luc and CCA1:luc, the maxima of the two signals differ by a large factor, which remains roughly constant over time while minima almost coincide. If the two traces were proportional to each other, the minima should be in the same ratio as the maxima. The only simple explanation for the systematic coincidence of minima is that they correspond to a zero protein level, and that there is the same bias on the two time series. This bias can therefore be removed by substracting the two traces, as is shown in Fig. 3, where the difference curves provide a better estimation of the actual protein profiles than the original traces. Importantly, the CCA1 protein level is predicted to touch zero very shortly near ZT9, while the TOC1 protein level appears to stay near zero for a large interval of time between ZT21 and ZT8 (Fig. 3). Note that this bias is not constant in time but varies slowly so that the different minima correspond to different luminescence levels.

FIG. 3: Translational fusion lines data recorded under 12:12 LD alternation. Time zero corresponds to dawn. On top panel (resp., bottom panel) the time evolution of photon count (in relative luminescence units) of CCA1 (resp., TOC1) translational fusion lines are drawn for two biological duplicates monitored at the same time in the same conditions (black and grey thin solid lines). Their difference is plotted as a thick black solid line.

While the existence of an experimental bias in the lu-minescence level is strongly supported by the observations above, we have currently no explanation for it. In the following, we will therefore adjust two versions of the experimental profiles : (i) the Fourier series F (t) determined directly from fitting p(t)F (t) to the time series, as described above, if there was no bias, and (ii) the same Fourier series with the floor level removed so that they touch zero, as with the difference curves in Fig. 3. More precisely, the target curves are then F (t) -F 0 where the F (t) are obtained from fitting p(t)F (t) and F 0 = min t F (t). This will allow us to assess the influence of the bias on model fitting. Again, using a more sophisticated model of the floor level bias could only lead to a better adjustment, which will not be needed, and would only make sense with a better understanding of the physical origin of the bias.

D. Measuring goodness of fit by passage times

Adjusting experimental data to a mathematical model requires a score function quantifying the discrepancy between experimental and simulated profiles, to be minimized over parameter space. One difficulty is that microarray and luminescent reporter data are very different in nature and in particular have different uncertainties so that it does make little sense to combine their adjustment errors. Therefore, we used passage times as a unifying measure of fitting errors, determining for each temporal profile the times of passage at 20% and 80% of the interval between the minimal and the maximal values. We then compared passage times for the experimental and numerical profiles, with the root mean square timing error serving as a goodness of fit. Another advantage of passage times is that they have a clear biological significance, with crossings of the 20% level bracketing the time interval with expression significantly above background level and the crossings with the 80% level bracketing the expression peak.

We denote by ∆ X 20↑ and ∆ X 20↓ (resp., ∆ X 80↑ and ∆ X 80↓ ) the time passage error in minutes for the concentration of species X (X = M C , P C , M T , P T ) at 20% (resp., 80%) of the interval between minimal and maximal values when increasing and decreasing. For each species X, one can therefore define a quadratic error function equal to:

Err(X) = ∆ X 20↑ 2 + ∆ X 80↑ 2 + ∆ X 80↓ 2 + ∆ X 20↓ 2 (2)
The first, and most natural, score function is a RMS error combining error for all species (mRNA and proteins):

S MP = 1 16 X∈MC ,PC ,MT ,PT Err(X) (3)
In some cases, profiles with similar passage times differ by their floor level. We therefore introduce a floor level error Φ X = min( X) max( X) -min(X * ) max(X * ) where X is the numerical profile for species X and X * is the experimental target for X. We choose a ponderation such that a floor level error Φ X = 0.05 is equivalent to a passage time error of 60 minutes. The score function of this second scheme is

S MP F = 1 √ 18 16 (S MP ) 2 + X∈PC ,PT (Φ X /Φ 0 ) 2 (4)
where Φ 0 = 0.05/60 = 1/1200 is the relative floor level error considered equivalent to one minute. Finally, microarray and luminescence data were not only adjusted simultaneously but also separately in order to assess the coherence between the two types of data. To this end, we use the following score functions:

S M = 1 8 X∈MC ,MT Err(X) (5)
S P = 1 10 X∈PC ,PT Err(X) + (Φ X /2Φ 0 ) 2 (6)
where the Φ X error term is scaled to 2Φ 0 so that it has the same relative importance in the total RMS error as in Eq. ( 4).

III. RESULTS: ADJUSTMENT OF MODEL TO RNA AND PROTEIN DATA

Figure 4 shows the numerical solutions of model ( 1) best adjusting the experimental data according to the score functions defined in Sec. II D, both without taking into account the floor level bias (Figs. 4(a)-(d)) or with floor levels removed (Figs. [START_REF] Hasty | [END_REF] I.

(e)-(h)). The corresponding parameter values are given in Table

A first remark is that although a set of four passage times is a very crude description of a temporal profile, it turns out that numerical solutions and experimental profiles that have close passage times generally follow each other closely all over the day. This is consistent with the hypothesis that experimental data are well described by a low-dimensional ODE set. Secondly, the quality of the adjustment is globally very good given the simplicity of the model. Clearly, Fig. 4 shows that Eqs. (1) can adjust simultaneously the protein profiles reconstructed from luminescence time series and the RNA profiles reconstructed from microarray data (score functions S MP and S MP F ). Especially, removing the protein profile floor levels so as to correct for the detected experimental bias leads to impressive agreement between numerical and target curves, which almost superimpose onto each other.

It is interesting to look more closely at the main discrepancies in Figs 4(a)-(d). When only passage times are taken into account (score function S MP ), the nonzero floor level of the TOC1 protein is not reproduced. This issue can be addressed by using a score function that takes into account floor level error (score function S MP F ). In this case, the floor level agreement is improved, at the cost of slightly advancing the protein peak and of inducing a small floor level for TOC1 mRNA. Although small, the latter is definitely much higher than the value that can be accurately determined from microarray data (which estimate the logarithm of mRNA concentration). Last, the numerical solution obtained by adjusting the protein profiles alone (score function S P ) predicts mRNA profiles rather poorly (Fig. 4), which could suggest that the two types of data are inconsistent.

To an uninformed mind, these minor discrepancies would most probably suggest shortcomings of model (1), not surprisingly given its caricatural simplicity, or the fact that microarray and luminescence data have been recorded in different experiments. However, these dis-crepancies are clearly linked to errors in floor level, which are naturally explained by the protein floor level bias discussed in Sec. II C and evidenced in Fig. 3.

Let us now consider more closely adjustment of model ( 1) to protein target profiles with floor levels removed, together with RNA target profiles, shown in Figs. 4(e)-(h). Goodness of fit is clearly improved compared to Figs 4(a)-(d), where it was already good. The simultaneous adjustment of the four target curves (score function S MP ) is achieved with an accuracy rarely obtained for a genetic circuit, with all numerical curves superimposing perfectly on the experimental ones, except that the numerical CCA1 protein profile rises slightly more slowly than the experimental one between dusk and CCA1 expression peak. It is quite striking to note that . the maximum value of protein profiles is 100 nM, the maximum value of CCA1 (resp., TOC1 ) mRNA profile is 10 nM (resp., 70 nM). The TOC1 and CCA1 mRNA maximum values are chosen in the same proportion as in microarray data. The third row of the table indicates whether the floor levels of luminescence data are removed (R) or not (NR). The last part of the table gives the degradation rate DX at the mean value X = (max(X) + min(X)) /2 ; Dx = δX KX / KX + X for each species. the two curves begin to separate precisely in the time interval where a transient window of CCA1 stabilization was predicted to occur in Ref. [2]. Remarkably, numerical profiles obtained by adjusting RNA data alone (score function S M ) reproduce all experimental data very well also. When protein data are adjusted separately (score function S P ), RNA profiles are much better reproduced than when the protein floor is not removed, with only the CCA1 mRNA profile showing noticeable discrepancies in certain parts of the day, but with the global shape of the profile preserved. Thus we can conclude that removing floor level significantly improves adjustment and restores coherence between microarray and luminescence data. Given the simplicity of the model, it is very unlikely that this may have occurred by chance. As we will discuss in the next section, it is moreover a remarkable finding that experimental data recorded with different techniques and in different experiments [1,15] can be adjusted simultaneously with such high accuracy by a minimal model. This certainly reveals an important property of the clock oscillator. The parameter sets in the last three columns of Table I, which correspond to profiles with protein floors removed, are generally more consistent between themselves than those in the first three columns, where the bias has not been corrected, and where important vari-ations can be seen (compare for example λ T across the first three columns). This reflects the fact that when the floor level is not removed, it is more difficult to predict mRNA profiles from adjusting protein profiles alone and vice versa. To make meaningful comparisons, it should however be kept in mind that the significant quantities are often not parameters themselves but combinations of them. When degradation is saturated, for example, the relevant quantity is the product δ X K X , and it may happen that δ X and K X fluctuate more between columns than their product. Therefore, we give at the bottom of Table I the effective degradation rates at a mean value of the concentration D X , which appear to be much more consistent than the individual degradation parameters. Another example is that although the values of µ t in the fourth and fifth columns are quite different, the minimum values taken by TOC1 transcription rate in both cases are in fact comparable (1.78 vs 1.45 nM/h), and result from different combinations of λ T , P C0 , and of the minimal value reached by the P C profile, which is very close to the repression threshold P C0 . This illustrates the general fact that although collective fitting can provide well-constrained predictions (here the time profiles), individual parameters may be poorly constrained, a feature observed in many systems biology models [42].

SMP

Some key ingredients of the dynamics can still be ex-tracted unambiguously from a careful examination of of K K PT so small that the concentrations of the respective molecular actors are well above them for almost all of the diurnal cycle. The signature of this behavior in the time profiles is a straight-line decay after the expression peak. A natural question is whether what appears as saturated degradation of a given actor may in fact result from an interaction with another actor not yet identified. In the case of TOC1, a post-transcriptional regulation acting after dusk has indeed been evidenced experimentally [START_REF] Djouani-Tahri | [END_REF].

Other key features are a maximum transcription rate significantly higher for TOC1 than for CCA1 and a small threshold of repression P C0 , comparable to the minimum value of the CCA1 protein level. These are related since the latter implies that TOC1 is repressed most of the time except in a very narrow time interval around the minimum of CCA1 protein level (this correlates with the narrow TOC1 mRNA peak). This must be compensated for by a high TOC1 transcription rate. Interestingly, the small value for P C0 can account for the experimentally observed relative inefficiency of antisense strategies against CCA1 [1]. If CCA1 level is reduced, even significantly, the small time interval where TOC1 is not repressed will extend only slightly, with the dynamical behavior at other times being mostly unchanged. Thus the global perturbation will remain limited and arrhythmia will be difficult to induce. In contrast to this, the activation threshold P T0 above which TOC1 activates CCA1 is relatively low so that activation occurs over a rather large time interval of more than 10 hours. Therefore, modifications of the activation threshold or TOC1 protein level will influence the dynamics for a long time and thus easily disrupt oscillations. This also explains why FRP is much more sensitive to overexpression of TOC1 induced by TOC1:luc insertion than of CCA1 induced by CCA1:luc [1].

Finally, we carried out data adjustment assuming freerunning periods of 23.8 and 25 hours, to verify that our results do not depend critically on FRP (Fig. 5, parameters given in Table II). In this case, we allowed day/night modulation of some parameters to ensure frequency locking with the diurnal cycle. The agreement remains very good, although it degrades noticeably compared to the case of a FRP of 24 hours. In particular, a phase shift of the TOC1 mRNA profile is induced. The parameter modulation remains very small (for a FRP of 25 hours, the repression threshold is modulated but remains below the minimum of the CCA1 temporal profile, so that the effect of the change is minimal). We furthermore observed that when the FRP was a freely adjustable parameter, it would systematically converge towards a value of 24 hours. All this confirms that there is no day/night parameter modulation in the TOC1 -CCA1 feedback loop and that coupling can only occur in specific time windows as proposed in [2].

IV. DISCUSSION AND CONCLUSION

In this work, we have carried out a careful and detailed analysis of time series characterizing temporal variations of expression of the two central genes of Ostreococcus circadian clock, TOC1 and CCA1 [1]. These time series had been previously obtained from microarray data and luminescent reporter data recorded in two different experiments [1,15]. From these time series, we have extracted periodic temporal profiles for the RNA and protein concentrations of the two genes, assumed to represent the mean circadian oscillations in individual cells so as to optimize comparison with numerical profiles, which are inherently periodic. In particular, approximating the periodic component of luminescence time series 2). The data processing also allowed us to detect a bias in the luminescence time series, which was confirmed by a direct comparison of individual signals from two genetically identical cell cultures (they were not proportional to each other). This bias manifests itself in a nonzero luminescence floor level. Taking it into account allowed us to show that both protein levels approach zero at some times of the day, which was essential for the quality of the subsequent model adjustment. This illustrates how important it is to exploit the redundancy of data by verifying that data that should provide the same information are indeed consistent. Thanks to the careful data processing, we could evidence an extraordinarily good agreement between a minimal model of a two-gene transcriptional feedback loop and the reconstructed concentration profiles. This model describes activation of CCA1 by TOC1 and repression of TOC1 by CCA1. It describes only regulated transcription of the two genes, translation and degradation, and comprises only differential equations, from which the time evolution of the two mRNA and the two proteins can be computed and compared to the target profiles. Therefore a biochemically detailed model taking into account compartmentalization or genomic insertions would have only served to adjust details without biological relevance, and could well have masked one of the main results which is the good adjustment by a free-running oscillator model. It is also quite remarkable that using only four points of the reconstructed profiles sufficed to constrain the numerical profiles to follow their targets throughout the day. Here also, adding more points would have only forced the model to adjust to unrelevant details without gaining information. This suggests that the complexity of the model and the constraining data should be carefully matched.

The excellent agreement found between the model and the data supports unambiguously the existence of a twogene oscillator at the core of Ostreococcus circadian clock. This not only builds a solid foundation on which future studies of this clock can rely but also provides what we believe is one of the clearest examples of a natural fewgene oscillator evidenced from experimental data. This is all the more important as the role of this oscillator in the circadian clock constrains it to be extremely robust to all kinds of fluctuations. Understanding the dynamical ingredients besides transcriptional regulation that underlie its robustness will certainly be of high interest for the study of genetic oscillators in general. Two remarkable features of the TOC1 -CCA1 oscillator have indeed emerged from our analysis. First, a strongly saturated degradation has been evidenced both for the CCA1 mRNA (already noted in our previous work [2]) and the TOC1 protein (detected in this work), which manifests itself by a straight-line decay at high concentrations after the expression peak. This behavior may be the signature of post-transcriptional and post-translation interactions, and thus compatible with the experimental observamechintroduce effective a negative feedback loop [10,[43][44][45] and the recent observation that this is a key ingredient to generate robust oscillations [8,46]. The second remarkable feature predicted by the model is that the TOC1 gene is repressed by CCA1 during most of the day, except during a short time interval located one or two hours before dusk, which is consistent with the very narrow peak of TOC1 mRNA expression observed in the experimental data. The small duration of TOC1 expression is compensated by a high transcription rate. One may wonder whether this design has an influence on the robustness to molecular fluctuations. In any case, it explains why functional studies of the oscillator showed that it was much more sensitive to pertubations in the TOC1 level than in the CCA1 level [1].

It has been recently been proposed that circadian clocks should not only be robust to fluctuations in molecule numbers [START_REF] Gonze | Proc. Nat. Acad. Sci. USA[END_REF][START_REF] Barkai | [END_REF] and to temperature variations [START_REF] Pittendrigh | Proc. Natl. Acad. Sci[END_REF][START_REF] Rensing | [END_REF] but also to fluctuations in the daylight intensity pattern, which is crucial to synchronize the circadian clock to the day/night cycle [2,23]. Our results demonstrate such robustness for the TOC1 -CCA1 oscillator in two different ways. First, the experimental data are reproduced accurately by a free-running oscillator model. This had already been noted by Thommen et al. [2], but the confirmation of this behavior using luminescence signals recorded every hour reinforces its plausibility significantly. As shown in [2], it can only be explained by assuming that coupling to light is scheduled so as to be active precisely when the oscillator does not respond to external perturbations. In this scheme, the oscillator is not affected by light in normal entrainment conditions, which naturally makes it blind to light fluctuations, and thus behaves as if it was free-running. When the oscillator drifts out of phase, light sensing occurs at a different time of its cycle, where it responds so as to recover the normal entrainment phase. The acute responses to light/dark and dark/light transitions that occur transiently in raw signals (Fig. 2) may be correspond to such resetting.

The second strong evidence for the robustness of the TOC1 -CCA1 oscillator comes from the fact that the profiles adjusted precisely and simultaneously by our simple model are reconstructed from time series recorded in different experiments by two different techniques, with totally different setups and in particular under very different lighting conditions. This strongly suggests that Ostreococcus clock is able to tick exactly in the same way in different conditions, and in particular in different levels of light, which is an obvious requisite for robustness to daylight fluctuations. If the clock can cope with randomly varying daylight profiles without being perturbed, it can certainly accomodate constantly low or constantly high daylight levels and deliver similar biochemical signals in both cases. The fact that the model adjusted from RNA profiles predicts protein profiles and vice versa (although with less precision in the latter case) clearly illustrates the consistency between the two experiments. Again, we stress that this remarkable finding may have remained masked without our careful data processing.

It is important to mention that both the presence of saturated degradation kinetics and the absence of effective coupling to light when entrained are strong predictions that can be tested experimentally.

While the impressive agreement between model and experiments obtained here clearly shows that a TOC1 -CCA1 oscillator underlies Ostreococcus clock, it should not make us forget that it is only a part of it. The freerunning average behavior in entrainment conditions, the acute responses at day/night and night/day transitions can only be explained if the TOC1 -CCA1 loop interacts with other molecular loops and feedback loops. Indeed it was shown in [2] that robustness to daylight fluctuations involves a precisely timed coupling activation. This necessarily requires additional feedback loops designed so as to generate the signal that will drive optimal coupling. Furthermore, interactions with other molecular actors may well be responsible for the saturated degradation of CCA1 mRNA and TOC1 protein.

We thus have now to identify these additional feedback loops as well as the light input pathways, which should help us to understand how the clock can synchronize in spite of fluctuations to day/night cycles of variable duration across the year. Adapting to different photoperiods and to different light levels are indeed unrelated evolutionary goals that must be simultaneously satisfied. In the end, we hope to build an accurate and comprehensive model of a simple and robust circadian clock. If the agreement between theory and experiment remains comparable to what was achieved in the present study, this may well provide us with new and deep insights into the design and function of circadian clocks. More generally, we believe our results promote Ostreococcus, whose low genomic redundancy [START_REF] Derelle | Proc. Nat. Acad. Sci. USA[END_REF] is probably crucial for allowing accurate quantitative approaches, as a very promising model for systems biology.

FIG. 1 :

 1 FIG.1: Construction of RNA target profiles as functions of Zeitgeber Time (ZT) describing phases within the dark/light cycle, with time ZT0 corresponds to dawn and time ZT12 to dusk. Interpolating curves going through data points (dots) are represented. (a) Interpolated curves from microarray data using cubic splines. (b) The exponentials of the interpolating curves were computed to obtain a smooth approximation of mRNA profiles. Large dots indicate passages through levels corresponding to 20% and 80% of maximum amplitude and serve as target points for adjustment.

FIG. 2 :

 2 FIG.2: Reconstruction of protein target profiles as functions of Zeitgeber Time (ZT) describing phases within the dark/light cycle, with time ZT0 corresponds to dawn and time ZT12 to dusk. (a) Luminescence time series (dashed lines) for individual wells are fitted by the product of a Fourier series of period 24 hours with 5 harmonics by a slowly varying polynomial function of degree 4 (solid lines). The slowly varying envelope is also shown with a dotted line. From top to bottom, the first (last) two panels correspond to two TOC1:luc (CCA1:luc) translational fusion lines with different insertions in the genome. The different Fourier series are normalized to have the same maximum and are averaged. (b) The average Fourier series (dashed lines, black: TOC1, grey: CCA1) provides a smooth approximation to the cloud of individual line time profiles renormalized using the slowly varying polynomial function and wrapped around 24 hours (thin solid lines). Large dots indicate passages through levels corresponding to 20% and 80% of maximum amplitude and serve as target points for adjustment.

FIG. 4 :

 4 FIG. 4: (Color figure) Adjustement of numerical solutions of the free-running model with FRP equal to 24 h (solid lines) to target curves (dash-dotted lines) using various score functions. Panels (a)-(d) in left column show adjustment of mRNA and protein profiles without bias correction, using score functions SMP (black solid line), SMP F (dark grey -blue thin solid line) and SP (light grey or red thin solid line). See text for the definition of score functions. Panels (e)-(h) in right column show adjustment of mRNA and protein profiles with protein floor level removed, using score functions SMP (black solid line), SP (light grey -red thin solid line) and SM (dark grey or blue thin solid line). (a), (e): CCA1 mRNA; (b), (f) CCA1 protein; (c), (g) TOC1 mRNA; (d), (h) TOC1 protein.

FIG. 5 :

 5 FIG. 5: (Color figure) Adjustement of numerical solutions of the light-dependent model (solid lines) to target curves with floor levels removed (dash-dotted line) using the SMP score function. Adjustment of CCA1 mRNA (a), CCA1 protein (b), TOC1 mRNA (c), TOC1 protein (d) without the floor level of proteins. The FRP of the autonomous oscillator is either of 25 h (black solid line) and 23.8 h (light grey or red solid line). Synchronization is obtained for the 25 h (resp., 23.8 h) FRP model by assuming that the parameter PC0 (resp., δP C ) take a different value on the day and night.

TABLE I :

 I Model parameter values. Optimal 24 hours

TABLE II :

 II Model parameter values. Optimal parameter values for adjustement of model to data using various score functions and assuming a free-running period of 25 or 23.8 hours. Parameters are rescaled as in TableI

	FRP	(h)	25	23.8
	score	(min)	27.8	34.3
	nC		1	1
	nT		2	2
	µC	(nM.h -1 ) 6.06 10 -4 4.73 10 -2
	λC	(nM.h -1 )	3.82	4.99
	PT O	(nM)	15.9	27.3
	βC	(h -1 )	2.71	2.34
	µT	(nM.h -1 )	22.7	6.07 10 -2
	λT	(nM.h -1 )	357	1130
	PC0	(nM)	5.94	6.15
	P night C0	(nM)	11.3
	βT	(h -1 )	0.805	0.809
	1/δM C (h)	0.305	0.296
	1/δP C (h)	2.28	2.50
	1/δ night P C	(h)		2.60
	1/δM T (h)	4.28 10 -2 4.51 10 -2
	1/δP T (h)	0.224	0.402
	KM C	(nM)	0.690	0.879
	KP C	(nM)	68.5	58.3
	KM T	(nM)	6.96	17.1
	KP T	(nM)	3.23	5.62

by a Fourier series allowed us to separate the dynamics at work throughout the day from fast transient processes activated only near light/dark and dark/light transitions, and probably involved in occasional resetting of the clock (compare raw data and target profile in Fig.

V. ACKNOWLEDGMENTS

This work has been supported by ANR grant 07BSYS004 to F-Y.B. and M.L., by CNRS interdisciplinary programme "Interface Physique, Biologie et Chimie: soutien à la prise de risque" to M.L., as well as by French Ministry of Higher Education and Research, Nord-Pas de Calais Regional Council and FEDER through Contrat de Projets Etat-Region (CPER) 2007-2013.