T T Nga Nguyen

Olivier Brun

Balakrishna J Prabhu

Joint downlink power control and channel allocation based on a partial view of future channel conditions

Keywords: Scheduling, utility maximization, average power constraint

We propose two downlink scheduling algorithms that take advantage of partial information on future channel conditions for improving the sum utility. The scheduling model allows for both power control and channel allocation. The objective of the scheduler is the long-term utility under an average power constraint. The two algorithms incorporate the channel predictions in their decisions.The STO1 algorithm computes the decision in each slot based on the means of future channel gains. Depending on the horizon considered, this can require solving a large-dimensional problem in each slot. The STO2 algorithm reduces the dimensionality by operating on two time-scales. On the slower scale it computes an estimation over a larger horizon, and in the faster scale of a slot, it computes the decision based on a shorter horizon. Numerical experiments with both fixed number of users as well as a dynamic number of users show that the two algorithms provide gains in utility compared to agnostic ones.

Introduction

Downlink scheduling and power control has been widely investigated for wireless networks [START_REF] Lee | Opportunistic power scheduling for dynamic multi-server wireless systems[END_REF][START_REF] Chiang | Power control in wireless cellular networks[END_REF][START_REF] Huang | Downlink scheduling and resource allocation for ofdm systems[END_REF][START_REF] Li | Network utility maximization over partially observable markovian channels[END_REF]. In each time-slot, one or multiple base-stations have to decide the users to serve and the transmit powers with the objective of maximizing the sum of the utilities of the users. Two different types of utilities can be defined: (i) opportunistic [START_REF] Lee | Opportunistic power scheduling for dynamic multi-server wireless systems[END_REF][START_REF] Chiang | Power control in wireless cellular networks[END_REF]; and (ii) long-term [START_REF] Huang | Downlink scheduling and resource allocation for ofdm systems[END_REF][START_REF] Subramanian | Joint scheduling and resource allocation in CDMA systems[END_REF]. In the opportunistic model, the utility function operates on the rate obtained in a time-slot whereas in the long-term model, the utility operates on the average rate obtained over an horizon.

The focus of this paper is on the long-term utility model. When transmit powers cannot be varied, the celebrated Proportional-Fair algorithm (see [START_REF] Kushner | Convergence of proportional-fair sharing algorithms under general conditions[END_REF] and references therein) is known to work well when the utilities are logarithmic functions of the total rate. This algorithm belongs to the class of gradient-based algorithms that choose the user that maximizes the marginal utility or the gradient of the utility function. In [START_REF] Huang | Downlink scheduling and resource allocation for ofdm systems[END_REF] the gradient-based solution was then extended to the setting in which joint power control and user scheduling is possible.

We revisit this problem in the context of vehicles which share their itineraries with the decision maker. With the availability of SINR maps in urban zones, the decisions can now be based also upon the future channel conditions of the users (or vehicles) [START_REF] Margolies | Exploiting mobility in proportional fair cellular scheduling: Measurements and algorithms[END_REF]. The future channel conditions are, however, not known exactly as the they vary randomly in time. The SINR maps are assumed to give the expected values of the channel gains on the routes takes by the users. With this additional information on the future expected channel gains, performance improvements can be expected compared to the setting when this information is not available.

We consider a downlink scheduling and power control problem for one base station with an average power constraint. The objective is the sum utility of the users, and the dependence of the utility function on the channel allocation and transmit power is similar to that in [START_REF] Huang | Downlink scheduling and resource allocation for ofdm systems[END_REF] with only one sub-channel and we do not have queues. In that respect, our model is a special case of that in [START_REF] Huang | Downlink scheduling and resource allocation for ofdm systems[END_REF]. However, the scheduling algorithms that we propose can be extended to the multiple sub-channel case though for this first paper we restrict ourselves to the simpler case of one sub-channel. There are two difference with the models studied previously. The first, and the main difference is that the base station is also aware of the mean channel gains in the future slots, and the second one is that we include an average power constraint over time in the optimization problem.

Contributions

We propose two heuristics for that use information on future channel conditions in order to improve the total utility. Both the heuristics were first proposed in [START_REF] Nguyen | An algorithm for improved proportional-fair utility for vehicular users[END_REF] in a setting without power control. The first one optimizes over future time-slots in every current slot whereas the second one reduces the complexity by optimizing over the future time-slots only once in a certain number of slots. The performance of these heuristics will be evaluated on three types of stochastic models for the channel gain processes. In the first, the channel gains are a stationary process with fixed mean channel gains (which can be different for different users) whereas in the second model, the mean channel gains are themselves varying on a slower time-scale to that of the channel gains themselves. Finally in the third model, the mean channel gains vary in every slot. For three models, these heuristics are shown to perform better than the other setting in which future information is not available.

Related work

Scheduling on the downlink but with fixed power examined in numerous papers [START_REF] Borst | Mobility-driven scheduling in wireless networks[END_REF][START_REF] Tan | Utility maximization resource allocation in wireless networks: Methods and algorithms[END_REF][START_REF] Chandur | Performance evaluation of scheduling algorithms for mobile wimax networks[END_REF]. For a logarithmic utility function, the Proportional-Fair (PF) scheduler has been known to be optimal for stationary channel conditions [START_REF] Kushner | Convergence of proportional-fair sharing algorithms under general conditions[END_REF] and its performance has been analyzed in both static and dynamic user scenarios [START_REF] Borst | User-level performance of channel-aware scheduling algorithms in wireless data networks[END_REF]. A restless bandit framework for network utility maximization for channel states modeled as partially observable Markov chains is proposed in [START_REF] Li | Network utility maximization over partially observable markovian channels[END_REF].

In [START_REF] Abou-Zeid | Long-term fairness in multi-cell networks using rate predictions[END_REF], future rates are assumed to be known accurately over the optimization horizon and improvement in data rates and fairness is compared with algorithms such as greedy (or max-rate) and equal share. Proportional-Fair algorithms with partial information on future channel conditions have been proposed in [START_REF] Bang | Channel predictive proportional fair scheduling[END_REF] (on short-time scales) and [START_REF] Margolies | Exploiting mobility in proportional fair cellular scheduling: Measurements and algorithms[END_REF] (on longer-time scales). The algorithm in [START_REF] Margolies | Exploiting mobility in proportional fair cellular scheduling: Measurements and algorithms[END_REF] is based on SINR maps that can be obtained from vehicles or users that were present earlier. They proposed a PF-like index algorithm that takes into account the future rate allocations using round-robin or a few other heuristics. In [START_REF] Nguyen | An algorithm for improved proportional-fair utility for vehicular users[END_REF], we proposed two heuristics that solve in each time-slot the utility maximization problem over a short-term horizon assuming that the means of the future rates are known over this short horizon.

The joint optimization of channel allocation and transmit power control has been investigated for different multiplexing schemes such as like CDMA [START_REF] Subramanian | Joint scheduling and resource allocation in CDMA systems[END_REF] and OFDM [START_REF] Huang | Downlink scheduling and resource allocation for ofdm systems[END_REF]. The proposed algorithms based their decisions on the current channel conditions and previous decisions. In the context of high-speed trains, [START_REF] Xu | Utility-based resource allocation in high-speed railway wireless networks[END_REF] solves the opportunistic utility maximization problem assuming all the future rates are known and with average power constraints.

This paper generalizes the algorithms in [START_REF] Nguyen | An algorithm for improved proportional-fair utility for vehicular users[END_REF] for joint channel allocation and transmit power control in order to maximize the total long-term utility. In addition to this, we also include an average power constraint on the base-station. These algorithms can be seen as being similar in spirit to model predictive control [START_REF] Mesbah | Stochastic model predictive control: An overview and perspectives for future research[END_REF].

Organization

The system model and the optimization problem is presented in Section 2 where we also explain the different channel-gain models. In Section 3, we describe the two algorithms that we propose as well as a baseline locally optimal algorithm which can be seen as gradient-based scheduling. These performance of these algorithms is then evaluated numerically in Section 4. Finally, we state the conclusions and future research direction in Section 5.

Problem formulation

Consider a base station with K mobile users in its coverage range. In time-slot t, the base station transmits to user i with power p i (t). Assume that user i has a channel gain of γ i (t) in slot t. These gains will be assumed to stochastic and independent between vehicles but not necessarily stationary for each user.

The received data rate for user i is computed according to the Shannon formula:

r i (t) = x i (t) log 1 + γ i (t)p i (t) x i (t) , (1)
where x i (t) is the fraction of the channel assigned to user i is slot t. User i gets a utility of U i (z) when it obtains an average rate of z. In order to keep the notation light, we will write r i (t) as a shorthand for r i (x, p, γ). Let S be the appropriate dimensional simplex. We shall the use the notation

[x i (t)] ∈ S (2)
to mean that the vector [x 1 (t), . . . , x K (t)] ∈ S, where S is the K-dimensional simplex.

The utility function will be assumed to concave and differentiable. A widely-used class of utility functions is that of the α-fair functions [START_REF] Mo | Fair end-to-end window-based congestion control[END_REF] that are parametrized by α ≥ 0:

U i (z) = z 1-α 1-α , α = 1; log(z), α = 1.
The special case of α = 1 is also known as the proportional-fair utility function.

The objective of the base station is to choose the power and the channel allocation so as to maximize the total utility of these K users over a horizon of T time slots. That is, the base station solves the optimization problem:

maximize K i=1 U i 1 T T t=1 r i (t) (OPT) subject to [x i (t)] ∈ S, ∀t; (3) 1 T t i p i (t) ≤ P ; (4)
i p i (t) ≤ P max ∀t. ( 5
)
Here P is the average transmit power budget available to the base station, and (4) is the average power constraint. This constraint also makes the problem different from that in [START_REF] Huang | Downlink scheduling and resource allocation for ofdm systems[END_REF] where there was no constraint on the average power.

Remark 1 (Short-term fairness).

A drawback of the utility function defined on the average rate is that if the scheduler knows that, for a particular user, the channel gain may be very high some time in the future then it might wait until this time to serve this user. This user may be starved of allocations in the short-term and the solution may be unfair to this user on short-time scales. One way to resolve the short-term unfairness is to introduce additional quality of service constraints such as requiring each user be scheduled at least once every given number of slots. This constraint was imposed in, for example, [START_REF] Margolies | Exploiting mobility in proportional fair cellular scheduling: Measurements and algorithms[END_REF]. We can also include this constraint in the optimization problem. For simplicity, we do not impose it. We believe the results will be similar as long as this constraint is not very restrictive.

Remark 2 (Fractional channel allocation). In (OPT) we have allowed for fractional channel allocations.

