L Nehaoua
email: lamri.nehaoua@ibisc.univ-evry.fr

Dalil Ichalal

H Arioui

S Mammar

L Fridman
email: lfridman@servidor.unam.mx

H L Nehaoua

Lean and Steering Motorcycle Dynamics Reconstruction : An Unknown Input HOSMO Approach

This paper deals with state estimation of Powered Two Wheeled (PTW) vehicle and robust reconstruction of related unknown inputs. For this purpose, we consider a unknown input high order sliding mode observer (UIHOSMO). First, a motorcycle dynamic model is derived using JourdainŠs principle. In a second time, we consider both the observation of the PTW dynamic states, the reconstruction of the lean dynamics (roll angle φ(t)) and the rider's torque applied on the handlebar. Finlay, several simulation cases are provided to illustrate the efficiency of the proposed observer.

I. INTRODUCTION

.

 to detect borderline cases of loss-of-control. These research are very few sustainable if they are not propped by a system estimating the dynamic states of the PTW.

, helping in states observation of motorcycle model and the reconstruction of rider's action.

Sliding mode based observer (SMO) for system states estimation in the presence and/or absence of unknown inputs has been the subject of several work from the control community. Today, one can notice that observation theory has matured and has succeed to deal with many technical issues where some restrictive conditions related to the observability and the reconstruction of unknown inputs were released even suppressed. In its simplest version, a SMO uses the same structure as a traditional Luenberger observer with an injection term related to the system's output. However, it is necessary that the system's unknown inputs and outputs satisfy the so-called matching condition [START_REF] Edwards | On development and applications of sliding mode observers[END_REF]. To overcome this limitation, system coordinates transformation are introduced and the use of sliding mode differentiator (SMD) for the auxiliary output generation is generalized [START_REF] Levant | Robust Exact Differentiation via Sliding Mode Technique[END_REF], [START_REF] Floquet | A canonical form for the design of unknown input sliding mode observers[END_REF]. This allows to use a powerful optimization technique, like linear matrix inequality (LMI), to offer a systematic design procedure of the observer gain [START_REF] Tan | An LMI Approach for Designing Sliding Mode Observers[END_REF]. Nevertheless, the chattering of these observers requires filtering [START_REF] Fridman | Observation and Identification Via High-Order Sliding Modes[END_REF].

To avoid filtering, the discontinuous output injection is replaced by a continuous super-twisting algorithm (STA) [START_REF] Davila | Second-order sliding-mode observer for mechanical systems[END_REF]. In this new version, the relative degree of the system's outputs with respect to the unknown inputs must be equal to the system order. This restriction is resolved by the introduction of the high-order sliding mode observers (HOSMO) [START_REF] Fridman | Observation of linear systems with unknown inputs via high-order sliding-modes[END_REF] based on the high-order robust exact sliding mode differentiator [START_REF] Levant | Higher-order sliding modes, differentiation and outputfeedback control[END_REF] where the notion of strong observability and strong detectability were presented [START_REF] Fridman | Observation of linear systems with unknown inputs via high-order sliding-modes[END_REF]. It remains at least that the outputs relative degree must exist which brings a novel restriction treated by the development of the concept of weakly observable subspaces detailed in [START_REF] Fridman | Higher-order sliding modes observation for linear systems with unknown inputs[END_REF].

This paper is organized as follows: a succinct problem statement is provided in section II. Sections III and IV are dedicated respectively to motorcycle modelling and the synthesis of the UIHOSMO. Simulation results are given is section V. The conclusion in section VI wrap up the paper.

II. PROBLEM STATEMENT

Our study concerns the identification of all relevant inputs and dynamic states helping in a next stage to quantify the risk of loss-of-control when cornering. Indeed, poor cornering is responsible for most PTW accidents (singlevehicle motorcycle crashes).

To perform a safe cornering, riders should consider : 1) the appropriate speed when starting the corner, 2) the road condition (under weak friction) and 3) weather conditions do not allow optimal visibility for driving.

Early warning systems are based generally on related work carried out for standard cars [START_REF] Slimi | Advanced Motorcycle-Infrastructure-Driver Roll Angle Profile for Loss Control Prevention[END_REF]. The goal is the synthesis of a function giving the maximum safe speed at which a vehicle can be kept on the road while moving at a constant speed on a circular section. This speed depends, among other factor, on the lateral friction which its computing involves all the dynamic states of the bike and a good modelling of the tireroad contact. This makes the success of such warning system strongly dependent on the availability of dynamic states of the motorcycle. This challenge constitutes the contribution of our ongoing work.

III. MOTORCYCLE DYNAMICS

A. Modeling Assumptions

The study of the dynamics of motorcycle vehicles highlights two main modes of motion: in-plane mode representing the motorcycle movements in its plane of symmetry including the longitudinal motion and that of suspensions and the out-of-plane mode which describes the lateral dynamics when cornering [START_REF] Sharp | The stability and control of motorcycles[END_REF], [START_REF] Sharp | A Motorcycle Model for Stability and Control Analysis[END_REF]. The last mode involves the roll inclination, the yaw rotation, the steering and the lateral motions of the bike. We consider here only the out-ofmode dynamics of the two-wheeled vehicle. The coupling between the two modes is materialized, when neccesery, by considering a variable longitudinal velocity that appears in the lateral dynamics. In this study, the motorcycle is represented as two linked bodies: the rear one G r which includes the chassis, engine and the rear wheel and the front body G f which represents the steering assembly and the front wheel. In that case, we recover the well known two bodies Sharp's model which allows to simulate 4 degrees of freedom (DOF): the lateral displacement v y , the yaw rotation ψ, the roll inclination ϕ and the handlebar steer angle δ with respect to the rider torque input τ r applied on the motorcycle's handlebar (Fig. 1).

B. Motorcycle Motion Equations

Let ℜ v be a reference frame attached to the motorcycle at point v in which the motion of the overall mechanical system will be expressed. With respect to the inertial reference frame ℜ o , ℜ v is referred by the velocity vector components v x , v y of point v, and the yaw rotation ψ around the ℜ o z-axis (Fig. 2). In ℜ v , the linear velocity vector of the two bodies are respectively defined by:

c r • v • c f • O • i o j o ℜ o ℜ v i v j v ψ
v oGr = v ov + ω oGr × r vGr (1) v oG f = v ov + ω oGr × r vB + ω oG f × r BG f
where ω oGr and ω oG f are the angular velocity vector of each body, expressed in ℜ v by:

ω oGr = ψk + φi (2)
ω oG f = ω oGr + δk δ
By differentiating (1) and ( 2), the linear and the angular acceleration of G r and G f can be written for i = r, f as following:

a oGi = ∂v oGi ∂ϑ θ + a R,Gi (3)
ǫ oGi = ∂ω oGi ∂ϑ θ + ǫ R,Gi
where ϑ = [v y , ψ, φ, δ] T denotes the vector of generalized velocities. From the Jourdain's principle of dynamics [START_REF] Rill | Simulation von Kraft-fahrzeugen[END_REF], the dynamic model of the motorcycle vehicle can be expressed by:

M θ = Q (4)
where the mass matrix M is symmetric and positive definite and obtained directly from the Jacobian matrices by:

M = i=r,f m i ∂v T oGi ∂ϑ ∂v oGi ∂ϑ + ∂ω T oGi ∂ϑ I i ∂ω oGi ∂ϑ ( 5
)
and Q is the vector of the generalized efforts, given by:

Q = Q δ + Q g + Q T + i=r,f ∂v T oGi ∂ϑ m i (g -a R,Gi ) - i=r,f ∂ω T oGi ∂ϑ (I i ǫ R,Gi + ω oGi × I i ω oGi ) (6)
In ( 4) and ( 5), I i is the inertia tensor matrix, Q δ includes the effect of the rider's torque and also the handlebar damping, Q g represents the gyroscopic effect resulting from the tires' spin, Q T includes the effect of the tire's sideslip force and g is the gravity force vector.

C. Tire/Road Interaction and Wheels Gyroscopic Effect

To describe the tire motion, a new reference frame ℜ T (c, i T , j T , k T ) is introduced at the contact point c of each wheel's tire. k T is the normal vector to the road surface along the F z force.

R • j ζ γ i T j T F y F z v oc α c • Fig. 3. Tire reference frame
The linear velocity vector at point location c i is expressed in ℜ v by:

v ocr = v ov + ω oGr × r vcr (7)
v oc f = v ov + ω oGr × r vc f + δe δ × r cηc f
From ( 7), the tires sideslip angles are defined by:

α r = -atan j v • v ocr i v • v ocr (8)
α f = -atan j v • v oc f i v • v oc f + δ cos ǫ
The equivalent tire effort at the center of each wheel is governed by:

F T = F y j T + F z k T (9)
where F y = F y (α, γ) is obtained from the sideslip angle α and the camber angle γ computed by sin γ = j ζ • k T . Consequently, the contribution Q T of the tire/road contact forces in the vector of the generalized efforts Q is given by:

Q T = i=f,r ∂v oci ∂ϑ T F T,i (10)
To compute the contribution of the gyroscopic effect Q g , let first write the spin velocity vector ω s,i equation of each tire in ℜ v reference frame:

ω osr = ω oGr + ξr j v (11)
ω os f = ω oGr + δk δ + ξf j ξ From (3), the contribution of the gyroscopic effect in the vector of the generalized efforts Q is :

Q g = - i=r,f i iy ∂ω T si ∂ϑ ǫ R,si (12)

D. Linearized Model

The motorcycle dynamic model ( 4) is linearized around the straight-running trim trajectory and can be expressed by the following state-space:

ẋv = A v x v + B v τ r (13)
Here, x v = [v y , ψ, φ, δ, ϕ, δ, F yr , F yf ] T denotes the state vector. A v is a time-varying matrix related to the forward velocity v x while B v is a time-invariant vector. F yr and F yf represent respectively the tires sideslip forces introduced in the state space representing the tire relaxation.

IV. STATES AND UNKNOWN INPUTS ESTIMATION

In this section, we aim to estimate the motorcycle states and reconstruct both roll angle ϕ and rider's torque τ r by using a HOSMO [START_REF] Fridman | Observation of linear systems with unknown inputs via high-order sliding-modes[END_REF].

At first, we recall some important definitions about observability and detectability of linear systems (for proofs see [START_REF] Fridman | Observation of linear systems with unknown inputs via high-order sliding-modes[END_REF], [START_REF] Isidori | Nonlinear Control Systems[END_REF]). Consider the following SISO system, where x ∈ R n and ζ ∈ R is the unknown input:

ẋ = Ax + Bu + Dζ (14) y = Cx Definition 4.1:
In the absence of unknown input (ζ = 0), system ( 14) is observable if and only if the observability matrix P such that:

P =     C . . . CA n-1     (15)
has the full rank. Otherwise, it is detectable if the system's invariant zeros are stable.

Definition 4.2: The relative degree of the output y with respect to the unknown input ζ is the number r such that:

CA j D = 0 j = 1, • • • , r -2 CA r-1 D = 0 Definition 4.3:
In the presence of an unknown input, system ( 14) is strongly observable if and only if the relative degree r satisfies: r = rank(P ). Otherwise, it is strongly detectable if and only if the relative degree r exists and ( 14) is minimum phase system. In that case, r ≤ rank(P ).

As previous, consider the MIMO system of the form [START_REF] Fridman | Higher-order sliding modes observation for linear systems with unknown inputs[END_REF], where x ∈ R n , y ∈ R m is the output vector and ζ ∈ R m is the unknown input vector.

Definition 4.4: In the absence of unknown inputs (ζ = 0), system ( 14) is observable if and only if the observability matrix P such that:

P =                 C 1 . . . C 1 A n-1 . . . C m . . . C m A n-1                 (16)
has the full rank, where

