Thorsten Schindler
email: thorsten.schindler@mytum.de

Binh Nguyen
email: nguyeb2@cs.rpi.edu

Jeff Trinkle

Understanding the difference between prox and complementarity formulations for simulation of systems with contact

To plan a robotic task involving intermittent contact, such as an assembly task, it is helpful to be able to simulate the task accurately and efficiently. In the past ten years, the prox formulation of the equations of motion has arisen as a competitive alternative to the well-known linear and nonlinear complementarity problem (LCP and NCP) formulations. In this paper, we compare these two formulations, showing through a set-based argument that the formulations are equivalent. Second, we provide simple examples to compare the most common approaches for solving these formulations. The prox formulation is solved by fixed-point iteration while the complementarity formulation is solved by a pivoting scheme, known as Lemke's algorithm. The well-known paradox of PAINLEVÉ is used in a case where two solutions exist to illustrate that the fixed-point scheme can fail while the pivoting scheme will succeed.

I. INTRODUCTION

In robotics, we are interested in the analysis of tasks involving intermittent contact, e.g., product assembly and optimal control of walking machines. The equations of motion of such kinds of applications can be written as a system of differential equations with set-valued constraints q = Y (q) u ,

M (q) u = h (q, u, t) + W (q) λ ,

(q, u, λ, t) ∈ N .

We refer to t as the time, to q as the configuration of the system and to u as its velocity. The elements of λ are rigid contact reactions e.g. between the robot's hand and tool or between feet of a walker and floor. Rigid means that these contact reactions are represented by constraints: the set N restricts the possible values such that there is no interpenetration and that friction forces maximize power dissipation at sliding contacts. As a consequence, the acceleration function u exhibits jump discontinuities over time. Further, impact equations have to be added in case of closing contacts, e.g. during stepping and grasping. The vector h contains all smooth external, internal and gyroscopic forces, which are functions of q, u and explicitly of t. It also includes compliant contacts. The symmetric and positive definite mass matrix M , the velocity kinematic map Y , and the constraint Jacobian W depend only on q.

In most physics engines used in robotic simulations, a linearized complementarity formulation has been used to represent the constraints N , and a pivoting algorithm has been used to solve the equations of motion directly. The prox formulation, first stated at the beginning of the 1990s [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF], is said to possess better numerical properties in practical simulation applications such as robotic manipulation tasks, the operation of continuously variable transmissions and electrical circuit breaker systems [START_REF] Pfeiffer | Mechanical system dynamics[END_REF], [START_REF] Transeth | 3-D snake robot motion: nonsmooth modeling, simulations, and experiments[END_REF], [START_REF] Acary | Nonsmooth modeling and simulation for switched circuits[END_REF]. We analyze the prox and complementarity formulations showing their mathematical equivalence. Further, we study the qualitative performance of typical solvers through two simple examples.

II. COMPLEMENTARITY FORMULATION

Joints of a robot system are represented by so-called bilateral constraints

B C := {(g B , λ B ) ∈ IR × IR | g B = 0} . (4)
Each such joint is physically enforced (i.e., g B is held to zero) by a contact force λ B . Note that the bilateral condition ( 4) is expressed at configuration level, but it can also be written at the velocity level when g B = 0, by replacing g B = 0 in (4) with ġB = 0. Similarly the condition can be written at the acceleration level when ġB = g B = 0, by replacing g B = 0 by gB = 0.

Unilateral contacts can form and detach intermittently. Such behavior is represented by SIGNORINI-FICHERAconditions [START_REF] Pfeiffer | Mechanical system dynamics[END_REF]