If the system imposes a binary constraint, that is only one user on one channel in any given slot, then these constraints can be imposed in OPT as well as in the algorithms we propose. In the experiments with a logarithmic utility function, we observed that allocations were mostly binary. So, we expect the qualitative conclusions will be valid whether allocations are binary or not.

Remark 3 (Maximum power constraint).

For conciseness, we shall not write the maximum power constraint explicitly in the optimization problems that we will define from now on. This constraint will be implicit and assumed to be applicable in all slots.

The current literature mostly solves (OPT) when the base station is aware of only the channel gains in the current slot 1 . For vehicles sharing their itineraries, partial information on the future channel conditions could also be available to the decision maker in the current time-slot. We assume that this partial information is in the form of the mean of the channel gains in the future slots. One would expect to improve the value of the objective when this information is incorporated in the decision-making.

Channel gains

The heuristics will be evaluated on three different stochastic models for channel gains: (i) the stationary model in which, for each vehicle i, γ i (t) are independent and identically distributed across t with mean γi ; (ii) non-stationary and slowly varying model in which, for each vehicle i, the mean of the channel gains varies on a slow time-scale. For example, for a time-slot of 2 ms, the means may vary every 100 slot or 200 ms; and (iii) a mobility model on a stretch of road, where cars can come and leave as shown in Fig. 2.1. For this model γ i (t) are non-stationary, and also the mean of γ i (t) changes every slot.

Partial information on future channel conditions will be shown to be helpful for the second model and the third model. Nevertheless, the first (or the stationary) model will still be useful to illustrate the benefits of being able to vary the allocated power under an average power constraint.

Algorithms

In this section, we present the three algorithms that will be evaluated in the numerical experiments. The first one is the standard gradient-based scheduling which will be taken as the baseline. The other two are the ones we propose in this paper.

Locally optimal algorithm

The locally optimal algorithm is a special case of the one in [START_REF] Huang | Downlink scheduling and resource allocation for ofdm systems[END_REF] for the scenario in which there is only one sub-channel and the power budget is P in each slot. We remark that the algorithm can be handle the multiple sub-channel case but that in this paper our focus is on the single sub-channel case. When there is no average-power constraint as in [START_REF] Huang | Downlink scheduling and resource allocation for ofdm systems[END_REF], this algorithm chooses maximizes the objective function computed without the knowledge of the future scheduling decisions. That is, in slot t, the scheduler solves maximize

K i=1 U i 1 t t t=1 r i (t) (LA) subject to [x i (t)] ∈ S, ∀t; (6)
i p i (t) ≤ P , ∀t. (7)
Here, the past rate r i (s) for s = 1, 2, ..., t -1 are known to the scheduler, so the decision variables are the channel allocations x i (t) and the transmit powers p i (t) in slot t. The maximum power constraint is set to P in order to make the comparison fair with algorithms which are subject to the average power constraint of P . For the logarithmic utility function, we shall call this the Proportion Fair (PF) solution.

Short-term Objective 1 (STO1)

We begin by defining a big-slot to be B consecutive time-slots. The first heuristic we propose uses future mean channel gains to improve the objective. It includes in its decisions the future channel allocations and transmit powers but only on the scale of big-slots. In each time-slot, the allocations are recomputed for the current time-slot as well as the future big-slots. The intuition behind this is to maximize in each slot the long-term objective based on the best information available. We explain its workings in the context of the non-stationary slowly varying channel model and the mobility model.

For the slowly varying model, B will be the number of time-slots during which mean channel gain remains constant whereas in for the mobility model B can be set by the system designer depending on how fast the mean channel gains vary. Let T = T /B be the number of big-slots in the horizon, and let τt ∈ {1, . . . , T } be the big-slot to which time-slot t belongs to, and let θ t be the number of slots remaining in big-slot τ t not including t, that is θ t = (τ t + 1)B -t.

We shall use the notation x for a quantity that is computed over a big-slot. For example pi (τ ) will denote the power used in all the slots inside big-slot τ . Similarly,

ri (τ ) = B xi (τ ) log 1 + γi (τ )p i (τ ) xi (τ ) (8)
is the total rate obtained by vehicle i in big-slot τ when it is served xi (τ ) fraction of time at a transmit power of pi (τ ). Note that here the rate is computed assuming that the channel gain is its mean value in big-slot τ . With slight abuse of notation,

ri (τ t ) = θ t xi (τ t ) log 1 + γi (τ t )p i (τ t ) xi (τ t ) , (9)
shall denote the total rate in the remaining slots in current big-slot τ t . Also, define

P t = T P - K i=1 t-1 s=1 p i (s)
to be the total remaining power available to the scheduler in slot t.

In each slot, STO1 maximizes

K i=1 U i   1 T   t-1 s=1 r i (s) + r i (t) + T τ =τt ri (τ )     (STO1)
subject to

[x i (t)] ∈ S; [x i (τ )] ∈ S, τ = τt , . . . , T ; (10)
K i=1   p i (t) + θ t pi (τ t ) + T τ =τt+1 B pi (τ )   ≤ P t . (11)
The variables in this problem are [p i (t)] and [p i (τ )], τ = τ t . . . T , and the corresponding channel allocations.

In [START_REF] Chandur | Performance evaluation of scheduling algorithms for mobile wimax networks[END_REF], the LHS is the total transmit power starting from the current slot which has to be less than the remaining power P t .

For the mobile model, instead of solving (STO1) over the whole horizon T , we solve it on a shorter horizon, which is equal to the maximum staying time of the users currently inside the system. This shorter time horizon can vary from one slot to another, and it explains the words 'short-term objective' in the name of the algorithm. The advantages of this is to reduce the computation time which can be helpful when the algorithm has to be executed every 1ms.

Short-term Objective 2 (STO2)

The STO1 algorithm recomputes in each time slot the optimal solution of the future big-slots. In STO2, we recompute the solution of the future big-slots only at the beginning of each big-slot. Inside a big-slot, we compute only the solution for the current slot assuming the solution for the future big-slots to be the same as that computed at the the start of the current big-slot. That is if t ≡ 1 (mod B), STO2 first maximizes

K i=1 U i   1 T   t-1 s=1 r i (s) + T τ =τt ri (τ )     (STO2-Big)
subject to

[x i (τ )] ∈ S, τ = τt , . . . , T (12)
B K i=1 T τ =τt pi (τ ) ≤ P t . (13)
The variables in this problem are [x i (τ )] and [p i (τ )], τ = τ t . . . T . Next, in each slot t, we compute the optimal allocation and transmit power assuming that the allocations and transmit powers in the future big-slots are those computed from solving (STO2-Big). In slot t, STO2 maximizes

K i=1 U i   1 T   t-1 s=1 r i (s) + r i (t) + T τ =τt ri (τ )     (STO2-Small) subject to [x i (t)] ∈ S, [x i (τ t )] ∈ S; ( 14
) K i=1 (p i (t) + θ t pi (τ t )) ≤ P t - K i=1 T τ =τt+1 B pi (τ ). ( 15
)
The variables in this problem are

[x i (t)], [p i (t)], [x i (τ t )] and [p i (τ t )].
As can be seen, in a slot the dimension of the problem is no bigger than the one for STO1. As in STO1, for the mobile model, STO2 solves (STO2-Big) and (STO2-Small) over a shorter horizion which is the maximum staying time of the users currently in the coverage range. The pseudo code for STO2 is shown in Algorithm 1.

Algorithm 1 The STO2 algorithm

t ← 1 while t ≤ T do if t ≡ 0 mod B then Solve (STO2-Big) and obtain x(τ ) and p(τ ).

Solve (STO2-Small) else

Solve (STO2-Small) end if end while

In terms of computational effort, compared to STO1, STO2 solves a lower dimensional problem in each slot except at the starting of every big-slot where the dimension is same as for STO1. Thus, one can except it to be faster but further from the optimal solution.

Numerical experiments

The numerical experiments were run in Python, and all optimization problems were solved using the python package CVXPY [START_REF] Diamond | CVXPY: A Python-embedded modeling language for convex optimization[END_REF] and the solver MOSEK. The results will be presented according to the channel-gain models. For the stationary and slowly varying channel models, we assume that the number of users is fixed. The third model will be evaluated in a dynamic setting in which users will arrive and leave the network.

In all the experiments, P is set to 15 and P max is set to 30. Whenever we show any performance measure of the optimal solution, it will be assumed to mean that the optimal is computed assuming all the future channel gains are known exactly.

Stationary channel

For the first experiment, we take K = 4, that is four users, and a logarithmic utility function for every user. The vector of means [γ i ] = [6.76, 5.45, 4.35, 1.31]. The channel gain in slot t for user i is generated as follows:

γ i (t) = γi A i (t; η) (16)
Here A i (t; η) is a sequence of i.i.d. uniform random variable in the range [1 -η, 1 + η]. We shall refer to η as the noise level. It is assumed that A i (t) and A j (t) are assumed to be independent for i = j. Varying η from 0 to 1 changes the variance of A i (t) from low to high. For η = 0, the channel gains become deterministic and known to the decision maker.

Remark 4. The method for generating γ i (t) need not be necessarily multiplicative as in [START_REF] Mesbah | Stochastic model predictive control: An overview and perspectives for future research[END_REF]. Our heuristics can be used as long as the means of the future channel gains are known. In this paper, we limit the numerical evaluation to the form in [START_REF] Mesbah | Stochastic model predictive control: An overview and perspectives for future research[END_REF].

Figure 4.1 shows the total utility as a function of the noise level η. The time horizon T was taken to be 500 with 5 big-slots, that is, one big-slot has B = 100. If the scheduling slot is 1 ms as in 4G [START_REF] Dahlman | 4G LTE-Advanced Pro and The Road to 5G[END_REF], then the scheduling horizon is of 500 ms. Five sample paths for channel gains were generated, and the plot shows the average of these 5 samples. The label PF is for the local optimization algorithm. The suffix FP attached to STO1 means that STO1 was run wilth a fixed power budget of P in each slot. As expected, allowing for an average transmit power constraint and using future information (even if it is just the mean channel gains) improves the utility. When η = 0, all the algorithms are equivalent and give the same utility since the channel gain is the same is every time slot and is known. Further, performance improvement is more when the noise variance is higher which is again to be expected.

Slowly varying channel

Next, we conduct experiments with the slowly varying channel model. It is assumed that the channel means are constant during B = 100 slots. The optimization horizon is T = 2000, that is there are 20 big-slots in each run. The mean channel gains were first determined for each big-slot. Within a big-slot, the channel gains were then generated using the same method as for the stationary case and given in [START_REF] Mesbah | Stochastic model predictive control: An overview and perspectives for future research[END_REF]. The number of users was again set to 4.

