Guillaume Andrieu

Christophe Alias

Laure Gonnord

Modular termination for general C Programs

Keywords: Static Analysis, Termination, Modularity, C programs Analyse statique, terminaison, modularité, programmes C

) and preliminary results shows the feasibility of the method.

Introduction

Although proving program termination is known to be undecidable, recent progress in program analysis made it possible to predict the termination of an always increasing class of sequential programs. The standard approach of the seminal paper [START_REF] Floyd | Assigning meaning to programs[END_REF] which consists in nding a function from the states of the program to some well-founded set, which strictly decreases at each program point, remains standard. In [START_REF] Alias | Multi-dimensional Rankings, Program Termination, and Complexity Bounds of Flowchart Programs[END_REF], we proposed a general algorithm to discover multidimensional ranking functions from owcharts programs, by means of the resolution of linear programing (LP) instances. However, our experiments showed two main problems to scale to larger C programs :

Although our C parser, namely C2fsm ( [START_REF] Feautrier | Accelerated Invariant Generation for C Programs with Aspic and C2fsm[END_REF]) was designed to handle a large part of C syntax, many syntactical variants are not handled. Most of them could be ignored to prove the termination.

Moreover, the size of the LP problems increases with the number of code lines, and thus quickly become intractable.

In this paper, we propose an eective modular interprocedural termination analysis, which relies on the previous method, but enables to analyse a broad range of C programs, and validate this approach on a large benchmark of the literature. Our method relies on classical methods from static analysis and compilation, such as slicing and summaries, but as far as we know, is the rst attempt in proving termination of a large set of C programs in a modular way.

The rest of the paper is organized as follows. In section 2, we introduce our motivation of a challenging program of the literature. In section 3 we quickly introduce our notation and theoretical foundations. In section 4 we describe our method, and we evaluate it on a large bench of middle-sized programs in section 6. We end with related works (section 7) and a conclusion.

RR n°8166 2 Motivating example