U C := {(g U , λ U ) ∈ IR × IR | 0 ≤ g U ⊥ λ U ≥ 0 } (5)
where g U denotes the normal distance and the symbol ⊥ implies orthogonality (i.e g U λ U = 0) to the contact force λ U . These conditions can also be formulated at the velocity or acceleration level as described for g B . The respective forcedisplacement laws are shown in Figs. 1(a For both bi-and unilateral constraints, dry friction can be considered. In order to establish COULOMB's law, the force of a single contact is decomposed in a positive component λ N ∈ {λ B , λ U } normal to the contact tangent plane and tangential components λ T in the tangent plane. When the relative velocity ġT at the contact point is nonzero, COULOMB's friction law is given as follows:

T C (λ N ) := { ( ġT , λ T ) ∈ IR 2 × IR 2 | ġT = 0 ⇒ λ T ≤ µ|λ N | , ġT = 0 ⇒ λ T = -ġT ġT µ|λ N | } (6)
where µ > 0 is the coefficient of friction. In the degenerate case defined by µ = 0 or λ N = 0, these conditions are not needed in the problem formulation. If the contact is sticking (i.e. ġT = 0) then in [START_REF] Ferris | Interfaces to PATH 3.0: design, implementation and usage[END_REF], ġT is replaced by gT to calculate forces. One can easily show that in the non-degenerate cases T C (λ N ) equals the complementarity formulation:

T σ (λ N ) := { ( ġT , λ T ) ∈ IR 2 × IR 2 , σ ∈ IR | 0 = µ|λ N | ġT + σλ T , 0 ≤ σ ⊥ µ|λ N | -λ T ≥ 0 } . (7)
In the planar case, COULOMB's law is depicted in Fig. 2.

λ T ġT +µ|λ N | -µ|λ N |
Fig. 2. Planar COULOMB friction. The vertical part of the graph reveals the set-valued nature of the constraint, i.e., when sliding velocity is zero, the friction force may take on any value in (-µ|λ N |, +µ|λ N |).

Finally, the complementarity formulation of a system of bodies in three dimensions with bilateral and unilateral constraints with COULOMB friction is given by (1)-( 3) and all of the equations and inequalities in (4), ( 5) and [START_REF] Acary | On the equivalence between complementarity systems, projected systems and unilateral differential inclusions[END_REF]. These equations represent the instantaneous equations of motions whose form is known as a nonlinear differential complementarity problem (DCP). Trajectories of the system constructed during simulation are obtained by solving a timesequence of nonlinear complementarity problems (NCPs). These can be solved reliably with the PATH solver [START_REF] Trinkle | On dynamic multi-rigidbody contact problems with Coulomb friction[END_REF], [START_REF] Ferris | Interfaces to PATH 3.0: design, implementation and usage[END_REF].

III. FORMULATION WITH PROX FUNCTION

A perhaps more appropriate formulation can be obtained using convex analysis. For x ∈ IR n and n ∈ IN, we define

C x 2 x 1 = prox C (x 1 ) prox C (x 2 )
Fig. 3. prox function for a convex set. the closest point, i.e. the proximal point, to a given convex set C ⊂ IR n :

prox C (x) = arg min x * ∈C x -x * . (8)
In Fig. 3, it is illustrated that the proximal point is the point itself when it is in C (cf. x 1 ) or the closest point on the boundary of C otherwise (cf. x 2 ). The contact conditions can be formulated as:

B P := { (g B , λ B ) ∈ IR × IR | f B (λ B , g B ) := λ B -prox CB (λ B -r g B ) = 0 } , (9a)
U P := { (g U , λ U ) ∈ IR × IR | f U (λ U , g U ) := λ U -prox CU (λ U -r g U ) = 0 } , (9b)
T P (λ N ) := { ( ġT , λ T ) ∈ IR 2 × IR 2 | f T (λ T , ġT ) := λ T -prox CT (λN ) (λ T -r ġT ) = 0 } (9c)
where the corresponding convex sets are specified by:

C B = IR , (10)
C U = {x ∈ IR | x ≥ 0} , (11)
C T (y) = {x ∈ IR 2 | x ≤ µ |y|} (12)
with y ∈ IR. The equations are given at the level of the lowest-order derivatives, configuration in normal directions N ∈ {B, U } and velocity in tangential directions. Higherorder derivative formulations can be obtained simply by increasing the number of dots over the g's. An independent auxiliary parameter r ∈ IR + appears in [START_REF] Cottle | The linear complementarity problem, ser. Classics in applied mathematics[END_REF]. This parameter is used to match units in the equations, but more importantly, its value strongly affects the stability and rate of convergence when finding a solution by fixed-point iteration. It is also important to notice that the solution set of ( 9) is nonsmooth. Typical spatial plots of f B , f U and f T are given in Fig. 4 for r = 1.0. The solution set is the intersection of the dark surface with the light plane. Notice the corner(s) for f U and f T . In the tangential cases, µ = 0.1 and λ N = 1 have been used. The intersection with the zero-constraint and the solution sets are reminiscent of the force laws from Figs. 1 and2. Mathematically this is proved in Section IV.

0 0 0 g B λ B (a) Bilateral 0 0 0 g U λ U (b) Unilateral 0 0 0 ġT λ T (c) Friction

IV. EQUIVALENCE BETWEEN COMPLEMENTARITY AND PROX

Complementarity and prox formulations represent the same constraints

B C = B P , U C = U P , T C = T P ( 13
)
as has been shown in [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF] using variational equations. Therefore, numerical solution strategies can be devised to make use of their advantages while avoiding their disadvantages. We show the equivalences for the case of the lowest-order derivatives starting basically from the set definitions like in [START_REF] Acary | On the equivalence between complementarity systems, projected systems and unilateral differential inclusions[END_REF], [START_REF] Förg | Mehrkörpersysteme mit mengenwertigen Kraftgesetzen : Theorie und Numerik, als Manuskript gedruckt[END_REF] but proceeding more directly.

Bilateral constraint

Starting with the prox formulation (g B , λ B ) ∈ B P , it is obvious that the argument (λ Br g B ) of prox CB is always an element of the corresponding convex set C B = IR. This yields:

0 = λ B -λ B + r g B ( 14
)
and finally g B = 0 as well as

(g B , λ B ) ∈ B C .
In the complementarity formulation, it is λ B ∈ IR and g B = 0. The relations both λ B = prox CB (λ Br g B ) and (g B , λ B ) ∈ B P are obvious.

Unilateral constraint

Starting with the prox formulation (g U , λ U ) ∈ U P , one has to distinguish two cases.

1

) λ U -r g U ∈ C U or λ U -r g U ≥ 0: This yields 0 = λ U -λ U + r g U and g U = 0 as well as λ U ≥ 0. 2) λ U -r g U / ∈ C U or λ U -r g U < 0: This yields λ U = prox CU (λ U -r g U ) = 0 and g U > 0.
This is a representation of the unilateral corner law by a complementarity formulation