In the first experiment, the mean channel gains are relatively homogeneous with an empirical average of the mean channel gains being [5.98, 5..55, 4.69, 5.30]. Figure 3a, plots the total utility as a function of the noise level while Fig. 3b shows the total transmit power in a slot as a function of time-slot. The data in the latter plot was obtained on a separate run with B = 50, η = 1, and only STO2 and OPT are shown so as to have a more readable figure.

We observe that all algorithms except PF are close to optimal almost throughout. Since the total transmit power in a slot is not far from P , STO1FP is almost as good as STO1. However, in this scenario prediction is still useful as PF is away from OPT even for η = 0, that is when there is no noise.

The mean channel gains in the second experiment are widely varying with the empirical average of the mean channel gains being [2.63, 10.3, 3.76, 0.009]. One user is in a very bad channel state, whereas another one has much better mean channel gains than the others. The plots for the total utility and total transmit power in a slot are shown in Fig. 4.2. Again, the total power was computed from a separate run.

This time we observe that both power control and channel prediction result in improvements. Due to one user being much worse than the others, the optimal transmit power varies much more than in the previous experiment and hits the maximum constraint quite often. Since the fixed-power version is inflexible in this respect, it performs worse for all noise levels.

Mobility Model

Consider a stretch of road of length 1 km covered by one base station, in which vehicles enter from the left and leave on the right (Fig. 2.1). To simplify for illustration, we assume they move with same velocity v = 25 m/s, but the algorithms presented in Sec. 3 do not depend on this assumption. So they stay in the coverage range of the base station in 40 seconds. New vehicles can enter the network only at the start of the big-slots. The length of a big-slot is set to 1 sec. The probability of a new arrival in a big-slot is set to p = 0.3. For each value of noise level (η) we simulate this network for 800 seconds. In this model, the channel gain of mobile user is non-stationary and varies every time slot (1 ms), STO1 and STO2 solves over an horizon of 40 seconds which is the staying time of users in the network. These two aspects make this model different from the two first models. Since the dimension of (OPT) is very large in this case, we do not compute the optimal solution in the experiments. Also, since STO1 solves a high dimensional problem in every slot compared to STO2, it is much slower when the horizon is large. So, we show only the performance of STO2 and PF. Fig. 4.3 illustrates the channel gain curve in the noise and no noise cases. Here, the rate is a function of users' position which is in fact a function of the distance to the base station. Let us take the left margin of the road be 0, and the right margin be 1000, then for position x ∈ [0, 1000] inside the coverage range, the channel gain is equal to f (x) = β(1 + κ exp(|500 -x|/σ). Here, β, κ, σ are adjustable parameters; in our experiment β = 0.01, κ = 40, σ = 200.

Fig. 6 illustrates the numerical results for the above mobile system. The total allocated power is shown only for a small interval of time. Again, we observe that channel prediction leads to a better total utility. Further, by observation, STO2 allocates more power when there are many users (than usual) close to the peak where channel gains are good with high probability.

Conclusions and future work

We proposed two heuristics for joint power control and channel allocation that exploit partial information future channel conditions to improve the utility. Even little information such as mean channel gains is sufficient to observe improvement compared to when no information is used.

This preliminary work opens several directions for future work. First, these algorithms can be generalized to multiple base-station networks. Second, the robustness of these algorithms with respect to errors in future information will be worth investigating. Third, dual-decomposition methods could be applied to get a STO2-like algorithm but with a penalty term instead of the average power constraint. For this relaxed problem, structural results similar to that for gradient-based scheduling could be hoped for. Finally, these scheduling problem can also be studied on the uplink, and similar algorithms can be investigated for the opportunistic utility model.

Figure 1 :

 1 Figure 1: Mobility model.

Figure 2 :

 2 Figure 2: Total utility as a function of noise level. Channel is stationary. Log utility function.

  Total allocated power per slot vs. time-slot

Figure 3 :

 3 Figure 3: Slowly varying channel. Log utility function. Experiment 1. η = 1.

Figure 4 :

 4 Figure 4: Slowly varying channel. Log utility function. Experiment 2. η = 1.

Figure 5 :

 5 Figure 5: Channel gain for two cases: with noise (η = 0.2) and without noise.

  Total allocated power vs. time (in second). η = 0.2.

Figure 6 :

 6 Figure 6: Mobility model. Log utility function.

We shall use slot and time-slot interchangeably to mean the decision making instants.

Acknowledgment

We gratefully acknowledge the funding received from Continental Digital Services France.


Dieter Fiems

Balakrishna Prabhu

Koen De Turck

Travel times, rational queueing and the macroscopic fundamental diagram of traffic flow

published or not. The documents may come

Travel times, rational queueing and the macroscopic fundamental diagram of traffic flow

Introduction

"How long will it take?" is one of the first questions that is asked before one goes on a journey. This question is often accompanied by apprehension for journeys through urban areas during rush hour when the risk of congestion and grid-locks is high. Can we save ourselves the anxiety by estimating average travel times just like we do with weather bulletins? This is one of the main questions that we seek to answer affirmatively in this paper.

The estimates for the travel times will be derived using concepts from queueing theory [START_REF] Cohen | The Single Server Queue[END_REF]. The analogy between customers in queues and cars on roads is a natural one that has been brought up in the modelling of various phenomena in road traffic networks [START_REF] Vandaele | A queueing based traffic flow model[END_REF][START_REF] Chen | A traffic breakdown model based on queueing theory[END_REF][START_REF] Boon | Applications of polling systems[END_REF]. The computation of the travel times will be based on two factors that have a strong influence on the time spent in a queue: (i) the arrival rate of customers (vehicles, in our case); and (ii) the speed with which the queue can be drained. The statistics of the number of vehicles arriving per hour can be obtained from measurements at various locations in a city (see, for example, [5]). The second factor, which is the average departure rate of vehicles, can be computed at two different spatial scales. At the microscopic level of a road segment, the relationship between the density of vehicles on a road segment and the average rate of departure of vehicles is given by what is known as the fundamental diagram in the transportation theory literature [START_REF] Kühne | Transportation Research Circular E-C149 (75 Years of the Fundamental Diagram for Traffic Flow Theory: Greenshields Symposium)[END_REF][START_REF] Lehmann | Distribution function properties and the fundamental diagram in kinetic traffic flow theory[END_REF][START_REF] Kerner | The Physics of Traffic -Empirical Freeway Pattern Features, Engineering Applications, and Theory[END_REF][START_REF] Seyfried | The fundamental diagram of pedestrian movement revisited[END_REF][START_REF] Ambühl | Empirical macroscopic fundamental diagrams: Insights from loop detector and floating car data[END_REF].

At a coarser scale, Daganzo [START_REF] Daganzo | Urban gridlock: Macroscopic modeling and mitigation approaches[END_REF] proposed a macroscopic fundamental diagram that relates the total outflow (vehicles per unit time) from a neighbourhood to the aggregate accumulation (vehicles per unit length of road) for a single-lane system. Figure 1 shows two macroscopic fundamental diagrams from literature, which were obtained from actual measurements in Yokohama [START_REF] Geroliminis | Existence of urban-scale macroscopic fundamental diagrams: Some experimental findings[END_REF] and Toulouse [START_REF] Buisson | Exploring the impact of homogeneity of traffic measurements on the existence of macroscopic fundamental diagrams[END_REF]. The macroscopic fundamental diagrams will serve as the departure rate curves for the computation of journey times in our queueing model.

Starting from the arrival rate statistics and the macroscopic fundamental diagram of a given area, we give a recipe based on queueing theory to derive the average journey time through that area at 1 any given time of the day. The aim of this study is two-fold: (i) estimation of the journey times per se; and (ii) prediction of the influence of public transport in reducing congestion. The first aim has an inherent interest, for example, to plan at what time to leave in order to arrive at a certain hour at the destination. The second aim is motivated by the park and ride schemes introduced by many cities. In these schemes, commuters can park their vehicles at a public transport hub (usually on the outskirts of a city) and then hop onto the public transport network to arrive at their destination inside the city. The journey times by public transport are usually easy to estimate, especially for metros and trams with dedicated lanes. Bus services too are now being provided with dedicated corridors to reduce the influence of traffic jams and to make the commuting times more predictable. If commuters could also obtain an estimate on how long it would take them to drive to their destination instead of taking public transport, then they could choose to switch to public transport if there is time to be gained. Assuming commuters are rational and seek to minimise their travel time, how many of them will choose to switch to public transport? What influence will this have on the congestion on the roads? The answers to these questions can provide guidelines for dimensioning public transport systems especially at peak hours.