C i , i = 1, • • • , m is the i-th row of the matrix C.
Definition 4.5: The relative degree of the output y with respect to the unknown input ζ is the vector r = [r 1 , • • • , r m ] such that:

C i A s D j = 0 i, j = 1, • • • , m, s = 1, • • • , r i -2 C i A ri-1 D j = 0 and: det     C 1 A r1-1 D 1 • • • C 1 A r1-1 D m . . . . . . . . . C m A rm-1 D 1 • • • C m A rm-1 D m     = 0 Definition 4.6:
In the presence of an unknown input, system ( 14) is strongly observable if and only if the total relative degree r T = r 1 + • • • + r m satisfies: r T = rank(P ). Otherwise, it is strongly detectable if and only if the relative degree vector r = [r 1 , • • • , r m ] exists, and ( 14) is minimum phase.

From definition (4.1), the motorcycle dynamics ( 13) is neither observable nor detectable. Indeed, for all v x in the allowable velocities range, the observability index is equal to 6 which is less than the system order (n = 8), in addition, the motorcycle dynamics has one unstable invariant zero which makes the motorcycle dynamics to be a non-minimum phase.

A. Estimation of the Roll Angle and Rider's Torque

In order to make the system observable, the motorcycle model ( 13) is rewritten such that the roll angle ϕ and the steering torque τ r are considered as unknown inputs. In fact, the unstable invariant zero is a direct consequence of the counter-steering phenomena generated by the motorcycle roll. In this case, the new system equation is written as:

ẋp = A p x p + D pp ζ (17)
y pp = δ ψ = C pp x p D pp = D p B p ζ T = ϕ τ r
where x p = [v y , ψ, φ, δ, δ, F yr , F yf ] T denotes the state vector and matrices A p , B p and D p are given by:

B p =              0.016 -0.006 -0.019 0.008 0 0 0              D p =              -26
            
From definition (4.5), the output y pp in (17) has a relative degree vector r = [2, 1] with respect to the unknown input vector In addition, system (17) has 3 stable invariant zeros for all v x in the allowable velocities range. It results from (4.6), that ( 17) is also strongly detectable. This definition implies that only r T = r 1 + r 2 system's states can be estimated exactly while the observation of the remaining states are asymptotically exact.

As before, to estimate the system's state x p and the unknown input vector ζ, it is necessary to separate the clean states from those contaminated by the unknown inputs. To achieve this, system ( 17) is transformed to a new coordinates system ξ p = T x p that the closed-loop system dynamics