(g U , λ U ) ∈ U C .
In a complementarity formulation (g U , λ U ) ∈ U C , there are two cases for the branches of the unilateral corner law.

1) g U = 0 with λ U ≥ 0:

This yields λ U -r g U = λ U ≥ 0 or λ U -r g U ∈ C U
and the validity of the prox formulation (g U , λ U ) ∈ U P in the specific case. 2) g U > 0 and λ U = 0:

This yields λ Ur g U < 0 or λ Ur g U / ∈ C U and the validity of the prox formulation (g U , λ U ) ∈ U P in the specific case.

COULOMB friction

Let µ ≥ 0 and λ N ∈ IR.

Starting with the prox formulation ( ġT , λ T ) ∈ T P (λ N ), one has to distinguish two cases.

1

) λ T -r ġT ∈ C T (λ N ) or λ T -r ġT ≤ µ |λ N |: This yields 0 = λ T -λ T + r ġT and ġT = 0. One has λ T ≤ µ |λ N | and finally ( ġT , λ T ) ∈ T C (λ N ). 2) λ T -r ġT / ∈ C T (λ N ) or λ T -r ġT > µ |λ N |: This yields: λ T = λ T -r ġT λ T -r ġT µ |λ N | = - r ġT λ T -r ġT -µ |λ N | k ġT ġT µ |λ N | . ( 15
)
If one were in the case ġT = 0, this would be a contradiction because

λ T > µ |λ N | ∧ λ T = λ T λ T µ |λ N | . ( 16
)
One has to discuss the case ġT = 0. Then, the contact is sliding and we must have

λ T -r ġT = -ġT ġT µ |λ N | -r ġT ġT ġT = µ |λ N | + r ġT > µ |λ N | (17)
and λ Tr ġT / ∈ C T (λ N ). Finally:

λ T = -ġT ġT µ |λ N | = - µ |λ N | + r ġT µ |λ N | + r ġT ġT ġT µ |λ N | = λ T -r ġT λ T -r ġT µ |λ N | = prox CT (λN ) (λ T -r ġT ) (18)
which means ( ġT , λ T ) ∈ T P (λ N ).

V. EXAMPLES

In this section, we provide several examples that serve to illustrate the connections between the LCP and prox formulations graphically and numerically.

A. Point mass on frictional plane

Is the prox formulation numerically most appropriate? The discussion of this question for a simple planar example is the goal of the following section. The planar restriction is not fundamental. It was imposed only to allow clear graphical illustrations of the solution process.

The instantaneous equations of motion of a point mass in contact with a rigid linear surface (see Fig. 5) can be written as the following system of differential equations:

qx = u x , (19)
qy = u y , (20)
m ux = h x + λ T , (21)
m uy = h y + λ U (22)
where h y is the component of the external force applied to the particle in the direction of the outward normal (the ydirection) of the surface and h x is the component tangential to the surface in the positive x-direction. Notice also that for this problem, the x-and y-coordinates of the particle are the displacement functions, q x = g T and q y = g U . To formulate the dynamics with contact constraints, one can augment the differential equations with contact conditions represented at the acceleration level; ( 5) and ( 6) for a LCP formulation or (9b) and (9c) for a prox formulation:

( uy , λ U ) ∈ U C/P , (23)
( ux , λ T ) ∈ T C/P (λ U ) (24)
where uy = gU and ux = gT .

Fixed-point iteration and prox formulation

Substituting ( 21) and ( 22) into the right-hand-sides of (9b) and (9c) (at the acceleration level) yields the fixed-point equations:

λ U = prox CU 1 - r m λ U - r m h y , (25)
λ T = prox CT (λU ) 1 - r m λ T - r m h x . (26)
The simultaneous solution of these equations yields the contact force and the acceleration of the particle. Notice that the first equation is completely independent of the second, but the second is coupled to the first through the set C T (λ U ).

Consider the problem of finding the normal component of the contact force for when the particle is in contact with the plane (q y = 0) with zero normal velocity (u y = 0) and with the external force pushing the particle against the plane. Fig. 6 illustrates the progress of the fixed-point iterations (dashed curve) with r in the range 0 < r < 2m, starting from an (infeasible) initial guess λ 0 U to the unique solution λ * U = -h y . Graphically, each iteration consists of two steps: first project the current solution estimate λ k U vertically to the prox function (dotted curve), then project horizontally from that point to the line with unit slope, f (λ U ) = λ U . In this picture, as r is varied, the slope of the non-constant part of the prox function changes while the intersection point (i.e. the solution) remains at λ * U = -h y . Notice that when r = m, the prox function becomes a horizontal line and only one iteration is needed to find the solution. The other case of interest for the normal contact force problem arises when the external force acts to separate the particle from the plane (i.e. h y > 0). In this case, the prox function and the line of unit slope have a unique intersection point at λ * U = 0, to which the iterates converge if 0 < r < 2m.

Next consider the tangential subproblem of the particleplane system where the particle is on the plane and sticking (i.e. g U = 0, ġU = 0 and ġT = 0) at the current time instant. Assume that the external force acting on the particle has its normal component towards the rigid surface (i.e. h y < 0). Further assume that the solution is already known. Then, the prox formulation is given as: In practice, the numerical solution proceeds by interleaving fixed-point iterations between the normal and tangential prox functions. For the particle-on-plane problem, the graph in Fig. 6 is unchanged, hence the normal force iterations are unchanged. However, the graph in Fig. 7 is dependent on the result of the preceding iteration for the normal force. The tangential iterations are guaranteed to converge to the correct value as the normal direction iterations converge, because the iterations for the tangential component converge from any starting guess.