Road traffic models have been investigated from different perspectives, including queueing theory [START_REF] Van Woensel | Modeling traffic flows with queueing networks: a review[END_REF], transportation research [START_REF] Daganzo | The cell transmission model -a dynamic representation of highway traffic consistent with the hydrodynamic theory[END_REF][START_REF] Daganzo | The cell transmission model 2. network traffic[END_REF][START_REF] Sumalee | Stochastic cell transmission model (SCTM): A stochastic dynamic traffic model for traffic state surveillance and assignment[END_REF], and statistical physics [START_REF] Chowdhury | Statistical physics of vehicular traffic and some related systems[END_REF][START_REF] Helbing | Traffic and related self-driven many-particle systems[END_REF][START_REF] Nagatani | The physics of traffic jams[END_REF]. Most queueing models divide a road into a number of sections and model each section as a queueing system, such that a stretch of road constitutes a queueing network [START_REF] Van Woensel | Modeling traffic flows with queueing networks: a review[END_REF]. Queueing models may be either stationary [START_REF] Heidemann | A queueing theory approach to speed-flow-density relationships[END_REF], in which case the traffic conditions do not change over time or non-stationary (transient) [START_REF] Boel | A compositional stochastic model for real time freeway traffic simulation[END_REF] when they do. In both cases, these queueing models are explanatory for the fundamental diagram of traffic theory. In transportation research, cell transmission models [START_REF] Daganzo | The cell transmission model -a dynamic representation of highway traffic consistent with the hydrodynamic theory[END_REF][START_REF] Daganzo | The cell transmission model 2. network traffic[END_REF][START_REF] Sumalee | Stochastic cell transmission model (SCTM): A stochastic dynamic traffic model for traffic state surveillance and assignment[END_REF] also divide the road network into sections (cells), and describe the evolution of the number of cars in the cells over time. While in queueing theory, each vehicle spends some random time in a cell, before moving to the next, the transition is deterministic in cell transmission models. In particular, the fundamental diagram is explained by noting that the traffic is constrained by the free space in the next cell, the number of cars in the preceding cell, and a bound on the flow from one cell to the next. The road traffic models from the statistical physics community can be roughly classified in hydrodynamic [START_REF] Lighthill | On kinematic waves ii. a theory of traffic flow on long crowded roads[END_REF][START_REF] Leibig | Pattern-formation characteristics of interacting kinematic waves[END_REF][START_REF] Lee | Steady-state solutions of hydrodynamic traffic models[END_REF][START_REF] Nagatani | Modified kdv equation for jamming transition in the continuum models of traffic[END_REF], gas-kinetic [START_REF] Anderson | On the statistical distribution function theory of traffic flow[END_REF][START_REF] Prigogine | Kinetic theory of vehicular traffic[END_REF][START_REF] Iannini | Kinetic theory of vehicular traffic[END_REF], car-following [START_REF] Helbing | Generalized force model of traffic dynamics[END_REF][START_REF] Shamoto | Car-following model with relative-velocity effect and its experimental verification[END_REF][START_REF] Kurtze | Development of kinks in car-following models[END_REF][START_REF] Jin | On the equivalence between continuum and car-following models of traffic flow[END_REF], and cellular automata [START_REF] Nagel | A cellular automaton model for freeway traffic[END_REF][START_REF] Maerivoet | Cellular automata models of road traffic[END_REF][START_REF] Schadschneider | Cellular automata models of highway traffic[END_REF][START_REF] Balouchi | Finite-size effects in the nagel-schreckenberg traffic model[END_REF][START_REF] Zeng | Research on three-phase traffic flow modeling based on interaction range[END_REF] models. They can be used to study the formation of traffic jams [START_REF] Kerner | Physics of traffic gridlock in a city[END_REF][START_REF] Orosz | Exciting traffic jams: Nonlinear phenomena behind traffic jam formation on highways[END_REF] and traffic waves [START_REF] Ray | Formation of density waves in traffic flow through intersecting roads[END_REF], for routing traffic based on real-time information [START_REF] Yokoya | Dynamics of traffic flow with real-time traffic information[END_REF], or for optimising traffic signals [START_REF] Arita | Signal optimization in urban transport: A totally asymmetric simple exclusion process with traffic lights[END_REF]. An in-depth treatment on these and many other models at both microscopic and macroscopic levels can be found in [START_REF] Kerner | The Physics of Traffic -Empirical Freeway Pattern Features, Engineering Applications, and Theory[END_REF][START_REF] Prigogine | Kinetic theory of vehicular traffic[END_REF][START_REF] Treiber | Traffic Flow Dynamics: Data, Models and Simulation[END_REF][START_REF] Schadschneider | Stochastic Transport in Complex Systems[END_REF]. Some traffic models have been combined with game-theoretic concepts to understand the influence of individual decision making. Examples include the computation of the price of anarchy in routing games [START_REF] Youn | Price of anarchy in transportation networks: Efficiency and optimality control[END_REF][START_REF] Rose | Price of anarchy on heterogeneous traffic-flow networks[END_REF], the study of the dilemma of switching lanes [START_REF] Yamauchi | Dilemma game structure observed in traffic flow at a 2-to-1 lane junction[END_REF][START_REF] Nakata | Dilemma game structure hidden in traffic flow at a bottleneck due to a 2 into 1 lane junction[END_REF], or the explanation of Braess' paradox [START_REF] Sousa | Braess paradox at the mesoscopic scale[END_REF]. In a spirit similar to that of statistical physics, models based on fluid queues have also been proposed for the dynamics of traffic flow in discrete time [START_REF] Vickrey | Congestion theory and transport investment[END_REF] and in continuous time [START_REF] Ban | Continuous-time point-queue models in dynamic network loading[END_REF][START_REF] Jin | Point queue models: A unified approach[END_REF].

Most of these references, including the literature on queueing and cell transmission models, are concerned with flows on fixed-length segments, and do not deal with the distribution of the journey times on the scale of cities or urban areas with the exception of [START_REF] Li | Bounding tandem queuing system performance with variational theory[END_REF] in which bounds have been obtained for two segments. Further, these works also assume First-in-First-Out (FIFO) traffic which is more suited for segments or routes in which the traffic dynamics after the departure time do not influence the journey times. In our model, vehicles arrive according to a Poisson process and have exponentially distributed journey lengths (distance). Different vehicles usually have different destinations within an urban area and it is natural assume that vehicles have different travel distances and that a vehicle that arrives later can leave earlier.

A problem related to ours that has been investigated by the transportation research community is the so-called morning commute problem. Assuming that drivers have to arrive at their common destination on or before their chosen time, the morning commute problem seeks to determine the departure times of the vehicles in order to minimise the waiting time at the destination. One can seek a centralised solution in which a central planner decides the departure times in order to minimise the total cost or a decentralised solution in which each vehicles chooses its departure time with the aim of minimising its own cost. Initiated by Vickery [START_REF] Vickrey | Congestion theory and transport investment[END_REF], this problem has been investigated in various settings: [START_REF] Newell | Traffic flow for the morning commute[END_REF] assumes a hydrodynamic model and a single route, [START_REF] Arnott | Departure time and route choice for the morning commute[END_REF] looks at multiple routes as well as at the impact of tolls, and more recently, [START_REF] Gonzales | Morning commute with competing modes and distributed demand: User equilibrium, system optimum, and pricing[END_REF] adds transit systems to the mix while [START_REF] Ma | The morning commute problem with ridesharing and dynamic parking charges[END_REF] includes ridesharing systems in their model. These papers again are based on assumptions of fixed segment-lengths and FIFO traffic.

Travel time analysis

Consider a neighbourhood of a city, the size of which depends upon the scale at which journey times need to be computed, and for which the arrival rate statistics and the macroscopic fundamental diagram are known. The neighbourhood could, in theory, be the entire city itself. We shall assume that new vehicles arrive in accordance to a Poisson process with time dependent rate λ(t), where this rate typically has a shape as in Figure 2.

The duration of a journey depends on the length of the journey as well as the speed during the journey. We assume that the lengths are independently and identically exponentially distributed with rate 1, taking some other value being equivalent to rescaling time. The speed of vehicles inside this area at a given time instant will depend upon the average density within the area in accordance with the macroscopic fundamental diagram (MFD). Moreover, the number of vehicles that either leave the area or reach their destination will also depend on the density, or equivalently, on the number of vehicles in the area.

Let Q(t) denote the number of vehicles (moving around) in the area at time t and let π n (t) = P[Q(t) = n] denote the probability that there are n vehicles in the area at time t. Let D(t) be the density of vehicles at time t, and let F (D(t)) be the average flow when the density is D(t); the function F is obtained from the MFD. Measurements have shown that the departure rate, µ, (by leaving the area or by reaching the destination) is proportional to the average flow: µ ∝ F (D(t)), see [START_REF] Geroliminis | Properties of a well-defined macroscopic fundamental diagram for urban traffic[END_REF]. Since the traffic density is proportional to the number of vehicles in the area, D(t) ∝ Q(t), we obtain the following relation between the departure rate and the number of vehicles:

µ Q(t) = CF (cQ(t)) ,
with c and C appropriate proportionality constants.

Markov model

With the above assumptions, Q(t) constitutes a Markov chain with state space N. When there are Q(t) = n vehicles in the area at time t, there is an arrival in (t, t + dt] with probability λ(t)dt + o(dt), and a departure with probability µ n dt + o(dt). Hence, the Kolmogorov backward equations read:

πn (t) = λ(t)π n-1 (t) + µ n+1 π n+1 (t) -(λ(t) + µ n )π n (t)
for n ∈ N + , and

π0 (t) = µ 1 π 1 (t) -λ(t)π 0 (t) .
Let W (τ ) denote the remaining travel time of a vehicle in the area at time τ . Since the arrival rate of vehicles depends on time, the travel time of a vehicle will not be stationary. Let w n (τ, t) = P[W (τ ) > t|Q(τ ) = n] denote the complementary distribution function of the remaining travel time W (τ ) of a vehicle in the area, conditioned on having Q(τ ) = n vehicles in the area at time τ . By conditioning on the number of vehicles at time τ + h, we can relate the remaining travel time of vehicles at time τ and τ + h. We obtain the following relation:

w n (τ, t + h) = w n+1 (τ + h, t)λ(t)h + w n (τ + h, t)(1 -λ(t)h -µ n h) + w n-1 (τ + h, t)µ n h 1 - 1 n + o(h) ,
for n ∈ N + . Sending h → 0, we further have

∂ ∂t w n (τ, t) - ∂ ∂τ w n (τ, t) = w n+1 (τ, t)λ(t) -w n (τ, t)(λ(t) + µ n ) + w n-1 (τ, t)µ n n -1 n , (1)
for n ∈ N + . Finally, the expected travel time of a vehicle that arrives at time τ is

W (τ ) = ∞ n=0 π n (t) ∞ 0 w n+1 (τ, t)dt ,
where we used the observation that the mean travel time of a vehicle arriving at time τ while there are n vehicles in the system equals the mean remaining travel time of a vehicle when there are n + 1 vehicles in the area.

Fluid limit

As the fundamental diagram translates to state-dependent service rates, there is no simple analytic solution for the infinite systems of differential equations for π n (t) and w n (τ, t). Moreover, numerically solving these equations is computationally demanding. We however can study the fluid limit of the system, that is, we consider a sequence of models Q N (t), the N th model having arrival rate λ(t)N and study q(t) = lim N →∞ Q N (t)N -1 . The parameter N is the scale at which we are making measurement.

As the scale becomes large, the number of vehicles also grows proportionally, and hence we have to divide Q by N to obtain a quantity which goes to a finite limit. This quantity can be thought of as the density of the vehicles. Thanks to the law of large numbers, the randomness related to the movement of individual vehicles is washed out by the fluid scaling and the distribution of the travel times concentrates on a single point. This reduces the complexity of the analysis since we are no longer analyzing individual cars but are dealing with a large number of cars which behave in a similar way. The fluid limit models the regime when the macroscopic area is much larger than the space occupied by a vehicle, and when the number of vehicles in the area is large, which is typically the case during rush hour. The traffic flow can then be thought of as a fluid with smooth dynamics.

Let µ(x) denote the departure rate in the fluid regime, that is,

µ(x) = lim N →∞ µ N x N .
Then, by writing the evolution of Q(t) as the difference of randomly time-changed Poisson processes and applying the functional strong law of large numbers for Poisson processes, we find that the limiting process q(t) adheres the following ordinary differential equation:

q (t) = λ(t) -µ(q(t)) .
In the limiting regime, there is no uncertainty regarding the number of vehicles in the system over time. Setting w n (τ, t) = w n-1 (τ, t) = w n+1 (τ, t) and n = q(t) in (1), we obtain