In this paper we detail our method on an implementation of the merge sort of an array. The code is taken from [START_REF] Colón | Practical methods for proving program termination[END_REF] and is depicted in Figure 1. For sake of lisibility, we drawed boxes around innerloops, and commented the end of outer loops. int main () { int n , i , j , k , l , t , h ,m, p , q , r ; Calling our tool suite (Figure 2) detailed in [START_REF] Alias | Multi-dimensional Rankings, Program Termination, and Complexity Bounds of Flowchart Programs[END_REF] and [START_REF] Feautrier | Accelerated Invariant Generation for C Programs with Aspic and C2fsm[END_REF] gives an intermediate automaton of 12 variables, 80 transitions and 9 locations. Aspic successfully computes the invariants but Rank fails with TWO_MUCH_VARIABLES error because the size of the underlying linear programming problem is too big. The whole computation time is 12 seconds on a Intel Core2 @ 1.60GHz. In the following, we will show how our method SToP , nally manages this code to prove its termination.

(q>0 && r >0) { i f ( a [ i ] < a [ j ] ) { a [ k ] = a [ i ] ; k = k+h ; i = i +1; q = q-1; } else { a [ k ] = a [ j ] ; k = k+h ; j = j -1; r = r
3 Preliminaries

Grammar for intraprocedural analysis

We dene in Figure 3 On this kind of programs, the problem of proving termination is undecidable, but like in [START_REF] Alias | Multi-dimensional Rankings, Program Termination, and Complexity Bounds of Flowchart Programs[END_REF] we will provide a sound (but not complete) analysis that is able to deal with a large bench of problems. Basically, the method described in this paper is based on the fact that while loops can be proven to be terminating quite independently of the other ones, provided that we keep an execution context precise enough, thanks to the notion of summary.

Summaries

Denition 3.1 Let C ∈ G be a statement and R C (x 0 , y 0 , x, y) the relationship between the values of the variables before and after C. C ′ ∈ G is said to be an abstraction of

C if R C ⊆ R C ′ .
It is obvious that if we manage to prove that C ′ terminates, then C terminates. Thanks to linear relation analysis, we are able to compute such abstractions of code behaviours. We again use Aspic ( [START_REF] Feautrier | Accelerated Invariant Generation for C Programs with Aspic and C2fsm[END_REF]) with the appropriate option to obtain a polyhedral over-approximation R C (x 0 , y 0 , x, y) of the behaviour of a given code C.

Denition 3.2 (Summary) Let C be a code and R C (x 0 , y 0 , . . . , x, y) be an over-approximation of the relation between initial variables x 0 , y 0 , . . . and nal variables x, y, . . .. Then the following code is called a summary of C :

.)) break

The summary of C is obviously an abstraction of C and thus can be used for proving the termination of C.

Inria 4 Intraprocedural modular termination

This section presents a termination procedure able to handle large programs. Our termination procedure is based on the following remarks.

The program terminates if and only if all the conditional loops (referred as while loops in the following) terminate. Hence, the termination analysis can focus on the smallest program slice which contains the while loops and the minimum subset of instructions to preserve their behaviour. This preprocessing is called termination-specic slicing and described in subsection 4.1. In general, the slicing reduces drastically the program size. However, this may not be sucient (see for example our motivating example) and the following analysis is generally required.

Basically, a while loop terminates if and only if (i) the body terminates for each iteration and (ii) the while loop itself terminates. This simple idea allows to decompose a big termination problem into several termination sub-problems, small enough to be processed successfully by WTC. The detail of the decomposition is presented in subsection 4.2.

Step 1. Termination-Specic Slicing

The termination-specic slicing step is a preprocessing which extracts the program subset relevant for proving termination. This subset must contains the while loops and the instructions required to preserve their behaviour. Consider the example given in gure 4. Clearly, for each variable used in the while condition, the denition sites must be kept (here, i=1 and i=i+1). But the variables read by these denitions sites must be kept as well, and so on. This process boils down to a backward traversal of the data dependencies, starting from the variables used in the while condition. Also, control dependencies must be taken into account. On the example, the while loop will be executed only when the then branch is chosen. In turn, this depends on the up variable, on which the same process must be applied. These are the basic ideas of program slicing, for which many variants has been developed [START_REF] Tip | A survey of program slicing techniques[END_REF] for various purposes as reverse engineering, program comprehension,. . . Our termination-specic slicing algorithm is depicted in Algorithm 1.

1 i f ( up == 0 ) { 2 i = 1 ; 3 while ( i <= n ) {
Do_Slice(function_body), the main function, selects the set of while statements of the current function_body and runs the main slicing function, Slice. Do_Slice is applied to each function body of the program.

Slice(stmt_set) proceeds the denition sites of stmt_set with a depth-rst traversal, ooding the connex part of stmt_set in the meaning of Definition_Sites. Notice that the function mark(stmt) adds the line of stmt to the slice. When stmt is a control structure (if, for,while,...), only the control lines are added (e.g. if, else, {, }), not the body.

RR n°8166

Definition_Sites(stmt) computes the set of statements whose stmt depends immediately, in the meaning of data-and control-dependence. More precisely, Definition_Sites computes the reaching denitions [START_REF] Aho | Compilers: Principles, Techniques, and Tools[END_REF] RD stmt (xi) of each variable xi read by stmt. The result is a set of assignments which can possibly dene the value of xi read by stmt. Also, the immediate compound control structure of stmt, Compound(stmt), is added if it exists.

This way control dependencies are taken into account.

Algorithm 1 Termination-Specic Slicing Algorithm 'd). Let us apply our slicing algorithm to the piece of code given gure 4. This code has been picked from the end of the motivating example. For the sake of the presentation, the while statement will be denoted by while and the if statement by if . Do_Slice call Slice with while_loop_set = { while }. Initially, the slice is empty, no lines are marked. Lines 3 and 6 are marked, then Definition_Sites(while) is run

and gives RD 3 (i) ∪ RD 3 (n) ∪ {if} = {i=1, i=i+1, if}.
Slice is then called on {i=1,i=i+1,if}. i=1 is processed rst. Line 2 is marked and Definition_Sites(i=1) is computed. i=1 does not read variables, thus we get RD 2 (∅) ∪ {if} = {if}. Slice is then called on { if }. Lines 1 and 7 are marked. Definition_Sites(if) gives RD 1 (up) ∪ ∅ = ∅. Indeed, up is not dened in this piece of code and if has no compound control structure. So, no additional recursive call is issued.

Going back on the recursive calls, Slice({ i=i+1 }) is called. Line 5 is marked, then Definition_Sites(i=i+1) is computed and gives RD 5 (i) ∪ {if} = {i=1, i=i+1, if}. Slice is then called on {i=1,i=i+1,if}. All these statements were already processed, so no Inria further recursive exploration is done. Do_Slice ends up with the slice dened by the marked lines {1,2,3,5,6,7}. This slice excludes a[i]=a[i+n] which, indeed, does not inuence the termination of the while loop.

4.2

Step 2. Scalable Termination This section presents a scalable termination procedure, which scatters the termination proof into several proof obligations, on small programs with a single conditional loop. These proof obligations will be small enough to be handled successfully with the termination method WTC, presented in [START_REF] Alias | Multi-dimensional Rankings, Program Termination, and Complexity Bounds of Flowchart Programs[END_REF].

Here, although WTC is a tool that can handle programs (with nested loops), we can see WTC as a termination engine that takes as input a single loop and an execution context, and tries to prove that this loops terminates under the given context. WTC can either answer OK, which means that the program actually terminates, or DK, which is inconclusive.

Basically, our algorithm is based on the structure of the program to analyse. As expressions have no side-eects, they terminate. Breaks and asserts terminate. A sequence of instructions terminates if each terminates (given their execution context). For tests, we have to prove the termination of each branch under their respective contexts. A for loop terminates if its body terminates (in our simple version of for loops). A (single) while loop terminates if we are able to exhibit a ranking function. This is done by a call to our tool WTC. To handle nested loops, we can replace the inner loop by a summary (3.2) of its behaviour, and call WTC on the outer loop with a lesser body.

Algorithm

The previous remarks are implemented in Algorithm 2. The main function Mod-ularTerm takes as input a statement and a structure that is able to give an over approximation of the context of each statement of the program under analysis. These invariants are computed using linear relation analysis [START_REF] Cousot | Automatic discovery of linear restraints among variables of a program[END_REF], implemented in the Aspic tool ([8]). ModularTerm proceeds the abstract syntactic tree with a depth-rst traversal, evaluating each sub-tree to the value OK if it terminates, or DK if the analysis does not succeed to answer. If a sub-tree is evaluated to DK, the analysis fails and stops with DK output.

Let us focus on while loops. To analyse loops, the algorithm : 1. First tries to prove the loop body to terminate. If the loop body cannot be proved to terminate, then the analysis fails without having to analyse the loop itself. 2. If the analysis succeeds, the termination of the loop is checked by applying the WTC method described in [START_REF] Alias | Multi-dimensional Rankings, Program Termination, and Complexity Bounds of Flowchart Programs[END_REF] on the loop itself, under its context. 3. Then the loop is replaced by its summary (section 3.2).The analysis continues with this last code. Notice that the subsequent calls to WTC will have to deal with simpler programs with a single loop, which reduces scalability issues. Notice that linear relation analysis expresses the invariants by a convex polyhedron, which over-approximate the actual possible values of the variables. As a consequence, the precondition of the loop needed to apply WTC (stored in pcinvs), as well as the summary of the loop are over-approximated. These approximations can make unprovable programs which could be proved successfully with WTC, on the cases where the ILP can be processed with a reasonable amount of resource (time and memory) (see Section 6).

Example (cont'd). Our termination procedure processes the loops 0, 1, 2, 5, 4 and 3 in this order (deep-rst search). For the sake of the presentation, we will restrict the discussion to the rst steps of the algorithm. contains while loops, which were all substituted by an abstraction of their eect on the variables.

The termination of loop5 comes from the two rst conditions at the beginning of loop5, whose global eect would be equivalent to the assignment m = m -2*min(p,m). As p ≥ 1 , m strictly decreases at each step of the loop and will end up with the value 0, ensuring the termination. A call to WTC ends with m as ranking function and gives a summary for loop5. In the same way, the process continues with the remaining loop4 and loop3, and ends up with the nal answer OK, meaning that the program always terminates.

Extensions

Interprocedural analysis

We add the notion of function calls and a main. Dealing with local variables, parameters, and function calls is rather classical. The main restrictions are that we avoid recursive functions and modications of the parameters and global variables during the function calls. To handle functions, we rst compute invariants for all program functions f : pcinvs f . Then, if the current statement is a function call x:=f(x1,...xn):

1. Compute the projection of the call context at the current control point (pcinvs caller (progpoint)) on the union of variables x 1 . . . x n and the global variables. This invariant is called callcontext.

2.

Call recursively M odularT erm on the function body whose context is modied by the call context : pcinvs callee ∩ callcontext.

Extensions to full C

For the moment, our tool is restricted to a subset of C and could not handle a general C program. However, it is always possible to extend the set of analysable programs with an appropriate set of normalizations, implemented as a preprocessing. We give thereafter two important normalizations which should enlarge drastically the class of analyzable programs.

For loops We assume the for loops to be Pascal-compliant: for i from expr1 to expr2 do BLOCK(i);. In C, the for loops can be much more general and used, for example, to iterate on linear data structures: for(ptr=first; ptr != NULL; ptr=next(ptr)) BLOCK(ptr). It is easy to design a preprocessing which detects non-Pascal-compliant for loops, and turn them into a while loop: ptr=first; while(ptr!=NULL) { BLOCK(ptr); ptr=next(ptr); } Pointers We assume the program to be pointer-free. In C, the use of pointers complicate the data-dependence analysis, hence the slicing (because of reaching denitions) and the invariant analysis. On the following program: p = &a; q = &a; *p=1; *q=2; x = a;

The reaching denition of a in x=a should be RD x=a (a) = *q=2. Also, the invariant of x=a should implies x = 2. To make these analysis possible, the interferences must be analyzed carefully, thanks to pointer analysis [START_REF] Amme | Data dependence analysis in programs with pointers[END_REF]. However, the data-dependences are overapproximated, thus this will aect the size of slices. As for the invariant computation, the obtained data-dependencies should be analyzed carefully to produce an integer interpreted automaton leading to correct invariants. This analysis is not trivial and is left for future work. and for which SToP now answers DK). We also give average execution times of each method.

Experimental results

We implemented our method as a driver over WTC, written in C++. The total number of LOC is 3000 for the driver itself, and 150 additional lines for statistics. The code intensively uses the source-to-source compiler infrastructure Rose ( [START_REF] Quinlan | Rose: Compiler support for object-oriented frameworks[END_REF]) : functions for parsing, constructing ow graphs, searching in the code structure, pretty printing and instrumenting C codes.

Related works

As far as we know, this is the rst tentative of proving termination of C programs in a modular way. Previous papers prove terminations of kernels, and other papers such as [START_REF] Bradley | Linear ranking with reachability[END_REF] use what they call C abstractors to make benchmarks but no general algorithm is given. Our algorithm is a driver built upon WTC, a previous toolchain to prove the termination of (little) C codes. To improve its precision, we could use (instead of WTC) other methods and tools to prove fragments of code :

The methods described in previous papers on termination such as [START_REF] Colón | Practical methods for proving program termination[END_REF], [START_REF] Podelski | A complete method for the synthesis of linear ranking functions[END_REF] or more specialised methods for other kinds of programs such as lists or trees ( [START_REF] Habermehl | Proving termination of tree manipulating programs[END_REF]). In these last cases, the numerical invariants that we precompute may not be precise enough, though.

Instead of testing if a given code fragment terminates under a given over-approximation of its context, we can use the method described in [START_REF] Cook | Proving conditional termination[END_REF] to compute (sucient) termination preconditions, and check if these preconditions are satised under the given context.

As proving termination is undecidable, all these methods may positive results when other ones timeout or give DK answer, thus they could be used as on-demand back-ends.

Inria 8 Conclusion and Future works

In this paper we have presented a framework for proving the termination of C programs in a modular way. We used WTC as a tool for proving termination of simple loops, but the method is fully adaptable to other termination engines. The preliminary experimental results show that they are relatively few regression cases comparing to the previous work published in [START_REF] Alias | Multi-dimensional Rankings, Program Termination, and Complexity Bounds of Flowchart Programs[END_REF], and we are able to deal with larger programs such as sipmamergesort in reasonable time.

Future work include more tests to improve the reliability of our tool suite with respect to C variants, and on much bigger programs such as linux drivers and image processing codes.

In particular, the size of the sub-programs generated by SToP expresses a trade-o execution time/precision which needs to be investiguated.

  int up ; /* r e a l l y b o o l e a n */ int a [2*n+1] up = 1; p = 1; loop4 : do{ // s o r t i n g a h

Figure 1 :Figure 2 :

 12 Figure 1: Our motivating example: iterative merge sort

  x0 = x; y0 = y; ... x = random(); y = random();... if(!R_C(x0,y0, ..., x, y, ..

Figure 4 :

 4 Figure 4: Example for slicing taken from our motivating example

  function ::= declaration list statement list funcall ::= var := functionname ( var1 ... varn )

  stmt = x := expr[x1,...,xn] then return RD stmt (x1) ∪ . . . ∪ RD stmt (xn) ∪ Compound(stmt) stmt = if cond[x1,...,xn] then S1 else S2 then return RD stmt (x1) ∪ . . . ∪ RD stmt (xn) ∪ Compound(stmt)

	1: function Do_Slice(function_body)
	2: 3:	while_loop_set = set of while statements Slice(while_loop_set)
		//At this point, all the statements belonging to the slice are marked
	4: end function 5: function Slice(stmt_set) 6: for all stmt ∈ stmt_set do 7: if is_not_marked(stmt) then
	8: 9: 10: 11: 12: end function mark(stmt) //Add stmt to the slice Slice(Definition_Sites(stmt)) end if end for 13: function Definition_Sites(stmt)
	14:
	15:	end if
	16:
	17: 18: 19: 20: 21: 22: end function end if if stmt = for i from expr1[x1,...,xn] to expr2[y1,...,yp] do S then return RD stmt (x1) ∪ . . . ∪ RD stmt (yp) ∪ Compound(stmt) end if if stmt = while cond[x1,...xn] do S then return RD stmt (x1) ∪ . . . ∪ RD stmt (xn) ∪ Compound(stmt) end if
		Example (cont

if if

Table 1 :

 1 WTC1 and WTC2 denote middle-sized benchmarks taken from our previous experiments in[START_REF] Alias | Multi-dimensional Rankings, Program Termination, and Complexity Bounds of Flowchart Programs[END_REF]. WTC2 consist on a bench of 8 sorting functions from the literature, and WTC1 are also classical codes including sipmamergesort,our motivating example. For each benchmark, we give its number of programs, the average number of lines of codes before(LoCA) and after (LoCB) slicing, Regressions denote the number of regressions (examples that were proven OK with WTC

	Benchmark #Progs LoCB LoCA Regressions WTC(s) SToP(s)
	WTC1	50	28	21	3	0.92	2.45
	WTC2*	6	55	39	0	3.66	3.52

Table 1

 1 experiments were done on a Dual Core 2Ghz. As we expected, the execution times on middle-sized examples are much larger with SToP than WTC, mainly because of the cost of slicing and producing intermediate les. In future versions of SToP , we will investigate the opportunity of calling WTC on larger subprograms (some nested loops can be handled with a unique call to WTC). For sorting functions, the method shows it pertinence in terms of timing results and in terms of precision (it handles sipmamergesort in 22 seconds, but we had to manage the very last step by hand by providing a coarser abstraction for loop5).

give an extract of the experiments we made with our tool SToP (the entire benchmark is available on http://compsys-tools.ens-lyon.fr/stop.

The Thus SToP is able to deal with larger programs than WTC and seems to scale well. More experiments remains nevertheless to be done.

RR n°8166

end if 39: end function After instrumenting the code with loops labels and calling C2fsm and Aspic on the resulting code, we nd the following execution contexts for the loops : loop0 :

WTC succeeds to prove the termination of loop0 under its context (its nds the ranking function q + r), and also gives the following approximation for its behaviour : {j ≤ j 0 ∧ i ≥ i 0 ∧ h = h 0 ∧ n = n 0 ∧ i + q = i 0 + q 0 ∧ j + r 0 = r + j 0 } whose summary is : i0 = i; j0 = j; q0 = q; r0 = r; h0 = h; n0 = n; if(!(<=j0 && i>= i0 && h=h0 && n=n0 && i+q=i0+q0 && j+r0=r+j0)) break;

Loop1 and loop2 are processed in the same way. At the end, the outer loop loop5 no longer


Aurélien Monot

Nicolas Navet

Bernard Bavoux

Cristian Maxim

Fine-grained Simulation in the Design of Automotive Communication Systems

Early in the design cycle, the two main approaches for verifying timing constraints and dimensioning automotive embedded networks are worst-case schedulability analysis and simulation. The first aim of the paper is to demonstrate that both provide complementary results and that, most often, none of them alone is sufficient. In this paper, we present a simulation approach accounting for the clock drifts that occur on the network nodes at runtime and evaluate the extent to which the results obtained with this approach are relevant for the designers in order to validate the performances of a CANbased communication system. One of the practical outcome of this study is to show that initial phasings between nodes, as well as the values of the clock drifts, do not significantly impact the frame response time distributions that can be observed on the long run.

Introduction

Context. Electronic architectures in the automotive domain are defined years in advance and their realtime properties need to be evaluated and validated early in the design process. Amongst numerous design concerns, engineers are driven by two primary objectives which are to ensure safety and to dimension the electronic architecture so as to obtain a good trade-off between cost and performances. In order to achieve these goals two approaches are mainly used when designing an automotive communication system: worst-case analysis and simulation. As CAN is the main communication protocol between ECUs (Electronic Control Units) in cars, CAN frame response times have been studied extensively by the means of analysis and simulation. If Worst Case Response Times (WCRT) were initially the main focus of these research studies in order to ensure the safety requirements, response time distributions were eventually looked into as they provide richer information on the real-time behavior of the system. Experiments suggest that WCRT are rare events and quantiles such as 99,99% might be more relevant when dimensioning the electronic architecture as they can be significantly lower and might still provide enough guarantees to meet the safety requirements for many of the frames. Thus, being able to obtain accurate response time distributions might lead to significant cost optimization on hardware.

Previous work.

The timing analysis of CAN has been rather extensively investigated in the past. Bounds on the worst case response times have first been provided in [START_REF] Tindell | Calculating Controller Area Network (CAN) message response times[END_REF] and [START_REF] Davis | Controller Area Network (CAN) schedulability analysis: Refuted, revisited and revised[END_REF]. Studies have integrated the limitations of hardware [START_REF] Khan | Schedulability Analysis of CAN with Non-abortable Transmission Requests[END_REF][START_REF] Davis | Controller area network (can) schedulability analysis with fifo queues[END_REF][START_REF] Khan | Integrating hardware limitations in CAN schedulability anal-ysis[END_REF], and the communication stack, as the first analyses usually overlooked them, and also considered the effect of aperiodic traffic on CAN frame response times [START_REF] Khan | Aperiodic traffic in response time analyses with adjustable safety level[END_REF] and the consequences of transient perturbations [START_REF] Navet | Worstcase deadline failure probability in real-time applications distributed over CAN (Controller Area Network)[END_REF]. Of course, methods to minimize the response times, or make them more predictable, were investigated too, e.g. in [START_REF] Grenier | Pushing the limits of CAN -scheduling frames with offsets provides a major performance boost[END_REF] and [START_REF] Nahas | Reducing message-length variations in resourceconstrained embedded systems implemented using the Controller Area Network (CAN) protocol[END_REF]. In the previous works, some studies provide probabilistic analysis, and specifically response time distributions as a consequence of 1) all the possible phasing configurations between ECUs [START_REF] Zeng | Stochastic Analysis of CAN-Based Real-time Automotive Systems[END_REF][START_REF]Using statistical methods to compute the probability distribution of message response time in Controller Area Network[END_REF] and 2) CAN bit-stuffing mechanism [START_REF] Nolte | Probabilistic worst-case response-time analysis for the Controller Area Network[END_REF]. However, though they are a subject of interest in the field of wireless sensor networks [START_REF] Biasi | Simulation of process control with WirelessHART networks subject to clock drift[END_REF], none of these works assumed the possibility of ECU clock drifts until [START_REF] Monot | Impact of clock drifts on CAN frame response time distributions[END_REF] which this paper builds upon. Clock drifts result from the fact that the clocks of the ECUs which drive the instantiations of the CAN frames do not exactly operate at the same frequency. Due to production tolerances, the oscillators are not exactly identical and their frequency may also change over time because of environmental factors such as the temperature. In practice, the clock drifts are unavoidable and cause the phasings between the communicating ECUs to vary continuously over time. As a direct consequence, the response times of frames vary continuously too.

Paper organization. The paper is organized as follows. Section 2 describes our approach to simulate CAN networks. Then, section 3 describes how simulation can be used along analysis to study the worst case scenario. In the ensuing section 4, we present how to use our simulation approach to derive and interpret frame response time distributions. Finally, section 5 concludes this paper.

Fine-grained simulation of a CAN Network

In this section, we describe our simulation model and present the simulation tool. Specifically, we detail the clock model and assumptions made on the CAN nodes as well the kind of results obtained from the simulation. Clock drifts are measured in "parts per million" or ppm, which express how slower or faster a clock is compared to a "perfect" clock. For instance, 1 ppm corresponds to a deviation of 1µs every second. Typical values provided by quartz manufacturers for an automotive usage are tolerances of +/-50 ppm for the fabrication, +/-5 ppm per year for aging and +/-150 ppm for temperature in their operating temperature range (i.e., -40°C/125°C). In practice, the observed values tend to be smaller, e.g. +/-20 ppm over the whole temperature range.

Modeling the CAN network

In this study, we choose a rather simple and widely applicable drift model, based on fixed deviations of clock speeds (positive or negative) with respect to the nominal bus speed. Here, the variation due to aging is assumed to be negligible for the simulation lengths that we are considering. Furthermore, we assume that the CAN nodes are operating at a constant temperature. For a given clock c driving a CAN node, its local time t c with respect to a global time t is simply as follows :

t c (t) = φ c + δ c • t
where φ c is the initial start time (the phasing) of the CAN node with regard to the bus time referential, and δ c is the constant drift value. For instance, a drift rate of +50ppm means that δ c = 1.000050. For this work, every CAN node j is assigned a clock defined by the tuple (φ j , δ j ).

Hypotheses on the CAN communication system

For this work, no specific hardware characterics, other that the clock drifts, are considered. At each communicating node, the messages are queued in priority order and there is a sufficient number of communication buffers such that the highest priority message is always considered for arbitration. Additionally, the copy time from the message queue to the communication buffers is assumed to be negligible.

The CAN frames are considered to be strictly periodic. A CAN message m i is described by (T i , O i , C i , P i ) where T i is its period, O i its initial offset, C i its transmission time and P i its priority. In the context of this study, the worst case of bit-stuffing is always considered for C i . Offsets serve to desynchronize the streams of frames sent by a node so as to reduce peak loads and improve the response times as typically as done in automotive networks today (see for instance [START_REF] Grenier | Pushing the limits of CAN -scheduling frames with offsets provides a major performance boost[END_REF]).

Since, the instantiations of the CAN messages are driven by their sending node clock, quantities T i and O i are measured by the local clock of the sending node. As a result, messages with equal periods but sent by different CAN nodes may have slightly different transmission periods with respect to the bus time reference. On the other hand, C i , the transmission time on the CAN bus, depends only of the frame payload and the bit stuffing, and is expressed in the nominal time referential. Simulating this system consists in computing the frame release dates and reproducing the CAN arbitration and transmission process.

Assessing frame response time distribution

Frame response times essentially depend on the frame characteristics and of the phasings between the communicating nodes. Upper-bounds on the response times, valid whatever the phasings between ECUs, can be computed with analysis tools such as NETCAR-Analyzer 1 and they correspond to the worst-case behavior of the network. Response times in the typical case can be better evaluated by simulation. The novelty of the simulation study performed in this paper is that clock drifts are explicitly modeled and their impact on the network performances are evaluated. The classical simulation approach is to derive statistics on a number of configurations corresponding to distinct initial configurations (i.e., different phasing between ECUs at the system startup). This is needed because a single simulation run only captures a limited set of response times for each frame, as the period of the system is small (typically a few seconds). Here we explore another approach which is to gather statistics on a single trajectory of the system simulated during a long time (i.e., one simulation run). As the phasings between the ECUs vary continuously due to the clock drifts, simulating a single trajectory is sufficient to observe a wide range of response times for each frame. Because there is a single simulation run to perform, this approach is more convenient from an experimental point of view. Besides, it is more realistic in the sense that clock drifts, which unavoidably occur at run-time, are accounted for. The aim of the paper is to evaluate the extent to which the results obtained with this approach is relevant for the designers in order to validate the performances of a CAN-based system.

Software toolset

The main software used in this study, RTaW-Sim [START_REF]RTAW-Sim: Fine-grained simulation of Controller Area Network with faultinjection capabilities[END_REF] is a fine-grained discrete-event CAN bus simulator providing the frame response time distributions. It is free for all uses and available for download. The granularity of the simulation is 1µs allowing thus to simulate accurately a CAN bus at the bit level even at the highest possible speed of 1 Mbit/s. RTaW-Sim has various other features such as a precise modeling of the communication buffers, and fault injection capabilities that will not be used here. Also, we are going to use its import feature from NETCAR-Analyzer that allows to replay the phasing scenarios leading to the worst case response times 2 . In this paper, we present the experimental results obtained by simulation on realistic sets of message.

The simulator is able to compute about 700 kiloevents per seconds on a 2 GhZ laptop computer, which means for instance that it can simulate 24 hours of communication on a 60% loaded CAN bus in about 20 minutes. When simulating a CAN bus, the usual approach (see §2.2) is to compute a few hyperperiods for as much phasing configurations possible, resulting often in tens of hours of simulation even using a rather coarse-grained time granularity.

Exploring worst case scenario

In order to ensure safety, design engineers have to look at the worst case response time of frames which is typically done by analysis. However, analysis techniques do not tell us for how long the worst cases last, nor how often these situations happen.

Finding the worst case scenario

Finding the worst case response times using analysis is a problem that can be solved with some COTS software tools such as Netcar-Analyzer [START_REF]NETCAR-Analyzer: Timing analysis and resource usage optimization for Controller Area Network[END_REF]. However, being able to find the worst case scenario through simulation alone proved in our experiments to be almost impossible on realistic-size network. Determining a worst case response scenario for a frame consists in identifying the specific phasings between sending ECUs leading to the longest possible response time.

The number of phasing configurations for a CAN network is ( T τ ) N with T being the hyperperiod of the set of periodic frames, τ the bit-time of the bus and N the number of ECUs. For a typical high-speed CAN network at 500 kbit/s made up of 10 ECUs sending a set of periodic frame having a hyperperiod of 2s, this amounts to 10 60 different phasing configurations. Amongst them, several phasing configurations could possibly lead to the same response time but the number of phasing configurations giving different response times is still very large. As a consequence of the clock drifts, numerous phasing configurations will be explored in a single simulation run but this still does not guarantee to find a phasing configuration triggering a worst case response time. Fortunately, we can use the results of an analysis tool to reproduce the phasing parameters leading to a worst case scenario for a specific frame. This will be done in the next subsection.

Reproducing a worst case scenario

In these first experiments, we will replay a configuration leading to the wort-case response time for a specific frame and examine what happens shortly after a worst case response time. In order to do that, we start a simulation with the phasing parameters between ECUs leading to the a frame WCRT that we obtain through analysis. Here, the message set used is the SAE benchmark at a speed of 250 kbit/s corresponding to an average bus load of 62,3%. The frame under study has identifier 49. The offsets for the frames were computed using the SOA algorithm from NETCAR-Analyzer. Simulation statistics are shown in figure 2. We simulated four configurations corresponding to different drift parameters. The reference is represented by the black and green curves giving the average and maximum response time of the frames with no drift. Since there is no drift, and simulation is larger than the hyperperiod of the system, these two curves are identical for the three simulation length. The other curves correspond to statistics obtained with random drift values in +/-50 ppm, +/-150 ppm (realistic drift values) and +/-1000 ppm (very pessimistic with respect to suppliers data sheet). The top graph shows that after 10s of simulation (5 hyperperiods), the average response times for all the configurations (except +/-1000 ppm) are very close to the maximum response times of the case without drift which is the sign that most response times are large. The middle graph shows that after 1mn, the average response time have significantly decreased. It is worth pointing out that configurations with higher drifts lead to lower average values because the system get away from the worst case phasing configuration quicker. Finally, the bottom graph shows that after 10mn the average response time curves stabilize around values that are significantly lower than with smaller simulation lengths, whatever the clock drifts. Considering that 150 ppm is a realistic drift bound value, these experiments suggest that the frames stop experiencing large average response times above 1mn. Interestingly, this means that clock drifts, which are usually seen as being detrimental because they reduce the predictability of the system, might also be beneficial in some cases: they help getting away from phasing configurations leading to large response times.

Maximum observed response time versus worst-case scenario

In this second set of experiments, we simulated during 4 hours (i.e., a long car trip) different CAN networks using 20 random different initial phasings. For a low priority frame, we recorded the worst observed values among the simulation runs for the response times, average response times and two high quantiles. For this set of experiments, we used the SAE benchmark, the message set introduced in [10] and a message set which was derived from an existing CAN set from PSA Peugeot-Citroën. Results are given in Figure 3. Interestingly, we were able to find the worst case scenario for the first network. Though, it should be noted that the number of ECUs is small, thus considerably lowering the number of possible phasing configurations. However, it should be noted that the average and high quantile values of the response time distributions are always noticeably lower than the worst case, especially for the third network where the num-ber of nodes is realistic in view of today's automotive networks.

CAN networks

Deriving frame response time distributions

In this section we study the extent to which the results of simulations with drifts relate to the expected behavior of the network.

In other words, we explore the impact of the simulation parameters on the resulting frame response time distributions to determine the relevance of simulating a single trajectory with clock drifts for a long time.

The purpose of having a long observation window is to be able to observe the overall response time statistics over long scenarios (a single trip or the car estimated total functioning time, e.g. 7000 hours).

The statistics obtained by simulation should indicate what kind of performance can be expected from the studied configuration. This should allow engineers to draw conclusions about the dimensioning of the CAN bus: its speed and set some characteristics of the frames (e.g., priority, offsets, periods).

Impact of the simulation length

Here we study how the statistics stabilize over time.

While doing preliminary simulations with large simulation length, we observed that frame response time statistics vary at the begining but eventually stabilize.

In fact, with drifts the system is still periodic but the period is however considerably larger. Supposing that the hyperperiod of the frames sent by a node i is T i , the hyperperiod of the arrival of the frames from this node becomes T i = T i • δ i , with δ i being the drift value for the sending node. If we do that for every node, we are still simulating a periodic system though the hyperperiod of the system became much higher (ppcm of the T i • δ i ). For instance, assuming hyperperiods equal to T for N sending nodes and drift value close to δ leads to a period that can be up to (T • δ) N which leads to approximately (2s • 1.00015) 10 ≈ 17mn for T = 2s, 10 nodes and drift values in the range of +/-150ppm. The scheduling of the CAN frames become thus periodic after an initial phase (assuming the load of the system is less than 100%). As a consequence, the distribution of response times corresponding to the periodic phase of the scheduling will be repeated as the simulation lasts, making the response times recorded during the initial phase negligible. This means that it is not needed to simulate the CAN networks more than for a certain time as the response time distributions will no longer change after some point. It should be noted that this only explains the fact that the response time distributions stabilize as the simulation goes on but this does not mean that the resulting distributions should be identical. In the next experiments we compare the distributions resulting from different initial offsets and initial drift values.

Impact of initial offsets between sending nodes

In these experiments, we study the overall behavior of the CAN frames response time statistics over long simulations. We provide here results for a simulated time of 24h. The message set is the one introduced in [START_REF] Zeng | Stochastic Analysis of CAN-Based Real-time Automotive Systems[END_REF] which is made up of 69 CAN frames sent by 6 different nodes over a 250Kbit/s CAN bus. The total load of this message set is 60.25%. The offsets for the frame were computed using the SOA algorithm from NETCAR-Analyzer [START_REF]NETCAR-Analyzer: Timing analysis and resource usage optimization for Controller Area Network[END_REF]. Figure 4 shows response times statistics for three configurations corresponding to the same clock drift parameters but different initial phasing configurations. The initial phasing and clock drift values were chosen randomly (clock drifts are in +/-150 ppm range). The first thing to point out is that the maximum observed response times (the green curves) are slightly different depending on the initial configuration. This is especially true for the lowest priority frames (right-hand side of the graphic). As expected, the minimum response time curves are all equal because the phasings between sending nodes changes in such a way that every frame will be sent without delay at least once during the 24h. What is more surprising is that the average response times are equal and that the 99% and the 99.9% quantiles are almost identical. This is a meaningful result because it suggests that, for a long simulation time, the initial phasings are not significant for the final response time distributions statistics (except of course for the maximum observed response times).

Impact of the drift values

The aim here is to evaluate the impact of the clock drift values on the resulting frame response time statistics. Our first experiments showed similar results as seen with the initial phasings and differences were hardly noticeable. However, since using different drift values effectively affects the bus load, it should affect the response time distributions too. As a consequence, we tried to use a larger value range for the drifts : +/-1000 ppm. We also picked values so as to impact the average load while staying within the tolerance bounds. The figure 5 shows the resulting cumulated frame response time distributions of the lowest priority frame from the SAE benchmark for 24h simulation runs corresponding to different drift configurations. The simulation length is chosen large enough so that response time statistics have stabilized. The purple and green curves correspond to configurations with random drift value respectively in the range +/-1ppm (almost no drift) and +/-1000 ppm (high possible drift range). The blue and red curve correspond to configurations with close but different drift values (so that stations do not keep the same synchronization) that are chosen near the limit of the +/-1000 ppm range (e.g., -1000 ppm, -999 ppm, -998 ppm, etc, for the red curve and +1000 ppm, +999 ppm, +998 ppm, etc, for the blue curve).

Case study

This last experiment aims to evaluate the impact of clock drifts with sets of parameters corresponding to realistic cases. Since initial phasings depend of the relative start-up times of the nodes which cannot be predicted and since clock drifts are small, we study the relevance of a single simulation run to assess the expected performances of a CAN network.

In order to do that, we compare the distances between the resulting distributions for different simulations with random initial phasings and random drift values within pre-established drift ranges (respectively +/-100 ppm and +/-1000 ppm for the two sets of experiments). Each set of experiments correspond to 10 simulations of 48h of the SAE benchmark. We study and compare the response time distributions of the lowest priority frame at different times in the simulation. It is shown in Figure 6 how the distances between these simulations evolve by computing the standard deviations of the different response time intervals and compare their average.

Here the response time range is subdivided into 100 intervals of equal size. The graphic clearly shows that the standard deviations decrease significantly over time which indicates that the response time distributions are getting closer to each other. Another observation is that this happens faster with larger clock drift values as the set of experiments with +/-1000 ppm give lower standard deviation values. In combination with the fact that each distribution stabilizes over time, these curves indicate that all distributions are getting close to some limit distribution as they stabilize whatever their starting point and clock drift values. It should be noted that the results are more clear-cut with higher priority frames for which standard deviations are even closer to zero. This suggests that, assuming realistic clock drifts values, we can capture the behavior of a CAN network over a long time with a single limit distribution. This also suggests that a single long simulation run is sufficient to assess the expected behavior of a CAN network over time. Additionally, other experiments have shown that higher clock drifts help to stabilize the statistics faster which means that it can reduce the simulation times to obtain long-term CAN frame response time distributions. However, we have seen in section 3 that this simulation approach is not appropriate to capture corner cases: i.e, worst case scenario and other rare and large response times.

Conclusion

In this paper we present a simulation approach for CAN networks taking into account the clock drifts affecting the sending nodes which causes response times to vary over time. In section 3 of this paper we have seen experimentally that simulation and analysis alone are not enough. On the one hand, these experiments suggest that simulation alone is not appropriate to find the worst case scenarios because they are too rare. On the other hand, worstcase analysis cannot help to quantify how rare these events are, nor how long they last, nor what the average (or any other relevant statistics) of the response times are. However, it is possible to find by analysis the phasing conditions, specific to each frame, that provoke its worst case response time. Then, using a simulation tool, it becomes possible to observe for how long this situation lasts and where the clock drifts lead from there. Such simulations also contribute to validate the results obtained from the analysis tool, which is needed because these tools are usually black boxes and, though progresses are steadily being made [START_REF] Khan | Schedulability Analysis of CAN with Non-abortable Transmission Requests[END_REF][START_REF] Davis | Controller area network (can) schedulability analysis with fifo queues[END_REF], have to make simplifications about the hardware and the communication stack [START_REF] Navet | Frame latency evaluation: when simulation and analysis alone are not enough[END_REF].

We have then established that the frame response time distributions stabilize over time and that the initial phasings between the communicating nodes do not seem to have a significant impact on the resulting distributions. However, other experiments have shown that clock drifts may have an impact on the frame response time distributions when their values are far from realistic parameters. A noteworthy outcome of the last study is that, for typical ECU clock drift values, the response time distributions converge whatever the initial phasings between the sending nodes. This suggests that a single simulation is enough to capture the relevant statistics.

As future work, we plan to experiment with variable clock drifts (due for instance to temperature variation) which are already implemented in the simulation tool. As a result of using variable clock drift, the system will no longer be periodic meaning that the stabilization of the response time distribution is no longer a given.

Figure 1 :

 1 Figure 1: Screenshot of RTAW-Sim. The largest window in the underground contains the frames properties. The two smallest windows show the simulation options and its progress. The last bottom right window displays simulation statistics after a run.

Figure 2 :

 2 Figure 2: Simulation statistics with initial phasings leading to worst case response time for the frame 49. From top to bottom, the graphics show the frame response time for the first 10s, 1mn and 10mn of simulation. The X-axis are the frames sorted by decreasing order of priority. The black curve is the maximum response times recorded without drift. The other curves are the average response times recorded without drift (green) and drift bounds resp. equal to 50 ppm (red), 150ppm (blue), and 1000 ppm (purple).

Figure 3 :

 3 Figure 3: Statistics for a low priority frame of each benchmark obtained for random initial offsets and a drift range of +/-150 ppm. The values in the 4 last rows are the highest values observed in the set of experiments done for the corresponding CAN network.

Figure 4 :

 4 Figure 4: Response time statistics of a CAN message set for three different initial phasings. The X-axis are the frames sorted by decreasing priorities ( i.e., increasing IDs) while the Y-axis shows the response times. The different curves correspond to different response time statistics. From the lowest to the highest curve: minimum, average, 99% quantile, 99.9% quantile and maximum response times.

Figure 5 :

 5 Figure 5: Cumulative response time distributions of frame 49 of the SAE benchmark for different drift values recorded for a simulation length of two days. The X-axis is the response times in µs. The different curves correspond to different clock drift configurations.We clearly observe that the resulting frame response time distributions are different which indicates that the clock drift values do have an impact on the response times. Indeed, choosing clock values close to +1000 ppm and -1000 ppm makes a difference of 0,2% of the initial load which may have an impact on the response times of the lower priority frames. Because the load is lower with -1000 ppm, the red curve is above the others until 2,5 ms and below the others afterward because the response time distribution contains more lower response times. The blue curve (i.e., -1000 ppm) is the exact opposite. However, the green and purple curve are close to each other as their average load are most likely close to

Figure 6 :

 6 Figure 6: Evolution over time of the standard deviation of the response time distributions for sets of 10 simulations with random initial phasings and drift values. The frame under study is the lowest priority frame of the SAE benchmark. The Y-axis correspond to the standard deviation value and the X-axis to the simulation time.

NETCAR-Analyzer is a worst-case response time analysis tool[START_REF]NETCAR-Analyzer: Timing analysis and resource usage optimization for Controller Area Network[END_REF] ©INRIA/INPL/RTaW

For each frame, the initial offset configuration leading to the worst case is usually different.

Acknowledgement

We would like to thank Rob Davis (University of York) who helped us to find the experimental set-up showing the impact of clock drift in §4.3.


Matthieu Roustit
email: matthieuroustit<mroustit@chu-grenoble.fr

Marcin Hellmann

Claire Cracowski

Sophie Blaise

Jean-Luc Cracowski

Keywords: Adult, Aged, Female, Humans, Male, Middle Aged, Phosphodiesterase 5 Inhibitors, adverse effects, pharmacology, Piperazines, adverse effects, pharmacology, Purines, adverse effects, pharmacology, Raynaud Disease, drug therapy, physiopathology, Regional Blood Flow, drug effects, Skin, blood supply, Skin Temperature, drug effects, Sulfones, adverse effects, pharmacology Raynaud s ', sildenafil, microcirculation, local cooling

est destinée au dépôt

Sildenafil increases digital skin blood flow during all

phases of local cooling in primary Raynaud's phenomenon.

Sildenafil increases digital skin blood flow during all phases of local cooling in primary Raynaud s phenomenon

Introduction

Raynaud s phenomenon (RP) is defined as episodic vasospasm of the extremities in response to cold or emotional stimuli, often ' accompanied by pain. RP can be primary (i.e. idiopathic) or secondary to a connective tissue disease. The pathophysiology of primary RP is multifactorial and complex, with both vascular and neural abnormalities (1, 2). Although the morphology of skin capillaries is normal in primary RP, available evidence suggests that skin microvascular function is impaired (3 5). By using a local cooling test we recently showed increased digital skin vasoconstriction in participants with primary RP compared to healthy controls [START_REF] Roustit | Impaired transient vasodilation and increased vasoconstriction to digital local cooling in primary Raynaud s phenomenon[END_REF].

Local cooling of the skin induces an initial vasoconstriction followed by a transient vasodilation and finally, a prolonged vasoconstriction [START_REF] Johnson | Local thermal control of the human cutaneous circulation[END_REF]. Cooling exerts a significant portion of its vasoconstrictor effects through inhibition of the nitric oxide (NO) system [START_REF] Johnson | Local thermal control of the human cutaneous circulation[END_REF]. Indeed, NO is a potent vasodilator as it generates cyclic guanosine-5-monophosphate (cGMP), which induces vascular smooth muscle relaxation leading to vasodilation. cGMP is hydrolyzed by phosphodiesterases, and particularly by cGMP-specific phosphodiesterase 5 (PDE5).

PDE5 inhibitors have been mostly studied in RP secondary to systemic sclerosis as continuous therapy [START_REF] Levien | Advances in the treatment of Raynaud s phenomenon[END_REF][START_REF] Lavega | Phosphodiesterase-5 inhibitors for the treatment of Raynaud s: a novel indication[END_REF]. In less severe primary RP, single dose PDE5 inhibitors (e.g. before exposure to cold) could be an interesting therapy in the prevention of vasospasm, thus avoiding continuous exposure to vasodilators.

Although local cooling cannot be used as a surrogate outcome for RP, it could be a safe and reproducible experimental test to explore the effect of drugs on the different phases of microvascular reactivity to local cooling in primary RP [START_REF] Roustit | Reproducibility of a local cooling test to assess microvascular function in human skin[END_REF].

We hypothesized that sildenafil, a PDE5 inhibitor, may reverse the exaggerated microvascular response to local cooling in primary RP. The primary objective of this study was to assess the effect of a single oral dose of sildenafil (50 mg and 100 mg) on cutaneous vascular conductance (CVC) in the fingers of patients with primary RP, while cooling locally at 15 C or 24 C. We also assessed the effect °°of sildenafil on resting skin blood flux, CVC and skin temperature, as well as on the perfusion gradient between the different phalanges and the hand. The safety of sildenafil in primary RP patients was also assessed.

Results

Population characteristics

The demographic and clinical characteristics of the study population are summarized in Table 1. Antinuclear antibodies were negative and the nailfold videocapillaroscopy pattern was normal in all subjects. One postmenopausal woman was taking hormonal substitutive treatment, and three of the premenopausal women were taking oral contraceptives. Among the premenopausal women, 4 were enrolled during the follicular phase and 4 during the luteal phase of the cycle.

Effect of sildenafil on resting cutaneous vascular conductance and skin temperature

Oral 50 mg sildenafil significantly increased skin temperature from 28 3.3 C before administration to 29.3 3.1 C 50 min later (P ± °± °= 0.02, Wilcoxon rank test). In the same way, 100 mg sildenafil significantly increased skin temperature from 28.9 3.3 C to 30.7 2.7 C ± °± °(P 0.008, Wilcoxon rank test). There was a significant difference in skin temperature over time between 50 mg sildenafil and 100 mg = sildenafil (P 0.048, two-way ANOVA) (Figure 1A).

=

Resting skin blood flux was 32.1 19.8 PU before and 32.3 25.1 PU 50 min after administration of 50 mg sildenafil (P 0.92,

± ± =
Wilcoxon rank test), and 44.9 30 PU before and 58.7 37.9 PU after 100 mg sildenafil (P 0.007, Wilcoxon rank test) (Figure 1C).

± ± =

Resting CVC was 4.1 2.7 mV/mm Hg before and 4 2.8 mV/mm Hg 50 min after administration of 50 mg sildenafil (P 0.34, Wilcoxon

± ± =
rank test). Alternatively, CVC was 4.8 3.9 mV/mm Hg before and 7 5.3 mV/mm Hg 50 min after administration of 100 mg sildenafil ± ± (P 0.001, Wilcoxon rank test). There was a significant difference in skin blood flux and CVC over time between 50 mg sildenafil and 100 = mg sildenafil (P 0.012 and P 0.035, respectively; two-way ANOVA) (Figure 1D).

= =

Effect of sildenafil on cutaneous vascular conductance during local cooling

Sildenafil at 100 mg, but not 50 mg, significantly lessened the 15 C local cooling-induced decrease in skin blood flux compared to °control (Table 2, Figure 2A). This difference is not uniquely due to the change in baseline, as the AUC (i.e. excluding baseline CVC)

5-35

were also significantly different (6856.7 5977.4 and 2248.7 1778.6, respectively; P 0.04).

± ± =

This difference is observed at baseline and for all phases of local cooling, i.e. initial vasoconstriction, transient vasodilation and prolonged vasoconstriction (Table 2). When data were expressed as a percent change from baseline however, there was no significant difference between sildenafil at 100 mg, 50 mg and control (Table 2).

Similar results were observed when the cooling temperature was set at 24 C (Table 2, Figure 2B). Again, the difference in skin blood °flux between 100 mg sildenafil and control was not uniquely due to the change in baseline (AUC were 5592.1 3434.2 and 2518.5

5-35

± ± 920.9, respectively; P 0.008). The difference between 100 mg sildenafil and control was significant at baseline and for all phases of local = cooling when expressed as CVC (Table 2), but not as a function of baseline.

Effect of sildenafil on the skin perfusion gradient

At each visit we observed a perfusion gradient between the distal, middle and proximal phalanges and the dorsum of the hand: 15.9 ± 13.4, 11.5 8.3, 10.2 6.9 and 7.9 4.5 mV/mm Hg respectively for control (P 0.008); 13 9.9, 10.6 7.1, 9.4 6.1 and 8.1 4.7

± ± ± = ± ± ± ±
mV/mm Hg respectively for 50 mg sildenafil (P 0.004); and 22.3 13.5, 18 13.4, 16.2 13.8 and 10.3 5.7 respectively for 100 mg

= ± ± ± ±
sildenafil (P 0.001). The difference between 100 mg sildenafil, but not 50 mg, and control is significant for each phalange and for the = hand (Figure 3).

Safety

None of the subjects complained of pain associated with the local cooling and it did not induce any symptoms of RP. No adverse event was reported during visit 1 (no treatment). At visit 2 (50 mg sildenafil), four patients experienced headache: 3 were mild and regressed spontaneously, the fourth one was moderate and ceased after taking 1 g acetaminophen. Mild flushing was observed in 6 patients. At visit 3 (100 mg sildenafil), the same patient complained of moderate headache and was given 1 g acetaminophen. Mild flushing was observed in 4 patients, one of them also reporting nasal congestion.

Mean arterial blood pressure before 50 mg sildenafil (visit 2) and 100 mg (visit 3) were 90 10.9 and 93 10.8 mm Hg, respectively.

± ±

The nadir were reached 37 and 43 min after sildenafil administration, and were 85.9 9.6 (P 0.07 vs baseline) and 83.5 14.6 mm Hg (P

± = ± =
0.01 vs baseline), respectively (Figure 1B).

Discussion

Our data show that a single dose of 100 mg sildenafil significantly increases cutaneous vascular conductance during the three phases of localized digital cooling in patients with primary RP. It also increased resting CVC and skin temperature. The effect was observed within 1 hour after sildenafil intake, at 15 C and 24 C. Finally, we observed a non significant trend towards increased CVC while cooling °°with Microvascular reactivity impairment in primary RP is not yet fully understood. A loss of CGRP containing nerve fibers in the digits has been suggested [START_REF] Bunker | Deficiency of calcitonin gene-related peptide in Raynaud s phenomenon[END_REF], as well as increased endothelin-1 dependent vasoconstriction (1). Another mechanism could involve postjunctional -adrenoreceptors, which are clustered distally, whereby their translocation from cytosol to cell surface is enhanced by

α 2c
cooling, thus leading to vasoconstriction (2). In vitro, cooling-induced -adrenergic constriction of arterioles isolated from patients with α 2

primary RP was increased compared to controls, which was suppressed by protein tyrosine kinase inhibitors [START_REF] Furspan | Increased tyrosine phosphorylation mediates the cooling-induced contraction and increased vascular reactivity of Raynaud s ' disease[END_REF].

Local cooling is a reproducible test to assess skin microvascular reactivity [START_REF] Roustit | Reproducibility of a local cooling test to assess microvascular function in human skin[END_REF]. Most of its vasoconstrictor effect depends on inhibition of the NO system and on adrenergic function [START_REF] Hodges | Adrenergic control of the human cutaneous circulation[END_REF], particularly through the translocation of -adrenoreceptors mediated by

α 2c
RhoA-Rho kinase [START_REF] Thompson-Torgerson | Cold-induced cutaneous vasoconstriction is mediated by Rho kinase in vivo in human skin[END_REF]. We recently showed that local cooling-induced vasoconstriction on the fingers of patients with primary RP was increased compared to healthy controls [START_REF] Roustit | Impaired transient vasodilation and increased vasoconstriction to digital local cooling in primary Raynaud s phenomenon[END_REF].

This study further demonstrates that sildenafil, a selective PDE5 inhibitor, increases digital skin blood flow both at baseline and during local cooling in primary RP. However, sildenafil did not restore the typical pattern of skin blood flow during local cooling described in healthy subjects, i.e. initial vasoconstriction followed by a transient vasodilation and finally, a prolonged vasoconstriction [START_REF] Johnson | Local thermal control of the human cutaneous circulation[END_REF]. Indeed, in most cases the transient vasodilation is blunted in patients with RP compared to controls, and partially restored by local anesthesia, suggesting an abnormal neural response in primary RP depending on cold sensitive nerves [START_REF] Roustit | Impaired transient vasodilation and increased vasoconstriction to digital local cooling in primary Raynaud s phenomenon[END_REF]. Sildenafil induced a shift towards higher flux values at all phases, explaining why its effect was not significant when data were expressed as a percentage of baseline.

Moreover, assessment of resting perfusion on the contralateral arm shows an increase in CVC on all phalanges and on the dorsum of the hand with sildenafil at 100 mg, but not 50 mg. All together, our data suggest that 100 mg sildenafil increases skin blood flow at baseline and in response to local cooling, one hour after oral administration. This delay was chosen according to previous work by our group showing that the effect of 100 mg oral sildenafil on the microcirculation was maximal one hour after administration [START_REF] Blaise | Oral sildenafil increases skin hyperaemia induced by iontophoresis of sodium nitroprusside in healthy volunteers[END_REF]. In that study, we showed that sildenafil at 100 mg, but not 50 mg, significantly increases forearm CVC [START_REF] Blaise | Oral sildenafil increases skin hyperaemia induced by iontophoresis of sodium nitroprusside in healthy volunteers[END_REF], which is consistent with the present findings.

Sildenafil was well tolerated; none of the participants experienced a serious adverse event. The mild adverse drug reactions observed in this study are similar to those observed with PDE5 inhibitors in the treatment of erectile dysfunction. The frequencies of headache and flushing were comparable to those previously reported [START_REF] Goldstein | Oral sildenafil in the treatment of erectile dysfunction. Sildenafil Study Group[END_REF]. The decrease in mean arterial pressure was maximal at around 40 min, and was 10 and 16 from baseline after administration of 50 mg sildenafil and 100 mg sildenafil, respectively.

-% -%

Several controlled studies have shown improved outcomes (i.e. RP duration, frequency of attacks, RCS or a visual analogue scale score) in secondary RP when treated with PDE5 inhibitors [START_REF] Fries | Sildenafil in the treatment of Raynaud s phenomenon resistant to vasodilatory therapy[END_REF][START_REF] Herrick | Modified-release sildenafil reduces Raynaud s phenomenon attack frequency in limited cutaneous systemic sclerosis[END_REF]. While these studies enrolled patients with secondary RP (particularly patients with systemic sclerosis) and assessed the continuous use of PDE5 inhibitors, little data is available for primary RP.

Friedman have studied the effect of another PDE5 inhibitor, tadalafil, on cold-induced skin vasoconstriction in 18 patients with et al primary RP and 2 with secondary RP [START_REF] Friedman | The effects of tadalafil on cold-induced vasoconstriction in patients with Raynaud s phenomenon[END_REF]. However, they did not show any significant change in skin blood flow at baseline and during local cooling between a single 10mg dose of tadalafil and placebo. The discrepancy between their findings and ours might be explained by an insufficient skin concentration of tadalafil after a 10 mg dose. Indeed, sildenafil shows a higher volume of distribution than tadalafil, and the maximal concentration obtained with 100 mg sildenafil is approximately that obtained with 20 mg tadalafil [START_REF] Wright | Comparison of phosphodiesterase type 5 (PDE5) inhibitors[END_REF]. This is also consistent with the lack of response we observed with 50 mg sildenafil.

Another explanation could be the difference in methods between the studies. Indeed, in the study of Friedman cooling was et al regional (whole hand) and flux was recorded on the finger pad [START_REF] Friedman | The effects of tadalafil on cold-induced vasoconstriction in patients with Raynaud s phenomenon[END_REF], whereas in the present work cooling and measurements were restricted to the dorsum of the finger. Importantly, in the study by Friedman , some blood flow measurements in response to cooling were et al paradoxically higher in patients with RP than in controls enrolled in preliminary studies [START_REF] Friedman | The effects of tadalafil on cold-induced vasoconstriction in patients with Raynaud s phenomenon[END_REF]. Therefore, the methods might be inappropriate to detect impaired microvascular reactivity to cooling in primary RP, thus making it more difficult to show a difference after tadalafil administration. The test we used in the present study allowed us to demonstrate microvascular dysfunction to local cooling in primary RP compared to healthy controls [START_REF] Roustit | Impaired transient vasodilation and increased vasoconstriction to digital local cooling in primary Raynaud s phenomenon[END_REF].

An open-label study has assessed the effect of vardenafil (10 mg bid) on skin blood flow in basal conditions and during exposure to cold, in 7 patients with primary RP and 33 with secondary RP [START_REF] Caglayan | Phosphodiesterase type 5 inhibition is a novel therapeutic option in Raynaud disease[END_REF]. Caglayan showed a significant increase in skin blood flow after et al vardenafil administration during exposure to cold at 4 C, but not at room temperature. The effect was visible after a single dose and after 2 °weeks [START_REF] Caglayan | Phosphodiesterase type 5 inhibition is a novel therapeutic option in Raynaud disease[END_REF]. Another small open-label study has reported significant improvement in digital temperature in response to a mild cold challenge (immersion in water at 15 C) after 50 mg sildenafil in 3 out of 5 patients with secondary RP [START_REF] Kumar | Thermographic and symptomatic effect of a single dose of sildenafil citrate on Raynaud s phenomenon in patients with systemic sclerosis: ' a potential treatment[END_REF]. We did not assess the effect of °local cooling on skin temperature (as the temperature was set at 24 C or 15 C), but we observed a dose-dependent effect of sildenafil on °°resting skin temperature and cutaneous blood flow.

In order to allow for the variability of baseline skin blood flow (i.e. after sildenafil administration but before cooling) we homogenized baseline skin temperature to 33 C, as previously described [START_REF] Minson | Nitric oxide and neurally mediated regulation of skin blood flow during local heating[END_REF]. Homogenizing skin temperature considerably decreased inter-individual ° variability (data not shown for clarity) although baseline CVC remained significantly increased by sildenafil. This increase in baseline may explain the trend towards greater vasoconstriction when data are expressed as a percentage of baseline while cooling at 15 C after 50 mg °sildenafil.

This is probably the consequence of a greater elevation in baseline CVC when homogenizing skin temperature at 33 C at 50 mg °than 100 mg. This could be explained by the potentiating effect of sildenafil on the vasodilation induced by local heating (which is partly NO-dependent). Indeed, after 50 mg sildenafil skin temperature increased from 29.3 to 33 C, whereas the difference was smaller (30.7 to °33 C) after 100 mg sildenafil. °Gender is another concern when studying microvascular function. The effect of the phase of the menstrual cycle or of oral contraception on microvascular function has been explored, but with conflicting results. Resting cutaneous blood flux and conductance are affected by sex, females having lower values than males [START_REF] Hodges | Influence of age, sex, and aerobic capacity on forearm and skin blood flow and vascular conductance[END_REF]. Finally, healthy females show greater vasoconstriction on local cooling than males, the response being more pronounced during the luteal phase than the follicular phase of the menstrual cycle [START_REF] Cankar | Gender differences in cutaneous laser doppler flow response to local direct and contralateral cooling[END_REF]. As our study design (7-day intervals between visits) implied differing hormonal status, we took care to include as many premenopausal women in the follicular phase as those in the luteal phase. We deliberately chose to shorten the inclusion period to limit seasonal differences in microvascular reactivity, which have been previously described in patients with primary RP [START_REF] Gardner-Medwin | Seasonal differences in finger skin temperature and microvascular blood flow in healthy men and women are exaggerated in women with primary Raynaud s phenomenon[END_REF][START_REF] Leppert | Seasonal variations in cyclic GMP response on whole-body cooling in women with primary Raynaud s phenomenon[END_REF] The main limitation of this study is that it was not randomized and double-blind. Indeed, we were unable to obtain a placebo for sildenafil from the manufacturer. We therefore chose objective measurements as the primary outcome (i.e. cutaneous vascular conductance). Moreover, given the lack of safety data on single doses of 50 or 100 mg in young and low body weight women, as well as an adverse drug event in a previous study by our group ( 14), we sequentially tested 50 mg sildenafil followed by 100 mg rather than randomize the visits.

In conclusion, this study shows that a single dose of 100 mg sildenafil significantly increases skin blood flow within one hour in patients with primary RP, both at baseline and during all phases of digital local cooling. Sildenafil administration did not induce any serious adverse drug event. This pilot pharmacology study suggests that 100 mg sildenafil could be used as required before exposure to " " cold in primary RP. This should be confirmed in a randomized double-blind controlled trial. Another interesting approach would be to determine in a laboratory study whether the higher skin blood flow decreases the pain or the frequency of vasospasm to whole body cooling in a laboratory study.

Methods

Study population

All the participants enrolled in this study were recruited through local newspaper advertisements and included between February 2011 and June 2011. All subjects were 18 years of age or older. Primary RP was diagnosed according to the criteria of LeRoy and Medsger [START_REF] Leroy | Raynaud s phenomenon: a proposal for classification[END_REF], and had to affect more than 2 fingers on a hand. Subjects taking calcium-channel blockers were instructed to stop medication one week before enrolment in the study. The Raynaud s Condition Score (RCS) was assessed as previously described [START_REF] Merkel | Measuring disease activity and functional status in patients with scleroderma and Raynaud s phenomenon[END_REF].

'

Non-inclusion criteria included pregnancy (urine pregnancy tests were performed at the beginning of each visit), cigarette smoking, any associated chronic disease, and an abnormal capillaroscopic pattern [START_REF] Cutolo | Raynaud s phenomenon and the role of capillaroscopy[END_REF]. Antinuclear autoantibodies were determined for all participants. In cases with positive antinuclear autoantibodies (>80 UI/mL), specific autoantibodies against topoisomerase I (Scl-70) or centromere-associated proteins were sought. Positive autoantibodies against topoisomerase I (Scl-70) or centromere-associated proteins were exclusion criteria.

The investigation conforms to the principles outlined in the Declaration of Helsinki. Grenoble Institutional Review Board (IRB n°6 705) approval was obtained and each subject gave written informed consent before participation.

Study design

This was an open-label pharmacology study. Three consecutive visits were planned for each volunteer, separated by 7 days / 3 at the +same time of the day. Upon arrival at the center subjects were placed in a temperature-controlled room (23 / 1 C). They remained supine + -°for the whole experiment.

Visit 1 included a clinical examination, followed by a local cooling test during which digital skin blood flux was continuously recorded. Blood pressure was also recorded continuously by using digital photoplethysmography (Nexfin monitor, Bmeye B.V., Amsterdam, The Netherlands). Before recording started, the arm was immobilized with a vacuum cushion to ensure stable positioning, as previously described [START_REF] Roustit | Impaired transient vasodilation and increased vasoconstriction to digital local cooling in primary Raynaud s phenomenon[END_REF]. The procedure was repeated at visit 2 and visit 3, one hour after the oral administration of sildenafil citrate (purchased from Pfizer France, Paris, France). Patients were being given 50 mg sildenafil on visit 2 and 100 mg sildenafil on visit 3 (Figure 4).

Local cooling and skin blood flow measurements

Cutaneous blood flow was assessed by LDF (Periflux System 5000, Perimed, J rf lla, Sweden). Two skin sites were equipped with ä ä custom-designed LDF cooling probes (Probe 415-317, Perimed, J rf lla, Sweden) [START_REF] Roustit | Reproducibility of a local cooling test to assess microvascular function in human skin[END_REF]. Risk analysis was performed and was fully ä ä compatible with human use. The skin sites were chosen on the dorsum of two fingers affected by RP between the index, the middle and the ring finger (left hand). Maricq color charts [START_REF] Maricq | Diagnosis of Raynaud s phenomenon assisted by color charts[END_REF] were used to confirm the diagnosis of RP and to specify its topography. When the three fingers were equally affected by RP, two of them were randomly chosen. After a 30-min acclimation at room temperature (23 C / 1), resting skin blood flux was recorded for 10 minutes. Skin temperature °+was then homogenized at 33 C and skin blood flux was recorded during the following 10 minutes to obtain baseline flux (visit 1). In visits °2 and 3, the 10-min baseline (skin temperature homogenized at 33 C) was done 50 min after taking sildenafil. Indeed, we previously °showed that the effect of sildenafil on microvascular reactivity was maximal 60 min after oral administration [START_REF] Blaise | Oral sildenafil increases skin hyperaemia induced by iontophoresis of sodium nitroprusside in healthy volunteers[END_REF]. After baseline, local cooling was maintained over 30 min at 15 C on one skin site and at 24 C on the other, simultaneously. Blood flux measurements were 50 °°min (visit 1) to 100 min (visits 2 and 3).

Perfusion gradient measurements

Skin perfusion of the dorsal surface of the hand (right hand) was recorded before cooling with laser Doppler imaging (PeriScan PIM 3, Perimed, J rf lla, Sweden). Regions of interest (ROI) were defined on the distal, middle and proximal phalanges of three fingers (index, ä ä middle and ring fingers), as well as on the dorsal surface of the hand (without the fingers). Mean CVC was calculated for all ROIs. Data from the 3 fingers were averaged.

Data analysis

Data were digitized, stored on a computer, and analyzed off-line with signal processing software (PeriSoft 2.5.5, Perimed, J rf lla, ä ä Sweden). Skin blood flow was expressed as cutaneous vascular conductance (CVC) in mV/mmHg (i.e. flux in millivolts divided by mean arterial pressure). Expressing data as conductance is a more physiological approach, as it takes into account differences and variations in blood pressure [START_REF] Leary | Regional vascular resistance vs. conductance: which index for baroreflex responses?[END_REF]. Then, a minute by minute analysis of CVC was performed to assess the kinetics of the response and data were expressed as area under the curve over a 5-min baseline and the 30-min cooling period (AUC , in mV.s/mm Hg) as the primary 0-35 outcome. We subsequently analyzed the 3 phases of the response as previously described: 1. initial vasoconstriction: CVC was averaged over 1 min around the lowest flux value within the first five minutes; 2. transient vasodilation: CVC was averaged over 1 min around the highest flux value within the first fifteen minutes; 3. late prolonged vasoconstriction: CVC was averaged over the last 3 minutes of the measurement [START_REF] Roustit | Impaired transient vasodilation and increased vasoconstriction to digital local cooling in primary Raynaud s phenomenon[END_REF]. In order to take into account baseline (BL) flux variations, data were expressed as a percentage change from baseline CVC as a secondary outcome.

Statistical analysis

Categorical data were reported as frequency and percentage, and continuous data as mean and standard deviation. They were analyzed by repeated measures ANOVA, and paired t tests for 2 2 comparisons. Mauchly s test of sphericity was used to assess equality of variance × ' in the data. When significant (i.e. inequality of variance cannot be excluded) a Greenhouse-Geisser adjustment was used. When the conditions of application of parametric tests were not respected, nonparametric tests were used (Friedman test, and Wilcoxon test for paired comparisons). Two-sided significance tests were used throughout. We considered p values <0.05 as significant, corrected by Bonferroni s method for multiple comparisons. Statistical analysis was performed with SPSS 13.0 for Windows (SPSS Inc, Chicago IL, ' USA). 4.1 (1.8) Quantitative data are expressed as mean (standard deviation). Qualitative data (female, thumb involvement and feet involvement) are expressed as number (percentage). BMI: body mass index.

Figure 1 Figure 2

 12 Figure 1Effect of a single dose of 50 mg (red line) and 100 mg (green line) of sildenafil given orally on resting skin temperature in a temperature-controlled room (A); systolic (SAP), mean (MAP) and diastolic (DAP) arterial blood pressure (B); skin blood flux (C) and cutaneous vascular conductance (CVC) (D). APU: arbitrary perfusion units. P<0.05 *

Figure 4

 4 Figure 4 Study day schema. LDI: laser Doppler imaging.

'

Table 1

 1 Demographic and clinical characteristics of the study population (N 15).

	=

Raynaud s Condition Score

[START_REF] Bunker | Deficiency of calcitonin gene-related peptide in Raynaud s phenomenon[END_REF]

[START_REF] Bunker | Deficiency of calcitonin gene-related peptide in Raynaud s phenomenon[END_REF]

[START_REF] Bunker | Deficiency of calcitonin gene-related peptide in Raynaud s phenomenon[END_REF]

[START_REF] Bunker | Deficiency of calcitonin gene-related peptide in Raynaud s phenomenon[END_REF]

[START_REF] Bunker | Deficiency of calcitonin gene-related peptide in Raynaud s phenomenon[END_REF]

[START_REF] Bunker | Deficiency of calcitonin gene-related peptide in Raynaud s phenomenon[END_REF]

[START_REF] Bunker | Deficiency of calcitonin gene-related peptide in Raynaud s phenomenon[END_REF]

[START_REF] Bunker | Deficiency of calcitonin gene-related peptide in Raynaud s phenomenon[END_REF]

[START_REF] Bunker | Deficiency of calcitonin gene-related peptide in Raynaud s phenomenon[END_REF]

Acknowledgements:

We thank Claire Sors and Florence Gaillard-Bigot for their assistance in data collection, as well as Dr Alison Foote for correction of English language usage.

Financial support

Association des Scl rodermiques de France; Groupe Fran ais de Recherche sur la Scl rodermie; D l gation la Recherche Clinique et l

Footnotes:

Conflicts of interest JLC and MR have received research grants from Pfizer for other studies.

References:

1 .