λ * U = -h y . . . . . . λ 0 U λ 1 U λ 2 U λ k U λ U λ ⋆ U -hyr/m -hy λ U = p r o x C U p r o x C U
λ T = prox CT (λU ) (λ T (1 -r/m) -h x r/m).
To summarize despite the trivial nature of this system, the analysis suggests that solving multibody dynamics timestepping subproblems by fixed-point iteration of the prox formulation requires one to be careful about the choice of the parameter r. At this point, it is not clear how to choose r to guarantee convergence for general problems. We conjecture that the choices of r yielding global convergence to a physically correct solution of the time-stepping subproblem are likely to be related to eigenvalues of the Delassus matrix of the contact configuration [START_REF] Förg | Mehrkörpersysteme mit mengenwertigen Kraftgesetzen : Theorie und Numerik, als Manuskript gedruckt[END_REF].

Complementarity formulations at the acceleration level

Let us replace λ T with the sum of its positive and negative parts: λ T = λ f 1λ f 2 where λ f 1 ≥ 0 is the friction force component acting in the positive x-direction and λ f 2 ≥ 0 acts in the (-x)-direction. Now, the equations of motion (21-24) become linear in the unknown forces, accelerations and in the auxiliary variable σ:

m ux = h x + λ f 1 -λ f 2 , (27)
m uy = h y + λ U , (28)
0 ≤ λ U ⊥ uy ≥ 0 , (29)
0 ≤ λ f 1 ⊥ ux + σ ≥ 0 , (30)
0 ≤ λ f 2 ⊥ -ux + σ ≥ 0 , (31)
0 ≤ σ ⊥ µλ U -λ f 1 -λ f 2 ≥ 0 . (32)
Also at a solution, we have σ = | ux | [START_REF] Trinkle | On dynamic multi-rigidbody contact problems with Coulomb friction[END_REF]. The above mixed LCP can be reformulated as a standard LCP of size four by eliminating ux and uy .

The numerical solver has to test pivots until a solution is found. In the worst case, this is 2 n pivots where n is the number of unknowns. However, the average time required by Lemke's algorithm is approximately cubic in n.

B. PAINLEVÉ's paradox

PAINLEVÉ's paradox refers to problems in multibody dynamics in which the solution of the instantaneous dynamics is not unique. One such "paradox" arises in the planar case of a rigid object translating in contact with a rigid horizontal surface. The number of solutions of the acceleration-level formulation depends on the values of the mass, moment of inertia, coefficient of friction, and the angle and distance from the contact point to the center of mass of the object (see Fig. 8). In the case considered here, a gravitational force is assumed to act downward through the center of mass perpendicular to the contact tangent, i.e. h y = -mg. Since m and g are positive, h y is negative. The lateral external force h x is zero.

The dynamic equations are:

qx = u x , (33)
qy = u y , (34)
qθ = u θ , (35)
m ux = λ T , (36)
m uy = λ U + h y , (37)
J uθ = L 2 (λ T sin(q θ ) -λ U cos(q θ )) ( 38
)
where J is the moment of inertia about an axis perpendicular to the plane of motion and through the mass center. L/2 is the distance from the contact point to the center of mass.

Let us assume that at the current time instant, the object is translating toward the left and leaning toward the right (i.e. u x < 0, u y = u θ = 0 and 0 < q θ < π/2). The direction of the friction force is determined from velocity kinematics to be in the positive x-direction, because the contact point is initially sliding. The unilateral contact constraint must still be evaluated at the acceleration level because u y = 0 and uy = 0. Thus, the contact constraints for this problem are:

(g U , λ U ) ∈ U C/P , (39)
( ġT , λ T ) ∈ T C/P (λ U ) . (40)
The kinematic quantities, g U and g T , and their relevant time derivatives are:

g T = q x - L 2 cos(q θ ) , (41)
g U = q y - L 2 sin(q θ ) , (42)
ġT = u x + L 2 u θ sin(q θ ) , (43)
gU = uy + L 2 u 2 θ sin(q θ )uθ cos(q θ ) .

Note that (40) with ux < 0 implies that λ T = µλ U . This allows us to reduce the formulation of the problem to an LCP of size one:

0 ≤ λ U ⊥ Gλ U + b ≥ 0 (45)
where G = 1 m + L 2 cos(q θ ) 4J

(cos(q θ )µ sin(q θ )) and b = -g. Note that for physically plausible values of m, J, L, q θ , and µ, G can be made negative. By choosing the sign of the gravitational acceleration, b can be made negative (the force acts downward) or positive (the force acts upward). For this problem, it is known that when G is strictly positive (Figs. 9 and 10), there is a unique solution for any b [START_REF] Cottle | The linear complementarity problem, ser. Classics in applied mathematics[END_REF]. Also if G is strictly negative then if b is strictly positive, there are two solutions (Fig. 12), and if b is strictly negative (Fig. 11), there is no solution.

Using the initial conditions, the prox formulation can be written as a single scalar fixed-point problem:

λ U -prox CU (λ U -r gU ) = 0. (46)
Figs. 9-12 show that there is a unique solution when G > 0 and either zero or two solutions when G < 0 by intersecting the prox function (dotted curve) with the line of unit slope. From the point of view of solution existence, the prox formulation completely agrees with complementarity theory. A difference appears when attempting to find a solution via the fixed-point iteration scheme. In the case with g and G negative, for any initial guess λ 0 U less than the second solution (λ * U ) 2 , the iterates will converge to the first solution (λ * U ) 1 = 0. Iterations beginning with λ 0 U > (λ * U ) 2 will diverge. While a pivoting algorithm solving the LCP formulation could terminate at either solution (or be modified to find both), the fixed-point scheme will only find (λ * U ) 1 . Another interesting difference is that in the case with g positive and G negative, fixed-point iterations are divergent from any starting guess. Thus, one sees that a pivoting method is superior in that it will be able to recognize the case of solution non-existence. Whereas, the fixed-point scheme will not be able to determine if divergence is due to a poor initial guess or solution non-existence.

VI. CONCLUSIONS AND FUTURE WORKS

In this paper, we have shown the equivalence of two popular formulations of nonsmooth dynamic systems. Prox formulations lend themselves to solution by fixed-point iteration or root-finding schemes for nonsmooth equations. Complementarity formulation can be approached in a number of ways including pivoting methods. While some have reported faster solution times with fixed-point schemes, the simple problems discussed here shows that fixed-point iterations can diverge when a solution exists and that this solution will be found by a pivoting scheme. However, fixed-point schemes are worth pursuing further to explore the exploitation of finegrained parallelism in the solution process. While regularization of complementarity formulations have allowed one to exploit parallel implementation of quadratic program solvers, they do so at the expense of the introduction of a-physical effects. Parallel implementations of fixed-point schemes will not require such model adjustments. For the analysis, general approaches will have to be worked out in detail and applied to use case examples.

  ) and 1(b).

Fig. 1 .

 1 Fig. 1. Force laws for bilateral and unilateral contacts.

Fig. 4 .

 4 Fig. 4. prox functions for contacts and friction with r = 1.0.

  λ T = µ |λ N |. Comparing the expressions in (15), one sees that k must equal 1. Therefore, the only solution is given by λ T = -ġT ġT µ |λ N | and ( ġT , λ T ) ∈ T C (λ N ). If one has a representation by a complementarity formulation ( ġT , λ T ) ∈ T C (λ N ), again there are two cases. 1) ġT = 0 and λ T ≤ µ |λ N |: This yields λ Tr ġT ∈ C T (λ N ) and trivially ( ġT , λ T ) ∈ T P (λ N ). 2) ġT = 0 and λ T = -ġT ġT µ |λ N |: This yields:

Fig. 5 .

 5 Fig. 5. Configuration of point mass in contact with a surface.

Fig. 6 .

 6 Fig.6. Progression of iterates for finding the normal contact force.

Fig. 7 Fig. 7 .

 77 Fig. 7. Progression of iterates for finding the tangential contact force.

Fig. 8 .

 8 Fig. 8. PAINLEVÉ problem with g > 0.

Fig. 9 .Fig. 10 .

 910 Fig. 9.Globally convergent unique solution: g, G > 0. Contact is maintained.

Fig. 11 .Fig. 12 .

 1112 Fig. 11. No solution exists: g > 0, G < 0. The acceleration-based model fails.


Marco Pizzolato

Timothé Boutelier

Rutger H J Fick

Rachid Deriche

Elucidating Dispersion Effects in Perfusion MRI by Means of Dispersion-Compliant Bases

Keywords: perfusion, dispersion, delay, deconvolution, bases

come   L'archive ouverte pluridisciplinaire

INTRODUCTION

The perfusion in the brain can be characterized with Dynamic Susceptibility Contrast Magnetic Resonance Imaging (DSC-MRI). A bolus of paramagnetic agent (PA) is injected into the subject's vascular system and the corresponding concentration time-curve is measured for each voxel. In standard practice, the concentration curve C ts (t) of a tissular voxel is deconvolved with that measured in an arterial area C a (t) to obtain the residue function R(t) that characterizes the local perfusion. From this, the hemodynamic parameters such as the blood flow (BF) and mean transit time (MTT) are computed. However, the bolus of PA may undergo dispersion after injection, causing the calculated residue function to change according to macro-vascular properties [START_REF] Calamante | Delay and dispersion effects in dynamic susceptibility contrast MRI: simulations using singular value decomposition[END_REF]. Dispersion can be present in healthy subjects [START_REF] Østergaard | Cerebral blood flow measurements by magnetic resonance imaging bolus tracking: comparison with [ 15 o]h 2 o positron emission tomography in humans[END_REF] but it can also be caused by specific pathological conditions, such as the presence of † The author expresses his thanks to Olea Medical and to the Provence-Alpes-Côte d'Azur Regional Council for providing grant and support. a steno-occlusive disease in the artery [START_REF] Calamante | Estimation of bolus dispersion effects in perfusion MRI using image-based computational fluid dynamics[END_REF]. Moreover, dispersion leads to an erroneous quantification of the true perfusion deficit [START_REF] Calamante | Delay and dispersion effects in dynamic susceptibility contrast MRI: simulations using singular value decomposition[END_REF]. Indeed, the dispersion of the injected PA bolus causes a modification in the shape of the actual arterial input, which induces a non-monotonic shape in the computed R(t) [START_REF] Calamante | Delay and dispersion effects in dynamic susceptibility contrast MRI: simulations using singular value decomposition[END_REF]. This effect is described by an additional unknown convolution kernel -a dispersion kernel -called the vascular transport function (VTF). Another issue in perfusion data processing is the possible presence of a time delay between the tissular concentration time-curve and the measured arterial input. This is a confounding factor for dispersion characterization [START_REF] Willats | Improved deconvolution of perfusion MRI data in the presence of bolus delay and dispersion[END_REF][START_REF] Mehndiratta | Modeling and correction of bolus dispersion effects in dynamic susceptibility contrast MRI: Dispersion correction with cpi in DSC-MRI[END_REF]. Assessment of dispersion is a challenging task but is crucial for decoupling perfusion and macro-vascular influences on the residue function estimated via deconvolution.