∂ ∂t w(τ, t) - ∂ ∂τ w(τ, t) = -w(τ, t) µ(q(τ )) q(τ ) . (2)
Solving for w(τ, t) yields w(τ, t) = e -t+τ τ µ(q(u)) q(u) du = e -t+τ 0 µ(q(u)) q(u) du e -τ 0 µ(q(u)) q(u) du .

The mean travel time of a vehicle arriving at time τ can be computed by integrating w(τ, t) over t, that is,

w(τ ) = ∞ 0 w(τ, t)dt . (3)
Substituting from (2) into (3) and using the fact that -∂ ∂t w(τ, t) is a probability density function whose integral is equal to 1 for all values of τ , we obtain the following differential equation for the mean travel time:

d dτ w(τ ) = w(τ ) µ(q(τ )) q(τ ) -1 , (4)
that has the solution

w(τ ) = w(0)e τ 0 µ(q(u)) q(u) du - τ 0 e τ s µ(q(u)) q(u) du ds, (5)
where

w(0) = ∞ 0 e -t 0 µ(q(u)) q(u) du dt. (6)
The above calculations can be used to predict the travel times, given the arrival rate curve and the macroscopic fundamental diagram. To illustrate the approach, Fig. 3 shows three different arrival curves λ(t), the corresponding occupancy q(t) and travel time w(t) on different panes. The MFD of Yokohama (see Fig. 1) was fitted on a function of the form µ(q) = aqe bq 2 +cq+d and this smooth function was used for all the calculations. The arrival curves were chosen such that the occupancy reached the range 20 -60% during rush hour, which corresponds to the most characteristic part of the MFD.

The total number of vehicles is the same for all arrival curves, the arrival curves only differing in how concentrated the vehicles arrive during rush hour. These results clearly demonstrate that travel times can significantly be reduced if rush hour is slightly stretched. Using the same arrival curves, Fig. 4 shows λ(t) as a function of q(t). The curves are traversed clockwise for increasing values of t. At a given value of q, the derivative of q (t) is much more negative when λ(t) is on the way down. Hence, q decreases more rapidly for the same change in λ, which results in a flatter lower part.

Traffic game

Other than the estimation of travel times of vehicles, ( 5) and ( 6) can also be used for dimensioning the capacity of public transport systems. Consider a scenario in which users have access to travel time estimates both for commuting by personal vehicles as well as by public transport. If the public transport system (PTS) offers a faster ride, a certain number of users may switch from personal vehicles to the PTS. Assuming users are rational, that is, they choose the mode of transport that minimises their travel times, we use concepts from game theory to compute the number (or the fraction) of users that will switch to the PTS. This can be used to estimate the load on the PTS and dimension its capacity accordingly.

For simplicity, assume that taking the PTS incurs a cost, which may depend upon the itinerary and the time of departure but not upon the congestion in the road network. This cost includes the travel time in the public transport system, the travel time to the access points of the public transport system, the inconvenience of using public transport, costs related to possible disruptions of the service, etc. We assume that this cost can be expressed in terms of travel times: the cost of the public transportation system is the maximal travel time one is willing to undergo.

Note that the assumption that the travel cost does not depend on congestion is natural in case the public transportation system has dedicated lines. This can be either completely separated from the road network like a metro network, or consist of separate lanes embedded in the road network.

User equilibrium

We first investigate the choice made by rational users. When the influence of the choice of one user has negligible impact on the travel times of the population, the solution of the game was proposed by Wardrop [START_REF] Wardrop | Some theoretical aspects of road traffic research[END_REF] in the context of choosing routes in road traffic networks. Known as the Wardrop equilibrium, it says that at the equilibrium choice of routes is such that the cost on each used route is the same and the cost on the unused ones is larger than that of the used ones. We shall apply the concept of Wardrop equilibrium to our time-dependent system to compute the fraction of users who switch to the PTS.

Let p(t) denote the fraction of users that prefer their vehicles over public transport, such that the arrival traffic equals λ(t; p) = λ(t)p(t) .

Moreover, let w(t; p) be the travel time of users arriving at time t that do not choose public transport. Following the arguments of the preceding section, we have q (t; p) = λ(t; p) -µ(q(t; p)) ,

w(τ ; p) = ∞ 0 e -t+τ 0 µ(q(u;p)) q(u;p) du e -τ 0 µ(q(u;p)) q(u;p) du dt .
For the computation of p, we use the same argument as that for a Wardrop equilibrium. For a given p, whenever w(t; p) < C(t), vehicles arriving at time t are better off not using the public transportation system. Hence, unless p(t) = 1, rational behaviour requires that more users opt for private transportation. Similarly, when w(t; p) > C(t), rationality requires that more users opt for public transportation. Summarising, the optimal proportion of users that choose private transportation, p • (t), should satisfy the following constraints:

     w(t; p • ) < C(t) for p • (t) = 1 , w(t; p • ) = C(t) for 0 < p • (t) < 1 , w(t; p • ) > C(t) for p • (t) = 0 . (7)
This means that at the user equilibrium, a fraction p • (t) of users will use private transportation whereas 1 -p • (t) will use the PTS. The function p • has no nice and easy analytical solution so we resort to the following iterative numerical algorithm for its computation. Given the arrival curve λ(t) and the fundamental diagram µ(x), we set p(t) = 1, and then update it according to p(t) ← min 1, p(t)e -α(ω(t;p)-C(t)) , such that p → p • . Here α is a parameter which determines the speed of convergence of the recursion. Small α corresponds to conservative updates, such that convergence is slow. In contrast, for large α, the update rule may never converge.

To illustrate the approach, Fig. 5 shows the rational traffic arrival curves λ(t; p * ), the corresponding occupancy q(t; p * ) and travel time w(t; p * ) on different panes for different values of the cost C as indicated. Note that for C = ∞ the PTS is not used, i.e. these curves represent the case that there is no PTS (and no game to be played). The arrival curve (before the game) corresponds to the solid arrival curve from Fig. 3, and the MFD of Yokohama (see Fig. 1) is again used for all calculations. One easily verifies that the travel times indeed do not exceed the fixed cost C of the public transport. Depending on the cost, the PTS is a viable option from the start of the rush hour (C = 20), from the peak of the rush hour (C = 30) or at the end of the rush hour (C = 40). Even when the PTS cost is high and only few people use the PTS, the impact on travel times is considerable. This is in line with the observations of Fig. 3, which showed that a small shift in traffic has a huge impact on travel times.

Social optimum and Price of Anarchy

To benchmark the cost of free transportation choice, we also calculate the socially optimal traffic mix. The total transportation cost K(t) at time t, assuming policy p equals,

K(t; p) = w(t; p)p(t) + (1 -p(t))C(t) .
The socially optimal policy is then given by,

p * (t) = arg inf p T2 T1 λ(t)K(t; p)dt (8)
where [T 1 , T 2 ] is the rush hour interval. The computation of the social optimum is not exactly an optimal control problem as one would guess at first. The difference is that w depends upon the Figure 6: Socially optimal arrival intensity λ(t; p * ), mean occupancy q(t; p * ) and mean travel time w(t; p * ) for different public transport costs as indicated, and for the Yokohama fundamental diagram. whole trajectory of q through w(0) (see [START_REF] Kühne | Transportation Research Circular E-C149 (75 Years of the Fundamental Diagram for Traffic Flow Theory: Greenshields Symposium)[END_REF]). This makes it difficult to define the dynamics and hence the Hamiltonian locally as is done in a classical optimal control problem formulation. Define x(t) = [q(t), w(t)]. We can then approximate (8) by the solution of following optimal control problem:

min p(t)∈[0,1] T2 T1 λ(t)K(t; p)dt , (9)
subject to

ẋ(t) = λ(t)p(t) -µ(q(t)), w(t) µ(q(t)) q(t) -1 , (10) q
(0) = q 0 , w(T 2 ) = w(∞), (11)
where q 0 = q(∞) = µ -1 (λ(∞)) and w(∞) = q 0 /λ(∞) are the steady-state values of q and w when there is no congestion. The inverse of µ is the one closer to 0 to reflect the lack of congestion. The constraint on w is taken on its final value to reduce the influence of the whole trajectory of q on this constraint. We solved (8) using an algorithm similar to gradient descent where the function p(t) is perturbed and the new function is accepted if it yields a lower social cost. The approximate problem (9) was solved using BOCOP which is a freely available solver [START_REF]BOCOP: an open source toolbox for optimal control[END_REF]. Both algorithms gave similar results. From the optimal control formulation, it can be seen that if λ(t) goes above a certain value, λ(t; p * ) becomes constant, the value of this constant being equal to the solution of (8) when q(t) does not depend on t. More precisely, for fixed cost C(t) = C, the constant is equal to µ(q) where q is the positive fixed point (which exists and is unique for the MFD at hand) of the functional equation, q = µ(q) C .

Recall that p • (t) is the solution of the set of equations [START_REF] Lehmann | Distribution function properties and the fundamental diagram in kinetic traffic flow theory[END_REF]. A measure of inefficiency of individual decision making is the Price of Anarchy [START_REF] Koutsoupias | Worst-case equilibria[END_REF] which is the ratio of the cost at the game equilibrium to the optimal social cost, PoA = T2 T1 λ(t)K(t; p • )dt T2 T1 λ(t)K(t; p * )dt . The socially optimal arrival curves λ(t, p * ) and the corresponding occupancy q(t, p * ) and travel times w(t, p * ) are depicted in Fig. 6. We use the same parameters as in Fig. 5, so that we can compare the rational and socially optimal solutions. It is readily observed that a considerably larger fraction of the traffic switches to the PTS in the socially optimal solution, both when the PTS cost is low (C = 20) and high (C = 40). To calculate the price of anarchy, we consider the cost in the interval [START_REF] Wardrop | Some theoretical aspects of road traffic research[END_REF]150] minutes. We find a 23.2%, 35.2% and 46.8% increase of the cost for C = 20, C = 30 and C = 40, respectively.

Conclusion

We have proposed a queueing-theory-based model for computing the time-dependent travel times during rush hour in a congested area. In contrast to classical queueing models, we account for the macroscopic fundamental diagram of traffic flow by including departure rates which depend on the queue content. As the traffic intensity is high during rush hour, we study the fluid limit of the queueing model at hand. Numerical experiments show that the travel times can be reduced significantly by reducing a small fraction of traffic on the roads. Public transport systems can thus be an effective means of reducing congestion in cities. Relying on game theory, we show how our model can be used to assess the impact of public transport systems on congestion and compare the game-theoretical and socially optimal traffic mixes.