(A p -C pp L pp , D pp , C pp ) is expressed by: ξ11 = ξ 12 (18
where ξ T p = [ξ T 1 , ξ T 2 ] ∈ R n and ξ T 1 = [ξ 11 , ξ 12 , ξ 12 ] ∈ R rT and ξ T 2 = [ξ 21 , • • • , ξ 24 ] ∈ R n-rT .
Next, the state observer is designed as:

ż = A p z + L pp (y pp -C pp z) (19) θ2 = A 21 ϑ 1 + A 22 ϑ 2 xp = z + T -1 ϑ
in which xp is the vector of estimated states, z ∈ R n and ϑ ∈ R n is given by (recall that n = 7, r = [2, 1] and r T = 3):

ϑ = ϑ 1 ϑ 2 R rT R n-rT =      v 1,1 v 1,2 v 2,1 ϑ 2      (20)
Once again, v i ∈ R rM +k+1 is the nonlinear part of the observer where i = 1, • • • , m and r M = max(r i ).

           
is a solution of the discontinuous vector differential equation by considering k = 1 as:

vi,1 = -3λ

1 4 |v i,1 -y p + C p z| 3 4 sign(v i,1 -y p + C p z) + v i,2
vi,2 = -2λ

1 3 |v i,2 -vi,1 | 2 3 sign(v i,3 -vi,1 ) + v i,3 (21)
vi,3 = -1.5λ

1 2 |v i,3 -v2 | 1 2 sign(v i,3 -vi,2 ) + v i,4 vi,4 = -1.1λsign(v i,4 -vi,3 )
Finally, the reconstruction of the unknown input vector is possible by using:

φ τr = D-1 v 1,3 v 2,2 - a 11,11
Remark 4.1: In the strong observability case, all system's states are exactly estimated and the coordinates transformation matrix ξ = T x is none than the observability matrix P . However, in the strong detectability case, only r system's states are exactly estimated where r is the relative degree. Consequently, it is simple to show that the r first lines of the coordinates transformation matrix T are the same of those of the observability matrix P . The expression of the remaining nr lines are explicitly described and proven in [START_REF] Ilchmann | Tracking with prescribed transient behaviour for nonlinear systems of known relative degree[END_REF].

V. SIMULATION RESULTS

In this section, the unknown input high order sliding mode observer is constructed for the presented motorcycle model. Some simulations and discussions are provided to illustrate the effectiveness and the ability of the UIHOSMO in estimating simultaneously the dynamic states and both roll angle and the applied torque by the rider on the handlebar. The observer is designed in such a way to estimate all the dynamic states and unknown inputs from only the knowledge of steering angle δ(t) and the yaw rate ψ(t). The parameters λ i , i = 1, 2 of the differentiator [START_REF] Slimi | Advanced Motorcycle-Infrastructure-Driver Roll Angle Profile for Loss Control Prevention[END_REF] are chosen as follows λ 1 = λ 2 = 10000. The Luenberger gain L pp of ( 25) is computed by pole placement at: -15, -30, -45, -60, -150, -165, -180, and the obtained matrix is:

L pp =              -254
            
The change of coordinates matrix T is given in equation ( 23). With these parameters and matrices, the UIHOSMO is implemented with initial conditions x(0) = [0.1745 0.1745 -0.2618 -0.0873 50 50 0] and φ(0) = 0. These initial conditions correspond to a early cornering case with longitudinal speed at v x (t) = 15m/s.

A. Initial Conditions

A first simulation is performed with the same initial conditions x(0) = x(0). For different initial conditions, let us consider the initial conditions of the system as x(0) = [0 0 0 0 0 0 0.1] and φ(0) = 0. The simulation results concerning the dynamics states and their estimations are depicted in figures 5 (top), 6 and 4. The lateral forces on each wheel (rear and front) are illustrated in figure 7. Finally, the unknown inputs (roll angle φ(t) and rider torque τ (t) are shown in figures 5 (bottom) and 8. One can conclude that the UIHOSMO provides satisfactory results. In this work, only state and unknown input estimation is considered, but it is possible to deal with the problem of oscillations by parameter adjustment of the observer (λ).

T =             0 0 0 0 1 0 0 0 -51, 4 0 1 -208, 44 0 0 0 1 0 0 0 0 0 0 6, 28 -1 0, 14 -4, 08 0 0 -1 -3, 82 0 -0.08 4, 3 0 0 0 -1.14 × 10 4 0 -275, 19 1.3 × 10 3 1 0 0 -8, 13 × 10 3 0 -194, 85 -4, 3 × 10 3 0 1             (23)

B. Simulation With Uncertain Parameters

In order to illustrate the performances of the proposed observer in the presence of modeling uncertainties, the observer is designed by using the motorcycle nominal model and applied to the system with 20% parameters variation with respect to nominal values. As shown in figures 9, 10 and 11 the observer provides accurate estimations for almost the state variables. The effect of the uncertainties can be seen only in the estimation of the lateral forces and rider torque. Nevertheless, the obtained estimations are acceptable. As To end the simulation part, let us consider the abilities of the proposed observer in noisy measurement case. A centered random noise in the range [-0.0001, 0.0001] is added to the measurements δ(t) and ψ(t). The obtained results are depicted in the figure 14. The states are correctly estimated expect the steering angle rate δ(t). This problem can be solved by reducing the value of the parameters λ i , i = 1, 2, for example, with the values λ 1 = λ 2 = 10, the obtained

0

VI. CONCLUSION

In this article, the problem of observer design for simultaneously estimating the dynamics states of a motorcycle, an unknown inputs (rider's action) and the lateral forces on each wheel (front and rear) is considered. For that purpose, an Unknown Inputs High Order Sliding Mode observer is proposed. A motorcycle model similar to the well-known Sharp model is derived using JourdainŠs principle. The observability of the initial model is not guaranteed, then a transformation of this last into a model with two inputs by considering the unobservable state φ(t) as an unknown input as well as the rider torque applied on the handlebar. The obtained model is then exploited to construct the Unknown Inputs High Order sliding mode observer. Some simulation results are provided in order to illustrate the efficiency of the High Order Sliding Mode Observer.

Fig. 1 .

 1 Fig. 1. Geometrical representation of the Sharp's motorcycle model

Fig. 2 .

 2 Fig. 2. Configuration of system axis ℜv w.r.t the inertial frame ℜo

Fig. 4 .Fig. 5 .

 45 Fig. 4. Actual state (solid blue line) and estimated state (dashed red line)

Fig. 6 .Fig. 7 .

 67 Fig.6. Actual state (solid blue line) and estimated states (dashed red line)

Fig. 8 .

 8 Fig. 8. Actual states (solid blue line) and estimated states (dashed red line)

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Actual states (solid blue line) and estimated states (dashed red line) in the presence of 20% variations in the parameters

Fig. 12 .Fig. 13 .Fig. 14 .Fig. 15 .

 12131415 Fig. 12. Longitudinal velocity with respect to time

a ij a 11 =

 11 m f + mr a 12 = m f k a 13 = mrh + m f j a 14 = m f e a 22 = m f k 2 + Irz + I f z cos 2 ǫ + I f x sin 2 ǫ a 23 = m f kj -Crxz + (I f z -I f x ) cos ǫ sin ǫ a 24 = m f ek + I f z cos ǫ a 33 = mrh 2 + m f j 2 + Irx + I f x cos 2 ǫ + I f z sin 2 ǫ a 34 = m f ej + I f z sin ǫ a 44 = m f e 2 + I f z parameters b ij b 12 = (m f + mr)vx b 22 = m f kvx b 23 = -(i f y /R f + i Ry /Rr)vx b 24 = -i f y /R f sin ǫvx b 32 = (m f j + mrh + i f y /R f + i Ry /Rr)vx b 34 = i f y /R f cos ǫvx b 35 = -(m f j + mrh)g b 36 = Z f ηm f eg b 42 = (m f e + i f y /R f sin ǫ)vx b 43 = -i f y /R f cos ǫvx b 45 = Z f ηm f eg b 46 = (Z f ηm f eg) sin ǫ

  ) ξ12 = a 11,11 ξ 11 + a 11,12 ξ 12 + a 11,13 ξ 13 + a 12,11 ξ 21 + • • • + a 12,14 ξ 24 + d 11 ϕ + d 12 τ r ξ13 = a 11,21 ξ 11 + a 11,22 ξ 12 + a 11,23 ξ 13 + a 12,21 ξ 21 + • • • + a 12,24 ξ 24 + d 21 ϕ + d 22 τ r ξ2 = A 21 ξ 1 + A 22 ξ 2

y T pp,new = ξ 11 ξ 13

  Each unknown input ζ i is bounded with |ζ i | ≤ ζ i,max and the (r M -r i +k) successive derivatives of ζ i are bounded by the same constant ζ ′ i,max , consequently the auxiliary variable v i

	Ap =	           	0 0 0 0 0 -64934.37 -45832.65	-0.899vx -0.015vx 0.028vx -0.163vx -2.373vx 0 31155.51 -42841.04	0.033vx 0.033vx -0.009vx 4.516vx 0 0 0	-0.057 0.016vx + 0.374 0.109 -0.059vx -0.025vx -15.611 1 0 5307.421	-1.862 -2.443 2.949 127.454 0 0 42580.46vx	0.018 -0.018 -0.021 0.035 0 -4.102vx 0	0.009 0.024 -0.012 -0.309 0 0 -4.101vx

  Av state matrix Av = M -1 E Bv input matrix Bv = M -1 [0, 0, 0, 1, 0, 0, 0]T , mr front and rear mass 30.65, 217.45 [kg] Z f front tire normal force -1005.3 [N] I f x , I f z front body inertia 1.23, 0.44 [kg.m 2 ] Irx, Irz, Crxz rear body inertia 31.18, 21.07, 1.73 [kg.m 2 ] i f y , iry wheels spin inertia 0.718, 1.051 [kg.m 2 ] R f , Rr wheels radius 0.304, Rr = 0.304 [m] σ f , σr tire's relaxation [m] a, b, h, e, f 0.948, 0.479, 0.615, 0.024, 0.028 [m]

		VII. APPENDIX
	motorcycle
	vx, vy:	longitudinal and lateral velocity
	ϕ, ψ, δ:	roll, yaw and steer rotations
	τr:	rider torque
	F yf , Fyr:	lateral force
	M :	motorcycle mass matrix
	E:	motorcycle generalized effort vector
	notations
	ẋ, ẍ:	derivatives of a variable x w.r.t time
	x:	estimate of a variable x
	x T :	transpose of vector or matrix x
	x f , xr	denotes front and rear
	motorcycle [15]
	m f η	pneumatic trail 0.116 [m]
	ǫ	castor angle 0.4715 [rad]
	C δ	handlebar damping [N.m -1 .s]
	g	gravity force 9.81 [N]
	l	vehicle wheelbase [m]


Deepesh Agarwal

Frédéric Cazals

Noël Malod-Dognin

Malod - Dognin

Stoichiometry Determination for Mass-spectrometry Data: the Interval Case

Keywords: Protein complex, mass spectrometry, denumerant, dynamic programming, linear diophantine equation, enumeration, output sensitive algorithm Complexes protéiques, spectrométrie de masse, dénumérant, programmation dynamique, équation diophantienne, énumération, algorithmes sensibles à la sortie

In structural proteomics, given the individual masses of a set of protein types and the exact mass of a protein complex, the exact stoichiometry determination problem (SD), also known as the money-change problem, consists of enumerating all the stoichiometries of these types which allow to recover the target mass. If the target mass suers from experimental uncertainties, the interval SD problem consists of nding all the stoichiometry vectors compatible with a target mass within an interval.

We make contributions in two directions. From an algorithmic standpoint, we make two modications to exact SD algorithms, to inherently address the interval SD problem. The rst modication concerns the classical tree-like exploration, resulting in algorithm DIOPHANTINE, which solves the interval SD problem using constant-memory space. The second modication concerns the classical dynamic programming based algorithm, resulting in algorithm DP++, which solves the interval SD problem in an output sensitive manner. From an applied perspective, we raise three points.

First, we show that DIOPHANTINE and DP++ yield an improvement from 3 to 4 orders of magnitude over state-of-the-art exact SD algorithms, for typical protein complexes facing uncertainties on the target mass in the range 0.1-1%. It is shown that this gain comes from the ability of DP++ and DIOPHANTINE to avoid redundant calculations occurring when solving the interval SD problem as a sequence of exact SD problems. Second, we show that DIOPHANTINE behaves like an outputsensitive algorithmespecially when the interval width increases, albeit such a property cannot be expected in general. Third, from a biological perspective, using a panel of biological complexes (eukaryotic translation factor, yeast exosome, 19S proteasome sub-unit, nuclear pore complex), we stress the importance of enumeration, even at a null noise level. We also propose a representation of all solutions of an interval SD problem using a directed acyclic graph stressing the prominent protein types.

The programs accompanying this paper are available from http://team.inria.fr/abs/

Mass spectrometry. Biology rests on multimeric assemblies, so that reconstructing models of such assemblies is a central problem in (computational) structural biology. While the structure determination of proteins and small complexes has increased at a fast pace thanks to structural proteomics projects [START_REF] Janin | structural genomics: winning the second half of the game[END_REF], that of larger complexes, which typically do not crystallize or are too big for nuclear magnetic resonance studies has remained elusive. Fortunately, mass spectrometry (MS), a technique initially developed in chemistry to infer the mass-to-charge ratio of molecules, is in the process of lling this gap, allowing to study complexes that vary in mass, size, solubility, and bound/unbound states [START_REF] Sharon | The role of mass spectrometry in structure elucidation of dynamic protein complexes[END_REF][START_REF] Benesch | Mass spectrometry: come of age for structural and dynamical biology[END_REF]. In a nutshell, MS rst requires isolating sample complexes from the cellular environment, a process typically resorting to Tandem Anity Purication [ea01]. Operating on these samples in solution, a soft ionization technique such as electrospray ionization or matrix-assisted laser desorption/ionization allows ionizing the proteins or protein assemblies, while preserving non covalent interactionsa key feature to study proteins and complexes [START_REF] Sharon | The role of mass spectrometry in structure elucidation of dynamic protein complexes[END_REF]. The ions produced are then passed to a mass selection instrument (time of ight analyzer, ions traps, orbitraps, . . . ), recording information based on the ration mass-to-charge of the ions.

In classical time-of-ight MS, these primary ions are sorted by a mass analyzer according to their mass over charge (m/z) ratio, and the relative abundance of the samples are measured by a detector. Unfortunately, for a solution containing a mixture of complexes, dierent complexes with identical m/z ratio may contribute to the same peak in the mass spectrum. To disentangle such cases, tandem mass spectrometry (tandem MS) can be used. In this rened protocol, the precursor ions undergo a collision induced dissociation, which consists of removing highly charged peripheral sub-units thanks to gas collisions, so that the product ions can indeed be sorted by m/z ratio.

MS can also be used to gain information on the sub-complexes of an assembly [THS + 08].

Denaturating solvent or high salt concentrations can be used to break the complex into pieces, which are then identied by MS or tandem MS. In extreme conditions, complete denaturation occurs, so that the masses of the individual components can be determined. In milder conditions, multiple overlapping complexes are generated, providing information on the building blocks.

Stoichiometry determination. Upon determining the masses (those of the complexes in a mixture, or those of the sub-complexes of an assembly), the stoichiometry determination (SD) problem consists of determining the stoichiometry of the protein types accounting for each target mass. Denoting p the number of protein types, w i and s i the mass and the stoichiometry of the ith protein type respectively, one wishes to reconstruct the target mass M as follows:

i=1,...,p s i w i = M.

(

) 1
This problem, which was rst addressed in [START_REF] Dromey | Calculation of elemental compositions from high resolution mass spectral data[END_REF] using an algorithm exploring all possible stoichiometry vectors in a tree-like fashion, is also known as the money exchange problem, for obvious reasons [START_REF] Böcker | The money changing problem revisited: Computing the frobenius number in time o(ka1) *[END_REF]. Mathematically and since the masses are integers (expressed say in Daltons 1 ), SD is related to the theory of integer partitions [START_REF] Comtet | Advanced Combinatorics: The art of nite and innite expansions[END_REF], to linear diophantine equations [START_REF] Smart | The algorithmic resolution of Diophantine equations[END_REF], and is also coupled to so-called knapsack and subset sum problems [START_REF] Pisinger | Where are the hard knapsack problems?[END_REF]. 1 The unied atomic mass unit, also known as Dalton, is approximately equal to the mass of one neutron or one proton.

Inria

Stoichiometry Determination for Mass-spectrometry Data 5 Uncertainties in the mass determination. Practically, two diculties are faced in mass spectrometry. First, because the proteins involved in dierent copies of a complex may have undergone dierent post-translational modications 2 [START_REF] Kelleher | Peer reviewed: Top-down proteomics[END_REF], there might be a discrepancy between the masses of the sub-units. Second, any mass experimentally determined is generally larger than the theoretical one, due to extra molecules (solvent and electro-spray buer molecules) sticking to the molecule or complex analyzed. The magnitude of this second source of error actually depends on the type of ionization technique used and on the type of sample.

Mass analyzers such as orbitraps, which were initially designed to handle small molecules, typically achieve relative mass error of 0.005% (50 part-per-million or ppm). Interestingly, it has been reported recently that upon modifying the instruments, in particular the ability to completely desolvate the samples [RDD + 12], intact macromolecular assemblies with molecular weights of the order 1 MDa were amenable to studies with the same precision. Yet, this tour de force has been achieved so far on very few macro-molecular complexes as follows, IgG antibody (146 kDa), bacteriophage HK97 (253 kDa), yeast 20S proteasome (730 kDa), E. coli GroEL (801 kDa).

On the other hand, classical Time of Flight (TOF) analyzers accommodate a much wider range of samples, and continue to be used for characterization of large macromolecular assemblies [START_REF] Zhang | Native mass spectrometry of photosynthetic pigment-protein complexes[END_REF]. Using such instruments, a relative mass error as high as ∼ 1% are typical for heterooligomeric complexes of up to 1 MDa [THS + 08, MR12].

The interval problem. To account for the uncertainties on the masses just discussed, let ε be a positive integer corresponding to the mass discrepancy to be accommodated. For the sake of conciseness, the relative mass error r is called the noise level in the sequel, the number ε being dened by ε = r × M . Instead of solving Eq. (1), we are instead interested in the following enumeration problem Report all stoichiometry vectors S = (s 1 , . . . , s p ) such that | i=1,...,p

s i w i -M |≤ ε. (2)
For the sake of conciseness, the previous problem uses symmetric errors terms ±ε with respect to the target mass M . However, an interval [M -ε 1 , M + ε 2 ] with ε 1 = ε 2 can be accommodated at the cost of trivial modications of our algorithms.

Contributions

To the best of our knowledge, the question of mass accuracy for the complex composition to be determined has only been touched upon indirectly. In [BLM + 08], the interval SD problem is solved by repeatedly calling the exact SD algorithm on each target mass in the interval. In [THS + 08], 100 hypothetical complexes from 6 to 14 sub-units have been created with masses in the range 10-50kD, and 100 dimers, trimers or tetramers were generated from these complexes.

With an error rate of 1%, it is observed that no tetramer has a unique stoichiometry vector.

In this context, we make two contributions. First, we present a constant memory space algorithm (DIOPHANTINE), and an output sensitive dynamic programming based algorithm (DP++), both inherently addressing the interval SD problem. For DIOPHANTINE, we also show that sorting the masses w i yields a tree of optimal size, an observation which had not been made previously.

Second, we present a detailed experimental study on various biological and synthetic datasets, which results in three conclusions. The rst observation is that for biological datasets, enumeration does matter even for a moderate (and sometimes even null) noise level. The second 2 The chemical modication of a protein after its bio-synthesis by the ribosome.

RR n°8101

nding is that our algorithms DIOPHANTINE and DP++ outperform state-of-the art dynamic programming based approaches by three to four orders of magnitude, for a noise level in the range 0.1% to 1%, which is typically faced in structural proteomics. Not surprisingly, we show that this improvement owes to the ability of our algorithms to avoid redundant calculations which are carried out when calling the exact dynamic programming based algorithm on each target mass in an interval. The last one is that DIOPHANTINE actually exhibits an output sensitive behavior. Thus and interestingly, we show that one of the very rst algorithms designed to solve the exact SD problem [DF80] can be modied not only to solve the interval SD problem, but also to outperform advanced strategies based on dynamic programming.

Theory and Algorithms

In section 2.1, we recall the rich background of knapsack and integer partition problems. In section 2.2, we present classical counting results, and derive bound on the number of solutions for the SD problem.

Denumerants, Unbounded Knapsack and Subset-sum Problems

In the following, we consider a vector W = {w 1 , . . . , w p }, of positive integers or real numbers. A stoichiometry vector is denoted S = {s 1 , . . . , s p }. The vector is called positive if s i > 0 for all i, and non-negative if s i ≥ 0 for all i.

Denumerants. Assume that the vector W contains integers. In combinatorics, the number of non-negative integer solutions to Eq. ( 2) is known as the denumerant D(M ) of the target mass. It has been known since Bell [Bel43] that if lcm(W) denotes the least common multiple of the w i s, then, for each b ∈ {0, . . . , lcm(W) -1}, if M = m lcm(W) + b, the denumerant is a polynomial in m of degree p -1, that is:

D(m lcm(W) + b) = i=0,...,k-1 c i m i . (3)
More generally, the number of solutions reads from the power series expansion of the generating function 1/Π i (1 -X wi ), but the diculty precisely relies in extracting the coecient of X n in that expansion. In analytic combinatorics [FS09, Chapter IV], a classical result obtained by singularity analysis states that if gcd(W) = 1, one asymptotically gets

D(M ) ∼ M p-1 (p -1)!w 1 . . . w p . (4)
Upper and lower bounds on the denumerant have also been obtained based on binomial coefcients solution of certain recurrence relations [START_REF] Agnarsson | On the sylvester denumerants for general restricted partitions[END_REF]. However, these bounds are not of real interest for the SD problem for two reasons: rst, the bounds concern a particular denumerant rather than the solutions corresponding to an interval, as specied by Eq. (2); second, as we shall see with experiments, the bounds returned are not tight.

Frobenius numbers. The asymptotic behavior given by Eq. ( 4) apparently contradicts the existence of (small) masses which cannot be decomposed, this issue being related to the so-called

Frobenius number and its generalizations. The Frobenius number g 0 (W) is dened as the largest integer which cannot be represented as a non-negative integer combination of the integers in W

Inria

Stoichiometry Determination for Mass-spectrometry Data 7 [START_REF] Ramírez Alfonsín | The Diophantine Frobenius Problem[END_REF]. While explicit formula are known up to M = 3, only upper and lower bounds are known in general [START_REF] Ramírez Alfonsín | The Diophantine Frobenius Problem[END_REF]. Computing the Frobenius number is a NP-hard problem [START_REF] Ramírez-Alfonsín | Complexity of the frobenius problem[END_REF] for which pseudo polynomial time algorithms exist, such as the one by Round-Robin [START_REF] Böcker | The money changing problem revisited: Computing the frobenius number in time o(ka1) *[END_REF]. (The running time of such an algorithm is polynomial with respect to the numerical values of the input data, but exponential in their bit-length.)

A related number is the positive Frobenius number g + 0 (W), namely namely the largest integer which does not admit any positive integer solution. As observed in [BDF + 10, Lemma 4], both numbers satisfy:

Either g + 0 (W) = g 0 (W) = 0, or g + 0 (W) = g 0 (W) + i=1,...,p w i . (5)
Unbounded knapsack (UKP) and subset-sum (SSP). Assume that the weights in W are real numbers coding the weights of p object types, and that each object also has a value v i ∈ R + . Given a target mass M , the unbounded knapsack problem (UKP) consists of nding for each type the integral quantity s i ≥ 0, such that the corresponding sum of values is maximum while the corresponding sum of weights does not exceed M . This is formally dened by the following integer programming model:

Maximize p i=1 s i v i , (6)
Under the constraint

p i=1 s i w i ≤ M. (7)
The special case where w i = v i holds is known as the subset-sum problem (SSP), and consists in nding the quantities s i such that the corresponding sum of weights is the closest to M , without exceeding it. Changing Eq. ( 7) into an equality allows subset-sum to solve the exact SD problem of Eq.(1).

UKP and SSP: algorithms and complexity. The optimization version of SSP is one of the rst problem proved to be NP-Hard [START_REF] Karp | Reducibility among combinatorial problems[END_REF], but it is also known to be pseudo-polynomial.

An example of such pseudo-polynomial time algorithms is well known Bellman recursion [START_REF] Bellman | Dynamic Programming[END_REF],

which solves UKP in O(M p) time. More recently, an output sensitive algorithm solving SSP has been developed [START_REF] Bocker | Ecient mass decomposition[END_REF]. Denoting w 1 the smallest mass, the algorithm has complexity O(p w 1 D(M )), and relies on a data structure of size O(p w 1 ). We shall use our implementation of this algorithm, called DECOMP, as main contender.

Other approaches have been developed to solve UKP, either based on dynamic programming, on branch and bound, or on a combination of both [START_REF] Poirriez | A hybrid algorithm for the unbounded knapsack problem[END_REF]. Interested readers are referred to [START_REF] Pferschy | Knapsack Problems[END_REF] for detailed review on the dierent knapsack problems and on the dierent techniques used to solve them.

UKP and SSP: on the Number of Solutions

Non-negative solutions. The number of solutions #UKP(p, W, M ) to an unbounded knapsack problem or of a subset-sum problem that is constrained by Eq. ( 7) was rst studied by Beged-Dov [START_REF] Beged-Dov | Lower and upper bounds for the number of lattice points in a simplex[END_REF], and the bounds were later rened by Padberg [START_REF] Padberg | A remark on an inequality for the number of lattice points in a simplex[END_REF] and Lambe [START_REF] Lambe | Bounds on the number of feasible solutions to a knapsack problem[END_REF].

Denote #SD(p, W, M, ε) the number of solutions to the SD problem dened by the vector W, and by the target mass M ± ε. The knowledge of this value is of interest in the context of

RR n°8101

Agarwal and Cazals and Malod-Dognin stoichiometry determination, to avoid generating a large number of solutions which would be impossible to analyze. From the exact number of solutions to an UKP problem, one gets

#SD(p, W, M, ε) = #UKP(p, W, M + ε) -#UKP(p, W, M -ε -1). (8)
However, the previous works [START_REF] Beged-Dov | Lower and upper bounds for the number of lattice points in a simplex[END_REF][START_REF] Padberg | A remark on an inequality for the number of lattice points in a simplex[END_REF][START_REF] Lambe | Bounds on the number of feasible solutions to a knapsack problem[END_REF] only provide bounds. If # counts an exact number of solutions, denote # and # a lower and an upper bound on #, respectively. Using upper and lower bounds on #UKP(p, W, M ), one gets the following lower bound

#SD(p, W, M, ε) ≥ #UKP(p, W, M + ε) -#UKP(p, W, M -ε -1), (9)
together with the following upper bound

#SD(p, W, M, ε) ≤ #UKP(p, W, M + ε) -#UKP(p, W, M -ε -1). ( 10
)
In a nearby vein, a fully polynomial-time approximation scheme (FPTAS) Positive solutions. In line with the denition of the shifted Frobenius number of Eq. ( 5), the question of estimating the number of positive solutions is also of interest, in particular for the biological cases where positive stoichiometries are prescribed for selected protein types. To estimate the number of such solutions, it is sucient to compute the previous bounds with a target mass reduced by the minimal imposed stoichiometries.

3

The Interval Stoichiometry Determination Problem

Tree Like Enumeration

The exact case. A pedestrian way to compute the denumerant consists of exhaustively trying all solutions [DF80, MHFH10, MR12]. To solve the interval SD problem of Eq. (2) using the same principle, assume that the stoichiometry vectors are built incrementallysay from left to right, and that the stoichiometry vector S has been computed up to index i (see also Fig. 1 ). Under this assumption, we dene the remaining mass to be accounted for by the p -i proteins whose stoichiometry has to be determined by

M - 0 [S] = M, and M - i [S] = M - j=1,...,i s j w j for i = 1, . . . , p -1, M - i [S] ≥ 0. (11)
Then, denoting a | b the fact that the integer a divides the integer b, the naive counting strategy for the denumerant consists of computing the following nested sum: The interval case. To solve Eq. (2) using the generalized sum of Eq. ( 12), observe that upon reaching the last sigma, the following condition holds:

M/w1 s1=0 M - 1 [S]/w2 s2=0 • • • M - p-2 [S]/
M -ε - j=1,...,p-1 s j w j ≤ s p w p ≤ M + ε - j=1,...,p-1 s j w j , (14)
or equivalently M -ε -j=1,...,p-1 s j w j

w p ≤ s p ≤ M + ε -j=1,...,p-1 s j w j w p (15)
Therefore, denoting #(•) the length of an integer interval, counting the number of stoichiometry vectors solving the SD problem is done by

(M +ε)/w1 s1=0 (M +ε) - 1 [S]/w2 s2=0 • • • (M +ε) - p-2 [S]/wp-1 sp-1=0 I (w p : s 1 , s 2 , s 3 ....s p-1 ) (16)
with

I (w p : s 1 , s 2 , s 3 ....s p-1 ) = # sp ( (M -ε) - p-1 [S]/w p ≤ s p ≤ (M + ε) - p-1 [S]/w p ) (17)
Note that when ε = 0, the equations (16-17) and (12-13) coincide.

Remark 1. For the sake of conciseness, given a stoichiometry vector S, let S = si∈S s i w i . Also, let the remaining mass associated to S be dened as m = M + ε -S . One has:

M -ε ≤ S ≤ M + ε ⇔ 0 ≤ S -M + ε ≤ 2ε ⇔ 0 ≤ m ≤ 2ε. (18)
These equivalent conditions characterize a valid stoichiometry vector.

Algorithm DIOPHANTINE.

w p = w min ≤ 2 × ε, (19)
all leaves are fertilethe length of the interval dened in Eq. ( 14) is 2ε so that is contains at least one solution for s p .

Regarding the memory footprint, because the stoichiometry vector S is passed by reference, a unique vector is used along the whole recursion tree, and we have:

Observation. 1. Algorithm DIOPHANTINE takes Θ(p) storage.

Remark 2 (On positive versus non-negative solutions.). As seen from the sigmas of Eq. ( 16), algorithm DIOPHANTINE generates all non-negative solutions. If minimal stoichiometries for the proteins are imposed, compliant solutions can be generated by starting the summations at these stoichiometriesin particular, staring at 1 rather than 0 yields positive solutions. Such solutions can also be generated by subtracting the mass corresponding to these constraints from the target mass, and seeking non-negative solutions.

Agarwal and Cazals and Malod-Dognin

Algorithm 1 Algorithm DIOPHANTINE. The algorithm stores the stoichiometry vectors in a set

Sol. The recursive function takes four arguments : the stoichiometry vector S under construction, passed by reference as indicated by the ampersand & (C++ convention); the remaining mass

m = M - i-1 [S]
of Eq. ( 11); the error threshold ε; and the index i ∈ [1, . . . , p] of the protein type to be processed by the current recursive call. The weights are sorted in the decreasing order. The initial call is DIOPHANTINE(S 0 , M +ε, ε, 1), where S 0 is a stoichiometry vector whose p entries are set to 0, M is the target mass, and ε is the allowed error. The sub-routine keep_value checks whether the stoichiometry j of the ith protein is admissibleas one may impose a lower and/or upper bound on the stoichiometry of that type.

DIOPHANTINE(Stoi_vector & S, unsigned m, unsigned ε, unsigned i)

// Stop the recursion if we are on the last protein type if i == p then q min = (m -2ε)/w p // p-th type: min stoichiometry q max = m/w p // p-th type: max stoichiometry for j = q min to q max do if keep_value(p, j) then

S[p] = j

Insert S into Sol // Recurse to set the remaining stoichiometries else quotient = m/w i // i-th type: max stoichiometry for j = 0 to quotient do if keep_value(i, j) then On the sortedness of weights. The cost of a particular enumeration problem is provided by the number of edges of the recursion tree explored by DIOPHANTINE, which actually depends on the sortedness of the vector W of weights:

S[i] = j //Set the i-th stoichiometry DIOPHANTINE(S, m -(j w i ), ε, i + 1) //
Observation. 2. The size of the recursion tree of algorithm DIOPHANTINE is minimized when the vector of weights is sorted in descending order.

Proof. The proof is divided into three steps. First, we prove that swapping two consecutive weights w i and w i+1 only changes the number of node at depth i. Second, that swapping w i and w i+1 results in a smaller number of node at depth i only if after reordering, w i ≥ w i+1 . Third, that iteratively applying this results leads to the optimal weight reordering w

1 ≥ w 2 ≥ • • • ≥ w p
Point 1: Consider the upper bound of Eq. ( 14). Without reordering, a node of the exploration tree at depth k < p is a solution to k i=1 s i ×w i ≤ M +ε, and the number of such node is equal to the corresponding number of such solutions (which we denote here by # k ). For depth p, # p is the number of solution as specied by Eq. ( 17). Thus, the number of tree node is

# 1 + # 2 + • • • + # p .
Let now swap w i and w i+1 . The numbers of nodes at depth 1 to i -1 do not depend on w i nor on w i + 1, and thus do not change. # i , the new number of nodes at depth i is now the number of solutions to i-1

j=1 (s j × w j ) + s i+1 × w i+1 ≤ M + ε, which is dierent from # i . The number of nodes at depth i + 1 is the number of solutions to i-1 j=1 (s j × w j ) + s i+1 × w i+1 + s i × w i ≤ M + ε,
which is the same as # i+1 , and the same holds for larger depths.

Point 2: For any branch dened by a sub-solution S = (s 1 , ..., s i-1 ), the number of nodes at

depth i is (M + ε) - i-1 [S]/w i + 1.
This number is smaller for larger values of w i , and is thus minimum if we reorder w i and w i+1 so that w i ≥ w i+1 Point 3: The optimal weights reordering is then by decreasing order of weights (i.e. when

w 1 ≥ w 2 ≥ • • • ≥ w p )
. This is proved easily by remarking that using any other weights reordering suppose that there exist a i for which w i < w i+1 , implying that the number of nodes in the corresponding exploration trees can be decreased by swapping protein types i and i + 1.

Enumeration Based on Dynamic Programming

The classical solution to the money changing problem based on dynamic programming is well known [START_REF] Bocker | Ecient mass decomposition[END_REF]. It uses a binary table stating whether the target mass is decomposable using a subset of proteins, this table being used by the backtracking algorithm reporting all solutions.

In the sequel, we show that a slight modication of the construction of the binary table based on Eq. (15) helps in accommodating the interval case. The enumeration algorithm using this table shall be denoted DP++, and is called as DP++(M + ε, p).

Assume that the weights w i are sorted by increasing value. Along the course of the backtracking, denote m the mass remaining from M + ε once a set of protein instances have been used. That is, given the stoichiometry vector S, one has remaining mass,

m = M + ε - si∈S s i w i (20)
We wish to build a binary table whose semantics is the following:

B[i, m] = 1 ⇔ ∃ a mass x in the interval [max(0, m -2 ), . . . , m]
such that x is decomposable over {w 1 , w 2 , ...w i }.

RR n°8101

Denote #s 1 the number of non-negative integer values satisfying (m -2ε)/w 1 ≤ s 1 ≤ m/w 1 . For i = 1, . . . , M + ε and i = 1, . . . , p, the binary table B[i, m] as follows:

For i = 1 : B[1, m] = 1 if #s 1 ≥ 0, and 0 otherwise (21) For i > 1 : B[i, m] = B[i -1, m] if m < a i B[i -1, m] ∨ B[i, m -a i ] otherwise. ( 22
)
DP++ consists of backtracking as usual using this binary table. DP++ is also output sensitive, a property inherited from the original algorithm. The number of backtracking steps for each solution being equal to (M + ε)/w 1 in the worst case, if D(M, ε) denotes the total number of solutions, the enumeration runs in time O((M + ε) D(M, ε)/w 1 ).

Solution Sketches to Represent a Solution Set

Because an interval SD problem may admit a large solution set D, the question arises to represent all such solutions in a compact way. For a solution S ∈ D, we dene the stripped solution S by S = {. . . , s i , . . . } with s i = s i -m i , and m i = min

S∈D s i . ( 23
)
Note that the stoichiometries m i dene the background of D. In addition, we dene the skeleton

B of a solution S by B = {. . . , b i , . . . } with b i = 1 if s i > 0, and 0 otherwise. ( 24
)
The set of ones in a skeleton vector identify the protein types W B involved in a stripped solution.

Since inclusion among such sets denes a partial order, we construct the following Hasse diagram, called the sketch of the solution set D:

A node corresponds to a set W B . Each such node is associated the solutions of D whose skeleton denes W B . The multiplicity µ of the node is the number of such solutions.

An edge links the nodes W B and W D i W B ⊂ W D , and there does not exist a node

W C such that W B ⊂ W C ⊂ W D .
Note that stripping a solution set D is mandatory for positive solutions of an interval SD problem or the sketch would involve a single node corresponding to all protein types.

Experiments: Implementation and Datasets

In this section, we sketch the software developed, and introduce several datasets used in our experiments, namely Data-Bio (10 instances), Data-Pseudo-Bio10 (200 instances), Data-Prime10 (200 instances).

Implementation sketch

Algorithms. Four algorithms were implemented in (generic) C++:

The algorithm of [START_REF] Böcker | The money changing problem revisited: Computing the frobenius number in time o(ka1) *[END_REF] to compute Frobenius numbers.

The FPTAS of [GKM + ] to estimate the number of solutions of any knapsack problem.

Inria

Stoichiometry Determination for Mass-spectrometry Data 13

The algorithms DIOPHANTINE and DP++ from section 3 to solve the interval SD problem.

The algorithm of [START_REF] Bocker | Ecient mass decomposition[END_REF] to solve SSP, called DECOMP. Recall that this algorithm has a preprocessing step consisting of computing the so-called extended residue table (ERT), followed by the backtracking step which decomposes a particular mass. To solve the interval SD problem, we compute the ERT once, and call the decomposition procedure for every mass in the interval.

The rst three are made available within the program addict, which can be downloaded form http://team.inria.fr/abs/addict/.

To avoid overows, unsigned integers coded on 64 rather than 32 bits were used. Such number types respectively accommodate integers up to 2 64 -1 ∼ 0.18 × 10 20 and 2 32 -1 ∼ 0.42 × 10 10 , and instances with more than 10 10 solutions are commonly facedsee Experiments. We note that on Linux operating systems, the type unsigned long long is always coded on 8 byte, be the system 32 or 64 bits, see e.g. http://en.cppreference.com/w/cpp/language/types.

Running times. Practically, the experiments were conducted on a DELL PRECISION T7400 computer equipped with 8Go of RAM, running Fedora Core 14. For all programs, a cut-time of three hours was set, in order to handle the case of enumeration problems yielding an astronomical number of solutions. For such cases, we estimated the number of solutions using the FPTAS strategy of [GKM + ]. In fact, whenever the 3 hours time limit was hit, we obtained a lower bound on the number of solutions by taking the maximum value of the number of solutions computed by DIOPHANTINE upon hitting the cut-time, and the lower bound computed from the FPTAS. In analyzing running times, for algorithms DECOMP and DP++, the total running was split into pre-processing (binary table building), and post-processing (enumerating the solutions), e.g.

t Tot DECOMP = t Pre DECOMP + t Post DECOMP .

Biological Dataset

In the following, we briey present the 10 biological complexes processed, whose masses span the range 321,274 5,276,467 Da (Table 1). This table lists experimental molecular weights, except for NPC (Y -ring and single spoke), where theoretical molecular weights were computed from the known stoichiometries. The reader is referred to the supplemental section 8.1 for more details.

Yeast 19S Proteasome lid.

Articial Datasets

Prime numbers instances: the Data-Prime10 dataset. The use of prime number weights prevents relative multiplicity of individual weights, and therefore avoids outright hardness of the instance stemming from so-called spanner sets [START_REF] Pisinger | Where are the hard knapsack problems?[END_REF]. To generate the kth synthetic complex, with k ∈ 1, . . . , 10, p prime numbers were chosen randomly from the list of prime numbers from 7,000 to 70,000, in accordance with the weights (in Daltons) of individual proteins involved in the biological complexes of Data-Bio. Moreover, 20 instances for each kth weight vector(W k ) were generated by using 20 target masses uniformly spaced in the range

[ w i , 2g + 0 (W k )]. (25)
These 200 instances are referred to as Data-Prime10-i-j.

Random biological instances: the Data-Pseudo-Bio10 dataset. To generate synthetic complexes using the masses of real proteins, we replaced the set of prime numbers by the masses of ∼ 6700 YEAST proteins retrieved from the non-redundant UNIPROT protein database, see http://www.uniprot.org/help/about. The theoretical molecular weight of proteins were calculated by adding up the weight of constituent amino acids in the proteins. (To be consistent with Data-Prime10, only masses beyond 7 kDa were retained.) As above, 10 vectors of masses were picked, with 20 target masses in the interval of Eq. ( 25), yielding 200 instances denoted Data-Pseudo-Bio10-i-j.

Experimental Results

The following comments are in order. First, in structural proteomics, the protein types involved in a given protein complex can generally be determined by mass spectrometry upon dismantling the complex into its individual sub-units, thanks to denaturing conditions. This explains why Data-Prime10, enumerated in this order. Third, three error levels typical in structural proteomics were used, namely 0%, 0.1% and 1% the value ε of Eq. (2) being set to the error level times the target mass.

Notation-wise, the median of a list L of values is denoted median(L). Also, the linear correlation coecient (i.e. Pearson) between two variables X and Y is denoted Corr(X, Y ).

Biological Examples: Enumeration Matters Even at Null Noise Level

The formula recalled in Eq. ( 4) shows an asymptotic polynomial growth of the number of solutions, yet, are situations with multiple solutions commonplace in biology? Table 1 answers this question positively, with multiple solutions even at null noise level in some cases. The case of the NPC is interesting. The case of one full spoke shows that the stoichiometry determination problem in this case is ill-posed, with an astronomical number of solutions. The same holds for the NPC-Y-ring system, which consists of 8 copies of the Y-complex, since a total of 788 solutions are obtained.

The last column of this table also shows that a large value M/g + 0 is a good hint at a large number of solutions.

Also, a solution sketch for positive solutions, as dened in section 3.3, is presented on Fig. 2 for the system EhATPase-sub-4. While on this simple example every node of the Hasse diagram corresponds to a single solution, other example yield more complex (and cluttered) Hasse diagrams, with numerous solutions per node. Interestingly, we also observed that for non-negative solutions and even a 0.1% noise level, all solution sketches of biological complexes pretty much involve all possible tuples (supplemental Figs. 17 and 18). That is, the Hasse diagram essentially contains all p k=1 p k nodes which are dened by Eq. 24.

Counting Solutions and Convergence to the Denumerant

Rationale. As noticed in section 2, counting solutions ahead of computing them is of interest to avoid generating a large number of (useless) solutions. We compared the accuracy of the estimations provided in [Pad71, Bd72, Lam74] against those of the FPTAS algorithm of [GKM + ].

One comment is also in order about hard instances, i.e. instances such that DIOPHANTINE does not terminate within a prescribed time limitthree hours in this work. According to [START_REF] Pisinger | Where are the hard knapsack problems?[END_REF],

two parameters lead to hard UKP and SSP instances: the magnitude of weights, and the linear redundancy between the weights. Intuitively, if a given weight is a linear combination of other weights, additional solutions can be generated. However, beyond the Frobenius number and in the asymptotic regime, the number of solutions of SSP is given by Eq. 4. Therefore, instead of using a target mass proportional to the sum of masses, as in [START_REF] Pisinger | Where are the hard knapsack problems?[END_REF], we systematically explored ranges of target masses expressed in units of Frobenius number.

Results: counting solutions. The biological examples just discussed show that estimating the number of solutions is important to detect ill posed problems.

Table 2 compares various strategies to estimate the number of solutions to a SD problem.

First, we plugged the upper and lower bounds provided by [START_REF] Padberg | A remark on an inequality for the number of lattice points in a simplex[END_REF], [START_REF] Beged-Dov | Lower and upper bounds for the number of lattice points in a simplex[END_REF], and [START_REF] Lambe | Bounds on the number of feasible solutions to a knapsack problem[END_REF] into Eqs.

(9-10). For each biological system, we retained the best bounds obtainedthose yielding the tighter interval (Table 2, columns 2-4). Second, we also computed estimates using the FPTAS

of [GKM + ] (

RR n°8101

Agarwal and Cazals and Malod-Dognin

On the one hand, the combinatorial bounds are not of real interest, since at least two orders of magnitude of dierence between #UKP(p, W, M ) and #UKP(p, W, M ), yielding a null lower bound. On the other hand, the bounds from the FPTAS are much tighter.

Results: convergence towards the denumerant. Due to the rapid rise of the number of solutions, an interesting problem is the speed of the convergence of the number of solutions of a SD problem to the denumerant.

To study this convergence, a toy dataset called Data-Quad4 consisting of the quadruplet {10, 15, 32, 48} from [START_REF] Shallit | Unbounded discrepancy in frobenius numbers[END_REF] was used. A total of 100 instances uniformly distributed in the range [ w i , 41g 0 ] was used, the jth instance being denoted Data-Quad4 -j. We also studied this convergence for a real biological system of Yeast Exosome complex. A total of 40 instances having target masses uniformly distributed in the range [ w i , 20g 0 ]. The results on this system and on a real biological system are presented on Figs. 3 and4, and show a rapid convergence to the denumerant. The ratio of denumerant to the number of non-negative solutions grows to 96% at 40g 0 starting from 29%, in case of Data-Quad4 whereas, for the real biological instance, this ratio grows to 70% at 10g 0 from its initial value of 12%. The Yeast Exosome instances having target mass greater than 11g 0 did not terminate within the time limit of 3h. Note that to compare with Denumerant, non-negative solutions are plotted.

Running Times: DECOMP, DIOPHANTINE, DP++

Rationale. To compare the three running times t Tot DIO , t Tot DECOMP , t Tot DP++ , we processed our three datasets at the noise levels 0%, 0.1%, and 1% (Fig. 5).

Results DP++ over DECOMP owes to the constants involved in the output-sensitive complexity, which are respectively of O(D(M, 0) M/w min ) and O(D(M, 0) p w min ). Practically, biological complexes are such that p ∼ 10, w min ∼ 10 kDa, and M ∼ 1 MDa. At 0.1% and 1% noise level, DECOMP is not competitive anymore, since it is outperformed by three to four orders of magnitude. Clearly, the wider the interval [M -ε, . . . , M + ε], the more challenging the instances for DECOMP, an issue precisely analyzed in section 5.5.

In presence of noise, the two competitors are thus DIOPHANTINE and DP++. At 0.1% and 1%

respectively, one has median({t Tot DP++ /t Tot DIO }) = 1.07, median({t Tot DP++ /t Tot DIO }) = 1.2. More interestingly, DIOPHANTINE is never twice slower than DP++, while in the worst case, DP++ is 17 (resp. 20) times slower than DIOPHANTINE at 0.1% (resp. 1%). Thus, DIOPHANTINE is more stable.

We also notice that asymptotically when the target mass increases, the ratio of running times converges to one, hinting at an output sensitive behavior of DIOPHANTINE recall that DP++ is output sensitive.

Output sensitivity of DIOPHANTINE

useless operations. The running times of DIOPHANTINE being comparable to those of the output sensitive algorithm DP++, we analysis the output sensitive behavior of DIOPHANTINE.

Results: running time versus number of non-negative solutions. We rst computed the Pearson correlation coecient between t Tot DIO and the number of solutions (Fig. 6). With values equal to 0.97, 0.89 and 0.97 for the three datasets at 0.1% noise level, and to 0.97, 0.89 and 0.97 at 1% noise level, this correlation is excellent.

Results: running time versus the recursion tree size. The correlation between t Tot DIO and the recursion tree size is even better (Fig. 7), with coecients of 0.99, 0.99 and 0.99 at 0.1% noise level, and of 0.99, 0.95 and 0.98 at 1% noise level, respectively for the three datasets.

We further plotted the proportion of fertile edges as a function of the target mass (Fig. 8).

Note that this proportion is 1 when the condition of Eq. ( 19) holds, which is the case for 0/410 and 338/410 instances at 0.1% and 1% respectively. As expected, this proportion increases with the error level: at 1%, all the instances which do not meet the condition of Eq. ( 19) have a ratio above 0.6, while no instance has such a high ratio at 0.1%.

RR n°8101

Agarwal and Cazals and Malod-Dognin

Decreasing Performances of DECOMP upon Increasing the Noise Level

Rationale. The larger the noise level, the better the performances of DP++ w.r.t. DECOMP. Also, algorithm DP++ is output sensitive for the interval SD problem, while algorithm DECOMP is output sensitive for the exact SD problem but possibly not for the interval SD problem. We wish to relate these two facts, using two ingredients: First, the decomposition of the running time into pre and post-processing, namely t Tot DP++ = t Pre DP++ + t Post DP++ , and likewise for DECOMP.

Second, the number of nodes explored by DP++ and DECOMP. For DP++, this number is the number of nodes in the backtracking procedure, as explained in section 3.2. For DECOMP, this number is equal to the sum of the number of explored nodes for target masses in the range [M -ε, M + ε].

Results: for algorithm DP++, the post-processing time t Post DP++ converges to total runtime t Tot DP++ when the number of solutions increases. We rst analyze the relative importance of the pre-processing and post-processing times for both algorithms. As shown by the scatter plots (Fig. 9), the post-processing step always dominates for DECOMP, while for DP++ it becomes dominant as the number of solutions increases. For the latter algorithm, the linear trend between t Post DP++ /t Pre DP++ and the number of nodes but also the number of solutions is conrmed by correlation coecients beyond 0.9, at both noise levels 0.1% and 1% (Table 3, and supplemental Figs. SI-11, SI-12 ). Interestingly, the correlation between the ratio of running times t Post DP++ /t Pre DP++ and the target mass is poor, with coecients of 0.32 and 0.11 at 0.1% and 1% noise level respectively.

Results: the decreasing performance of DECOMP owes to a redundant post-processing. 27) with N = 203 3 and with N = 112 4 instances respec- tively at 0.1% and 1% noise levels is that the overall correlation is spoiled by the instances for which the pre-processing time still dominates in a run of DP++. To measure the incidence of the convergence of t Post DP++ to t Tot DP++ when the number of solution increases, we therefore resort to a sequence of correlation coecients, the set of instances used to investigate the correlation corresponding to the cases where the ratio t Post DP++ /t Pre DP++ is larger than some threshold. More precisely, given a set of N instances, we compute N correlation coecients as follows:

Based on the running times of algorithm DP++ we compute

α i = t Post DP++ /t Pre DP++ ( 28
)
for each instance, and sort the N instances by increasing α i value.

For each index i = 1, . . . , N -1, let S >i : the set of instances such that α j > α i .

(

) 29
We compute the Pearson correlation coecient Corr(

t Tot DECOMP t Tot DP++ , #nodes DECOMP #nodes DP++ ) on the set S >i .
The plot of these coecients is presented on Fig. 10. If one omits the last 10 coecients which are inherently unstable since they involve less than 10 points, the values obtained for the coecient of Eq. ( 27) now reach 0.9. Two points are noticeable. First, at 0.1% noise level, the correlation rises rapidly as a function of the α i value. As scatter plots show (supplemental Fig. 19), the decreasing section for the trailing 30 instances or so owes to the paucity of the dataset. Second, at 1% noise level, the correlation rises monotonically over the range of α i values explored. This indicates that the increase in number of nodes of DECOMP translates directly on the running time t Tot DECOMP , which is expected since when α i increases, the post-processing time which depends on the number of nodes converges to the total running time. Again, the ratio #nodes DECOMP /#nodes DP++ directly measures the redundancy of the work carried out by algorithm DECOMP with respect to the output sensitive algorithm DP++, this redundancy having a linear incidence on the running time.

Discussion. These observations are actually consistent with the complexities of the pre-processing steps and the number of nodes generated. Indeed, the pre-processing step in the case of DECOMP is independent of the target mass of the instance (O(p w min )), while in case of DP++ it is linearly dependent on the target mass (O((M + )p)).

As for the number of nodes generated during the post-processing step, the degraded performances of DECOMP for the interval SD problem is actually expected. Indeed, while a given stoichiometry vector is generated exactly once for any interval SD problem (the solutions to two dierent target masses are dierent), redundant calculations are possibly carried out by dierent runs of DECOMP for consecutive target masses. Ideally, one would measure the redundancy of calculations carried out to generate two dierent solutions is precisely measured by the longest common sub-sequence between the stoichiometry vectors of these two solutions. More plainly, we use the number of steps required by the backtracking procedure.

Conclusion and Outlook

Mass spectrometry is playing an increasing role to investigate large macro-molecular assemblies, thus complementing more classical techniques such as X ray crystallography or nuclear magnetic resonance, which are better suited for smaller complexes. In this context, the stoichiometry determination (SD) problem, which consists of determining how many copies of each sub-unit are required to account for the observed mass, is the rst one to be addressed, before investigating the geometry of the contacts between these sub-units. In this work, we develop a constant memory space enumeration algorithm (DIOPHANTINE), and an output sensitive dynamic programming based algorithm (DP++), which are shown to outperform state-of-the-art SD algorithms by several (three to four) orders of magnitude. These two algorithms exhibit comparable performances for typical noise levels in the range 0.1% to 1%, which is remarkable since DIOPHANTINE does not use any pre-processing and has constant memory footprint. Both algorithms performed to satisfaction on all the biological cases processed, which are the most challenging ones currently investigated in structural biology. It is noticeable that enumeration does matter, since systems with numerous solutions are observed at a moderate (sometimes even null) noise level. Also, the coupling between DIOPHANTINE (or DP++) and the fully polynomial-time approximation scheme (FPTAS) of [GKM + ] provides a powerful way to get a lower bound on the number of solutions, and thus to identify ill-posed SD problemswhich admit a too large number of solutions.

From a theoretical standpoint, several outstanding questions remain open. The rst one is related to the output sensitivity of DIOPHANTINE, which is observed in practice, but cannot be guaranteed in general. Getting ner insights on this phenomenon may allow developing more ecient enumeration schemesby pruning sterile portions of the recursion tree. The second one is reminiscent from phase transitions, a well known phenomenon for selected hard optimization problemssuch as 3-SAT for example. While the SD problem is about enumeration rather than optimization, linear redundancies within protein masses (one can trade a set of proteins against another set) have a direct impact on the rise of the number of solutions. A challenge is therefore to gain insights on the role of linear relationships and the total number of solutions. The third one is the asymptotic analysis of the number of solution of an interval stoichiometry determination problem. While the asymptotics of the denumerant D(M, 0) is known, working out the behavior or D(M, ε) in conjunction with the distribution of solutions per node of the Hasse of the solution diagram sketch would provide valuable information on solutions as a function of the input masses and the noise level. Artwork

Algorithms

Figure 1 The recursion tree explored by algorithm DIOPHANTINE. A leaf is called fertile if Eq. ( 17) has at least one solution, and sterile otherwise. Similarly, any edge of the tree is termed fertile if at least one fertile leaf is found downstream the tree, and sterile otherwise. The size of the recursion tree is dened as its number of edges.

Protein of mass w 2

RR n°8101

Figure 2 The 20 solutions for EhATPase-sub-4 at 0.1% noise level, see also Table 1 The Hasse diagram presents the sketch of the solution set, as dened in section 3.3, while the background presents the min stoichiometry of each type present in any solution. Each node of the Hasse diagram reads as follows: the left hand side presents the number of solutions. Each row on the right hand side reads as follows: (top) list of protein types (middle, bottom) min and max stoichiometries of all solutions using these types. In the block for Background information, rst row corresponds to the index of protein types and the second row corresponds to the minimum stoichiometry of these protein types across all the positive solutions. Note in particular that an arrow represents the inclusion between the protein types of two nodes. 2 Estimating the number of positive solutions for ve biological systems.

1 1 1 1 2 2 2 3 4 4 1 1 7 7 5 1 1 1 2 5 5 3 2 2 1 2 2 2 9 1 1 1 3 5 5 5 1 1 1 1 1 1 2 1 1 3 1 1 5 1 1 6 1 1 1 1 2 2 2 1 1 3 1 1 7 1 1 1 1 2 2 2 3 3 6 1 1 1 1 1 1 3 3 3 4 2 2 1 1 1 1 3 2 2 4 1 1 6 1 1 1 1 3 3 3 2 2 6 1 1 1 1 2 2 3 1 1 8 1 1 1 1 5 5 4 1 1 5 1 1 1 2 3 3 3 1 1 4 2 2 1 2 1 1 3 1 1 4 1 1 7 1 1 1 2 1 1 3 1 1 5 3 3 1 2 3 3 4 1 1 6 1 1 1 2 2 2 5 1 1 7 1 1 1 2 1 1 5 1 1 8 1 1 1 3 1 1 4 1 1 8 1 1 BACKGROUND 1 1 2 1
Columns 2 to 4: upper and lower bounds obtained by plugging the upper and lower bounds of [START_REF] Padberg | A remark on an inequality for the number of lattice points in a simplex[END_REF], [START_REF] Beged-Dov | Lower and upper bounds for the number of lattice points in a simplex[END_REF], and [START_REF] Lambe | Bounds on the number of feasible solutions to a knapsack problem[END_REF] into Eqs. (9-10). Columns 5 to 7: upper and lower bounds on the number of positive solutions obtained using the FPTAS of [GKM + ], with δ = 0.1.

Complex # M ± 1% M ± 1% solutions #SD #SD #SD #SD COP9 1 0 8 1 1 Y19-lid 1 0 10 1 1 Y-exosome 149 0 486 101 183 EIF3 1 0 14 1 1 NPC-1-Spoke Unknown 0
1.06 × 10 22 2.71 × 10 17 3.44 × 10 17

RR n°8101

Figure 3 Number of non-negative solutions vs. asymptototic behavior of the denumerant (Eq. ( 4)) for the quadruplet W = {10, 15, 32, 48}. The target mass is expressed in units of the Frobenius number.

Note that Y-axis is drawn with a logarithmic scale. (Left) Noise level of 0.1% (Right) Noise level of 1%. Note that Y-axes are drawn with a logarithmic scale. Figure 10 The parameterized correlation between the ratios of total running times t Tot DECOMP /t Tot DP++ and #nodes DECOMP /#nodes DP++ increases upon stepwise removal of in- stances. (Top) At an x-coordinate = k, the Pearson correlation value corresponds to the N -k instances of the set S >i , see Eq. ( 29). (Bottom) The x-axis features the sorted α i values. respectively.

An electrospray mass spectrum is recorded for rotary ATPases from E. hirae (EhATPase) and E. thermus (TtATPase). Each of these complexes have nine dierent protein types but, the membrane embedded rotor for EhATPase is larger because each K type contains four transmembrane helices as compared to two transmembrane helices in corresponding L type in TtATPase.

MS experiment data was retrieved from [START_REF] Morgner | Massign: An assignment strategy for maximising information from the mass spectra of heterogeneous protein assemblies[END_REF]. EhATPase was undergone controlled disassembly resulting in formation of sub-complexes in gas phase and solution phase using collision induced dissociation (CID) and partial denaturation by manipulating the ionic strength, respectively. We choose four sub-complexes 2,3,4 and 5 formed in the solution phase, and assumed that each subcomplex is composed of all the 9 types of proteins as does the intact complex.

Measured molecular weights of subcomplexes and percentage error in measurement are shown in table 5. While we are not aware of mass spectrometry experiments on the whole NPC, a number of sub-systems have been studied in detail. One of them is an heptameric sub-complex known as the Y-complex [FMPS + 12, DDC12], which participates to the formation of the scaold of the NPC, as one nds one copy of the Y-complex per half-spoke.

In the sequel, we simulated mass spectrometry data for two complexes. The rst one is a complete spoke, which contains a total of 57 proteins instances of 30 dierent types. The second

Inria

Stoichiometry Determination for Mass-spectrometry Data 39 one is the Y-ring complex, containing 8 copies of the Y-complexthat is 56 proteins of 7 types.

For each complex, noise levels of 0%, 0.1% and 1% were applied to the exact mass of the complex computed from the masses of the individual proteins.

Plots Corresponding to the

  Proteomics, Integer Partitions and Knapsack Problems

  Determination for Biological Examples: Enumeration Matters

  Solutions and Convergence to the DenumerantTable

Figure 4 Figure 5

 45 Figure 4 Number of non-negative solutions vs. asymptototic behavior of the denumerant (Eq. (4)) for Yeast Exosome. The target mass is expressed in units of the Frobenius number. Instances with a target mass beyond 11g 0 did not terminate within the time limit. Note that Yaxis is drawn with a logarithmic scale.

Figure 7 Figure 8

 78 Figure 7 DIOPHANTINE: running time t Tot DIO as a function of the recursion tree size. (Left)Noise level of 0.1% (Right) Noise level of 1%. Note that Y-axes are drawn with a logarithmic scale.

  Figure 9 DECOMP [Top] and DP++[Bottom]: post-processing time (i.e. t Post ) as a func- tion of the pre-processing time (i.e. t Pre ). Note that Y-axes are drawn with a logarithmic scale.

Figure 14

 14 Figure 14 Pairwise plots among ratio of runtimes and ratio of number of nodes for DECOMP over DP++ (Columns) The left and right columns respectively correspond to noise levels 0.1% and 1%. (Top row) Ratio of number of nodes explored vs. Ratio of total run time for DECOMP over DP++ (Middle row) Ratio of number of nodes explored vs. Ratio of post-processing time (backtracking) for DECOMP over DP++ (Bottom row) Ratio of number of nodes explored vs. Ratio of pre-processing time for DECOMP over DP++

  Recursion

Inria

Stoichiometry Determination for Mass-spectrometry Data 11

  Nuclear Pore Complex (NPC). The NPC, which is the largest protein assembly known to date in the eukaryotic cell[START_REF] Wente | The Nuclear Pore Complex and Nuclear Transport[END_REF][START_REF] D'angelo | Structure, Dynamics and Function of Nuclear Pore Complexes[END_REF], is a protein assembly anchored in the nuclear envelope, regulating the nucleo-cytoplasmic transport. It is composed of ∼ 30 distinct proteins types each present in multiple copies. It has eight-fold radial symmetry and consists of eight spokes, each containing 57 protein instances. One particular complex, the Y -complex,

	Agarwal and Cazals and Malod-Dognin
	Yeast Exosome. The exosome is a 3'-5' exonuclease complex involved in RNA processing and
	degradation. The yeast exosome is composed of 10 dierent protein types with unit stoichiometry
	[HDT + 06].
	Rotary ATPases. Rotary ATPases are membrane associated molecular machines involved in
	energy conversion by coupling ATP hydrolysis (or synthesis) with proton (or Na	+ ) translocation
	across biological membranes. Practically, we investigate the intact TtATPase, and four sub-
	complexes of EhATPase.
	cap involving a peripheral lid, composed of 9 distinct protein types each with unit stoichiometry
	[STA + 06].
	COP9 Signalosome. The COP9 Signalosome is a multi-functional complex primarily involved
	in ubiquitin mediated proteolysis linked to diverse cellular activities such as signal transduction,
	cell cycle progression and transcriptional regulation. It is composed of 8 protein protein types
	each with unit stoichiometry.
	Eukaryotic Translation factor EIF3. Eukaryotic initiation factors (eIF) are proteins in-
	volved in the initiation phase of the eukaryotic translation. They form a complex with the 40S
	ribosomal subunit, initiating the ribosomal scanning of mRNA. Among them, EIF3 consists of
	13 dierent protein types each with unit stoichiometry.
	RR n°8101

Proteasomes are protein assemblies involved in the elimination of damaged or misfolded proteins, and the degradation of short-lived regulatory proteins. The most common form of proteasome is the 26S, which involves two ltering caps (the 19S), each Yeast involving 7 protein types with unit stoichiometry, is present in two sets of 8 copies. Each such set presumably forms one ring (the Y -ring) involving 56 protein instances [FMPS + 12, DDC12].

Table 1

 1 is about positive solutionsrather than non-negative solutions. But all the remaining experiments have been conducted with non negative solutions. Second, when three series of numbers are reported, they refer to the three datasets Data-Bio, Data-Pseudo-Bio10 and

Inria

Stoichiometry Determination for Mass-spectrometry Data 15

Table 2

 2

, columns 5-7).

  . While DIOPHANTINE and DP++ completed all instances within the imparted time,

	regardless of the noise level, on a per-dataset basis, DECOMP nished on 10/10, 200/200, 200/200
	cases at zero noise level, 9/10, 104/200, 93/200 cases at 0.1% noise level, and 8/10, 56/200,
	51/200 cases at 1% noise level. At null noise level, one gets median({t Tot

DECOMP /t Tot DIO }) = 1.4, median({t Tot DIO /t Tot DP++ }) = 13.66, median({t Tot DECOMP /t Tot DP++ }) = 16.02. That is, DP++ is faster than DECOMP and DIOPHANTINE by about one order of magnitude. It should be noticed in particular that the better performance of

  Tot DECOMP , and the convergence of t Post DP++ to t Tot DP++ establish that the increased total running time t Tot DECOMP linearly depends on the number of nodes visited by DECOMP. In particular, taking as reference#nodes DP++ since algorithm DP++ is output sensitive, the ratio #nodes DECOMP /#nodes DP++ measures the redundant work of DECOMP corresponding to the successive calls to solve the individual SD problems.

	Stoichiometry Determination for Mass-spectrometry Data	19
	in computing the coecient of Eq. (
	correlation coecient	Corr(	t Post DECOMP t Post DP++	,	#nodes DECOMP #nodes DP++
							DECOMP to
	t Since the previous line of argument holds asymptotically, namely when t Post DP++ converges to t Tot DP++ , we also investigate directly
		Corr(	t Tot DECOMP t Tot DP++	,	#nodes DECOMP #nodes DP++	).	(27)
	With values of 0.48 and 0.35 at noise levels 0.1% and 1%, this correlation is weak. The problem
							Inria

Since algorithm DP++ is output sensitive, we take it as a yardstick and postulate that that if the number of nodes explored by DECOMP incurs an x-fold increase with respect to the number of nodes explored by DP++, so should the running-times. Phrased dierently, we compute the ).

(26)

As seen from Table

4

, this rst correlation is strong, with values of 0.75 and 0.79 at noise levels of 0.1% and 1% respectively. On removing two outliers at each of the noise level, the coecients improve to 0.86 and 0.93, respectively (Supplemental Fig.

. The linear relationship between the two terms of the correlation dened by Eq. (

26

), together with the convergence of t Post

Table 1

 1 Biological instances can lead to multiple positive solutions, even with null noise level. For the NPC-1-spoke, lower and upper bounds were obtained using UKP-FPTAS

	of [GKM	+ ] with an error factor of δ = 0.1.
	Complex COP9 Y-19S-lid EhATPase-sub-5 Y-exosome EhATPase-sub-4 EhATPase-sub-3 EhATPase-sub-2 TtATPase EIF3 NPC-Y-ring NPC-1-spoke	#sol, 0% noise #sol, 0.1% noise #sol, 1% noise 0 1 1 0 0 1 0 4 39 0 13 149 0 20 190 0 74 700 0 224 2,213 21 24,487 246,242 0 0 1 788 6,900,664 69,042,257 [1.72 × 10 16 [1.72 × 10 16 [2.71 × 10 17 , 1.73 × 10 16 ] , 4.77 × 10 16 ] , 3.44 × 10 17 ]	M (Da) 321,274 ± 35 376,151 ± 369 921,712 g + 0 961,855 387,356 ± 230 682,901 402,708 ± 68 649,185 424,441 ± 148 682,901 461,674 ± 324 682,901 500,178 ± 294 682,901 659,202 ± 131 607,304 797,999 ± 180 1,257,629 4,603,280 2,282,543 5,276,467 3,169,210	#sol(1% noise)/2ε 3.12 × 10 -4 1.87 × 10 -4 5.03 × 10 -3 1.85 × 10 -2 2.38 × 10 -2 8.10 × 10 -2 2.32 × 10 -1 18.7 8.05 × 10 -5 750 ≥ 2.57 × 10 12	+ 0 M/g 0.33 0.41 0.57 0.62 0.62 0.68 0.73 1.08 0.63 2.02 1.66

Table 3

 3 DP++: post-processing time (i.e. t Post DP++ ) dominates the pre-processing time (i.e. t Pre ) for instances with large #Nodes and large #Solutions. Numbers provided are Pearson correlation coecient at 0.1% and 1% noise level, respectively.

	#N odes #Solutions T arget M ass #Solutions T arget mass
	N oise level 0.1% 1% 0.1% 1% t Post DP++ t Pre DP++ 0.93 0.95 0.93 0.91	0.1% 1% 0.32 0.11 0.94 0.91 0.1% 1%

Table 4

 4 DP++ and DECOMP: increase in t Post DECOMP over t Post DP++ is due to relative size of trees (#nodes) for both the algorithms. Pairwise Pearson correlation coecients to compare DECOMP and DP++ at two noise levels. (See also gs. SI-13 SI-16).

		Noise Level 0.1% #nodes DECOMP #nodes DP++ t Tot DECOMP t Tot DP++	t Pre DECOMP t Pre DP++	t Post DECOMP t Post DP++
	error width	0.37	0.74 -0.21 0.16
	#nodes DECOMP #nodes DP++ t Tot DECOMP t Tot DP++ t Pre DECOMP t Pre DP++		0.48 0.46 0.75 -0.01 0.43 0.49
		Noise Level 1% #nodes DECOMP #nodes DP++ t Tot DECOMP t Tot DP++	t Pre DECOMP t Pre DP++	t Post DECOMP t Post DP++
	error width	0.51	0.66 -0.14 0.58
	#nodes DECOMP #nodes DP++ t Tot DECOMP t Tot DP++ t Pre DECOMP t Pre DP++		0.35 0.29 0.79 -0.13 0.73 0.12

  ) and of two regulatory caps ltering the entry of the proteins (the 19S). The 19S subcomplex itself subdivides into two other subcomplexes, the base that binds directly to the 20S core particle and a peripheral lid. The latter is composed of 9 dierent protein types with a single copy for each type. It is involved in recognition of the polyubiquitin chain of the substrate and followed by it's deubiquitination. This substrate is then unfolded and translocated to the core particle for further degradation [STA + 06].For yeast, the measured mass of the intact 19S Proteasome lid complex is 376,151 ± 369 Da.Summing the theoretical weights of the protein types amounts to 374,576 Da, i.e. the error in the measurement is 0.42%. There are primarily two types of ATPases, F-type and V-type complexes each having two domains F 0 , F 1 and V 0 , V 1 . The membrane embedded domains F 0 /V 0 mediate proton (or Na + ) translocation and F 1 /V 1 are involved in ATP production or consumption

	Stoichiometry Determination for Mass-spectrometry Data	37 Agarwal and Cazals and Malod-Dognin
	8	Supporting Information homologous to ribonucleases e.g. RNase PH and RNase II and others, e.g. S1, KH, PINc and
		HRDC [ACL	+ 02].
	8.1 Detailed Description of Biological Complexes The measured mass of the intact Yeast exosome complex is 397,860 ± 99 Da. Summing the
	theoretical weights of the protein types amounts to 397,881 Da, i.e. error w.r.t. to the theoretical Yeast 19S Proteasome lid. Proteasomes are protein assemblies involved in elimination of weights is 0.005%.
	0.95 damaged or misfolded proteins and the degradation of short-lived regulatory proteins. They are 1
	0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 Summing the theoretical weights of the protein types amounts to 321,270 Da, i.e. the error in 0 20 40 60 80 100 120 140 160 180 200 Pearson correlation coff. Compliment of Instance set size 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 10 20 30 40 50 60 70 Pearson correlation coff. Compliment of Instance set size Pearson correlation coff. 0.5 0.6 0.7 0.8 0.9 1 the measurement is 0.0012% [SMBE + 09]. Pearson correlation coff. found in all eukaryotic cells and archaea, and also in selected bacteria. The most common form Rotary ATPases. Rotary ATPases are membrane associated molecular machines involved in 80 90 of the proteasome is 26S, named after its sedimentation coecientexpressed in Svedberg. energy conversion by coupling ATP hydrolysis (or synthesis) with proton (or Na + ) translocation The 26S Proteasome, consists of one core particle corresponding to the degradation chamber coexpression in E.coli reveals the molecular weight for an intact complex to be 321,274 ± 35 Da. (the 20SAn Electrospray MS experiment conducted on human COP9 signalosome reconstituted by across biological membranes.	100
		0.55
				0.4
	0.45 0.5 Eukaryotic Translation factor EIF3. Eukaryotic translation initiation factors (EIFs) are the 0.3
	0 proteins involved in assembling of elongation competent 80S ribosome to initiate the translation 2 4 6 8 10 12 14 16 18 20 Alpha Value 0 0.5 1 1.5 2 2.5 Alpha Value	3	3.5	4
	Noise Level: 0.1% process. There are atleast nine EIFs involved in the initiation process. They carry out their Noise Level: 1%
	function in two steps: formation of 48S complex with established codon-anticodon base pairing
	in the P-site of 40S ribosomal subunits, and the joining of 48S complex with 60S subunits
	[JHP10].
		Among them, EIF3 binds to the 40S subunit of the ribosome to initiate protein synthesis
	followed by recruitment of messenger RNAs. It also promotes attachment of 43S complexes
	to mRNA. It has 13 dierent protein types with a single copy for each type. These 13 pro-
	tein types have been unambiguously identied when MS/MS spectra, and were scanned against
	UniProt/Swiss-Prot and NCBInr database using the MASCOT search engine [DFZ	+ 07]. The
	measured mass of the intact Yeast EIF3 complex is 797,999 ± 180 Da. Summing the theoretical
	weights of the protein types amounts to 793,558 Da, i.e. error w.r.t. to the theoretical weights
	is 0.56% [ZSF	+ 08].
	Yeast Exosome. The exosome is a 3'-5' exonuclease complex involved in RNA processing and
	degradation. In eukaryotes it is present in the cytoplasm and nucleolus and therefore reacts with
	dierent substrates in respective compartments. It is composed of 10 dierent protein types each
	with unit stoichiometry [HDT + 06]. The yeast exosome complex is known to contain domains
	RR n°8101			Inria

COP9 Signalosome. The COP9 Signalosome is a multifunctional complex primarily involved in ubiquitin mediated proteolysis linked to diverse cellular activities such as signal transduction, cell cycle progression and transcriptional regulation. It is composed of 8 protein types with a single copy for each type and shares remarkable homology with Yeast 19S proteasome lid complex.

Table 5

 5 Measured molecular weights (in Da) and percentage error in measurement w.r.t. weighted sum of the theoretical masses of individual protein types complying with their known stoichiometries. Refer to g. 6 of[START_REF] Morgner | Massign: An assignment strategy for maximising information from the mass spectra of heterogeneous protein assemblies[END_REF] and TableS1Aand S1B of supplementary information Yeast Nuclear Pore Complex (NPC). The NPC is a protein assembly anchored in the nuclear envelope, regulating the nucleo-cytoplasmic transport. It is composed of ∼ 30 distinct proteins types each present in multiple copies, and is the largest protein assembly known to date in the eukaryotic cell[START_REF] Wente | The Nuclear Pore Complex and Nuclear Transport[END_REF][START_REF] D'angelo | Structure, Dynamics and Function of Nuclear Pore Complexes[END_REF]. It has eight-fold radial symmetry and consists of eight spokes, each with a cytoplasmic and a nuclear side.

	of [ZMB	+ 11]
	Complex	Stoichiometries	Measured Mass (Da)	s i w t i	%error in measurement
	EhATPase-sub-2	A 3 B 3 CDE 2 F 2 G	500,178 ± 294	499,131	0.21
	EhATPase-sub-3	A 3 B 3 DE 2 F 2 G	461,674 ± 324	460,968	0.15
	EhATPase-sub-4	A 3 B 3 DEF G	424,441 ± 148	423,813	0.15
	EhATPase-sub-5	A 3 B 3 DG	387,356 ± 230	386,658	0.18
	TtATPase	A 3 B 3 CDE 2 G 2 F IL 12	659,202 ± 131	657,979	0.19

Table 3

 3 Figure 11 DP++: The Dominance of t Post DP++ over t Pre DP++ increases for instances having large tree size (#nodes) and large #solutions. [Top] t Post DP++ /t Pre DP++ as a function of number of nodes and [Bottom] Number of solutions, at 0.1% and 1% noise level for three datasets.

à l'algorithme DP++, lequel est sensible à la sortie. D'un point de vue appliqué, nous présentons trois résultats. Tout d'abord, nous montrons que DIOPHANTINE et DP++ sont plus rapides de 3 à 4 ordres de grandeur que l'algorithme état de l'art résolvant le problème exact, pour des complexes protéiques typiques présentant des incertitudes de l'ordre de 0.1% a 1%. Nous montrons que ce gain tient à l'évitement de calculs redondants lorsque l'algorithme exact est appellé sur chacune des masses d'un intervalle. D'autre part, nous montrons que DIOPHANTINE a un comportement d'autant plus sensible à la sortie que la taille de l'intervalle augmente. Enn et d'un point de vue biologique, en étudiant plusieurs complexes protéiques (facteur de traduction eucaryote, exosome de la levure, sous-unité 19S du protéasome, pore nucléaire), nous soulignons l'importance des solutions multiples, même à un niveau d'incertitude nul. Nous introduisons également une représentation de toutes les solutions d'un problème par intervalle, mettant en evidence les types de protéines les plus utilisés.

RR n°8101

Rationale. In the exact case, DIOPHANTINE can only be output sensitive for a target mass beyond the Frobenius number, since any call with a non representable mass necessarily results in Inria

number of instances for which both the algorithms terminated within 3h with 0.1% noise level. Note that DP++ terminated for all

instances.4 number of instances for which both the algorithms terminated within 3h with 1% noise level. Note that DP++ terminated for all 410 instances. RR n°8101

RR n°8101RESEARCH CENTRE SOPHIA ANTIPOLIS -MÉDITERRANÉE 2004 route des Lucioles -BP 93 06902 Sophia Antipolis Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Agarwal and Cazals and Malod-Dognin