Dispersion effects have been tackled by few techniques. One allows to classify voxels where either dispersion or delay is present [START_REF] Willats | Improved deconvolution of perfusion MRI data in the presence of bolus delay and dispersion[END_REF]. The VM+VTF and CPI+VTF [START_REF] Mehndiratta | Modeling and correction of bolus dispersion effects in dynamic susceptibility contrast MRI: Dispersion correction with cpi in DSC-MRI[END_REF] aim at quantifying perfusion and dispersion separately, but assume a model for the VTF together with the prior information required for the optimization procedure. The CPI 0 [START_REF] Mehndiratta | Modeling and correction of bolus dispersion effects in dynamic susceptibility contrast MRI: Dispersion correction with cpi in DSC-MRI[END_REF] retains a "modelfree" nature but fixes the initial value of the residue function to 0, thus preventing an optimal solution in the absence of dispersion. The block-circulant Singular Value Decomposition (oSVD) [START_REF] Wu | Tracer arrival timinginsensitive technique for estimating flow in MR perfusionweighted imaging using singular value decomposition with a blockcirculant deconvolution matrix[END_REF], despite not explicitly aiming at dispersion characterization, allows the reconstruction of dispersed residue functions with no assumptions, and performs comparably to other methods for relative cerebral blood flow estimation in the presence of dispersion [START_REF] Mehndiratta | Modeling and correction of bolus dispersion effects in dynamic susceptibility contrast MRI: Dispersion correction with cpi in DSC-MRI[END_REF].

We propose to characterize bolus dispersion in DSC-MRI data without any a priori assumption on the residue function and vascular transport function. We perform deconvolution with dispersion-compliant bases [START_REF] Pizzolato | Perfusion MRI deconvolution with delay estimation and non-negativity constraints[END_REF] obtaining a representation of the residue function which allows the decoupling of time-delay and dispersion effects. We derive an analytic representation of the dispersed residue function, to be used for in silico experiments that validate the reliability of the adopted approach against oSVD [START_REF] Wu | Tracer arrival timinginsensitive technique for estimating flow in MR perfusionweighted imaging using singular value decomposition with a blockcirculant deconvolution matrix[END_REF] at different levels of dispersion. We introduce the use of shape parameters, the dispersion time and index, to detect presence and amount of dispersion in the shape of the residue function. We show the effectiveness of the approach with synthetic experiments and on stroke MRI data.

DISPERSED RESIDUE FUNCTION

According to the indicator-dilution theory [START_REF] Meier | On the theory of the indicator-dilution method for measurement of blood flow and volume[END_REF], the tissular concentration time-curve C ts (t) in a voxel is expressed as the convolution between the arterial input C a (t) and the unknown residue function R(t), which expresses at each time t the residual amount of PA in the voxel

C ts (t) = BF • C a ⊗ R(t) ( 1
)
where BF is the unknown cerebral blood flow parameter.

In presence of dispersion, the true arterial input is the result of the convolution of C a (t) with the vascular transport function V T F (t), which represents the probability density of the transit times t from the arterial measurement location to the voxel where C ts (t) is measured [START_REF] Calamante | Delay and dispersion effects in dynamic susceptibility contrast MRI: simulations using singular value decomposition[END_REF]. Therefore eq. ( 1) becomes

C ts (t) = BF • [C a ⊗ V T F (t)] ⊗ R(t) = BF • C a ⊗ [V T F ⊗ R(t)] (t) = BF • C a ⊗ R d (t) (2)
where R d (t) is the dispersed residue function resulting from the use of the associative property.

In order to perform in silico experiments, we derive a formulation of R d (t) = V T F ⊗ R(t). Among the several analytic formulations of R(t), the bi-exponential shows the best in vivo fitting performance in normal and Diffusion Positive tissue [START_REF] Mehndiratta | Modeling the residue function in DSC-MRI simulations: Analytical approximation to in vivo data[END_REF]. It includes a fast and a slow flowing capillary components

R bi-exp (t) = f • e -τ F t + (1 -f ) • e -τ S t (3)
where τ F and τ S are the fast and slow time-rates respectively, and f is the relative weight of the fast component. We specialize R d (t) for the bi-exponential model of eq. ( 3) in the case of a mono-exponential dispersion kernel V T F (t) = βe -βt [START_REF] Calamante | Delay and dispersion effects in dynamic susceptibility contrast MRI: simulations using singular value decomposition[END_REF][START_REF] Calamante | Estimation of bolus dispersion effects in perfusion MRI using image-based computational fluid dynamics[END_REF][START_REF] Willats | Improved deconvolution of perfusion MRI data in the presence of bolus delay and dispersion[END_REF] obtaining

R d (t) = -β (β -τ F )(β -τ S ) [(f (τ F -τ S ) + β -τ F )e -βt + f (τ S -β)e -τ F t + (f (β -τ F ) -β + τ F )e -τ S t ] (4)
which holds for β > 0 and, similarly to what observed for the case of the mono-exponential residue function [START_REF] Calamante | Delay and dispersion effects in dynamic susceptibility contrast MRI: simulations using singular value decomposition[END_REF], depends on the original mean transit time

M T T = f /τ F +(1-f )/τ S .
The dispersion amount induced by V T F (t) on the calculated residue function can be expressed by the vascular mean transit time, which in the case of eq. ( 4) is

M T T v = β -1 [1].
Its influence is illustrated in Fig. 1, where the top image shows that dispersion effects increase with M T T v , being a pure biexponential decay for M T T v = 0, whereas the bottom image shows that for a specific value of M T T v > 0 dispersion effects increase also with M T T . We propose to quantify the deviation of the residue function from a purely decaying shape as the difference between the integrals of the decreasing and increasing parts of the curve, normalized by the total area

τ M T T τ δ MTT v τ MTT=cost MTT v =cost>0 0 t (s) 0 t (s) R(t) R(t) t max
δ = ∞ tmax R(t)dt - tmax τ R(t)dt ∞ τ R(t)dt (5)
where t max is the time of maximum and τ the time-delay as shown in Fig. 1. The dispersion index indicates the absence of dispersion when δ = 1 or its presence when δ < 1. In addition we note that Fig. 1 illustrates that the difference between t max and τ , henceforth called dispersion time τ δ = t max -τ , increases with M T T and M T T v . We then propose to calculate dispersion time τ δ and index δ as shape descriptors of the residue function sensitive to the presence and amount of dispersion.

The next section describes the deconvolution strategy to calculate these parameters and the experimental setup to test dispersion characterization reliability in silico and in vivo.

METHOD AND EXPERIMENTAL RESULTS

The detection of dispersion is carried out by deconvolving the measured C a (t) and C ts (t) while representing the residue function by means of dispersion-compliant bases [START_REF] Pizzolato | Perfusion MRI deconvolution with delay estimation and non-negativity constraints[END_REF]

R(t) = Θ(t -τ ) N n=1 [a n + b n (t -τ )]e -αn(t-τ ) (6)
where N is the maximum basis order and Θ(t) is the Heaviside step function assuming Θ(0) = 1. The deconvolution problem is solved linearly for a n , b n subject to a n , b n ≥ 0 to avoid negative solutions, with α n estimated a priori and τ determined via grid search [START_REF] Pizzolato | Perfusion MRI deconvolution with delay estimation and non-negativity constraints[END_REF].

We perform in silico experiments to test the reliability of the adopted deconvolution strategy against oSVD [START_REF] Wu | Tracer arrival timinginsensitive technique for estimating flow in MR perfusionweighted imaging using singular value decomposition with a blockcirculant deconvolution matrix[END_REF]. The measured arterial input C a (t) is simulated in range [0, 90]s with time step ∆t = 1s as a gamma-variate function [START_REF] Calamante | Delay and dispersion effects in dynamic susceptibility contrast MRI: simulations using singular value decomposition[END_REF]. The tissular concentration C ts (t) is obtained via eq. ( 2) using the derived dispersed kernel in eq. ( 4) when β > 0. Gaussian noise with zero mean is added to obtain concentration timecurves as previously reported [START_REF] Pizzolato | Perfusion MRI deconvolution with delay estimation and non-negativity constraints[END_REF] with SN R = 50 [START_REF] Mehndiratta | Modeling and correction of bolus dispersion effects in dynamic susceptibility contrast MRI: Dispersion correction with cpi in DSC-MRI[END_REF]. Deconvolution is performed with the dispersion-compliant bases (N = 20) searching for τ in range [-5, 15]s [START_REF] Pizzolato | Perfusion MRI deconvolution with delay estimation and non-negativity constraints[END_REF] and with oSVD [START_REF] Wu | Tracer arrival timinginsensitive technique for estimating flow in MR perfusionweighted imaging using singular value decomposition with a blockcirculant deconvolution matrix[END_REF]. The effective blood flow BF * is obtained as the maximum value of the estimated residue function R(t), the blood volume BV as the ratio between the time integrals of C ts (t) and C a (t), and the mean transit time as the ratio M T T * = BV /BF * . Results are studied as dispersion effects increase for increasing values of M T T v , assuming no dispersion for M T T v = 0, i.e., data is generated with eq. ( 3). We simulate the concentration curves 100 times for each combination of dispersion, M T T v ∈ [0 : 1 : 10]s, delay τ ∈ [-5 : 1 : 5]s, original value of blood flow, BF ∈ [20 : 10 : 60]ml/100g/min with blood volume BV ≈ 2% and original mean transit time M T T ≈ 4s (τ F = 0.34, τ S = 0.025, f = 0.97). Results for τ = 2s and original BF = 30ml/100g/min are shown in the left column of Fig. 2 and mean relative or absolute errors among all the tested combinations are shown in the right column. We note that deconvolution via dispersion-compliant bases results in more accurate estimates than oSVD. The mean relative errors on BF * estimates are lower than 20% and uncorrelated with the underlying amount of dispersion (M T T v ). As in the case of oSVD, the performance in estimating M T T * ameliorates as dispersion increases with an improved accuracy. The mean absolute precision on t max always falls below the sampling time ∆t = 1s and below 2s for τ δ , as depicted from the green bars in Fig. 2, allowing for dispersion time estimation.