As the current work focuses on the fluid limit, only average travel times are accounted for when comparing public and private transport. One can argue that other travel time metrics like the travel time variance or the complete travel time distribution should be accounted for as well. For example, one can think of a scenario where users opt for public transport if the risk/probability of overly long travel times exceeds a certain value. The tools for approximating such metrics (diffusion limits, control of diffusions) are different and considerably more complicated and therefore considered subject for further study.

Figure 1 :

 1 Figure 1: Macroscopic fundamental diagram

Figure 2 :

 2 Figure 2: A two-way arrival rate curve (see, e.g., [5] which has statistics from the Vancouver area in Canada)

Figure 3 :

 3 Figure 3: Arrival intensity λ(t), mean occupancy q(t) and mean travel time w(t) for different arrival intensity curves, and for the Yokohama fundamental diagram.

Figure 4 :

 4 Figure 4: λ(t) vs. q(t) (increasing t corresponds to clockwise motion along the curve

Figure 5 :

 5 Figure5: Rational arrival intensity λ(t; p • ), mean occupancy q(t; p • ) and mean travel time w(t; p • ) for different public transport costs, and for the Yokohama fundamental diagram.


Baudouin Denis De Senneville

Nora Frulio
email: nora.frulio@chu-bordeaux.fr

Hervé Laumonier
email: herve.laumonier@gmail.com

Cécile Salut
email: cecile.salut@chu-bordeaux.fr

Luc Lafitte
email: luclafitte@gmail.com

Hervé Trillaud
email: herve.trillaud@chu-bordeaux.fr

Liver contrast-enhanced sonography: Computer-assisted differentiation between focal nodular hyperplasia and inflammatory hepatocellular adenoma by reference to microbubble transport patterns

Keywords: Ultrasound Imaging, Adenoma, Perfusion imaging, Computer-Assisted Diagnosis, Retrospective studies Ultrasound CEUS: contrast-enhanced ultrasound FNH: Focal nodular hyperplasias HCA: Hepatocellular adenomas I-HCA: Inflammatory hepatocellular adenoma CT: Computed tomography MRI: Magnetic resonance imaging T: Tesla

Objective: A new computer tool is proposed to distinguish between focal nodular hyperplasia (FNH) and an inflammatory hepatocellular adenoma (I-HCA) using contrast-enhanced ultrasound (CEUS). The new method was compared with the usual qualitative analysis.

Methods:

The proposed tool embeds an "optical flow" algorithm, designed to mimic the human visual perception of object transport in image series, to quantitatively analyse apparent microbubble transport parameters visible on CEUS. Qualitative (visual) and quantitative (computer-assisted) CEUS data were compared in a cohort of adult patients with either FNH or I-HCA based on pathological and radiological results. For quantitative analysis, several computer-assisted classification models were tested and subjected to cross-validation. The accuracies, area under the receiver-operating characteristic curve (AUROC), sensitivity and specificity, positive predictive values (PPVs), negative predictive values (NPVs), false predictive rate (FPRs) and false negative rate (FNRs) were recorded.

Results: Forty-six patients with FNH (n = 29) or I-HCA (n = 17) with 47 tumors (one patient with 2 I-HCA) were analysed. The qualitative diagnostic parameters were: accuracy = 93.6%, AUROC=0.94, sensitivity = 94.4%, specificity = 93.1%, PPV = 89.5% and NPV = 96.4%, FPR = 6.9%, FNR = 5.6%. The quantitative diagnostic parameters were: accuracy = 95.9%, AUROC = 0.97, sensitivity = 93.4%, specificity = 97.6%, PPV = 95.3%, and NPV = 96.7%, FPR = 2.4%, FNR = 6.6%.

Conclusions: Microbubble transport patterns evident on CEUS are valuable diagnostic indicators.

Machine-learning algorithms analysing such data facilitate the diagnosis of FNH and I-HCA tumours.

Key Points:

-Distinguishing between focal nodular hyperplasia and an inflammatory hepatocellular adenoma using dynamic contrast-enhanced ultrasound is sometimes difficult.

-Microbubble transport patterns evident on contrast-enhanced sonography are valuable diagnostic indicators.

-Machine-learning algorithms analysing microbubble transport patterns facilitate the diagnosis of FNH and I-HCA.

-The technique offers a potential future means for accurately characterizing other hepatic lesions, potentially obviating the need for biopsy or surgical resection.

GB: Gigabit RAM: Random access memory

INTRODUCTION

Benign hepatocellular tumours are rare, constituting 10% of all hepatic tumours [START_REF] Cherqui | Management of focal nodular hyperplasia and hepatocellular adenoma in young women: a series of 41 patients with clinical, radiological, and pathological correlations[END_REF]. Two large groups of benign hepatocellular tumours can be distinguished: reactive regenerative lesions (focal nodular hyperplasias [FNHs]), and tumoural lesions (hepatocellular adenomas [HCAs]). Both lesions are most common in young females [START_REF] Cherqui | Management of focal nodular hyperplasia and hepatocellular adenoma in young women: a series of 41 patients with clinical, radiological, and pathological correlations[END_REF]. Diagnostic imaging is essential to guide treatment decisions, which range from no treatment to surgical resection or confirmatory biopsy.

Traditionally, multiphase computed tomography (CT) or magnetic resonance imaging (MRI) has been used for detailed evaluation of hepatic lesions. However, the high-level radiation associated with multiphase CT and the limited accessibility of MRI have rendered dynamic contrast agent-enhanced ultrasound (CEUS) an attractive, safe, non-invasive, accurate, and economic tool for evaluating hepatic lesions (2)- [START_REF] Onofrio | Contrast Enhanced Ultrasound of Focal Liver Lesions[END_REF]. Although the appearance is not always typical in some cases, both FNH and HCA demonstrate typical, reproducible, arterial phase enhancement patterns on CEUS in most cases. The diagnostic criteria for FNH are a hyper-enhancing lesion in the arterial phase with rapid centrifugal filling from a central vessel, and radial vascular branches (the "spoke and wheel" sign) (2) [START_REF] Baranes | Imaging benign hepatocellular tumors: atypical forms and diagnostic traps[END_REF] and also sustained enhancement in portal and late phase [START_REF] Bartolotta | Hepatic focal nodular hyperplasia:contrast-enhanced ultrasound findings with emphasis on lesion size, depth and liver echogenicity[END_REF]. HCAs constitute a heterogeneous group of tumours exhibiting multiple histological subtypes (inflammatory, with FNH1A or catenin gene mutations, or unclassified) [START_REF] Bioulac-Sage | Hepatocellular adenoma subtype classification using molecular markers and immunohistochemistry[END_REF]. On CEUS, HCAs are hyper-enhancing in the arterial phase; the enhancement pattern commences peripherally and exhibits rapid centripetal filling; this pattern is characteristic of 86-90% of all inflammatory HCAs (I-HCAs). Other HCA subtypes exhibit iso-vascularity or moderate hyper-vascularity, with mixed filling patterns in the arterial phase (2) [START_REF] Laumonier | Role of contrast-enhanced sonography in differentiation of subtypes of hepatocellular adenoma: correlation with MRI findings[END_REF]. In clinical practice, it is essential to distinguish FNH from adenoma to ensure appropriate management. Confirmed FNHs are managed conservatively (with regular follow-up); HCAs require cessation of oral contraceptive use, (commonly) biopsy, and either surgery or (at least) follow-up imaging. I-HCA show the most important hypervascularity and 10-15% of I-HCA are also found to be β-catenin activated with a risk for malignant transformation.

Distinguishing between FNH and I-HCA using CEUS is sometimes difficult because both lesions evidence hyper-enhancement during the arterial phase and it can be challenging to qualitatively differentiate centrifugal from centripetal tumour filling, particularly for larger nodules. Computerassisted methods are thus required for quantitative spatiotemporal assessment of organ perfusion. Such techniques must be faster and more reproducible than visual analysis, and must lack learning curves. Efforts have been made to quantify enhancement parameters in vascular compartments as indicators of several pathological conditions [START_REF] Tranquart | Perfusion quantification in contrast-enhanced ultrasound (CEUS)-ready for research projects and routine clinical use[END_REF][START_REF] Rognin | Parametric imaging of dynamic vascular patterns of focal liver lesions in contrast-enhanced ultra-sound[END_REF][START_REF] Dietrich | An EFSUMB introduction into dynamic contrast-enhanced ultrasound (DCE-US) for quantification of tumour perfusion[END_REF][START_REF] Strouthos | Indicator dilution models for the quantification of microvascular blood flow with bolus administration of ultrasound contrast agents[END_REF][START_REF] Mischi | Angiogenesis imaging by spatiotemporal analysis of ultrasound contrast agent dispersion kinetics[END_REF]. In particular, transport equations have been recently derived to estimate microbubble velocity at the time of bolus contrast arrival [START_REF] De Senneville | Development of a fluid dynamic model for quantitative contrast-enhanced ultrasound imaging[END_REF]. In practice, a so-called "optical flow" algorithm is employed to mimic the human visual perception of microbubble transport in CEUS [START_REF] Corpetti | Dense estimation of fluid flows[END_REF][START_REF] Horn | Determining optical flow[END_REF][START_REF] Zachiu | An improved optical flow tracking technique for real-time MR-guided beam therapies in moving organs[END_REF]. Here, we use this approach to quantitatively distinguish between FNH and I-HCA. We quantify divergence (sources and sinks), curling (shearing), amplitudes, and convergence towards the centre of tumour (centrifugal/centripetal nature) in dense transport fields [START_REF] Corpetti | Dense estimation of fluid flows[END_REF]; these are very simple indicators of displacement vector directions, orientations, and magnitudes. In turn, these serve as inputs to a binary FNH/I-HCA classifier.

The purpose was to compare, as a preliminary study, the original concept of computer-assisted method with the usual qualitative analysis for the diagnosis of two benign hepatocellular tumours (FNH and I-HCA) with hypervascularity during the arterial phase of the CEUS.

MATERIALS AND METHODS

Study design and population

In this retrospective single-centre study conducted from July 2005 to July 2018, we identified images from patients who underwent CEUS and were (otherwise) definitively diagnosed with FNH or I-HCA. We included I-HCA patients who had been histologically diagnosed [START_REF] Onofrio | Contrast Enhanced Ultrasound of Focal Liver Lesions[END_REF] and FNH patients diagnosed based on commonly accepted MRI criteria (3), imaging follow-up, or histologically. All MRIs were performed using a 1.5-T machine running a published imaging protocol [START_REF] Laumonier | Hepatocellular adenomas: magnetic resonance imaging features as a function of molecular pathological classification[END_REF] [START_REF] Van Aalten | Hepatocellular adenomas: correlation of MR imaging findings with pathologic subtype classification[END_REF]. The study adhered to all local regulations and data protection agency recommendations (the National Commission on Informatics and Liberty (CNIL) dictates). Patients have been informed for the use of their data anonymously.

Demographic characteristics

We enrolled 46 patients (Table 1) with the inclusion criteria, 29 had FNH and 17 I-HCA (18 I-HCA tumours were analysed because one patient had two tumours). Of 29 FNH patients, 23 (79%) were female and the median age was 44 (21-61) years; of 17 I-HCA patients, 16 (94%) were female and the median age was 40.5 (21-66) years. The median diameters of FNH and I-HCA lesions were respectively 2.9 (3-10) and 6.9 (3.4-12) cm. Histological data of the 18 I-HCA tumours were available for 15 surgical specimens and 3 percutaneous biopsies Histological data on 7/29 FNH tumours (24%) were available (percutaneous biopsy, six samples; one surgical sample); imaging follow-up data were available for 15/22 patients without histological diagnosis (68%) with a median follow-up of 12 (4-84) months ; CEUS was performed using Sequoia (n = 37), S2000 (n = 4), and S3000 (n = 5) instruments.

Histological analysis

Histological samples were obtained by biopsy or during surgical resection; for ethical reasons, no samples were taken purely for the purpose of this study; clinical indications were required. All analyses were performed as previously described (8)(9) [START_REF] Laumonier | Hepatocellular adenomas: magnetic resonance imaging features as a function of molecular pathological classification[END_REF], in the same laboratory.

CEUS protocol

CEUS was performed by abdominal radiologists who had 5-10 years of experience. Each patient received a bolus injection of ultrasound contrast agent (SonoVue, Bracco). Contrast-enhanced sequences were obtained using dedicated, low mechanical index (MI) contrast-imaging software (MI < 0.2) employing one of three US machines (Sequoia, S2000 and S3000; a Siemens Medical Solution instrument featuring Cadence Contrast Pulse Sequencing [CPS]; and a Convex Array 4C1-S probe). Standard pre-settings were used, but it was possible to adjust settings for individual patients. SonoVue was injected intravenously as a bolus of 2.4 mL via a 20-gauge cannula into the antecubital vein, followed by flushing with 5 mL normal saline. Digital cine clips showing dynamic contrast enhancement within the lesion and surrounding liver tissue were continuously recorded, commencing 5 s before SonoVue injection and covering the arterial (10-45 s post-injection), portal (60-90 s), and late (120-150 s) phases. Injection was repeated using the same dose (2.4 mL SonoVue) if the data were of poor quality. All sequences were digitally stored. Intra-tumoural vascular geometry and lesional enhancement patterns were evaluated.

CEUS analysis of lesional type

Qualitative visual analysis

Data were reviewed in consensus by two abdominal radiologists blinded to pathological and MRI diagnoses. Each lesion was classified using pre-defined criteria for FNH and I-HCA. For FNH, these were hyper-enhancement in the arterial phase, with rapid centrifugal filling; (usually) an obvious central vessel and radial vascular branches (especially in larger lesions; the "spoke and wheel" sign); and iso-or hyper-enhancement in the portal and venous phases, without washout. For I-HCA, the criteria were hyper-enhancement in the arterial phase, frequently accompanied by rapid centripetal filling; no radial vascular structure; and iso-or hyper-enhancement in the portal and venous phases, without washout (3) (9) [START_REF] Quaia | Characterization of focal liver lesions with contrast-specific US modes and a sulfur hexafluoride-filled microbubble contrast agent: diagnostic performance and confidence[END_REF].

Computer-assisted quantitative analysis using a transport equation model

Microbubble transport fields in lesions were estimated (using a transport equation) on a pixel-bypixel basis employing the so-called "optical flow" process [START_REF] De Senneville | Development of a fluid dynamic model for quantitative contrast-enhanced ultrasound imaging[END_REF]. The "optical flow" problem has long been studied by vision scientists in efforts to analyse general visual motion in images of a moving target [START_REF] Corpetti | Dense estimation of fluid flows[END_REF] [START_REF] Horn | Determining optical flow[END_REF]. For each lesion, the absolute changes in four image-based displacement indicators were calculated: (i) the divergence δ (reflecting the presence of sources and sinks); (ii) the curl ρ (reflecting local vortices); (iii) the amplitude γ (reflecting the magnitude of apparent displacement); and, (iv) the centripetal nature τ (reflecting the flow field convergence towards the centre of tumour). The analysis was restricted to a region of interest, manually drawn on a high contrast CEUS image, encompassing the tumour. The analytical window size was fixed at 2 s commencing at the bolus arrival time, and thus covered the filling phase. The reader is referred to Appendix A for additional information on numerical resolution and implementation. All computerassisted analyses were blinded to pathological data.

Statistical analysis

The accuracies, area under the ROC curve (AUROC), sensitivity, specificity, positive predictive values (PPVs), negative predictive values (NPVs), false predictive rate (FPRs) and false negative rate (FNRs) of qualitative and quantitative analyses were recorded (we considered the diagnostic of an adenoma as a "positive case" in the scope of this study).

For quantitative analyses, using one of the four microbubble displacement indicators (δ, ρ, γ or τ) as an input feature, we developed machine learning models to differentiate between FNH and I-HCA.

For this binary classification task, the following four machine learning algorithms were applied using the commercial software Matlab (©1994-2019 The MathWorks, Inc.)/"Statistics and Machine Learning" toolbox: random forest (RF), k-nearest neighbour (KNN), support vector machine (SVM), and logistic regression (LR). Default hyperparameters in Matlab implementations were employed. We refer the interested reader to [START_REF] Kohavi | A Study of Cross-Validation and Bootstrap for Accuracy Estimation and Model Selection[END_REF] [START_REF] Cantor | Determining the Area under the ROC Curve for a Binary Diagnostic Test[END_REF] for additional information about abovementioned computer-assisted classification algorithms. We evaluated the diagnostic performances through self-validation (the complete 47-tumours set was used for both train and test samples) and through 10-fold-stratified cross-validation (the 47-tumours set was randomly partitioned into complementary 90%-training and 10%-test subsets). The cross-validation steps were repeated 100 times with shuffling of the folds and test metric averages calculated. We also compared the medians and interquartile ranges of all four indicators using the unpaired Mann-Whitney U-test. A p-value < 0.025 was considered to reflect statistical significance.

RESULTS

Qualitative CEUS analysis

FNH and I-HCA were correctly identified via qualitative CEUS in 27/29 and 17/18 tumours, respectively (accuracy = 93.6%, AUROC=0.94, sensitivity = 94.4%, specificity = 93.1%, PPV = 89.5%, NPV = 96.4%, FPR = 6.9%, FNR = 5.6%; Table 2, first row).

Quantitative CEUS analysis

Figures 1 and2 show typical microbubble transport fields as revealed by dynamic contrast imaging; one clip (Fig. 1) is from an FNH patient and the other (Fig. 2) from an I-HCA patient. Of the four tested transport indicators, divergence and centripetal indicators differed most significantly between the two populations (Mann-Whitney test, p-value = 2 × 10 -4 for divergence, and 1 × 10 -7 for centripetal indicator) (Figs. 3 and4). The centripetal indicator served as a valuable binary classifier in all tested machine learning systems (Table 2). In particular, using the naïve Bayes classifier applied on the centripetal indicator, the diagnostic parameters were: accuracy = 95.9%, AUROC = 0.97, sensitivity = 93.4%, specificity = 97.6%, PPV = 95.3%, NPV = 96.7%, FPR = 2.4%, FNR = 6.6% (in average over the 100 cross-validation steps, FNH and I-HCA were thus correctly identified in 28.3/29 and 16.8/18 tumours, respectively).

DISCUSSION

We show that the dense transport fields provide valuable kinetic information in CEUS time series; the results are more accurate than those of qualitative visual analysis. Using the qualitative analysis, one false negative case and two false positive cases were to deplore. Concerning the false negative case, the filling direction was difficult to determine visually. Concerning the two false positive cases, one tumour (44 mm) presented two feeding pedicles, and one (22 mm) underwent a too fast filling. For both FNH tumours it was also difficult to appreciate visually the centrifugal filling. A quantitative approach delivers reproducible results and minimises operator dependency, as visual interpretation of CEUS images lacks a learning curve when the process is automated. Our method deals with the intrinsic variations in spatio-temporal greys that are inevitable during dynamic imaging. This allows numerical access to visual perceptions of microbubble trajectories. We used four simple indicators (δ, ρ, γ, and τ) of transport field direction/orientation and amplitude. The amplitude and curl indicators were not useful (Figs. 4b,c), whereas the divergence and centripetal indicators were (Fig. 4a,d). Best results were obtained using the indicator τ which best fits the initial centrifugal/centripetal tumour filling hypothesis (Figs. 4d). For its part, the divergence operator gave decent results. In theory, the divergence of any vector field is positive for sources (centrifugal trajectories) and negative for sinks (centripetal trajectories). The divergence operators were positive for both FNH and I-HCA data; bolus arrival manifested as one or several sources of microbubbles. However, for I-HCA lesions, the divergence operator was modulated by centripetal filling, whereas the divergence operator was enhanced by centrifugal filling in FNH patients.

Using our quantitative method, diagnosis is near-instantaneous once the region of interest (encompassing the tumour) is delineated. Although the duration of the temporal window for the analysis must be sufficient to cover the filling phase, 2 s was adequate; this is a great advantage, eliminating all long-term bias imparted by probe motion, and respiratory and other motion artefacts [START_REF] De Luca | Estimation of large-scale organ motion in B-mode ultrasound image sequences: A survey[END_REF] [START_REF] Pratikakis | Robust multiscale deformable registration of 3d ultrasound images[END_REF].

Several limitations of our work must be mentioned, particularly the small sample size. This was a single-centre retrospective study lacking an external validation cohort. Considering that only two categories of focal liver lesions were examined (FNH and I-HCA), an inherent overestimation of both qualitative and quantitative analysis has to be taken into account. Also, mean tumour diameter was significantly smaller in the FNH group, associated with recruitment bias: only patients with histological diagnoses obtained after surgical resection or via percutaneous biopsy were included in the I-HCA group. However, in our centre, when an I-HCA tumour is identified using MRI [START_REF] Bise | New MRI features improve subtype classification of hepatocellular adenoma[END_REF] or CEUS, a pathological analysis is performed only when the tumour diameter is > 3 cm. Thus, the I-HCA group featured only tumours that met this criterion, unlike the FNH group, for which tumours of all diameters (including small tumours) were evaluated. Also, the fact that any US artefacts can intrinsically be interpreted as "false" motions by the transport equation constitutes a major source of uncertainty. This may bias the microbubble estimations in transport fields, in turn affecting all four image-based indicators. This is also of concern when brief US "shadow" artefacts develop in obese patients (one of our cohort was obese and was constantly miss-classified by our quantitative approach due to poor image quality). Similarly, in-plane and/or out-of-plane organ motion within the image field-of-view must be no more than moderate. Please note that when organ motions are large or complex, it is possible to "pop" microbubbles on-line to virtually repeat the imaging session. Alternatively, image post-processing strategies may be valuable (24)(28) (Appendix A).

Finally, manually drawn masks encompassing lesions must exclude adjacent feeding arterials; otherwise, the estimated displacement is likely to be calculated from the border to the centre of the tumour, compromising FNH diagnosis.

In conclusion, this proof-of-concept study indicates that microbubble displacements evident on CEUS can be used to efficiently diagnose FNH/I-HCA lesions. Machine learning allows for computer-assisted diagnoses. In the future, we will optimise the model ( 28 Quantitative indicators are shown with standard deviations and 95% confidence intervals in parentheses. images, which are not prone to contrast enhancement, are used to this end. We estimated translational displacements restricted to the binary mask Γ. We used a gradient-driven descent algorithm maximising the inter-correlation coefficients. This translation was used to compensate for displacement of imaged tissues on CEUS images prior to microbubble transport estimation employing Eq. ( 2).

Pixel-wise analysis of dense flow fields

We next calculated a pixel-wise understanding of flow directions/orientations and amplitude, as follows:

Maps of sources and sinks: Sources and sinks in the transport were analysed using the divergence operator. Mathematically, the divergence of a two-dimensional vector is expressed as: div

The final, discrete divergence operator employed in numerical implementation was:

(div ) i,j,t (4)
where denotes the pixel coordinates and t the frame acquisition time. We emphasise that positive and negative values are associated with pixels located near sources and sinks, respectively.

Vortex maps:

Local vortices in the estimated transport vectors were analysed with the aid of the curl operator. Mathematically, the curl of is expressed as: curl

The resulting discrete curl operator is:

(curl ) i,j,t [START_REF] Onofrio | Contrast Enhanced Ultrasound of Focal Liver Lesions[END_REF] Amplitude maps: The amplitudes of estimated transport vectors were calculated as Euclidian distances. Mathematically, the magnitude of is expressed by: ( 7)

Centripetal indicator maps: The convergence of estimated transport vectors towards the centre of tumour were calculated with the aid of the scalar product. The cosine of the angle formed by two vectors and is expressed by: (

In our study, is any vector of the estimated flow field, and has the same origin as , but the extremity located at the gravity centre of tumour (i.e, the centre of mass of the binary mask Γ). That way, the cosine of the angle lies in intervals [-1,0] and [0,1] for centrifugal and centripetal , respectively.

Quantitative analysis of dense flow fields

As described above, we created four sets of pixel-wise maps (of divergence, curl, amplitude and convergence towards the centre of tumour). Each set was then simplified to a single parameter as follows. The spatiotemporal averages of each map were individually computed under a mask defining the imaged tissue (i.e., Γ) within the relevant time window. The duration spanned by that window is termed ΔT below and commenced at the bolus arrival time t 0 . Spatiotemporal averaging was weighted by the grey level intensity in CEUS image I; thus, the values for scenarios exhibiting identical microbubble transport behaviours were identical irrespective of the numbers of pixels evidencing microbubbles. The divergence and curl operators were termed δ and ρ, respectively. We measured the absolute value of curl; thus, the direction of vortex rotation didn't affected analysis.

The centripetal indicator was termed τ.

(9) (10) [START_REF] Rognin | Parametric imaging of dynamic vascular patterns of focal liver lesions in contrast-enhanced ultra-sound[END_REF]

Fig. 1 .

 1 Fig. 1. Typical results obtained when evaluating an FNH lesion. Data obtained at different CEUS

Fig. 2 .

 2 Fig. 2. Typical results from a patient with an I-HCA lesion. The manually drawn mask

Fig. 3 .

 3 Fig. 3. Boxplots of indicators of the four dense transport fields (divergence δ [a], curl ρ [b],

Fig. 4 .

 4 Fig. 4. ROC curves obtained using the four quantitative indicators (divergence δ [a], curl ρ [b],

Table 1 .

 1 ), enrol larger patient cohorts, include other lesional features, and study other pathologies. Demographic characteristics. Median values of age and tumor size are shown with standard deviations and minimum-maximum intervals in parentheses. Statistical comparison of age and tumor size between FNH and I-HCA groups were performed using an unpaired Mann-Whitney U-test (last column), a p-value < 0.025 was considered to reflect statistical significance.

	Characteristic	FNH	I-HCA	Total	Statistical significance (p-value)
	Age	44 ± 11 (21-61)	40.5 ± 11 (21-66)	42 ± 11 (21-66)	No (0.16)
	Gender (F/M)	23/6	16/1	39/7	-
	Tumor size	29 ± 16 (13-100)	60.5 ± 29 (34-126)	36 ± 28 (13-126)	Yes (10 -6 )
	Histological data (available/not available)	7/22	16/1	23/23	-

Random Forest 80.2 ± 16.2 (79.2-81.3) 0.70 ± 0.28 (0.69-0.72) 61.6 ± 38.6 (59.1-64.0) 91.9 ± 20.9 (90.6-93.3) 71.4 ± 40.7 (68.8-74.0) 80.8 ± 20.0 (79.5-82.1) Centripetal indicator (γ) Logistic Regression

		95.7 ± 10.1	0.97 ± 0.09	93.5 ± 19.6	97.1 ± 11.2	94.9 ± 17.7	96.8 ± 10.3
		(95.1-96.3)	(0.96-0.97)	(92.3-94.8)	(96.4-97.8)	(93.8-96.0)	(96.1-97.4)
	Support Vector	95.8 ± 10.0	0.97 ± 0.09	92.8 ± 21.3	97.5 ± 9.4	94.6 ± 19.1	96.7 ± 9.9
	Machines	(95.1-96.4)	(0.96-0.97)	(91.4-94.1)	(97.0-98.1)	(93.4-95.8)	(96.1-97.3)
	Naive Bayes	95.9 ± 9.8 (95.3-96.5)	0.97 ± 0.08 (0.96-0.97)	93.4 ± 19.8 (92.1-94.6)	97.6 ± 10.0 (97.0-98.2)	95.3 ± 17.5 (94.2-96.4)	96.7 ± 10.4 (96.0-97.3)
	Random Forest	91.8 ± 11.9 (91.0-92.5)	0.92 ± 0.13 (0.91-0.92)	83.4 ± 28.0 (81.6-85.2)	97.0 ± 10.5 (96.3-97.7)	91.5 ± 24.3 (89.9-93.0)	92.4 ± 12.6 (91.6-93.2)

Table 2 .

 2 Diagnostic performances of the various classifiers. Qualitative (i.e., visual) and quantitative scores are given; the latter were derived via evaluation of divergence δ, curl ρ, amplitude γ, and centripetal indicator τ (after 10-fold cross-validation) by various machine-learning algorithms. AUROC: area under the ROC curve; PPV: positive predictive value; NPV: negative predictive value. Accuracies, sensitivities, specificities, PPVs, and NPVs are shown in percentages.

APPENDIX 1: Estimation of apparent microbubble displacement during CEUS This appendix deals with numerical implementation of the algorithm estimating microbubble displacement during CEUS. For each CEUS clip, the lesion was first manually delineated on a hyper-enhanced image. A binary mask was then constructed (this is termed Γ below). We then proceeded as follows:

Estimation of CEUS dense transport fields

We used the transport equation to estimate the apparent microbubble transport field ( ), as suggested in ( 15): [START_REF] Cherqui | Management of focal nodular hyperplasia and hepatocellular adenoma in young women: a series of 41 patients with clinical, radiological, and pathological correlations[END_REF] where I denotes the grey level intensity on CEUS images and I t the partial temporal derivative of I.

Practically, the desired transport field V was estimated between two points in time (t and t + δt, respectively) on a pixel-by-pixel basis using the so-called "optical flow" process [START_REF] Horn | Determining optical flow[END_REF]. The algorithm yields the displacement between two images when the following function is optimised

where is the image-coordinate domain, the estimated pixel-wise transport vector components, and the spatial location. All slices were re-sampled via bi-cubic interpolation to obtain a common isotropic, in-plane, 0.25 × 0.25-mm 2 pixel representation. A spatial low-pass filter was then applied (29) (the cut-off frequency was the proportional pixel fraction of the original image, divided by 16, as suggested by ( 15)) to mitigate the impact of US speckles on the transport equation. An in-house developed, freely available, software provided 2D transport fields using the optical flow metric of Eq. (2) (http://bsenneville.free.fr/RealTITracker/). The reader is referred to [START_REF] De Senneville | Development of a fluid dynamic model for quantitative contrast-enhanced ultrasound imaging[END_REF] for additional details on the numerical implementation of Eq. ( 2).

Note that possible periodic, spontaneous, and drift displacements of tissue must be initially compensated for (24)(25), because they may change image intensities over time; Eq. 1 would (erroneously) attribute such changes to microbubble transport. As proposed in [START_REF] De Senneville | Development of a fluid dynamic model for quantitative contrast-enhanced ultrasound imaging[END_REF], B-mode [START_REF] Dietrich | An EFSUMB introduction into dynamic contrast-enhanced ultrasound (DCE-US) for quantification of tumour perfusion[END_REF] with:

(13)

Determination of the temporal window

The temporal window of analysis was of duration ΔT and commenced at the bolus arrival time t 0 ;

this was determined individually for each patient. To this end, we used a published time intensity curve (TIC) widely employed to determine time constants [START_REF] Tranquart | Perfusion quantification in contrast-enhanced ultrasound (CEUS)-ready for research projects and routine clinical use[END_REF]. The average US image intensity over Γ (termed ) was analysed as a function of time using a two-compartment model, as

where I ∞ is the asymptotic US signal enhancement, and k the uptake rate. I ∞ , t 0 , and k were computed using the Levenberg-Marquardt least-square fit [START_REF] Marquardt | An algorithm for least squares estimation on nonlinear parameters[END_REF] employing all images of the US sequence. The use of a simple two-compartment model was motivated by the fact that only the rise step was screened. The goodness-of-fit was considered acceptable when the Pearson correlation coefficient (r 2 ) was > 0.95. In such cases, the t 0 values chosen earlier served as the start times for the temporal windows.

Hardware and implementation

Our test platform was an Intel 2.5 GHz i7 workstation (eight cores) with 32 GB of RAM. The implementation was performed in C++ and parallelised via multi-threading.

The English in this document has been checked by at least two professional editors, both native speakers of English. For a certificate, please see: