Guillaume Allibert

Dinuka Abeywardena

Moses Bangura

Robert Mahony

Estimating Body-Fixed Frame Velocity and Attitude from Inertial Measurements for a Quadrotor Vehicle

A key requirement for effective control of quadrotor vehicles is estimation of both attitude and linear velocity. Recent work has demonstrated that it is possible to measure horizontal velocities of a quadrotor vehicle from strap-down accelerometers along with a system model. In this paper we extend this to full body-fixed-frame velocity measurement by exploiting recent work in aerodynamic modeling of rotor performance and measurements of mechanical power supplied to the rotor hub. We use these measurements in a combined attitude and velocity nonlinear observer design to jointly estimate attitude and bodyfixed-frame linear velocity. Almost global asymptotic stability of the resulting system is demonstrated using Lyapunov analysis of the resulting error system. In the current work, we ignore bias and leave it for future work. The performance of the observer is verified by simulation results.

I. INTRODUCTION

Inertial measurement systems have a long history of being used for estimating the orientation of aerospace vehicles such as spacecrafts, aircrafts and missiles [START_REF] Crassidis | Survey of nonlinear attitude estimation methods[END_REF]. Most orientation estimators are based on the principle of observing known vectorial directions. For micro-aerial vehicles such as quadrotors, the two most commonly used vectorial measurements are the gravity and Earth's magnetic fields. In earlier work, Mahony et. al [START_REF] Mahony | Nonlinear complementary filters on the special orthogonal group[END_REF] showed that with just one vectorial measurement along with angular velocity can be used to estimate roll and pitch angles and gyroscope biases. This approach has been applied extensively in practice for attitude estimation in quadrotor vehicles using the accelerometers as a gravity estimate, even though the assumptions do not hold exactly in this case [START_REF] Martin | The true role of accelerometer feedback in quadrotor control[END_REF]. In particular, if the inertial accelerations of the vehicle are significant, then they need to be subtracted from the accelerometer measurement to obtain the true gravity vector measurement. Global Positioning Systems (GPS) can be employed for this purpose [START_REF] Hua | Attitude estimation for accelerated vehicles using gps/ins measurements[END_REF] but high rate absolute position and velocity measurements are a luxury that most mobile robots operating in urban and indoor environments lack. Recent investigations into the dynamics of quadrotor vehicles have revealed that the horizontal drag force affecting a quadrotor vehicle is proportional to the horizontal body frame velocity [START_REF] Bangura | Aerodynamic power control for multirotor aerial vehicles[END_REF], [START_REF] Bangura | Nonlinear dynamic modeling for high performance control of a quadrotor[END_REF], [START_REF] Mahony | Multirotor aerial vehicles: Modeling, estimation, and control of quadrotor[END_REF] and that this force can be effectively measured using a triad of body 1 Guillaume Allibert is with I3S-CNRS, Université de Nice Sophia Antipolis, France. Allibert@i3s.unice.fr 2

Dinuka

Abeywardena is with Centre for Autonomous Systems, University of Technology, Sydney, Australia.

Dinuka.Abeywardena@ieee.org 3 Moses Bangura and Robert Mahony are with Research School of Engineering and Information Sciences, Australian National University, Canberra, Australia. Moses.Bangura@anu.edu.au, Robert.Mahony@anu.edu.au mounted accelerometers [START_REF] Martin | The true role of accelerometer feedback in quadrotor control[END_REF]. As a natural extension of this, other authors have demonstrated that these accelerometer measurements can be combined with body mounted gyroscope measurements to obtain an estimate of the horizontal velocity of the quadrotor vehicle, along with its attitude [START_REF] Abeywardena | Improved state estimation in quadrotors mavs[END_REF], [START_REF] Leishman | Quadrotors and accelerometers: state estimation with an improved dynamic model[END_REF].

In this paper we consider the problem of designing a nonlinear observer for the full body-fixed-frame velocity along with the attitude of a quadrotor vehicle. The paper makes two principal contributions: Firstly, we show how aerodynamic modeling of the rotors can be used along with the accelerometer measurements and measurements of rotor speed and motor shaft torque available in some electronic speed controllers, to measure the full body-fixedframe velocity of a quadrotor vehicle. This measurement depends directly on accelerometer measurements and is always subject to high level of noise. To address this, we propose a nonlinear observer in the spirit of Hua [START_REF] Hua | Attitude estimation for accelerated vehicles using gps/ins measurements[END_REF], but posed in the body-fixed-frame. The second principal contribution of the paper is to use a matrix decomposition of the general R 3×3 matrix that appears in this observer design in order to identify the orthogonal part of the observer directly. This provides a direct estimation of the attitude of the vehicle without requiring a second stage orientation observer as was necessary in prior literature. The performance of the proposed observer is verified by simulation.

The paper has four sections in addition to the present introduction. In Section II we introduce the model and some mathematical definitions. In Section III we use momentum theory to show how the full body-fixed-frame velocity can be estimated from only an inertial measurement unit (IMU) and a suitable electronic speed controller (ESC). Section IV presents the proposed observer and provides a stability proof for its performance while Section V presents some simulation results to verify the performance of the proposed observer.

II. BACKGROUND

This section introduces the kinematic model underlying the observer problem considered and defines some Lie theory terminology that will be used in the sequel.

Consider a quadrotor vehicle as a rigid airframe with four independent rotors. A body-fixed frame denoted {B}, is attached to the vehicle, and an inertial frame {A} is fixed to the ground. Following the usual convention in aerospace, the z-axis of the inertial frame is pointing downwards and the body-fixed frame corresponds to the inertial frame when the vehicle is in hover and oriented to the north. Let m denote the mass of the vehicle. Let V ∈ {B} ≡ R 3×1 denote the linear velocity and Ω ∈ {B} ≡ R 3×1 denote the angular velocity of the body frame {B} with respect to the inertial frame {A} expressed in {B}. Let R ∈ SO(3) denote the rotation matrix representing the orientation of the body-fixed frame {B} with respect to the inertial frame {A}. The gravitational acceleration expressed in the inertial frame {A} is given by g -→ e 3 where -→ e 3 = [0; 0; 1] is the unit vector in the z-axis.

Since the mass m of the airframe of a quadrotor is a scalar constant, we can write a scaled version of the linear dynamics of the vehicle as

V = -Ω × V + gR ⊤-→ e 3 - 1 m T - → e 3 - 1 m H,
where T ∈ R denotes the vertical thrust and H ∈ R 3×1 the horizontal drag force generated by the rotors. For slow moving quadrotors, H ⊤-→ e 3 = 0. More details on the thrust and drag terms are provided in Section III where we discuss the aerodynamics of a quadrotor in detail. In this paper, we will not be concerned with the attitude dynamics and there is no need to model them.

The second order linear kinematics of a quadrotor are given by

V = -Ω × V + gR ⊤-→ e 3 + a, (1)
where a ∈ {B} is the specific acceleration of {B} with respect to {A} expressed in {B}. That is a is the sum of all the exogenous accelerations applied to the rigid-body.

Coriolis and gravitational accelerations are associated with the internal dynamics of the rigid-body and are not modeled by the specific acceleration. Clearly, the second order kinematics and the scaled dynamics are closely related. However, we draw an important distinction associated with the nature of the exogenous signals. In particular, the second order kinematics incorporate the specific acceleration a, which in our case can be measured using an inertial measurement unit (IMU), while the dynamics incorporate the forces T and H that we will model using an aerodynamic analysis. Trivially, one has

a = - 1 m T - → e 3 - 1 m H. (2)
Given that a is measured, then this equation provides direct measurements of T and H. In practice, this measurement should be treated with care due to high noise levels in a as well as possible bias effects. However, the information provided by this correspondence lies at the core of a number of recent works for velocity and attitude estimation for quadrotors [START_REF] Abeywardena | Improved state estimation in quadrotors mavs[END_REF], [START_REF] Leishman | Quadrotors and accelerometers: state estimation with an improved dynamic model[END_REF], [START_REF] Mahony | Multirotor aerial vehicles: Modeling, estimation, and control of quadrotor[END_REF], [START_REF] Martin | The true role of accelerometer feedback in quadrotor control[END_REF].

The system model that we consider consists of the second order linear kinematics along with the first order attitude kinematics

V = -Ω × V + gR ⊤-→ e 3 + a, (3a) Ṙ = RΩ × , (3b)
where the linear operator (.) × maps any vector in R 3×1 to its corresponding skew-symmetric matrix in so(3) such that x × y is equal to the cross product x × y for all x, y ∈ R 3×1 .

Define |.| to be the Euclidian norm in R 3×1 and ||.|| to be the Frobenius norm in R 3×3 . That is for all A ∈ R 3×3 , ||A|| = tr(A T A), where tr(.) is the trace.

The set SO(3) denotes the special orthogonal Lie-group and so(3) denotes its Lie-algebra, the set of skew-symmetric matrices. Define U(3) to be the set of all invertible upper triangular including entries along the diagonal, matrices in R 3×3 . It is easily verified that U(3) is a Lie-group with Lie algebra u(3) given by the set of all upper triangular matrices: that is including upper triangular matrices that are not invertible. Any invertible matrix A has a unique decomposition into an orthogonal and upper triangular part

A = U Q, Q ∈ SO(3), U ∈ U(3).
The factorisation is locally unique, and globally unique if the diagonal elements of U are required to be positive.

Observe that the Lie-algebra's so(3) and u(3) span R 3×3 in a natural manner under the usual vector space addition of matrices. Let P as (anti-symmetric) and P ut (upper-triangular) denote the associated projections onto so(3) and u(3) respectively. For any matrix A ∈ R 3×3 given by

A =   a 11 a
  ∈ u(3).
It is easily verified that P 2 ut = P ut and P 2 as = P as . Moreover, P ut (A) + P as (A) = A, for any matrix A. It is straightforward to verify that P as (P ut (A)) = 0 and P ut (P as (A)) = 0.

III. AERODYNAMIC POWER AND ESTIMATION OF V

This section presents a lumped model for the aerodynamics of a quadrotor. The approach is based on recent work in measuring aerodynamic power for quadrotor applications [START_REF] Bangura | Aerodynamic power control for multirotor aerial vehicles[END_REF]. The model is based on using simplified momentum theory to relate thrust, horizontal drag force and aerodynamic power and combining these relationships with (2).

A. Horizontal drag force

There are several recent works [START_REF] Abeywardena | Improved state estimation in quadrotors mavs[END_REF], [START_REF] Bangura | Nonlinear dynamic modeling for high performance control of a quadrotor[END_REF], [START_REF] Leishman | Quadrotors and accelerometers: state estimation with an improved dynamic model[END_REF], [START_REF] Mahony | Multirotor aerial vehicles: Modeling, estimation, and control of quadrotor[END_REF], [START_REF] Martin | The true role of accelerometer feedback in quadrotor control[END_REF] that have discussed and exploited models for horizontal drag force in estimating attitude and velocity for quadrotor and other aerial vehicles. The horizontal force associated with induced drag, blade flapping, and translational drag [START_REF] Bangura | Nonlinear dynamic modeling for high performance control of a quadrotor[END_REF] all manifest as bilinear functions of horizontal velocity and thrust magnitude

H = -T K r V ∈ {B}, (4)
where T ∈ R is the thrust,

K r =   c 0 0 0 c 0 0 0 0   ,
and c > 0 is a lumped parameter that models the combined linear drag coefficients. In contrast, parasitic drag depends quadratically on vehicle velocity and since we consider quadrotors operating in near hover conditions, its contribution is negligible.

B. Vertical thrust force

Modeling the thrust or vertical force of a rotor is complex because of the dependence of thrust on the induced velocity v i of air passing through the rotor; that is air velocity at the rotor plane that is induced by action of the rotor in generating thrust and separate from the relative wind experienced due to the motion of the vehicle through the air.

We will assume that the aerodynamic conditions across all four rotors are similar and there is an 'average' rotor speed ̟ that can be used to model the aerodynamic performance of the quadrotor using momentum theory. We will also assume a single average induced velocity v i for all four rotors. We will use vector notation for the induced velocity v i ∈ {B}. However, following the lead of classical helicopter books [START_REF] Leishman | Principles of helicopter aerodynamics[END_REF], only the vertical component of the induced velocity is considered in the aerodynamic model. That is

v i = v i z - → e 3 ,
where v i z ∈ R is a scalar. We will assume that all rotors have the same flapping angles and generate the same aerodynamic drag forces. These assumptions are reasonable if the quadrotor is flying in near hover conditions.

Consider a control volume for a single rotor shown in Figure 1. Momentum theory for a single rotor yields [START_REF] Leishman | Principles of helicopter aerodynamics[END_REF] T k = 2ρAv i z U, where T k denotes the thrust associated with rotor k, ρ is the density of air, A is the area of a rotor disc, and

U = |v i -V | (5)
is the magnitude of the total wind velocity, including both vehicle velocity and induced velocity of the air through the rotor. Note that the wind due to the vehicle motion is -V , the negative of the vehicle motion. Let c 1 = 8ρA total ∈ R, then the total thrust generated by all four rotors is modeled by

T = c 1 v i z U. (6)
Unlike the classical helicopter analysis found in books such as Leishman [START_REF] Leishman | Principles of helicopter aerodynamics[END_REF], where only the steady-state performance in the separate lift or drag directions is required, we will provide a combined model that includes transient effect for power in the rotor and both thrust and drag in a single model. The aerodynamic power P a supplied to the air consists of two components: vertical P T and horizontal P H power, associated with providing thrust and overcoming drag respectively.

P a = P T + P H , (7)
The power associated with motion of the air in the vertical direction is [START_REF] Leishman | Principles of helicopter aerodynamics[END_REF] P

T = T (v i z -V z ), (8)
where T v i z is the induced power (required to operate the rotor and overcome induce aerodynamic forces) and -T V z is the power to climb. Note that due to the choice of axis, when V z < 0 the vehicle is ascending and -T V z is positive as expected. In the horizontal plane one has

P H = -H ⊤ V > 0, (9)
where H ∈ {B} is a vector quantity that lies uniquely in the horizontal plane and hence only the horizontal velocity contributes to the power P H .

To account for losses such as tip loss and rotational wake, there is an an efficiency factor between the mechanical power supplied and the aerodynamic power at steady state. We refer to this efficiency as the figure of merit (FoM) which is a number between 0 and 1 and is obtained from static thrust tests [START_REF] Bangura | Nonlinear dynamic modeling for high performance control of a quadrotor[END_REF]. It follows that the mechanical power P m , supplied to the rotor hub factors into

P m = P r + 1 FoM P a , (10)
where the power supplied to accelerate or decelerate the rotor is P r = I̟ ̟. If the torque on the rotor shaft is τ , then the mechanical power supplied to the rotor hub is given by P m = τ ̟. The torque τ can be measured by identifying motor parameters and measuring current. The rotor speed ̟ is measured in the operation of all Electronic Speed Controllers (ESCs) and a complementary filter for estimating ̟ is provided in [START_REF] Bangura | Aerodynamic power control for multirotor aerial vehicles[END_REF]. It follows that P r and P m can be measured, and once FoM is identified then P a can be treated as a measured variable. Counting equations, one has two constraints from (4) (note that although this is a compound equation with three components, only the two associated with the horizontal rotor plane contain meaningful information), and one constraint from both ( 6) and ( 8) in six unknowns V x , V y , V z , v i z , H, and T . However, adding the constraint measurement (2) that has three additional constraints ensures that this system of equations can be resolved algebraically.

Given that a is measured, then T , V x and V y can be solved for directly using ( 2) and ( 4). This is the same approach used in [START_REF] Abeywardena | Improved state estimation in quadrotors mavs[END_REF], [START_REF] Bangura | Nonlinear dynamic modeling for high performance control of a quadrotor[END_REF], [START_REF] Leishman | Quadrotors and accelerometers: state estimation with an improved dynamic model[END_REF], [START_REF] Mahony | Multirotor aerial vehicles: Modeling, estimation, and control of quadrotor[END_REF], [START_REF] Martin | The true role of accelerometer feedback in quadrotor control[END_REF] to measure horizontal velocity of the vehicle. From ( 4) and ( 9) and using T , V x and V y then P H is computed. It follows that using ( 7) and ( 10) then P T can be computed. With P T known, then ( 8) is used to compute

(v i z -V z ).
Finally, this is used to compute the total wind velocity U and substituted into [START_REF] Hua | Attitude estimation for accelerated vehicles using gps/ins measurements[END_REF]

to compute v i = v i z - → e 3 .
Once v i z is computed the z-axis velocity V z is available by substituting back into the known value of (v i z -V z ). The key outcome of this algebraic process is the computation of a measured value of the body-fixed frame velocity V = (V x , V y , V z ). We emphasise that this "measurement" is independent of orientation of the vehicle, unlike the approximations made in prior work [START_REF] Abeywardena | Improved state estimation in quadrotors mavs[END_REF], [START_REF] Bangura | Nonlinear dynamic modeling for high performance control of a quadrotor[END_REF], [START_REF] Leishman | Quadrotors and accelerometers: state estimation with an improved dynamic model[END_REF], [START_REF] Mahony | Multirotor aerial vehicles: Modeling, estimation, and control of quadrotor[END_REF], [START_REF] Martin | The true role of accelerometer feedback in quadrotor control[END_REF]. The assumptions on the aerodynamic conditions still limit applications to those where the vehicle is in hover or near hover conditions. However, the additional identification of z-axis velocity provides the potential for a more advanced observer design than developed in prior works.

IV. OBSERVER DESIGN

In this section, a nonlinear observer for estimating the body-fixed-frame velocity for a quadrotor is proposed. The approach is based on deterministic Lyapunov design techniques and draws heavily from approach pioneered in [START_REF] Hua | Attitude estimation for accelerated vehicles using gps/ins measurements[END_REF].

For a quadrotor vehicle with an inertial measurement unit (IMU), the angular velocity Ω and specific acceleration a are available. If the same vehicle is equipped with electronic speed controllers (ESCs) that measure both rotor speed and combined torque to all four motors then the mechanical power supplied to the rotor hubs can be measured. As we have seen in Section III, this is sufficient to provide a measurement of the body-fixed-frame velocity V ∈ {B} of the vehicle. If the IMU is also equipped with magnetometers, then these can be used to provide an additional measurement

µ = R ⊤ μ ∈ {B}
for μ ∈ {A} the inertial "known" magnetic field. In practice, the magnetic field measurement is often corrupted by onboard magnetic fields and cannot be used for attitude estimation. For this reason we will initially develop the proposed observer in the case where the magnetic field is not available.

The goal of the observer design is to provide estimates V ∈ R 3 and R ∈ SO(3) of the body-fixed-frame velocity and attitude of the quadrotor. We will distinguish between the true velocity V (t) and the measured velocity V (t) in the following theorem, for the sake of clarity, although the deterministic stability analysis is based on the relationship

V (t) = V (t).
Assumption 1 Assume that the trajectory of the quadrotor is sufficiently smooth such that Ω(t), Ω(t), V (t) and V (t) are bounded signals.

Theorem 1 Consider system (3) with Ω, a and V measured. Consider the observer

V = -Ω × V + gX ⊤-→ e 3 + a -k 1 ∆ 1 , V (0) = V (0) (11a)
Ẋ = U R + U Ṙ, X(0) = U (0) R(0) (11b)
where k 1 > 0 is a scalar gain and X = U R is the uppertriangular orthogonal decomposition of a general matrix X ∈ R 3×3 . The dynamics of U and R are given by

U = -k 2 gU P ut (U -1-→ e 3 ( V -V ) ⊤ R⊤ ), U (0) = I 3 (12a) Ṙ = RΩ × + k 2 R∆ 2 , R(0) = I 3 (12b) ∆ 1 = V -V (12c) ∆ 2 = -g R⊤ P as (U -1-→ e 3 ( V -V ) ⊤ R⊤ ) R, ( 12d
)
where k 2 > 0 is a scalar gain. Suppose that Assumption 1 is satisfied, then for almost all initial conditions, the estimate

V (t) → V (t) and R⊤-→ e 3 → R ⊤-→ e 3 .
Proof: Define a velocity error

Ṽ = V -V.
The time derivative of Ṽ is given by

V = V - V = -Ω × Ṽ + (X -R) ⊤ g - → e 3 -k 1 ∆ 1 . (13)
Define a candidate Lyapunov function

L := 1 2 Ṽ ⊤ Ṽ + 1 2k 2 ||X -R|| 2 .
The time derivative of L satisfies

L = Ṽ ⊤ (X -R) ⊤ g - → e 3 -k 1 Ṽ ⊤ ∆ 1 + 1 k 2 tr((X -R) ⊤ ( U R + k 2 X∆ 2 )) (14) = tr( Ṽ ⊤ (X -R) ⊤ g - → e 3 ) A + 1 k 2 tr((X -R) ⊤ ( U R + k 2 X∆ 2 )) B -k 1 Ṽ ⊤ ∆ 1 . C ( 15
)
Since the rotation matrix R is unknown, the innovation terms ∆ 2 and U must be chosen properly in order to cancel the expressions A and B in the previous equation while ∆ 1 is chosen to make L negative. Considering just terms A and B from above, one has

tr( Ṽ ⊤ (X -R) ⊤ g - → e 3 ) + 1 k 2 tr((X -R) ⊤ ( U R + k 2 X∆ 2 )) = tr((X -R) ⊤ (g - → e 3 Ṽ ⊤ + 1 k 2 U R + X∆ 2 )) = 1 k 2 tr((X -R) ⊤ (k 2 g - → e 3 Ṽ ⊤ R⊤ + U + k 2 U R∆ 2 R⊤ ) R) = 1 k 2 tr((X -R) ⊤ U Γ R), (16)
where the matrix Γ ∈ R 3×3 is given by

Γ = k 2 gU -1-→ e 3 Ṽ ⊤ R⊤ + U -1 U + k 2 R∆ 2 R⊤ .
Since X can be an arbitrary matrix, then the only choice to cancel terms A and B from (15) for all possible X is to set Γ ≡ 0. Using the projectors defined in Section II, one can split Γ into its components

Γ = P as (Γ) + P ut (Γ), (17)
with

P as (Γ) = k 2 gP as (U -1-→ e 3 Ṽ ⊤ R⊤ ) + k 2 R∆ 2 R⊤ P ut (Γ) = k gP lt (U -1-→ e 3 Ṽ ⊤ R⊤ ) + U -1 U since ∆ 2 = -∆ ⊤ 2
, P as (U -1 U ) = 0 and P ut (k 2 R∆ 2 R⊤ ) = 0. Substituting (12d) for ∆ 2 and (12a) for U , it is easily verified that P as (Γ) = 0 = P lt (Γ) and hence Γ ≡ 0. Substituting (12c), one obtains

L = -k 1 Ṽ ⊤ Ṽ ≤ 0. (18)
Since the time derivative of L is semi-negative definite and L is positive definite then Ṽ and X are bounded. In view of (13) and Assumption 1, one deduces that V is bounded and it follows that L is also bounded. This is sufficient to ensure that L is uniformly continuous along trajectories of the system. Applying Barbalat's lemma [START_REF] Khalil | Nonlinear systems[END_REF] ensures the convergence of L → 0 that implies the convergence of Ṽ to 0.

The same procedure is performed to prove that V is bounded (since Ω, V , Ẋ and Ṙ are bounded) and consequently to demonstrate the uniform continuity of V . Barbalat's lemma ensures the convergence of V to 0 which implies from (13) the convergence of

X ⊤ g - → e 3 to R ⊤ g - → e 3 .
Finally, substituting X = U R, one sees that R⊤ U ⊤-→ e 3 → R ⊤-→ e 3 . Taking norms, one has |U ⊤-→ e 3 | → 1, that is the (3, 3) entry of U converges to unity and U ⊤-→ e 3 → -→ e 3 . It follows that R⊤-→ e 3 → R ⊤-→ e 3 by continuity.

Since R⊤-→ e 3 → R ⊤-→ e 3 , then this observer will effectively identify the pitch and roll components of the attitude of a quadrotor. The remaining degree of freedom in the orientation, corresponding to yaw rotation around the -→ e 3 axis is unobservable without additional measurements. If magnetometer measurements are available, then the filter can be augmented with these measurements to obtain full estimation of the attitude.

V. SIMULATION RESULTS

In this section, we illustrate through simulation results, the performance and robustness of the proposed observer. Simulations are performed for a model of a quadrotor aerial vehicle. The complete observer scheme is shown in Figure 2. The trajectories of both V (t), a(t) and R(t) are generated assuming that the vehicle performs circular motion in the inertial frame combined with a sinusoidal changes of altitude. The vehicle's linear velocity, given in the inertial frame is v(t) = [0.5 sin(α 1 t); 0.5 sin(α 1 t + π 2 ); sin(α 2 t)] ∈ {A} with α 1 = 1 2 and α 2 = 1 3 . For the observer given by ( 11), the gains are k 1 = 10 and k 2 = 1. If R x (θ), R y (θ), R z (θ) are the rotations about the x, y, z-axis respectively by θ, then the initial conditions are R(0) = R x (10)R y (-20)R z (30) and U (0) = [.5 0 0; .2 .3 0; .4 .5 .2], ensuring that R ∈ SO(3) and U is an upper triangular matrix. In the following, we assume that we can measure P a = F oM (P m -P r ). In terms of practical implementation, [START_REF] Bangura | Aerodynamic power control for multirotor aerial vehicles[END_REF] showed that P a can be measured from the current and ̟ along with a complementary filter which estimates ̟ onboard an electronic speed controller board [START_REF] Team | The AutoQuad ESC32 -a yet unseen electronic speed controller[END_REF]. We also consider that the onboard IMU gives the specific acceleration and angular velocity measurements. In order to have a more realistic simulation, a random signal with a constant power spectral density is added to all inputs of the observer. The range of this signal is in accordance with the sensors used. For the accelerometer measurements, a high frequency noise of 200% the magnitude of typical accelerations experienced by such vehicles is added.

The noise in Ω of relative magnitude 20% is also added. For the P a used, the experimental results obtained in [START_REF] Bangura | Aerodynamic power control for multirotor aerial vehicles[END_REF] showed that the estimated values have relative noise that are under 10%. Hence 10% noise was added to P a . The performance of the proposed observer is shown in Figures 3 and4. The plots show the convergence and the performance of the proposed observer for the three body velocities. Note that despite the magnitude of noise in the input as well as in V , the outputs of the observer are smooth and track the true real body velocities. Finally, Figure 4 shows the convergence of R⊤ e 3 to R ⊤ e 3 and the convergence of U (3, 3) to 1.

Fig. 1 .

 1 Fig. 1. Momentum theory control volume and rotor hub showing the various powers for one of the rotors. In the control volume, we have the axis definition for {B} which may not always align with the rotor tip path plane. Vc = -Vz is the climb velocity described in the helicopter literature.

Fig. 2 .

 2 Fig. 2. Observer scheme

  12 a 13 a 21 a 22 a 23 a 31 a 32 a 33

			  ,
	then one has
	P as (A) =	  a 21 0 -a 21 -a 31 0 -a 32 a 31 a 32 0	  ∈ so(3),
	and
	P ut (A) =	a 11 a 12 + a 21 a 13 + a 31  0 a 22 a 23 + a 32  0 0 a 33

ACKNOWLEDGMENT

This research was supported by the Australian Research Council through Discovery Grant DP120100316 "Integrated High-Performance Control of Aerial Robots in Dynamic Environments".

Fig. 3. Results of the simulated V ∈ {B}. With V being the true velocity, we can see that despite the 200% noise in accelerometer measurements that is reflected in V determined algebraically as in Section III, the estimated velocities V track quite well the true velocity V . This is indicated by the small difference ǫ between estimated V and true V .

Time [s]

Fig. 4. Attitude error and scaling factor. We can see that the norm of the attitude error between estimated ( R⊤ e 3 ) and true (R ⊤ e 3 ) converges to zero in less than 0.5s which implies R⊤ e 3 → R ⊤ e 3 . This performance is further reaffirmed by U e 3 → 1.


Carinne Ndjembo Ezougou

Fathia Ben-Rached

David K Moss

Jing-Wen Lin

Sally Black

Ellen Knuepfer

Judith L Green

Shahid M Khan

Amitabha Mukhopadhyay

Chris J Janse

Isabelle Coppens

´le `ne Yera

Anthony A Holder
email: aholder@nimr.mrc.ac.uk

Gordon Langsley
email: gordon.langsley@inserm.fr

Myristoylated Rab GTPase That Is Targeted to the Parasite's Plasma and Food Vacuole Membranes

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

Introduction

Plasmodium falciparum is a parasite that infects red blood cells (RBC) and causes severe human malaria, a disease that kills approximately 750,000 people per year, mostly children living in tropical Africa [1]. The malaria parasite spends much of its life cycle inside RBC, cells that provide it with an abundant food supply in the form of haemoglobin and a degree of protection from the host's immune system. For survival inside the RBC the parasite has to both import nutrients and to export metabolic waste products and like other eukaryotes it uses Rab GTPases to regulate vesicular trafficking [START_REF] De Castro | Identification of a family of Rab G-proteins in Plasmodium falciparum and a detailed characterisation of pfrab6[END_REF] [3] [START_REF] Langsley | Comparative genomics of the Rab protein family in Apicomplexan parasites[END_REF].

Rabs are molecular switches belonging to the Ras-superfamily that are conserved from yeast to humans and regulate in space and time the budding and fusion of intracellular vesicles from donor to acceptor membranes [START_REF] Zerial | Rab proteins as membrane organizers[END_REF]. Rabs vary in size from 20 to 29 kDa and were first identified in yeast and named Ypt proteins. Eleven Ypts were characterized in Saccharomyces cerevisiae and seven in Schizosaccharomyces pombe [START_REF] Novick | Regulation of vesicular traffic by a GTP-binding protein on the cytoplasmic surface of secretory vesicles in yeast[END_REF], whereas higher eukaryotes may possess more than 60 different Rabs [START_REF] Grosshans | Rabs and their effectors: achieving specificity in membrane traffic[END_REF]. Rab proteins exist in GTP-bound and GDP-bound states with switching from one state to another regulated by GTP Exchange Factors and GTPase Activating Proteins [START_REF] Pfeffer | Rab GTPases: specifying and deciphering organelle identity and function[END_REF] [START_REF] Pfeffer | Vesicle tethering factors united[END_REF]. GDP-bound Rabs are considered to be inactive and are associated with Rab GDP Dissociation Inhibitors (rabGDIs), whereas GTP-bound Rabs are considered to be active and able to recruit and interact with their effectors, which is a key specific Rab-mediated function [START_REF] Stenmark | The Rab GTPase family[END_REF]. Rabs from different species display a significant degree of sequence conservation in the GTPbinding domain, but they possess a variable N-terminal sequence responsible for recruitment of effectors and a hyper-variable Cterminal domain. The C-terminal domain confers the specific intracellular location for each Rab and has a motif containing one or two cysteines (CAAX, CC, CXC, CCX, CCXX) necessary for the isoprenylation by geranylgeranyl-transferases that is required for their association with vesicle membranes [START_REF] Chavrier | Hypervariable C-terminal domain of rab proteins acts as a targeting signal[END_REF].

Mammalian cells encode three Rab5 isoforms (Rab5A, Rab5B and Rab5C) that share a high level of sequence identity and are involved in homotypic and heterotypic fusion of early and late endosomes [START_REF] Hirota | A role for Rab5 activity in the biogenesis of endosomal and lysosomal compartments[END_REF] [START_REF] Gorvel | rab5 controls early endosome fusion in vitro[END_REF]. In spite of this homology different Rab5 isoforms can regulate separate functions, for example Rab5A and Rab5B, but not Rab5C, regulate the transfer of epidermal growth factor receptor from early to late endosomes [START_REF] Chen | Rab5 isoforms differentially regulate the trafficking and degradation of epidermal growth factor receptors[END_REF]. siRNA knockdown of individual Rab5 isoforms demonstrated that collectively they organise early endosomes, late endosomes and lysosomes and can be considered master organisers of the endocytic system [START_REF] Zeigerer | Rab5 is necessary for the biogenesis of the endolysosomal system in vivo[END_REF]. However, lineage-specific expansions of the gene family combined with functional diversification have contributed to species-specific variations in membrane trafficking. A classic example is the plant-specific Rab5-like GTPase called ARA6 that mediates trafficking from endosomes to the plasma membrane of Arabidopsis thaliana, where it contributes to the regulation of the plant's response to saline stress [START_REF] Ebine | A membrane trafficking pathway regulated by the plant-specific RAB GTPase ARA6[END_REF][START_REF] Ebine | Endosomal trafficking pathway regulated by ARA6, a RAB5 GTPase unique to plants[END_REF][START_REF] Ueda | Ara6, a plant-unique novel type Rab GTPase, functions in the endocytic pathway of Arabidopsis thaliana[END_REF].

Multiple reciprocal BLAST analyses gave rise to the P. falciparum Rab (PfRab) nomenclature and identified a family of 11 PfRabs [START_REF] Quevillon | The Plasmodium falciparum family of Rab GTPases[END_REF]. This family contains three Rab5 isoforms with PfRab5A (PF3D7_0211200) belonging to a restricted orthology group (OG4_36791) present only in Apicomplexa parasites known to invade erythrocytes (Plasmodia, Babesia and Theileria) [START_REF] Rached | Construction of a Plasmodium falciparum Rab-interactome identifies CK1 and PKA as Rab-effector kinases in malaria parasites[END_REF]. However, PfRab5A has an insertion of 30 amino acids between the RabF1 and RabF2 motifs (Rab-effector binding motifs) [START_REF] Quevillon | The Plasmodium falciparum family of Rab GTPases[END_REF] that is not observed in Rab5A of Theileria and Babesia, suggesting that some PfRab5A effectors might be Plasmodium-specific. Putative PfRab5A-effectors might be involved in the uptake of haemoglobin, as GFP-tagged PfRab5A decorates vesicles containing haemoglobin [START_REF] Elliott | Four distinct pathways of hemoglobin uptake in the malaria parasite Plasmodium falciparum[END_REF]. The PfRab5B orthology group (OG4_18709) is found more broadly; the C-terminus of PfRab5B (PF3D7_1310600) has no lipid-modification motif [START_REF] Quevillon | The Plasmodium falciparum family of Rab GTPases[END_REF][START_REF] Rached | Construction of a Plasmodium falciparum Rab-interactome identifies CK1 and PKA as Rab-effector kinases in malaria parasites[END_REF], like human Rab8 and Rab23 [START_REF] Casey | Protein prenyltransferases[END_REF] and ARA6 [START_REF] Ueda | Ara6, a plant-unique novel type Rab GTPase, functions in the endocytic pathway of Arabidopsis thaliana[END_REF]. The PfRab5C (PF3D7_0106800) orthology group (OG4_10168) is the largest [START_REF] Rached | Construction of a Plasmodium falciparum Rab-interactome identifies CK1 and PKA as Rab-effector kinases in malaria parasites[END_REF]. PfRab5C looks like a classical Rab5, but comparison of a 3Dmodel of PfRab5C with the known 3D-structure of mouse Rab5C showed that these two Rab5s present different amino acids at their effector-interaction surfaces, implying that PfRab5C has the potential to recruit parasite-specific effectors [START_REF] Quevillon | The Plasmodium falciparum family of Rab GTPases[END_REF].

As S. cerevisiae possesses 11 Ypts/Rabs like P. falciparum, the two ''Rabomes'' were compared and a putative potential function was attributed to each PfRab [START_REF] Rached | Construction of a Plasmodium falciparum Rab-interactome identifies CK1 and PKA as Rab-effector kinases in malaria parasites[END_REF]. In addition, potential PfRab functions have been inferred from over-expression of wild type and dominant-negative (GDP-on) mutants of Toxoplasma gondii (Tg) Rabs [START_REF] Agop-Nersesian | Rab11A-controlled assembly of the inner membrane complex is required for completion of apicomplexan cytokinesis[END_REF] [24] [START_REF] Kremer | An overexpression screen of Toxoplasma gondii Rab-GTPases reveals distinct transport routes to the micronemes[END_REF]. Plasmodium versus Toxoplasma comparison is only valid for true orthologues in the two Apicomplexa and is therefore not applicable to PfRab5A, but some potential Plasmodium-specific functions have been inferred from studies on the over-expression of wild type TgRab5B and a FKBPmyc-TgRab5C fusion [START_REF] Kremer | An overexpression screen of Toxoplasma gondii Rab-GTPases reveals distinct transport routes to the micronemes[END_REF]. Over-expressed TgRab5B was found to concentrate at endosomal-like compartments and to a lesser extent at the surface of T. gondii, whereas over-expression of TgRab5C variants blocked secretion of specific cargo to a subset of micronemes at the apex of the parasite [START_REF] Kremer | An overexpression screen of Toxoplasma gondii Rab-GTPases reveals distinct transport routes to the micronemes[END_REF]. Over-expression of both TgRab5B and TgRab5C FKBPmyc fusions was deleterious for T. gondii in culture.

The PfRab-interactome predicted that Casein Kinase 1 (PfCK1; PF3D7_1136500.1) is a specific PfRab5B-interacting protein and indeed, PfRab5B physically interacts with PfCK1 in vitro [START_REF] Rached | Construction of a Plasmodium falciparum Rab-interactome identifies CK1 and PKA as Rab-effector kinases in malaria parasites[END_REF]. The subcellular distribution of PfCK1 has not yet been described [START_REF] Barik | Identification, cloning, and mutational analysis of the casein kinase 1 cDNA of the malaria parasite, Plasmodium falciparum. Stage-specific expression of the gene[END_REF], but in T. gondii CK1 is found both in the cytosol and at the plasma membrane [START_REF] Donald | Characterization of two T. gondii CK1 isoforms[END_REF]. We decided to test therefore, whether in vivo PfRab5B is competent to traffic PfCK1 to different subcellular locations in P. falciparum. Moreover, as PfRab5B lacks the necessary motif for isoprenylation, yet possesses recognisable myristoylation and palmitoylation motifs [START_REF] Quevillon | The Plasmodium falciparum family of Rab GTPases[END_REF], we examined directly whether the N-terminal glycine residue of PfRab5B is a bona fide substrate for P. falciparum N-myristoyltransferase (NMT) [START_REF] Gunaratne | Characterization of N-myristoyltransferase from Plasmodium falciparum[END_REF] in vitro and whether in vivo its myristoylation is involved in targeting PfRab5B to specific subcellular membranes. Based on our results we propose that myristoylation (and palmitoylation) of PfRab5B allows it to traffic specific cargo to the parasite's plasma and food vacuole membranes.

Results

The N-terminal glycine residue in PfRab5B is required for N-myristoylation

To test the validity of the myristoylation motif prediction we produced recombinant PfRab5B protein and variants in which either the glycine or the glycine and cysteine residues at the Nterminus of PfRab5B (MGCSS) were changed to alanine. These His-tagged full-length recombinant wild type and variant PfRab5B proteins were purified and used as substrates in an N-myristoylation assay using recombinant P. falciparum N-myristoyl-transferase (PfNMT) [START_REF] Gunaratne | Characterization of N-myristoyltransferase from Plasmodium falciparum[END_REF] and tritiated myristoyl-CoA [START_REF] Rees-Channer | Dual acylation of the 45 kDa gliding-associated protein (GAP45) in Plasmodium falciparum merozoites[END_REF] (see Figure 1). Wild type and G2A-variant recombinant PfARF protein (PF3D7_1020900) were used as controls [30] [29]. The results showed that PfRab5B is a substrate for PfNMT (although the efficiency of myristoyl transfer was less than that for PfARF) and the presence of a glycine at position two is required for PfRab5B myristoylation in vitro, which occurs via peptide bond formation with the a-amino group of the N-terminal glycine. As expected, modification of both the myristoylation and the palmitoylation site (G2AC3A) had no additional effect. PfRab5B myristoylation in the parasite is essential for its targeting to membranes

We generated a transgenic P. falciparum line that expresses GFP fused to the N-terminal 28 amino acids of PfRab5B (see Table S1 in File S1). An episomally maintained plasmid was used for transfection, which expresses the fusion protein under the control of the msp3 promoter [START_REF] Ridzuan | Subcellular location, phosphorylation and assembly into the motor complex of GAP45 during Plasmodium falciparum schizont development[END_REF] [START_REF] Knuepfer | RON12, a novel Plasmodium-specific rhoptry neck protein important for parasite proliferation[END_REF]. Asexual blood stages of this transgenic line were analysed by fluorescence microscopy. GFPfluorescence was consistent with a location of the fusion protein in association with the parasite's plasma membrane (Figure 2, top left hand panel). As a positive control, we used a second transgenic P. falciparum line that expresses GFP fused to 29 amino acids from the N-terminus of glideosome associated protein (GAP)45 (Table S1 in File S1). GAP45 is known to be associated with the plasma membrane at its N-terminus [START_REF] Ridzuan | Subcellular location, phosphorylation and assembly into the motor complex of GAP45 during Plasmodium falciparum schizont development[END_REF]. In asexual blood stages of this transgenic parasite GFP was observed at the plasma membrane (Figure 2, top right hand panel). In transgenic parasites that expressed GFP fused to G2A variants of both PfRab5B and GAP45 the association of GFP with the plasma membrane was lost and we observed a diffuse GFP-staining throughout the cytoplasm of the parasites (Figure 2, second row panels). This cytoplasmic GFP-staining was also observed in blood stages expressing GFPfusion proteins in which additional substitutions had been introduced (Table S1 in File S1): the C3A change in PfRab5B or C5A in GAP45 to give the double G+C variants (Figure 2, third row panels). The C to A variants alone allowed some membrane association (Figure 2, bottom row). Clearly, the N-terminal 28 or 29 amino acids in combination with myristoylation of G2 of both PfRab5B and GAP45, respectively, are sufficient to target the GFP-fusion proteins to the parasite's plasma membrane.

PfRab5B-GFP and GAP45-GFP fusion proteins are embedded in membranes of P. falciparum parasites

Schizonts of transgenic parasites expressing GFP fused to Nterminal sequences of PfRab5B and GAP45 (28 or 29 amino acids, respectively) were purified and submitted to subcellular fractionation using hypotonic and high salt buffers to release cytosolic proteins, followed by high pH carbonate extraction to fractionate peripheral membrane (carbonate soluble) and integral membrane (carbonate insoluble) proteins. The proteins in the four different fractions were resolved by SDS-PAGE and analysed by western blotting with an anti-GFP antibody (Figure 3). When the wild type PfRab5B N-terminal sequence was fused to GFP (PfRab5B 28 -GFP) the protein was found in all subcellular fractions (to differing degrees) indicating that not all of this fusion protein was associated with membranes despite the presence of both G2 and C3 residues. In contrast, both the GFP-fusion proteins containing the G2A and C3A variants were located solely in the hypotonic lysis soluble fraction, as was the variant with both substitutions (Figure 3, track 1). The GFP protein fused to the wild type GAP45 N-terminal sequence (GAP45 29 -GFP) was located exclusively in the carbonate-insoluble fraction (track 4) typical for a membrane-associated protein. Similar to the PfRab5B 28 G2A variant of the GFP fusion protein, GAP45 29 G2A was found exclusively in the soluble fraction consistent with myristoylation alone being necessary and sufficient for membrane association. In contrast, GAP45 29 C5A was distributed across the different fractions, suggesting that palmitoylation of the N-terminus of GAP45 contributes to its membrane association, but it is neither necessary, nor sufficient.

Endogenous PfRab5B is associated with both the food vacuole and plasma membranes

We generated antibodies to PfRab5B and verified their specificity by western blot (see Figure S1). These specific PfRab5B-specific antibodies were then used to determine the subcellular location of PfRab5B by immuno-electron microscopy. In Figure 4 micrographs are shown of two different stages of intraerythrocytic development, mature schizonts (Panel A) and a trophozoite (Panel B). Gold-labelled anti-PfRab5B antibodies decorated the parasite plasma membrane (PPM, white arrows) and were also observed associated with the food vacuole (FV) and neighbouring vesicles, as shown for the trophozoite in panel B. The distribution of endogenous PfRab5B was therefore similar to that of PfRab5B-GFP suggesting that the plasma-and food vacuole membrane associations are the result of myristoylation and/or palmitoylation at the N-terminus of PfRab5B.

Endogenous PfRab5B colocalises with merozoite surface protein (MSP)1 and PfCK1, but not with haemoglobin GFP-tagged PfRab5A has been described as associated with vesicles harbouring haemoglobin imported by the parasite from the erythrocyte cytosol [START_REF] Elliott | Four distinct pathways of hemoglobin uptake in the malaria parasite Plasmodium falciparum[END_REF]. The association of PfRab5B with the food vacuole membrane and adjacent vesicles led us, therefore, to look at an association with haemoglobin-containing structures using indirect immunofluorescence. Specific anti-PfRab5A antibodies (see Figure S1) were used as a positive control. In contrast to PfRab5A (r = 0.724; n = 3) we could find no evidence for colocalisation between PfRab5B and haemoglobin (HBA1, r = 0.081; n = 3) in trophozoites, the developmental stage that imports haemoglobin (Figure 5, panel A).

Following merozoite invasion, a C-terminal fragment of MSP1 is associated with the developing food vacuole membrane [START_REF] Dluzewski | Formation of the food vacuole in Plasmodium falciparum: a potential role for the 19 kDa fragment of merozoite surface protein 1 (MSP1(19))[END_REF], a localisation that is similar to that observed for PfRab5B, as shown in Figure 4. We confirmed therefore that antibodies to MSP1 and PfRab5B significantly (r = 0.849; n = 3) decorate the same structure (Figure 5, panel B).

As PfRab5B was found to physically interact with PfCK1 in vitro [START_REF] Rached | Construction of a Plasmodium falciparum Rab-interactome identifies CK1 and PKA as Rab-effector kinases in malaria parasites[END_REF] we looked for evidence that PfRab5B-positive vesicles might traffic the kinase to different subcellular locations. Depending on the developmental stage of the parasite (indicated by the number and size of DAPI-stained nuclei; Figure 5, panel C) the degree of association between PfRab5B and PfCK1 varied (r = 0.615; n = 3 in late schizonts and r = 0.475; n = 3 in young schizonts). Clearly, PfRab5B-positive vesicles can traffic more than one type of cargo, however the failure to observe any association of PfRab5B with haemoglobin suggests that PfRab5B and PfRab5A perform separate trafficking functions.

Plasmodium Rab5B is essential for asexual blood stage parasite development

Over-expression of individual Rab5s and/or Rab5 GDP-on variants in transgenic cells has the potential to generate phenotypes due to sequestering effector proteins that are shared by the different Rab5 isoforms, and therefore the use of such 'overexpressing' transgenic cells may hamper discrimination between the different functions of the Rab isoforms. Therefore, we decided to analyse the function of Plasmodium rab5b by attempting to switch off its expression using targeted gene-deletion by reverse genetics. We targeted the orthologous rab5b gene of the rodent malaria parasite P. berghei for gene deletion using standardized transfection methods [START_REF] Janse | High-efficiency transfection and drug selection of genetically transformed blood stages of the rodent malaria parasite Plasmodium berghei[END_REF]. In multiple experiments (Figure S2 in File S1) we were unable to select for mutants lacking the Pbrab5b gene indicating that PbRab5B has an essential function during asexual blood stage growth and multiplication. Moreover, these results also indicate that PbRab5A and PbRab5C cannot compensate for the function of PbRab5B.

Discussion

P. falciparum has a small family of 11 PfRabs and all but PfRab5B possess classical geranyl-geranylation motifs at the Cterminus. This isoprenyl modification is present in 10 of the PfRabs presumably anchoring these GTPases to membrane. In contrast, we have shown here that PfRab5B is myristoylated allowing membrane attachment at the N-terminal glycine. PfRab5B peptides have recently been detected in a global in vivo myristome of P. falciparum [START_REF] Wright | Validation of N-myristoyltransferase as an antimalarial drug target using an integrated chemical biology approach[END_REF]. Our inability to knock out the rab5b gene of P. berghei indicates that Plasmodium Rab5B performs functions essential for parasite asexual blood stage development. Potential functions may include trafficking MSP1-tagged vesicles from the plasma membrane to form the food vacuole membrane following RBC invasion [START_REF] Dluzewski | Formation of the food vacuole in Plasmodium falciparum: a potential role for the 19 kDa fragment of merozoite surface protein 1 (MSP1(19))[END_REF] and PfCK1 to the plasma membrane. However, in contrast to PfRab5A [START_REF] Elliott | Four distinct pathways of hemoglobin uptake in the malaria parasite Plasmodium falciparum[END_REF] we never observed PfRab5B associated with haemoglobin-positive struc-tures indicating that these two Plasmodium-specific Rab5 isoforms perform separate functions.

The presence of the wild type PfRab5B-GFP fusion protein in different subcellular fractions suggests that endogenous PfRab5B associates with varying affinity to membranes. However, as the wild type PfRab5B-GFP fusion protein was also found in the soluble fraction the results suggest that not all of the endogenous PfRab5B is myristoylated and/or palmitoylated at a given time. Substitution of the N-terminal glycine with alanine abolished the association of PfRab5B-GFP fusion protein with membrane, clearly demonstrating that myristoylation of the PfRab5B Nterminal amino acid sequence is necessary and sufficient to attach it to membranes. Although we did not directly examine the palmitoylation status of PfRab5B substituting a cysteine with an alanine ablated membrane-targeting of the GFP-fusion protein as judged by solubility analysis, although some apparent association with the plasma membrane was detected by live cell imaging. Palmitoylation of endogenous PfRab5B therefore, seems to contribute to its subcellular location and indeed, a recent analysis of the P. falciparum palmitoylation proteome identified PfRab5B, as a component [START_REF] Jones | Analysis of protein palmitoylation reveals a pervasive role in Plasmodium development and pathogenesis[END_REF][START_REF] Jones | Getting stuck in: protein palmitoylation in Plasmodium[END_REF]. Protein palmitoylation is reversible and dynamic and can result in proteins targeted to specific membrane domains [START_REF] Blaskovic | What does S-palmitoylation do to membrane proteins?[END_REF].

PfRab5B resembles ARA6 a Rab5-like GTPase considered unique to plants such as A. thaliana [START_REF] Ebine | A membrane trafficking pathway regulated by the plant-specific RAB GTPase ARA6[END_REF][START_REF] Ebine | Endosomal trafficking pathway regulated by ARA6, a RAB5 GTPase unique to plants[END_REF][START_REF] Ueda | Ara6, a plant-unique novel type Rab GTPase, functions in the endocytic pathway of Arabidopsis thaliana[END_REF]. When ARA6 was fused to GFP and expressed in protoplasts prepared from cultured Arabidopsis cells, both punctate and plasma membrane localisation was observed [START_REF] Ueda | Ara6, a plant-unique novel type Rab GTPase, functions in the endocytic pathway of Arabidopsis thaliana[END_REF]. Interestingly, when the N+3 cysteine was changed to serine leaving the N+2 glycine intact ARA6 became uniquely myristoylated and the GFP-fusion longer localised to endosomes, indicating that additional palmitoylation is necessary for endosome targeting [START_REF] Ueda | Ara6, a plant-unique novel type Rab GTPase, functions in the endocytic pathway of Arabidopsis thaliana[END_REF]. When the N-terminal glycine was changed to an alanine ARA6-GFP was found in the cytosol demonstrating that myristoylation is essential for membraneassociation, as we have found for PfRab5B here. In addition, the C-terminal domain of Ara6 also plays a role in determining its subcellular localisation [START_REF] Ueda | Ara6, a plant-unique novel type Rab GTPase, functions in the endocytic pathway of Arabidopsis thaliana[END_REF] and it is possible that the C-terminus of PfRab5B plays a similar role. We did not observe PfRab5B 28 -GFP associated with the food vacuole membrane (Fig. 2), unlike endogenous PfRab5B (Fig. 4), and as only the first N-terminal 28 amino acids of PfRab5B are present in the GFP-chimera this implies that the C-terminus of PfRab5B might contribute to its food vacuole location. However, unlike PfRab5B, palmitoylation of ARA6 is not required for membrane association per se, only for endosome-specific targeting. It is possible that additional targeting information is encoded in the N-terminal amino acid sequence of myristoylated and palmitoylated proteins, since the subcellular fractionation of PfRab5B-GFP and GAP45-GFP fusion proteins was not identical following replacement of the palmitoylated cysteine residue. For PfRab5B with the C3A substitution the fusion protein was found only in the soluble fraction, in contrast to GAP45 with the C5A substitution, for which the fusion protein was distributed across all subcellular fractions. The N-terminal sequence of GAP45 fused to GFP has ten basic and six acidic residues, whereas the N-terminal sequence of PfRab5B fused to GFP has only four basic and one acidic residue, so the GAP45 Nterminal sequence is much more highly charged and this greater charge might contribute to greater affinity for membrane phospholipid head groups.

ARA6 acts in an endosomal trafficking pathway by modulating the assembly of a distinct SNARE complex that includes VAMP727 [START_REF] Ebine | A membrane trafficking pathway regulated by the plant-specific RAB GTPase ARA6[END_REF]. A family of P. falciparum SNAREs has been described [START_REF] Ayong | Identification of Plasmodium falciparum family of SNAREs[END_REF] [40] and PfVAMP8 (PF3D7_1303200.1) displays the highest identity (5.6e-35) to VAMP727 suggesting that perhaps PfRab5B regulates trafficking via assembly of PfVAMP8 complex. Whether a P. falciparum PfRab5B/PfVAMP8 complex is involved in regulating a stress response, as is the ARA6/VAM727 complex in Arabidopsis [START_REF] Ebine | Endosomal trafficking pathway regulated by ARA6, a RAB5 GTPase unique to plants[END_REF] is an interesting and open question. Unlike PfRab5A we never observed PfRab5B associated with haemoglobin containing vesicles, implying that haemoglobin uptake would likely not involve the PfVAMP8 complex. Clearly, in malaria causing parasites in spite of having a small Rab family of only 11 GTPases, PfRab5A, PfRab5B and PfRab5C isoforms regulate independent and crucial functions.

Materials and Methods

Attempts to generate P. berghei rab5b gene deletion mutants

Female Swiss OF1 mice (6-8 weeks old; Charles River/Janvier) were used. All mouse experiments received approval from the Animal Experiments Committee at the Medical Centre of Leiden University (DEC 07171; DEC 10099). The Dutch Experiments on Animal Act is established under European guidelines (EU directive no. 86/609/EEC regarding the Protection of Animals used for Experimental and Other Scientific Purposes). In order to create P. berghei rab5b (PBANKA_140910) gene deletion mutants, a DNA construct pL1709 was generated that would target Pbrab5b by double crossover homologous recombination (Figure S1), as previously described [START_REF] Janse | High-efficiency transfection and drug selection of genetically transformed blood stages of the rodent malaria parasite Plasmodium berghei[END_REF]. The targeting plasmid was made to target both the 59 and 39 regions of Pbrab5b gene and contains the drug selectable marker human dhfr (hdhfr) (See Table S2 in File S1 for details of the primers). Prior to transfection, the plasmid was linearized by digestion with enzymes HindIII and EcoRI. The P. berghei ANKA reference reporter parasite line PbGFP-Luc con of the ANKA strain of P. berghei was used for two independent transfection experiments (See RMgm-29 in www.pberghei.eu for details of the PbGFP-Luc con line). All experimental procedures including animal work and transfection and selection protocols to generate gene knock-out parasites were performed as described [START_REF] Janse | High-efficiency transfection and drug selection of genetically transformed blood stages of the rodent malaria parasite Plasmodium berghei[END_REF]. Substitution of glycine at position two with alanine (PfRab5B 28 G2A) results in loss of all membrane association, as does substitution of cysteine at position three (PfRab5B 28 C3A). A known myristoylated and palmitoylated protein, GAP45 29 -GFP is found exclusively in membrane fractions and a G2A substitution results in an entirely soluble protein (GAP45 29 G2A). A substitution of cysteine at position five causes a reduction in the proportion of protein associated with membranes (GAP45 29 C5A). Double substitutions of glycine and cysteine result in a cytosolic protein for both PfRab5B 28 and GAP45 29 GFP fusion proteins. HL: hypotonic lysis buffer supernatant, HS: high salt buffer supernatant, CS: carbonate buffer supernatant, CI: carbonate buffer insoluble. doi:10.1371/journal.pone.0087695.g003

Expression constructs and purification of recombinant proteins

To make C-terminal-tagged PfRab5A-hexaHis and PfRab5B-hexaHis constructs for expression in Escherichia coli, the CDS of Pfrab5a (PF3D7_0211200) and Pfrab5b (PF3D7_1310600) were PCR amplified from P. falciparum 3D7 cDNA and cloned into BamHI/XhoI sites of the E. coli expression vector pET-21d (Novagen, EMD4Biosciences). The sequences of all expression constructs were verified.

Pfrab5a and Pfrab5b constructs were transformed into BL21-CodonPlus (DE3)-RIL strain (Stratagene). BL21 cells were grown at 37uC in LB containing 34 mg/ml chloramphenicol and 100 mg/ml ampicillin to an absorbance of approximately 0.6 at 600 nm. Protein expression was induced by adding 0.2 mM IPTG, and the cultures were incubated overnight at 20uC. Cells were harvested by centrifugation at 5000 g for 20 min and resuspended in 6 M urea buffer supplemented with protease inhibitor cocktail (Roche), then stored at 280uC. The His-tagged proteins were purified on Ni-NTA agarose (Qiagen).

To generate Pfrab5B G2A , Pfrab5B C3A and Pfrab5B G2AC3A mutations, PCR site-directed mutagenesis was performed using Pfrab5b cloned into pET-21d vector as a template. Thereafter, PCR products were incubated at 37uC with DpnI enzyme for an hour. The sequences of all expression constructs were verified.

Specific antibodies

Specific PfRab5A antibodies were raised in rabbits and specific anti-PfRab5B antibodies were raised in rats (Eurogentec) to fulllength His-tagged recombinant PfRab5A and PfRab5B (see above). Specific anti-PfCK1 antibodies were a gift of Dominique Dorin-Semblat and Christian Doerig and anti-MSP1 19 antibodies have been described [START_REF] Holder | The carboxy-terminus of merozoite surface protein 1: structure, specific antibodies and immunity to malaria[END_REF].

In vitro N-myristoylation assay

The in vitro myristoylation assay was performed essentially as described previously [START_REF] Rees-Channer | Dual acylation of the 45 kDa gliding-associated protein (GAP45) in Plasmodium falciparum merozoites[END_REF] using bacterially expressed recombinant PfNMT, recombinant substrate protein and [ 3 H]myristoyl-CoA. Recombinant PfARF and G2A PfARF were used as controls. 0.5 mg PfNMT was incubated with either 5 mg PfRab5B or 2.5 mg ARF proteins and myristoyl-CoA (0.007 MBq) for 2 h at room temperature in reaction buffer (30 mM Tris-Cl, pH 6.8, 0.5 mM EGTA, 0.45 mM EDTA, 4.5 mM 2-Mercaptoethanol, 1% Triton Immuno-electron microscopy of schizonts (A) and a large trophozoite (B) of P. falciparum labelled with specific anti-PfRab5B antibodies, revealing the presence of gold particles both on the food vacuole (FV) and the parasite plasma membrane (PPM, white arrows) in trophozoites and exclusively on parasite plasma membrane for late stages. P, parasite; r, rhoptry. Scale bars, 150 nm. doi:10.1371/journal.pone.0087695.g004 X-100). Fractions from each reaction were resolved by SDS-PAGE under reducing conditions on precast 12% NuPAGE acrylamide gels. The gels were treated for fluorography, dried and exposed to X-ray film to detect the tritium label.

Parasite culture

P. falciparum 3D7 was maintained in human O+ erythrocytes using RPMI 1640 medium supplemented with Albumax at 3% hematocrit in gassed (90% nitrogen, 5% oxygen, 5% carbon dioxide) in flasks at 37uC. Parasites were synchronized using Percoll gradient centrifugation.

Plasmid constructs and transfection of P. falciparum

To generate GFP-fusion protein expressing parasites we used primers to amplify the sequence encoding the N-terminal fragments of PfRab5B and GAP45, and inserting nucleotide changes to introduce alanine instead of glycine or cysteine as indicated in Table S1 in File S1. The PCR reactions were purified using a QIAGEN quick purification kit and digested with SacII and AvrII. The digested inserts were gel extracted then ligated with digested pHH3 vector distal to the msp3 promoter [START_REF] Ridzuan | Subcellular location, phosphorylation and assembly into the motor complex of GAP45 during Plasmodium falciparum schizont development[END_REF]. Resulting plasmids were transformed into DH5a and the sequence of each construct was established using DNA sequencing (Geneservice). 100 mg of sterile plasmid DNA in 100 mM Tris-HCl, 10 mM EDTA, pH 8.0 buffer was transfected by electroporation into parasites at 10% ring stage parasitaemia and culture was continued in the presence of blasticidin (2.5 mg/ml).

Subcellular fractionation of schizonts and western blotting

The method described previously [START_REF] Knuepfer | RON12, a novel Plasmodium-specific rhoptry neck protein important for parasite proliferation[END_REF] was used with minor modifications; all incubations were carried out at 4uC and the buffers contained complete protease inhibitor cocktail (Roche). Parasite pellets were solubilised in 10 times the volume of hypotonic lysis buffer (10 mM Tris, 5 mM EDTA, pH 8.0) followed by centrifugation at 100,000 g for 30 min at 4uC and recovery of the supernatant. The pellet was washed once by repetition of the process and subsequently extracted with high-salt buffer (10 mM Tris, 5 mM EDTA, 500 mM NaCl, pH 7.5), followed by centrifugation. The supernatant was retained and the pellet was extracted again in carbonate buffer (100 mM sodium carbonate, pH 11.0). Following a final centrifugation the three supernatants and the remaining carbonate-insoluble material were each mixed with SDS sample loading buffer and subjected to SDS-PAGE and western blotting with anti-GFP antibodies. Samples were heated at 95uC for 5 min prior to being separated by SDS-PAGE on a pre-cast 12% Bis-Tris NuPAGE polyacrylamide gel (Invitrogen) and transferred to nitrocellulose membrane according to standard protocols. Blots were blocked in a solution of 5% w/v of non-fat milk powder and 0.2% Tween20 in PBS for 1 h then probed with a mouse anti-GFP antibody (Roche) diluted in blocking solution for 1 h. Bound antibody was detected by incubation with a horseradish peroxidase conjugated anti-mouse secondary antibody (Biorad) and visualised using enhanced chemiluminescence western blotting detection reagents (GE Healthcare).

Fluorescence Microscopy

Live synchronized parasite populations expressing GFP-tagged proteins were examined by epifluorescence microscopy on an Axioplan 2 microscope (Zeiss), equipped with a Plan Apochromat 1006/1.4 oil immersion objective and an AxioCam HRc camera [START_REF] Ridzuan | Subcellular location, phosphorylation and assembly into the motor complex of GAP45 during Plasmodium falciparum schizont development[END_REF]. Parasites were labelled with DAPI DNA stain, then a cell suspension was placed on a slide and overlaid with a Vaseline rimmed cover slip. Parasitized erythrocytes were viewed live, and the dual colour fluorescence images were captured using Axiovision 4.6.3 software and edited using Adobe Photoshop CS4.

For indirect immunofluorescence smears of red blood cells infected with P. falciparum 3D7 strain were fixed for 5 min using 100% cold methanol. In the case of haemoglobin labelling, cells were pre-treated with saponin prior to fixation. Cells were washed with PBS, and then permeabilized with 0.1% Triton X100 in PBS for 5 min. After washing with PBS, slides were blocked in PBS containing 3%(w/v) BSA for 1 h at room temperature. The slides were incubated successively for 1 h with primary antibodies, depending on the co-labelling desired: mouse anti-haemoglobin (1:500, Sigma), rat anti-PfRab5B (1:200), rabbit anti-PfRab5A (1:500), mouse anti-MSP-1 19 (1:500), and rabbit anti-PfCK1 (1:500). The slides were washed four times and incubated with secondary antibodies depending on the primary antibodies used: AlexaFluor 594 anti-rat IgG antibody (1:3000, Molecular Probes), AlexaFluor 488 anti-mouse IgG antibody (1:2000, Molecular Probes), AlexaFluor 488 anti-rabbit IgG antibody (1:5000, Molecular Probes), and with DAPI (1 mg/ml). Then, the slides were washed and mounted (Dako) and examined under a microscope (Leica DMI 6000, 6100 objective, NA 1.4 oil) with a cooled charge-coupled device camera (Micromax). Z-stack images were acquired with MetaMorph (Universal Imaging) and de-convoluted with Huygens (SVI). Images were analyzed and processed with ImageJ (NIH) and Photoshop (Adobe Systems Inc). For the merge, the ImageJ co-localization plug-in was used. Pearson's coefficient of co-localization was attributed with the ImageJ JACoP plug-in [START_REF] Bolte | A guided tour into subcellular colocalization analysis in light microscopy[END_REF].

Immuno-electron microscopy P. falciparum-infected red blood cells (mixed development stages) were fixed and sectioned as described [START_REF] Tomlins | Plasmodium falciparum ATG8 implicated in both autophagy and apicoplast formation[END_REF]. The sections were immuno-labelled with anti-PfRab5B antibody at 1/10 dilution in PBS/1% fish skin gelatin, then incubated with anti-IgG antibodies, followed directly by 10 nm protein A-gold particles before examination with a Philips CM120 Electron Microscope (Eindhoven, the Netherlands) at 80 kV.

Supporting Information

File S1 Supporting Information contains Supplementary Tables S1 andS2 and legends for Supplementary Figures S1 andS2. Table S1: N-terminal sequences of PfRab5B, GAP45, and their variants used in the study. Under ''Construct name'' are listed all the different wild type and mutant GFP chimera expression plasmids used in the study. The subscript number corresponds to the number of amino acids fused to GFP. Under ''Sequence fused to GFP'' are the different wild type and mutant amino acids (underlined) sequences fused to the Nterminus of GFP. Table S2: Primers used for generation of the gene-deletion construct and for genotyping parasites. Listed in the number of each primer, its sequence with the restriction site underlined, its description and the name of the restriction enzyme. Figure S1: Specificity of anti-PfRab5A and PfRab5B antibodies. Different amounts of purified His-tagged PfRabs recombinant proteins were separated using a 15% polyacrylamide gel and transferred to a nitrocellulose membrane. (A) The blot was first incubated with the rabbit anti-PfRab5A (1:500) and then incubated with an anti-rabbit peroxidase-conjugated secondary antibody (1:15000, Sigma Aldrich). The lower panel shows protein loading by Ponceau S staining. (B) The blot was first incubated with the rat anti-PfRab5B antibody (1:1000) and then with an anti-rat peroxidase-conjugated secondary antibody (1:4000, Sigma Aldrich). The lower panel shows protein loading by Ponceau S staining. Each anti-PfRab5 antibody specifically reacted only with its corresponding recombinant protein.

Figure 1 .

 1 Figure 1. Recombinant PfRab5B is N-myristoylated by PfNMT in vitro. PfRab5B and the G2A and G2AC3A variants were incubated with [ 3 H]-myristoyl CoA in the presence of N-myristoyl transferase. Recombinant ARF and ARFG2A were treated in the same way, acting as controls. N-myristoylation was detected by the incorporation of radiolabel into the substrate following SDS-PAGE and fluorography and the upper band at 50 kDa is due to label binding to PfNMT. Both PfRab5B and ARF migrated with mobility slightly faster than the 25 kDa molecular mass marker. No incorporation was detected if the Nterminal glycine residue of either protein was replaced with alanine (G2A). doi:10.1371/journal.pone.0087695.g001

Figure 2 .

 2 Figure 2. Localisation of PfRab5B-GFP fusion proteins in P. falciparum. Parasites were transfected with constructs expressing the N-terminal 28 amino acids of PfRab5B or 29 amino acids of PfGAP45 fused to GFP under a schizont stage-specific promoter (msp3 59UTR region). Localisation of PfRab5B 28 -GFP (A) and GAP45 29 -GFP fusions (B) as well as myristoylation (G2A) and palmitoylation-deficient fusions (C3A in PfRab5B 28 -GFP and C5A in GAP45 29 -GFP) were investigated. The first image in a series corresponds to GFP fluorescence, the second a merge of GFP fluorescence with nuclear DAPI stain, and the third a merge of GFP fluorescence, DAPI and bright field images. Size bars are 2 mm. doi:10.1371/journal.pone.0087695.g002

Figure 3 .

 3 Figure 3. Subcellular fractionation of PfRab5B and GAP45 GFP chimeras. When fused to GFP, the first 28 amino acids of PfRab5B confer a partial membrane association to the fusion protein. PfRab5B 28 -GFP is found to varying degrees in the four fractions studied.Substitution of glycine at position two with alanine (PfRab5B 28 G2A) results in loss of all membrane association, as does substitution of cysteine at position three (PfRab5B 28 C3A). A known myristoylated and palmitoylated protein, GAP45 29 -GFP is found exclusively in membrane fractions and a G2A substitution results in an entirely soluble protein (GAP45 29 G2A). A substitution of cysteine at position five causes a reduction in the proportion of protein associated with membranes (GAP45 29 C5A). Double substitutions of glycine and cysteine result in a cytosolic protein for both PfRab5B 28 and GAP45 29 GFP fusion proteins. HL: hypotonic lysis buffer supernatant, HS: high salt buffer supernatant, CS: carbonate buffer supernatant, CI: carbonate buffer insoluble. doi:10.1371/journal.pone.0087695.g003

Figure 4 .

 4 Figure 4. Ultrastructural detection of PfRab5b in P. falciparum-infected RBC. Immuno-electron microscopy of schizonts (A) and a large trophozoite (B) of P. falciparum labelled with specific anti-PfRab5B antibodies, revealing the presence of gold particles both on the food vacuole (FV) and the parasite plasma membrane (PPM, white arrows) in trophozoites and exclusively on parasite plasma membrane for late stages. P, parasite; r, rhoptry. Scale bars, 150 nm. doi:10.1371/journal.pone.0087695.g004

Figure 5 .

 5 Figure 5. PfRab5B colocalises with PfMSP1 and PfCK1, but not with haemoglobin. (A) PfRab5A colocalises with haemoglobin (HBA1) containing vesicles (r = 0.724), unlike PfRab5B (r = 0.081; n = 3). (B) PfRab5B colocalises to differing degrees (r = 0.849; n = 3) with the C-terminal 19 kDa fragment of PfMSP1 on structures close to the food vacuole and the parasite nucleus shown in blue by DAPI staining. (C) PfRab5B colocalises with PfCK1 on intracellular structures (r = 0.615; n = 3) and at the parasite plasma membrane (r = 0.475; n = 3). Areas of colocalistaion are shown in white and used to calculate Pearson's r coefficients. Scale bars, 2 mm. doi:10.1371/journal.pone.0087695.g005

  Figure S2: Unsuccessful attempts to generate P. berghei rab5b (PBANKA_140910) genedeletion mutants. (A) Schematic representation of the genedeletion construct used for targeting the rab5b gene for deletion and the expected gene locus before and after disruption. The construct that has hdhfr as a drug selectable marker (SM, black) is designed to disrupt the open reading frame (ORF) of the Pbrab5b genes by double crossover homologous recombination. The expected genomic integration of the construct into the P. berghei genome is indicated, and the size and location of the Pbrab5b targeting regions (hatched boxes) are shown in relation to the Pbrab5b gene ORF. The targeting regions are indicated as +/2 bp distance from the putative start codon. The location and name of the primers used for diagnostic PCR are shown. (B) Diagnostic PCR of genomic DNA of parasites selected after transfection with gene-deletion construct pL1709 (see A) showing that rab5b ORF was not disrupted in the selected parasites. Two independent transfection experiments were performed and the parasites that survived drug selection with pyrimethamine contained both the selectable marker and the intact ORF. The following primers were used: 59 integration (59): 6909/3189; 39 integration: (3) 4592/ 6910; amplification of the hdhfr cassette (M): 307C/3187; ORF (O): 6911/6912. Genomic DNA of wild type P. berghei parasites (wt) was used as control. Two faint non-specific bands are amplified with the 39-primers. (DOCX)

PLOS ONE | www.plosone.org

February 2014 | Volume 9 | Issue 2 | e87695

Acknowledgments

We thank Dominique Dorin-Semblat and Christian Doerig for the gift of anti-PfCK1 specific antibodies. We also thank Kimberley Zichichi from the Microscopy Facility at Yale University for her technical electron microscopy competence. We gratefully acknowledge the Institute Cochin's Imaging Facility and PlasmoDB (http://plasmodb.org/plasmo/) and thank these community resources.

FBR was supported by a fellowship from the Fondation de France and CNE by an AUF fellowship. GL acknowledges INSERM, the CNRS, the Labex ParaFrap ANR-11-LABX-0024 and CEFIPRA (grant No. 3303-3) for support. The study was funded in part by the UK Medical Research Council (file reference number U117532067) and the European Union through the EviMalaR Network of Excellence Health-2009-2.3.2-1-242095. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Author Contributions

Conceived and designed the experiments: GL AAH AM CJJ. Performed the experiments: CNE FBR DKM JWL SB EK JLG SMK IC HY. Analyzed the data: GL AAH AM CJJ SMK IC HY. Wrote the paper: GL AAH. Rab5B his-tagged expression construct, CNE, FBR, HY. Rab5B Specific antibodies and imaging, CNE FBR HY. Rab5B-GFP constructs and imaging, EK JLG. NMT assays: DKM SB. Immuno-ems IC. Mouse experiments JWL SMK CJJ.


Olivier Brüls
email: o.bruls@ulg.ac.be@www.ltas-mms.ulg.ac.be

Vincent Acary
email: vincent.acary@inria.fr

Alberto Cardona
email: acardona@intec.unl.edu.ar

A PENETRATION-FREE NONSMOOTH DYNAMICS METHOD FOR FRICTIONLESS CONTACT/IMPACT PROBLEMS

Keywords: Nonsmooth Contact Dynamics, Flexible Multibody Systems, Time Integration, Generalized-α Method, Time-Stepping Schemes, Index Reduction

This paper studies numerical algorithms for the simulation of mechanical systems including rigid and flexible bodies, kinematic joints and frictionless contact conditions. The condition of impenetrability of the bodies in contact is expressed as a unilateral constraint, with the consequence that impacts and/or instantaneous changes in the velocities may arise in the dynamic response.

When dynamic contacts are analyzed between elastic solids and structures, the gap velocity is necessarily discontinuous otherwise a non-physical penetration of the bodies would occur. When contacts between rigid bodies are studied, e.g., in multibody systems, impulsive reaction forces can also occur leading to an instantaneous change in the linear and angular momenta of each body. Standard schemes, such as the Newmark, HHT or generalized-α methods, are not consistent in these cases since the numerical response may artificially generate energy when a contact occurs. In order to analyze impact phenomena, nonsmooth time integration methods have been proposed in the literature and can be classified into two main groups, namely, event-driven schemes and time-stepping schemes. Event-driven schemes are accurate for the free flight smooth motions and are especially suitable for small multi-body systems with a limited number of events, but they become inefficient if frequent transitions occur in a short time. Time-stepping methods, such as Moreau-Jean scheme, have been proven to be convergent and robust even for a large number of events, and are extensively applied as the solution to nonsmooth system models.

A fundamental property of the Moreau-Jean scheme is that the unilateral constraints are imposed at velocity level. The consequence is that some penetration can be observed in the numerical solution, which may not be physically acceptable. In order to prevent such penetration problems, this paper presents an algorithm which enforces the constraint not only at velocity level, so as to inherit good consistency and stability properties, but also at position level. For that purpose, the Gear-Gupta-Leimkuhler (GGL) approach, which was initially developed for the stabilization of index-3 differential-algebraic equations, is generalized to systems with unilateral constraints following a similar idea as in [START_REF] Acary | Projected event-capturing time-stepping schemes for nonsmooth mechanical systems with unilateral contact and Coulomb's friction[END_REF][START_REF] Schoeder | Discussion on the Gear-Gupta-Leimkuhler method for impacting mechanical systems[END_REF]. However, the present work differs from [START_REF] Acary | Projected event-capturing time-stepping schemes for nonsmooth mechanical systems with unilateral contact and Coulomb's friction[END_REF][START_REF] Schoeder | Discussion on the Gear-Gupta-Leimkuhler method for impacting mechanical systems[END_REF] by the formulation of specific complementarity conditions at position and velocity levels, as well as by the construction of an implicit algorithm which solves the dynamic equilibrium and the complementarity conditions in a monolithic way using a semi-smooth Newton process at each time step.

In order to improve the accuracy in the smooth part of motion, we also follow the strategy proposed by Chen et al [START_REF] Q.-Z. Chen | A nonsmooth generalized-α scheme for flexible multibody systems with unilateral constraints[END_REF], which involves a partitioning of the generalized forces in the dynamic equilibrium into smooth terms and nonsmooth (impulsive) terms. The generalized-α method is used to integrate the smooth terms whereas all impulsive terms are treated using an Euler implicit integration scheme to ensure consistency. This means that different integration formulae are used for the different contributions to the equation of motion. Nevertheless, the advantage of this approach is that the numerical dissipation of the generalized-α method is significantly smaller than the numerical dissipation of the Euler implicit scheme, so that the energy behavior is strongly improved compared to the Moreau-Jean scheme, especially for mechanical systems exhibiting both impacts and structural vibrations. In some sense, in this scheme, the high numerical dissipation of the Euler implicit scheme is only acting locally at the contact region and during the contact time, and not on the full system during the whole trajectory.

The numerical behaviour of the proposed method is studied and compared to other approaches for a number of numerical examples. It is shown that the formulation offers a unified and valid approach for the description of dynamic contact conditions between rigid bodies as well as between flexible bodies.