We then applied the compared techniques on 256 × 256 × 15 stroke MRI data, ∆t = 1.5s. The techniques' setups are identical to those of in silico experiments. We compute the hemodynamic parameters BF * , BV and M T T * and remove large vessels and cerebrospinal fluid (CSF) from the results by setting voxels with BV > 3% to zero. The data reveals an infarcted region in the right hemisphere and iso-perfused tissue. We report the maps of t max in Fig. 3 where oSVD results confirm the overestimation tendency shown in the last row of Fig. 2, particularly in the infarcted region. Moreover, our approach accounts for both the contribution of the timedelay τ and dispersion time τ δ which are shown in the top row of Fig. 4. Indeed, the dispersion-compliant bases encompass delayed and dispersed shapes allowing to distinguish between the two different contributions. Although the infarcted region is distinguishable in both the t max and τ maps, the proposed dispersion time τ δ map improves the contrast between the inside and the outside of the infarcted region itself. The map of the dispersion index proposed in eq. ( 5) is shown in the bottom left corner of Fig. 4. Non-dispersed regions are colored in red (δ = 1), whereas green and blue regions denote moderate and high dispersion respectively. We distinguish two areas within the infarcted region: an upper area mainly characterized by the absence of dispersion, and a lower one dominated by high dispersion (blue color). Moderately dispersed shapes dominate the iso-perfused region which may be caused by natural dispersion effects in the data [START_REF] Østergaard | Cerebral blood flow measurements by magnetic resonance imaging bolus tracking: comparison with [ 15 o]h 2 o positron emission tomography in humans[END_REF] and noise effects. The dispersion index δ is a relative shape parameter since it characterizes the residue function regardless of its actual size. In order to understand the actual spread of the dispersed residue function, it is interesting to relate its δ value with its measure of the aperture given by the M T T * . We show the correlation term ρ 1-δ,M T T = (1 -δ)M T T * /M T T * max in the bottom right corner of Fig. 4. As expected this map highly correlates with the above dispersion time map τ δ , [r ≈ 0.95, p < 0.001]. Indeed, as previously described and illustrated in Fig. 1, the dispersion time depends on both dispersion effects -M T T v and its influence on δ -and the mean transit time. We further note how both the τ δ and ρ 1-δ,M T T maps are insensitive to regions where no dispersion is detected -the red areas in the δ map -being non-zero only in complementary regions.

Overall our approach allows low relative error of BF * , M T T * estimates, sub-second precision for t max , and reliable estimation of delay τ and dispersion time τ δ regardless the amount of dispersion. Furthermore, the proposed index δ provides new insights on dispersion, and its correlation with M T T * is consistent with τ δ .

CONCLUSION

Dispersion is a phenomenon that can be present in perfusionweighted data and its characterization is fundamental to assess the reliability of perfusion-related estimates while potentially revealing pathological conditions [START_REF] Willats | Improved deconvolution of perfusion MRI data in the presence of bolus delay and dispersion[END_REF]. In this work we propose a way of characterizing dispersion effects in vivo without making any assumption on the dispersion kernel. This is achieved by introducing shape parameters, such as the dispersion time and index, that are calculated on the residue function estimated via deconvolution by means of dispersion-compliant bases. The approach outperforms oSVD and yields reliable estimates specially when the underlying residue function has a dispersed shape, such as in the case of the derived formulation. Results on stroke MRI data confirm in vivo the influence of the mean transit time on the dispersed shape, showing consistency with the estimated dispersion time and index. The dispersion time is useful in the delimitation of the infarcted area whereas the dispersion index gives new insights on the shape of the residue function. The proposed strategy allows to decouple dispersion effects from the confounding effect of the delay. This opens for dispersion characterization, thus providing a better understanding of the tissue perfusion and vascular dynamic.

Fig. 1 .

 1 Fig.1. Influence of dispersion on the calculated residue function and time-delay τ between C ts (t) and C a (t). Top: the maximum decreases, and the time of the maximum t max and dispersion time τ δ increase with the vascular mean transit time M T T v . Bottom: t max , τ δ increase with M T T (curves are normalized).

Fig. 2 .

 2 Fig. 2. Estimated ( * effective) hemodynamic parameters for τ = 2s, BF = 30 ml/100g/min (left column), and overall errors (right column) as dispersion increases according to the vascular mean transit time M T T v (SN R = 50).

Fig. 3 .Fig. 4 .

 34 Fig. 3. The time of the maximum t max (s) of the residue function R(t) calculated via oSVD, and dispersion-compliant bases: these reduce overestimation as shown in Fig. 2.τ τ δ