Chi-Dung Phung

Felix Silva

Michele Nogueira

Stefano Secci

MPTCP Robustness Against Large-Scale Man-in-the-Middle Attacks

Keywords: MPTCP, Man-in-the-Middle attacks, communication robustness

Multipath communications at the Internet scale have been a myth for a long time, with no actual protocol being deployed at large scale. Recently, the Multipath Transmission Control Protocol (MPTCP) extension was standardized and is undergoing rapid adoption in many different use-cases, from mobile to fixed access networks, from data-centers to core networks. Among its major benefits -i.e., reliability thanks to backup path rerouting, throughput increase thanks to link aggregation, and confidentiality being more difficult to intercept a full connection -the latter has attracted lower attention. How effective would be to use MPTCP, or an equivalent multipath transport layer protocol, to exploit multiple Internet-scale paths and decrease the probability of Man-in-the-Middle (MITM) attacks is a question which we try to answer. By analyzing the Autonomous System (AS) level graph, we identify which countries and regions show a higher level of robustness against MITM AS-level attacks, for example due to core cable tapping or route hijacking practices.

Introduction

The Multipath Transmission Control Protocol (MPTCP) [START_REF] Ford | Tcp extensions for multipath operation with multiple addresses[END_REF] is an extension of TCP to concurrently use multiple network paths for a given connection. Among many proposals to support these features at the transport layer, it is considered as the one having attracted the largest interest and deployment [START_REF] Bonaventure | Use cases and operational experience with multipath tcp[END_REF]. One of the main reasons for this success is the incremental deployability adopted in its design, with the required signaling transparently reusing existing features of the TCP options.

MPTCP employs multiple 'subflows' to route traffic from a source to a destination in an IP network via different network interfaces and/or TCP ports at the transmitting and/or receiving endpoints. Subflow IP traffic can then be routed independently in the network segment. However, besides the usage of multiple network interfaces at the source or destination, the presence of flow-level load-balancers sensible to port numbers, or multipath proxies aware of the network topology [START_REF] Behaghel | Extensions for networkassisted mptcp deployment models[END_REF] can differentiate the route followed by the subflow packets.

MPTCP is being adopted by major operating systems; it is already happening for Apple OSX and IOS, where it is used for some applications. Its integration in the mainstream Linux kernel is expected for the upcoming versions [START_REF] Raiciu | How hard can it be? designing and implementing a deployable multipath tcp[END_REF]. Among the motivations pushed forward in support of MPTCP, there are [START_REF] Li | Multipath transmission for the internet: A survey[END_REF]: (i) bandwidth aggregation, i.e., the increased network bandwidth offered to a connection; (ii) connection reliability, i.e., the possibility to use an alternative path in case of failure along the primary path or at the primary network interface level; (iii) communication confidentiality, i.e., the decreased ability for a Man-in-the-Middle (MITM) attacker to intercept all the traffic of a same connection.

While the first two aspects above have been largely explored in the last decade, the latter was marginally studied to date. In this paper, we report the results of an extensive measurement campaign aimed at assessing the degree of confidentiality one can expect using MPTCP. In particular, we focus on confidentiality from large-scale, i.e., Autonomous System (AS) level, MITM interception, i.e., looking at the empirical probability that a single connection can be intercepted by an organization or an attacker able to capture all the traffic going through an AS on a given direction (most of Internet communications being asymmetric). Such attacks can happen either by remote access to routing devices of an AS or even by Border Gateway Protocol (BGP) route hijacking. In our analysis, we focus on the case of MPTCP-capable source devices using two edge providers, analyzing measurement results on a geographical basis to identify which countries and regions MPTCP may grant higher confidentiality with respect to large-scale MITM threats.

An important assumption of our analysis is that the MPTCP scheduler behavior of endpoints or multipath converters can be tuned so that it does not only look for throughput maximization, but also for path diversity exploitation for increased confidentiality, as investigated in [START_REF] Coudron | Enhancing buffer dimensioning for multipath TCP[END_REF]. Solutions offering programmability of the MPTCP scheduler are making surface, as notably [START_REF] Frömmgen | A programming model for application-defined multipath TCP scheduling[END_REF][START_REF] Corbillon | Cross-layer scheduler for video streaming over MPTCP[END_REF].

It is worth noting that, despite we refer to MPTCP as our reference multipath transport-layer protocol, our study can apply as well to other functionally equivalent protocols, such as for instance multipath QUIC (Quick User Datagram Protocol Internet Connections) [START_REF] Coninck | Multipath QUIC: design and evaluation[END_REF].

The paper is organized as follows. Section 2 gives a background on MPTCP and related security concerns. In Section 3, we describe our measurement methodology. Section 4 presents the results, different application scopes of this work are discussed in Section 5, and in Section 6 we conclude the paper.

Background

In this section we provide the necessary background on the MultiPath TCP (MPTCP) protocol and on Internet-scale Man-In-The-Middle (MITM) attacks.

MultiPath TCP (MPTCP)

MPTCP extends TCP and allows fragmenting a data flow from a single connection into multiple paths (subflows TCP) [START_REF] Ford | Tcp extensions for multipath operation with multiple addresses[END_REF][START_REF] Peng | Multipath TCP algorithms: theory and design[END_REF], as illustrated in Figure 1. At the application layer, a connection appears as a normal TCP connection. At the network layer, each subflow looks like a regular TCP flow whose segments carry in their header a new type of TCP option [START_REF] Ford | Tcp extensions for multipath operation with multiple addresses[END_REF]. The protocol improves the performance offered by a single flow and makes the connection more reliable using concurrent and redundant paths. The initial TCP connection handshake carries an option, the MP_CAPABLE option, to enable MPTCP capability discovery and subflow creation. The handshake can carry additional information, such as a cryptographic key employed to authenticate the end-hosts and set up new subflows [START_REF] Ford | Tcp extensions for multipath operation with multiple addresses[END_REF]. The establishment of additional subflow may employ also a token and random numbers (nonces), to prevent replay attacks on the authentication method.

Further, an additional address identifier may be employed to identify the source IP address of a packet. Hence, even if the IP header has been changed by a middlebox (e.g. NATs, firewalls), end-hosts can identify an address without any doubt or ambiguity.

MPTCP can overcome some weaknesses inherent to TCP, achieving (i) a greater throughput, (ii) higher reliability, and (iii) higher confidentiality.

Indeed, a multipath connection can improve the throughput aggregating bandwidth over different paths by concurrent data transmission across all available paths. Moreover, a multipath connection can quickly overcome one path failure by sending data to another available path, increasing the data delivery reliability [START_REF] Raiciu | Improving datacenter performance and robustness with multipath TCP[END_REF]. Finally, fragmenting data flow across different subflows makes complete connection interception difficult because attackers would need to capture the transmitted content through all the subflows to build the content.

Therefore, MPTCP can provide a greater level of confidentiality than a regular TCP transmission if the subflows of a connection are routed along disjoint paths: the higher the level of disjointedness, the higher the confidentiality guarantee, and furthermore the higher the level of robustness against such attacks. The goal of this paper is to precisely quantify the level of robustness in use-cases where MPTCP is adopted not (only) to improve communication performance or reliability, but (also) to improve confidentiality. When addressing this aspect, router-level path disjointedness can be considered as being too weak in particular against AS-level traffic capturing and route hijacking. This is the reason why we focus instead on a larger scale of path disjointness, i.e., AS-level path disjointedness, which do make sense in practical scenarios as elaborated here after. Running an analysis on an even larger scale than AS-level scale (e.g., regional or country level) would likely be either infeasible or not sufficiently realistic.

Internet MITM Attacks

In Internet-scale communications, MITM attacks can happen when the attacker gains access to all the traffic transiting through an AS, or at least a portion of it that is enough to reconstruct the transmitted data. In practice, it can be possible by optical layer or BGP route hijacking MITM attacks.

At the optical layer, an attacker is able to split cables by using fiber optical taps, as described in [START_REF] Witcher | Extensions for Network-Assisted MPTCP Deployment Models[END_REF], with a low probability of being detected if peculiar strategies are adopted as explained in [START_REF] Furdek | Physical-layer attacks in transparent optical networks, optical communications systems[END_REF][START_REF] Shaneman | Optical network security: Technical analysis of fiber tapping mechanisms and methods for detection and prevention[END_REF]. Moreover, one can intercept the traffic by exploiting coupling and out-of-the-fiber light propagation phenomena [START_REF] White | An analysis of coupling attacks in highspeed fiber optic networks[END_REF], despite the fact that this is particularly challenging when performing wavelength-division-multiplexing.

At the BGP layer, MITM attacks exploit the natural way BGP works, stealthily hijacking Internet routes to modify or capture the traffic before it reaches the destination. BGP-based MITM attacks have been quite deeply studied for about twenty years; in a recent survey [START_REF] Conti | A survey of man in the middle attacks[END_REF] we have a detailed description of such attacks, their effects as well as mitigation and defense strategies.

This type of attack gained special attention in 2008, when a major provider in central Asia hijacked Youtube traffic to apply local policies. In the same year, a practical BGP MITM attack was demonstrated during the DefCon hacking conference [START_REF] Pilosov | Stealing the internet -an internet-scale man in the middle attack[END_REF]: authors successfully intercepted traffic bound for the conference network and redirected it to a system they controlled before routing it back to DefCon. A recent notable attack happened in 2014, attackers injected BGP routes to redirect traffic from Bitcoin miner nodes to a compromised host [START_REF] Greenberg | Hacker redirects traffic from 19 internet providers to steal bitcoins[END_REF]; it was estimated that at least $83,000 worth of Bitcoins, Dogecoins, HoboNickels, and Worldcoins were stolen over a period of four months. More recently, in 2017 all traffic heading to Visa, MasterCard and other service providers was hijacked for a short period of few minutes [START_REF] Toonk | Bgpstream and the curious case of as12389[END_REF].

The cost of such BGP incidents could be even more than what have been reported. Notable ones are documented in [START_REF] Goldberg | Why is it taking so long to secure internet routing?[END_REF][START_REF] Sankar | B-secure: A dynamic reputation system for identifying anomalous BGP paths[END_REF]; often they are not reported because they cannot be always detectable, they have limited scope, last for a short time etc.

At the transport layer, the advent of MPTCP raised new security specification questions and challenges [START_REF] Bagnulo | Threat analysis for tcp extensions for multipath operation with multiple addresses[END_REF][START_REF] Bagnulo | Analysis of residual threats and possible fixes for multipath tcp (mptcp)[END_REF]. In [START_REF] Kim | Efficient design for secure multipath tcp against eavesdropper in initial handshake[END_REF], cryptography based solutions are proposed against eavesdropping. The authors in [START_REF] Bagnulo | Analysis of residual threats and possible fixes for multipath tcp (mptcp)[END_REF] present an analysis of residual threats in the MPTCP signaling and propose some fixes. Recently, an extension of MPTCP to secure multipath communications was proposed in [START_REF] Jadin | Securing multipath TCP: design & implementation[END_REF], to offer authentication and encryption mechanisms not only to the connection but also to single TCP options. This prevents different types of MITM attacks where an attacker could force all the traffic to be sent only over the path under his control by hijacking the traffic and erasing the In general, most of the works at the state of the art aim at either investigating security threats for MPTCP or proposing solutions for them. It is worth mentioning the rising interest in using MPTCP to further enhance confidentiality when using Internet over-the-top Virtual Private Networks (VPN) services such as ToR and OnionCat [START_REF] Hackdopi | Hackdopi[END_REF]: MPTCP is used in the upstream direction from the client to many gateways accessible via the VPN, on the way to the server, thus increasing the confidentiality level of the connection. Nevertheless, such practices can have a gain which can be hard to assess: how can you ensure the upstream source-destination traffic does follow disjoint paths, hence decreasing MITM efficiency, if not at the routerlevel, at the AS level? In this paper, for the first time at the state of the art to the best of our knowledge, and going beyond the preliminary study presented in [START_REF] Nguyen | Can MPTCP secure internet communications from man-in-the-middle attacks?[END_REF], we attempt to provide a response to such questions.

Methodology

In this section, we first give a description on the datasets used for constructing a representative AS-level graph of Internet, the basis for our analysis. Then, we describe our approach for computing the number of valid vertex-disjoint paths between two arbitrary nodes over the constructed graph.

Finally, we detail how we evaluate path diversity at different geographical scopes. The datasets we employed as well as our scripts are given in [START_REF]LIP6-MPTCP[END_REF] for the sake of reproducibility. Figure 2 gives a schematic illustration of the different blocks of our measurement methodology.

It is worth noting that our methodologies imply that there is a way for a single MPTCP connection to have access to the network path diversity, by means of ad-hoc signaling or specific APIs. Solutions exist in this direction, as described in [START_REF] Kabbani | Flowbender: Flow-level adaptive routing for improved latency and throughput in datacenter networks[END_REF][START_REF] Coudron | Boosting cloud communications through a crosslayer multipath protocol architecture[END_REF].

Graph construction

We extract 2015 AS-level BGP-derived routing data from [START_REF]AS-level Topology Archive[END_REF], and couple it with the inter-AS relationship data [START_REF]AS-level Topology Archive[END_REF] (i.e., indicating which AS is provider, client or peer in an inter-AS link) 2 . The result is a new dataset containing all the AS links along with their frequency of occurrence and relationship type.

We choose this approach because, comparing with other resources [START_REF] Caida | Archipelago (ark) measurement infrastructure[END_REF] [32],

the topological data from [START_REF]AS-level Topology Archive[END_REF] revealed to be more reliable and able to capture a broadened view of the Internet topology. Indeed, it integrates data not only from Routeviews [START_REF] Routeviews | University of oregon route view projects[END_REF], but also from other resources such as RIPE RIS [START_REF]Ripe routing information service[END_REF]. It is worth noting that the alternative traceroute-based approach employed in [START_REF] Caida | Archipelago (ark) measurement infrastructure[END_REF] has known issues [START_REF] Oliveira | The (in)completeness of the observed internet as-level structure[END_REF] when converting router-level paths into AS-level.

Employing measurements over a long period allows us to capture interdomain connection dynamics as well as inter-AS economic relationships. For instance, in a one month period, only 85% of inter-AS links appear more than 20 days, the remaining links with lower frequency of occurrence being those used for backup operations or during BGP convergence periods. For the sake of consistency, we removed these unstable links.

Path diversity computation

In order to have a measure of the path diversity, we need to enumerate all the paths connecting two nodes over a graph that satisfy given routing properties. This problem is often referred to as policy-compliant path diversity computation in the literature [START_REF] Erlebach | Connectivity measures for internet topologies on the level of autonomous systems[END_REF][START_REF] Klöti | Policycompliant path diversity and bisection bandwidth[END_REF]. The common approach [START_REF] Erlebach | Connectivity measures for internet topologies on the level of autonomous systems[END_REF] to this problem is to convert the original graph into a type-of-relationship (ToR) graph [START_REF] Di Battista | Computing the types of the relationships between autonomous systems[END_REF], i.e., a directed graph in which (i) the relationship between two adjacent vertexes is expressed via the direction of the edge connecting them, then (ii) maximizing the total number of vertex-disjoint paths between nodes in this graph. However, the time-complexity experienced in such methods is relatively high hence intractable for a graph as big as the AS graph.

In order to better scale, we introduce a novel path search algorithm leveraging the scale-free characteristics [START_REF] Albert | Statistical mechanics of complex networks[END_REF] of the input AS graph (i.e., a graph with relatively few hubs capturing the majority of the paths) to optimize the execution time. In such a scale-free graph, the diameter (i.e., the length of the longest path among all the shortest paths) is not too high. Thus, the average path length (measured in number of AS hops) connecting any pair of nodes in the AS-level graph of Internet is around 5 as of today [START_REF] Kühne | Update on as path lengths over time[END_REF] (note that it is a bit lower with IPv6).

Searching for paths in a scale-free graph, i.e., a graph with a large minority of hub nodes connecting the rest of the nodes, is a problem of controllable complexity when adopting breadth or depth-first search algorithms with a limited depth; indeed, fixing a limited depth to a graph search, and that for a scale-free graph that has a limited diameter, strongly decreases the number of explored branches in the graph exploration 3 . From the constructed AS graph G, the breadth-first search algorithm we describe in Alg. 1 can be applied to discover all the policy-compliant paths between two nodes s and d, in a reasonable time.

Alg. 1 works as follows: (i) starting from the origin s, the algorithm explores every adjacent node n of s. (ii) A queue P is introduced to keep track of the explored paths; initially, it includes all the paths from s to n.

(iii) Following these paths, the algorithm continues discovering the adjacent nodes to look for destination d. (iv) For a path p dequeued from P , the last node n is extracted, all of its neighbors are checked in sequence to determine the valid next hops towards d. (v) Once a neighbor is determined as valid, link to that neighbor will be added into the current path forming a new valid path toward destination. This new explored path is then enqueued into P for the next discovering phase. (vi) A node is considered as valid once the path through it does not violate the valley-free routing property [START_REF] Gao | On inferring autonomous system relationships in the internet[END_REF] 4 ; we express such policy-compliant path (i.e., a path that complies with the valley-free routing policy), using the following regular expression c2p * p2p?p2c * [START_REF] Klöti | Policycompliant path diversity and bisection bandwidth[END_REF] in which c2p, p2p and p2c denote the relationship between interconnected nodes (where ? means that you can have one or none p2p link).

3 breadth-first search explore first all the neighbors of a node, and then explore deeper in the graph; depth-first search, instead, explore first in depth starting from a given neigbor, and proceeds to the next neighbour only when the exploration in depth from the first one has terminated. 4 A valley-free path is defined as a path that does not cross more than one peering agreements, which are agreements over which two ASes exchange only routes towards respective customers, which is justified by the fact that peering agreement are meant to be free-of-charge for both ASes It is worth noting that, within G links are labeled according to their inferred relationship. For example, assuming that n 1 , n 2 , n 3 are neighbors of node s, in which s is customer ('c') of n 1 , provider ('p') of n 2 and peer with n 3 ; the links (s, n 1 ), (s, n 2 ), and (s, n 3 ) are labeled as 'c2p', 'p2c' and 'p2p', respectively. With these labels, the regular expression for policy-compliant path then could be leveraged to determine the validity of next hop toward the destination. For instance, taking the customer-type neighbors among the neighbors of s (i.e., n 2 ), and looking at their neighbors x in turn, those (n 2 , x) links are not validated if they are either c2p or p2p because a customer is not expected to grant transit towards its other provider(s) to one among its providers, and a customer is not expected to give access to its peer(s) to its provider(s). By checking the labels of links along the explored path, the validity of next hops can be determined. Once a valid path is discovered, it is enqueued into P for the next discovering phase. The same exploration and validation processes are repeated for all the paths in P until reaching destination d or the path length goes over a given threshold τ .

The path validation executes at run-time to ensure that non-compliant paths are detected at the early stage, thus avoiding wasting time exploring invalid paths. By reducing the number of paths needed to be explored in the following phases, the search space is continuously optimized. Moreover, a proper choice of τ not only limits the time and space complexity, but can also avoid selecting long paths to be avoided in practice.

As a result of the path search algorithm, policy-compliant paths between two endpoints may share common nodes. To get the final set of vertexdisjoint paths, we run a simple off-line filtering linear algorithm to capture the shortest disjoint paths. Since the original list of valid paths turned out to be quite small most of the time and already sorted, the complexity of such a filtering operation is negligible.

Source-destination pairs

Within the constructed AS-level graph, multipath connections could be simulated by simply attaching end-hosts as virtual nodes into AS nodes of the original graph. For instance, a multi-homed device can be emulated by adding a new node, then linking it with at least two AS nodes. The connection from that node to any other virtual nodes forms a multipath transportlayer communication. Our approach for emulating multipath communication can therefore be simply referred to as a process of source-destination pair selection. In the following, we define the target set of AS nodes which we

if v / ∈ V isitedN odes then for n ∈ v.N eighborSet do if n / ∈ V isitedN odes and (label(v,n)='p2c' or label(v,n)='p2p') then for x ∈ n.N eighborSet do if label(n,x)='c2p' or label(n,x)='p2p' then g.RemoveEdge(n,x) end end end N ewP ath ←-list(path) N ewP ath.append(n) if n = d then
ValidPathSet.append(NewPath) end if length(N ewP ath) = τ + 1 then break queue.append(NewPath) end VisitedNode.add(v) end end consider for attaching the end hosts. A simulation process is then described in details explaining which communication scenarios are covered in our study.

The current Internet ecosystem is composed of more than 70 thousand ASes, out of which the large majority are stub ASes, i.e., ASes that are only origin or destination ASes. About 13% are Tier-3 or small Tier-2 ASes, we arbitrary define in this paper as those appearing at most in the third from last position and at least penultimate position in BGP AS paths; we refer to such ASes as 'edge provider' ASes, which can be considered as a representative set of national Internet Service Provider (ISPs). Such ASes are often referred to as 'eyeball' ASes. In this paper, an edge provider AS is not a stub AS, but is rather expected to be a regional or national ISP, most of the time (rare are the cases where an international/intercontinental ISP gives Internet access to end-users).

Rather than taking into account all possible communications, we target the connections among hosts at the edges, performing connections using multiple sub-flows such as done with MPTCP. Considering connections between hosts in different countries, we precisely address the MITM robustness of Internet connections crossing multiple ASes. To precisely determine which communications to cover in our study, we define a target set of sourcedestination pairs that address, in a reasonable yet arbitrary way, the communications that may be more sensitive to communication privacy. Our choice of source-destination pairs is as follows:

• the source is interconnected to two edge providers in a country.

• the destination is not multi-homed, i.e., it is reachable via a single ISP, the one given by the best BGP path from each source edge provider, and belongs to an AS at another country than the one of the source. Figure 3 illustrates an example of how we simulate multipath communications accordingly the above policy. For each two arbitrary edge provider ASes in a same country, one source is created (i.e., a dual-homed source).

For each edge provider in another country, one destination is paired with the source. Such a pair dual-homed source -single-homed destination defines the two endpoints of a multipath communication. Listing all pairs, i.e., combining a given source with every destination, all possible (international) communications of a dual-homed host can be covered.

Besides reducing the number of pairs to a reasonable and treatable number (requiring about one week of computation), it is worth noting that, in such a way, we consider communication in a single direction: from source to destination. That is, under such a path election strategy, we cover the case when a multi-homed device uploads to a single-homed server, as well as the case when a single-homed device downloads contents from multi-homed servers.

The scenarios that are not covered in our study include: (i) multi-homed devices downloading from single-homed server; (ii) single-homed devices uploading contents to multi-homed servers; (iii) a multi-homed device communicating with another multi-homed device. A dual analysis, quite expensive computationally, covering these additional cases may be performed as well in future works.

MiTM robustness metric aggregations

The ability to split traffic over different paths allows multipath protocol to realize MiTM attacks more difficult. Thus, the chance for an attacker to capture all the traffic sent by a source is reduced in proportion to the number of disjoint paths between the source and the destination. Path diversity is therefore a proper indicator to evaluate the MiTM robustness of a multipath communication.

Rather than considering the robustness against MiTM attacks of every connection individually, we are more interested in evaluating such a robustness at the end-host level, thus measuring the degree of robustness offered by a multipath-capable source device to secure its data sending over the Internet.

Source-specific MiTM robustness metric

With regard to the aforementioned approach for source-destination pair selection, we define the source-specific MiTM robustness metric as the aver-age number of disjoint paths over all the destination edge providers that are in a different country than the source. Such a metric can be considered as a level of unlikelihood that a MiTM attack takes place for that source configuration; the higher the value of the robustness metric, the more difficult it is for an attacker to capture traffic from that source.

Source country-specific MiTM robustness metric

Aggregating results from all the sources within a given country we can obtain a source country-specific MiTM robustness metric. Such a definition allows us to characterize the robustness level offered by different source countries to multipath communications.

Country-level source-destination MiTM robustness metric

As another way to aggregate the MiTM robustness metric computation, we also study a country-level source-destination based aggregation, i.e., leading to a robustness metric for a pair of source and destination countries.

Given a source (a pair of edge providers in a country) and a destination country, its MiTM robustness metric is defined as the average number of disjoint paths from the source over all edge providers belonging to the destination country.

Country-pair MiTM robustness metric

By grouping together the results from all the sources within a source country, we can define the country-pair MiTM robustness metric for the corresponding pair of countries.

Metric computation

Let us more precisely characterize the aforementioned source-destination pair selection process with respect to the two MiTM robustness metric aggregations we study in the following, i.e., the source country-level one and the country-pair one. We segment the set of edge providers, E, in countryspecific subsets, E c , where c denotes a country in the set of countries C, i.e., E = c∈C E c . We employ the AS-to-country mapping given by the CIDR Report [START_REF] Bates | Cidr report[END_REF]. Let us indicate with E c the restriction to a specific country c ∈ C. Overall, for a given country c, the number of source-destination pairs is therefore equal to the number of pairs of edge providers for the given country multiplied by the number of edge providers of other countries, i.e.

|E c| × (|E c| -1) 2 × c =c |E c | (1)
For a given source and destination countries, s and d respectively, the number of source-destination pairs connecting them is equal to:

|E s | × (|E s | -1) 2 × |E d | (2)
Doing so, we target a lower bound, pessimistic analysis, since we only take into consideration international communications and we suppose the destination is not multi-homed. The filter we set on the destination enumeration allows us to target communications that may need a higher level of confidentiality due to their international connotation. Moreover, in this way we also avoid a huge bias potentially due to the fact that a large majority of the AS paths available at the national level are not visible in backbone BGP routing tables such as the Routeviews ones (typically because of Internet exchange points, as recently shown in [START_REF] Ager | Anatomy of a large european IXP[END_REF]). We believe having a lower bound stand is more appropriate than an upper bound one, while allowing us to scientifically qualify the value of the relative trends.

Results

We report the results obtained for a set of 147 countries, i.e., those countries from the United Nations statistics [START_REF]Standard country or area codes for statistical use (m49)[END_REF] that appear to have at least two distinct edge providers officially based in the country; this automatically excludes Greenland territories, very small city-state countries, many African countries and Indonesia. The geographical coverage is given in Figure 5. In the following sections, we present the statistics for two different MiTM robustness metric aggregations, the country source specific one and the country pair one.

Source country aggregation

Let us recall the measurement approach for source country-specific MiTM robustness analysis:

• For each country, we generate all possible dual-homed sources, i.e., all possible pairs of edge providers.

• For each such source configuration, we compute the number of disjoint paths to each destination. For each edge provider that is a different country than the source country, one destination is generated.

• For a given source, we compute its corresponding robustness metric by taking the average of the number of disjoint paths over all the destinations.

• For each country, a series of MITM robustness metrics is hence generated, one for each source.

We characterize the resulting series using boxplot distributions (using a 0.1% outliers threshold). We overlay over the boxplots the average of the corresponding series with a red square, order them with increasing averages5 from left to right. We report the results in Figure 4, and with a geographical view in Figure 5. We express three different viewpoints:

• device view (Figure 4a): the MITM robustness is computed with the source node integrated in the AS graph as an 'artificial' node, i.e., the path search algorithm finds the number of AS-disjoint paths from this source node toward the destination. It provides therefore a device view; obviously, in this view the upper bound of the robustness is 2, i.e., the number of edge providers used by the source.

• edge provider view (Figure 4b): the MITM robustness is computed counting the number of disjoint paths from the first and the second edge provider, then decreased by those paths that share an AS hop.

Taking into account such a view, we assume that additional AS paths can be made available to MPTCP subflows acting at the edge providers level, e.g., by forms of flow path steering and load-balancing.

• differential view (Figure 4c): the differential robustness results, i.e., the edge provider view robustness minus the device view robustness, computed for each source configuration individually. This view more precisely quantifies the gain achievable for MPTCP communications when inter-AS load-balancing is enabled at the edge providers. The above viewpoints also reflect different levels of trust on the providers.

That is, while the edge provider view assumes MITM attacks do not happen at the source and destination edge providers (i.e., there is a high level of trust on those providers), the device view assumes that attacks can happen at the source edge providers, hence revealing a low level of trust in source direct providers. As a general assessment, Figure 4 shows a distribution to be interpreted.

For example, one could consider 1.5 as the rough threshold above which the likelihood of MiTM is to be considered low, and conversely high if lower than 1.5. Only about 5% of the countries show good chances of being robust against MITM from a device viewpoint, while looking at the maximum instead of the average and median values one could speculate that careful choice of the edge providers could make the MiTM likelihood low for a majority of the countries. From an edge provider viewpoint, this ratio grows to roughly 60%, and higher than 90% looking at the maximum, that is if the edge provider choice can be influenced by confidentiality concerns.

Moreover, the average number of paths connecting a dual-homed node to international destinations has a significant variance depending on the origin country. The average robustness ranges from 1 (and less) to 1.6 from a device viewpoint, and from 1 (and less) to 2.5 from an edge provider viewpoint.

It is worth noting that the reason why some minimum, and even average values, are below 1, is the partial view over the Internet topology and the incompleteness of inter-AS relationship inference; in fact, these factors make some destinations unreachable (counted as 0 path), but we left the 0 values in the series to also give an index of the level of topology incompleteness for different countries. In any case, the boxplot median is a metric robust against such outliers to look at.

In addition, observing the distributions in Figure 4, we can also remark that:

• Within a country, a high inter-quantile range indicates that the path diversity strongly depends on how the two upstream edge providers are selected for the source.

• The gap between the min and max robustness is another interesting fitness metric to observe. Some countries maintain a small gap (below 1) while others have a very big gap (up to 2). In other words, the deployment of multipath transport-layer communications for securing international communications in some countries can statistically yield a much better result than in other countries, where this gap is smaller.

Particularly interesting is the case of Angola (AO), Venezuela (VE) and Namibia (NA), with small robustness gaps, which may be correlated to the presence of inter-continental cables landing in or close to the country [START_REF]Greg's cable map[END_REF].

• The median is mostly higher than the average in the device view, and lower than the average in the edge provider view. This is essentially due to outliers, counted in the average and not in the median.

• From the edge provider viewpoint, the maximum value is higher than 2 in the most of the countries, suggesting that with a proper choice of trusted source providers, one can adopt multipath communications to statistically expect high confidentiality for its communications. Particularly alerting are the cases of Uzbekistan (UZ), Nepal (NP) and Lebanon (LB), with quite low maximum values.

• From the device viewpoint, in most of the cases the maximum robustness is not higher than 1.6, both averages and medians are quite far from the desirable target of 2. Hence, without the support of inter-AS load-balancing at source providers, path diversity from a dual-homed node is reduced significantly, indicating a non negligible probability of paths joining on the way to the destination.

• Considering the differential robustness, we can remark that among the countries that have the lowest device view MITM robustness, those that could most benefit from inter-AS load-balancing practices are Mongolia (MN), Pakistan (PK) and Korea (KR). However, the majority of those countries with low robustness do not improve much the situation going from the device view to the edge provider view.

Looking at macro geographical regions, many European countries seem to grant better security than countries in other regions. In order to look at continental characteristics, the plots in Figure 6 show the boxplot results (with 1% outliers) aggregated on a macro-region basis (a and c, sub-continental level) and on a relative position basis (b and d, in terms of seacoast and inland borders). We can remark that:

• Western Europe appears to be the best off, followed by Northern Europe and Northern America. In almost 50% of Western Europe countries there can be 2 disjoint paths from the source edge providers to Internet destinations.

• Central Asia shows the worst robustness, followed by Australia and New Zealand; the reasons are likely network centralization practices and geographical isolation. It is interesting to notice the relevant gap between Central and South-Eastern/Western Asia.

• within Europe, Western countries do offer a better diversity over Northern countries, and especially over Eastern and Southern countries. with a small range of variation and a high median value show the best result.

• A high variance is recorded at Southern Asia, Northern Europe and Sub-Saharan Africa, which indicates high differences among the countries within these areas.

• We could not find a strong correlation between the relative continental position, and the robustness metric, yet a positive correlation exists, with countries at the boundaries of oceans, with inter-continental cable landing and that are sea-oriented (most of the border on the coast) that offer higher robustness than fully internal and continental-oriented ones.

Source-destination country pair aggregation

As we may notice, the MiTM robustness level of a multipath communication could be affected not only by the country where the communication starts but also by the choice of upstream providers at that country. Besides that, within a source country, the robustness level for different destination countries can significantly vary. To evaluate this latter aspect further, we perform a source-destination country pair aggregation.

Over the set of 147 countries, we evaluate the robustness metric for 1547 directional country-to-country communication pairs in which the MiTM robustness metric for one pair is computed as follows:

• For a given source country, we generate all possible dual-homed sources, i.e., all possible pairs of edge providers.

• For each such source configuration, we compute the number of disjoint paths to each edge providers located in the destination country.

• For a given source, we take the average of the number of disjoint paths over all the destinations to get its source-destination based MiTM robustness metric.

• For a given source-destination country pair, a series of MiTM robustness metrics, one for each source, is therefore created.

In Figure 7, we report the CDF of the average MiTM robustness, for all the 1547 pairs. The high range of variation (between 0.4 and 6) shows us the big robustness gap between pairs. Only 20% of the country pairs show an average of two or higher. For the remaining pairs, approximately 73% of them have the average range from 1 to 2. The remaining 7% are country pairs with very low robustness, below one; besides the specific context related to a country pair, a factor behind such bad performance can be the already discussed topology view incompleteness.

To better understand the impact caused by different destinations, we further characterize the top 147 and bottom 147 pair in the CDF distribution, i.e., roughly the top 10% and the bottom 10% cases. The results are presented to, e.g., Montenegro and Republic of Congo is extremely low.

Considering 1 and 2 as the thresholds for very low (zero) and high (sufficient) robustness, respectively, a source-destination pair can be classified as:

(1) highly robust against MiTM if it has the average robustness level of at least 2, and (2) weak against the MiTM once maintaining the average of 1 or lower.

We visualize the country-to-country communications in these two classes by mapping them into a geographical map in Figure 9. To avoid too many lines, we first group countries with respect to their subregion, then converting these country-to-country connections into the corresponding subregion-tosubregion connections. Finally, the subregional connections are expressed using lines with different opacity reflecting the portion of country-to-country communications between subregions having the MiTM robustness level less than or equal to 1 as in Figure 9a, and equal to or higher than 2 as in Figure 9b.

In Figure 9a, we only show the connections between subregions when there are more than 30% of the country-to-country communications with a robustness metric of at most one. For subregion pairs with less than 30% of their country-to-country communications having a robustness metric lower than one, the connection lines are hidden. In other words, the lines point out the subregions where the deployment of MPTCP cannot offer any protection against large-scale MiTM attacks. As presented in the map, the area of Central Asia and Melanesia are the two subregions having the worst performance, most of their MPTCP communications with other subregions are classified as zero-robust. Thus, most of the subregions could not be benefit from the deployment of MPTCP to secure their communications with Central Asia.

In the sub-regional view of the high robustness group presented in Figure 9b, we show the connection lines between sub-regions with more than 50% of the country-to-country communications having robustness level of 2 or higher. In such a view, Micronesia and then Western Europe are the two areas that outperform the others in term of MiTM robustness. As depicted in the plot, except for a few low connected regions, like Central Asia, Caribbean and Northern Africa, etc., most of the multipath communications from and to Micronesia can profit from a high level of robustness. It is worth noting that in the region of Micronesia, Guam is the only country covered by our study. The high robustness result captured for communications from and to this region is therefore directly related to the highly connected network infrastructure of Guam being a crucial node in the Internet cable network [START_REF] Starosielski | Critical nodes, cultural networks: Re-mapping guam's cable infrastructure[END_REF].

(a) regions with more than 30% of country-to-country communications having at most one path (b) regions with more than 50% of country-to-country communications having at least two paths

Practical aspects

We focused our study on MPTCP-based communications. More precisely, it covers the following cases:

• MPTCP capable endpoints: both source and destination, client and server (or vice versa), are MPTCP capable, and the MPTCP communication is not filtered by middle-boxes. As argued in Section 3.3, the multi-homed endpoint can be either the server or the client.

• MPTCP proxied endpoints: at least one endpoint is not MPTCP capable, but the TCP communications are handled by MPTCP proxies, converting TCP packets into MPTCP packets and vice versa, as explained in [START_REF] Behaghel | Extensions for networkassisted mptcp deployment models[END_REF][START_REF] Boucadair | Network-assisted mptcp: Use cases, deployment scenarios and operational considerations[END_REF], possibly routed via Internet disjoint paths as proposed in [START_REF] Coudron | Crosslayer cooperation to boost multipath TCP performance in cloud networks[END_REF][START_REF] Benchaïb | Transparent cloud access performance augmentation via an MPTCP-LISP connection proxy[END_REF]. The multipath conversion proxies can sit at endpoint premises (customer premises equipment for the client, hypervisor or middle-box at the server) or at the edge provider level borders.

Besides MPTCP-based communications, other protocols offering Internetscale multipath, connection flow-level load-balancing could also be covered by our study. The following protocols are either not deployed, or they have only undergone a limited deployment at the Internet scale so far; they are:

• SCTP : the Stream-Control-Protocol (SCTP) [START_REF] Tuexen | Stream control transmission protocol (sctp) chunk flags registration[END_REF] is another multipath transport protocol absolving the same function as MPTCP, but less deployed than MPTCP due to the limited retrocompatibility.

• LISP : the Locator/Identifier Separation Protocol (LISP) [START_REF] Fuller | Locator/id separation protocol (lisp) map-server interface[END_REF] is able to perform inter-AS inbound load-balancing by means of encapsulation, routing locator mapping, and appropriate traffic engineering (TE) policy configuration. LISP primary scope is the edge provider one, hence results with the edge provider view are readily applicable. Furthermore, deployment of LISP as an intra-AS TE tool can also allow us to perform inter-AS multipath on the outbound direction as proposed in [START_REF] Misseri | Internet routing diversity for stub networks with a map-and-encap scheme[END_REF].

• MultiPath BGP : in BGP, the routing decision process only allows us to take one route per network prefix. The selected path can be inefficient in terms of global routing. Recently, forms of Multipath BGP were discussed in standardization fora, but finally not standardized; however, some recommendations have been published [START_REF] Lange | Issues in revising bgp-4[END_REF], and implemented by some vendors (see, e.g., [START_REF] Juniper | Configuring bgp to select multiple bgp paths[END_REF] and [START_REF]Bgp best path selection algorithm. Cisco guide, Troubleshooting tech notes[END_REF]). Such multipath mode can be adopted at the edge provider scope to enable load-balancing at the egress direction. Despite the study [START_REF] Elena | Characterisation of as-level path deviations and multipath in internet routing[END_REF] on core routing tables reports that in 2010 multipath BGP was practically not used, speculations report that it is used by major cloud providers.

The above protocols are a selection of those protocol communication contexts where load-balancing can affect the AS-path selection. There are also other load-balancing protocols which can potentially influence the egress AS selection as well, as for instance in data-center environments. In the case of MPTCP communications, these protocols, operated at the edge provider view, are able to perform inter-AS load-balancing in such a way that the path diversity exposed in our edge provider view can be made available to MPTCP devices, hence giving them the full potential of MPTCP in terms of communication confidentiality and robustness against MITM attacks.

Finally, additional multipath transport-layer protocols are making surface, as for example the already mentioned multipath extension to the QUIC protocol [START_REF] Coninck | Multipath QUIC: design and evaluation[END_REF], nicknamed MPQUIC. As MPTCP that authenticates its options to avoid interference from one path to the others as already discussed, MPQUIC also has a similar protection by default, because every control information in QUIC is authenticated.

Conclusion

We explored in this paper how Internet path diversity could be exploited by means of multipath transport-layer protocols such as MPTCP, when looking at increased security against man-in-the-middle attacks. We focused on such attacks acting at the autonomous system level, and at the robustness of multipath communications in what appear as a reasonable configuration where at least one endpoint is multi-homed with two edge providers.

We reported extensive, specific and aggregated results for most of the world countries and regions, looking at macro trends that could inspire further research in the area. Results show that, statistically speaking, a multipath protocols such as MPTCP does not help in guaranteeing robustness against MiTM attacks hence high confidentiality, unless (i) the choice of the edge provider is carefully taken, or (ii) one can rely on inter-AS loadbalancing features offered implicitly or explicitly by edge providers. Some continental regions are strongly more robust than others, and there seems to be a positive correlation with inter-continental cable landing proximity.

Moreover, the results show that there are countries surprisingly less well connected than one could think of, such as Northern America countries, and countries that are more obviously less robust against such attacks due to network centralization practices.

It is worth mentioning that the methodology we propose to measure MiTM robustness could be instrumental also for other types of analysis. For instance, having a high MiTM robustness may also represent an increased sensibility to distributed denial of service attacks (DDOS), as the set of possible sources not sharing a network bottleneck can be expected to increase with the AS-level path disjointness. This aspect may be object of further work.

Application

Figure 1 :

 1 Figure 1: Multipath TCP Connection: Overview

Figure 2 :

 2 Figure 2: Adopted measurement methodology

Figure 3 :

 3 Figure 3: Representation of the source-destination pair selection process.

Figure 4 :

 4 Figure 4: MITM robustness distribution for 147 countries.

Figure 5 :

 5 Figure 5: Countries covered with corresponding MiTM robustness distribution

Figure 6 :

 6 Figure 6: MITM robustness metric with continental subregion grouping.

Figure 7 :

 7 Figure 7: CDF of average MiTM robustness for 1547 pairs of source-destination country

FigureFigure 8 :

 8 Figure8areports the MiTM robustness metric distribution for the top 147 country pairs. The high inter-quartile range (IQR) with a pair highlights the strong impact caused by edge providers choice at the source to the robustness metric. Besides that, there are also some source countries, such as Morocco (MA), Madagascar (MG), Gibraltar (GI), Guam (GU), Jersey (JE), Namibia (NA), Liechtenstein (LI) and Belize (BZ), that suffer from the presence of only one edge provider pair; these countries result in pairs with a collapsed robustness point in the box. In addition, within these top 147 pairs, there are some destinations, like Namibia (NA), Guam (GU) and Belize (BZ), that appear to show high sensibility to the destination choice on the MiTM robustness.In Figure8b, we report the results for the bottom 147 country pairs. The majority of them have Montenegro (ME) as the destination. The second destination is Republic of Congo (CG). That highlights again the impact of destination choice on the MiTM robustness level. Unlike the top 10% case, we see a small inter quartile range (IQR) for most of the pairs, showing that even a careful choice on the edge providers at the source country cannot improve much the level of robustness for such connections. In other words, regardless of the origin country as well as the choice of source edge providers, the possibility of employing MPTCP to secure the communications destined

Figure 9 :

 9 Figure 9: Regional view of the source-destination based MiTM robustness

A preliminary version of the content of this paper was presented in[START_REF] Nguyen | Can MPTCP secure internet communications from man-in-the-middle attacks?[END_REF].

The inter-AS relationship data from[START_REF]AS-level Topology Archive[END_REF] is extracted monthly from the Cyclops database[START_REF] Cyclops | Cyclops[END_REF], which combines BGP data with Internet eXchange Point (IXP) data and adopts inference techniques proposed in[START_REF] Oliveira | The (in)completeness of the observed internet as-level structure[END_REF].

Average values do include outliers.

Acknowledgement

We acknowledge the support and relevant effort of Ho Dac Duy Nguyen in the preparation of the data processing and methodology. This work was partially funded by the French Investissement d'Avenir FED4PMR project.

We thank Xenofontas Dimitropoulos for his valuable feedback.


Karim Makki

Bhushan Borotikar

Marc Garetier

Douraied Sylvain Brochard

Ben Salem

François Rousseau

Ben Douraied

Salem

In vivo Ankle Joint Kinematics from Dynamic Magnetic Resonance Imaging using a Registration-based Framework

Keywords: dynamic MRI, ankle, motion estimation, high-resolution reconstruction, intensity-based registration, locally rigid transformation

In vivo ankle joint kinematics from dynamic magnetic resonance imaging using a registration-based framework

Introduction

To better understand the biomechanics of the pediatric ankle joint, it is crucial to establish in vivo normative joint biomechanics [START_REF] Borotikar | Dynamic MRI to quantify musculoskeletal motion : A systematic review of concurrent validity and reliability, and perspectives for evaluation of musculoskeletal disorders[END_REF] before focusing on pathomechanics studies. This approach would likely contribute to a better long-term follow-up for childhood disabilities such as cerebral palsy.

Conventional MRI techniques have been used for accurate diagnosis of the ankle and foot disorders given the complexity of their anatomy [START_REF] Bae | New techniques in MR imaging of the ankle and foot[END_REF]. Dynamic MRI was originally developed to study cardiovascular functions towards the end of the 1980s [START_REF] Pettigrew | Dynamic cardiac MR imaging. Techniques and applications[END_REF]. Since then, MRI acquisition and reconstruction techniques have improved by leaps and bounds with current sequences such as Fast-PC MRI, ultrafast MRI, and Fast Field Echo FFE. These techniques have been successfully employed to quantify the knee joint kinematics [START_REF] Sheehan | Kinematic determinants of anterior knee pain in cerebral palsy : a case-control study[END_REF][START_REF] Borotikar | A methodology to accurately quantify patellofemoral cartilage contact kinematics by combining 3D image shape registration and cine-PC MRI velocity data[END_REF], and to analyze in vivo 3D musculoskeletal dynamics [START_REF] Behnam | Assessing the accuracy and precision of musculoskeletal motion tracking using cine-PC MRI on a 3.0 T platform[END_REF][START_REF] Sheehan | 3d musculoskeletal kinematics using dynamic mri[END_REF], non-invasively. Sequences such as cine-PC MRI [START_REF] Pelc | Phase contrast cine magnetic resonance imaging[END_REF][START_REF] Sheehan | Using cine phase contrast magnetic resonance imaging to non-invasively study in vivo knee dynamics[END_REF] or fast-PC MRI [START_REF] Sheehan | Quantitative MR measures of three-dimensional patellar kinematics as a research and diagnostic tool[END_REF][START_REF] Rebmann | Precise 3D skeletal kinematics using fast phase contrast magnetic resonance imaging[END_REF] can provide in vivo joint velocity field measurements. However, these data require long acquisition times with a high number of repeated cycles. Sequences such as ultrafast MRI [START_REF] Clarke | A non-invasive, 3D, dynamic MRI method for measuring muscle moment arms in vivo : demonstration in the human ankle joint and achilles tendon[END_REF] allow faster scan speed by exploiting spatio-temporal redundancy and can acquire images during a single cycle. However, they require slow motion of the joint making the movement quasi-static. Thus, both these types of sequences are not appropriate to acquire dynamic in vivo data for children especially with musculoskeletal disorders. In this work, we use real-time FFE sequences [START_REF] Schaeffter | Projection reconstruction balanced fast field echo for interactive realtime cardiac imaging[END_REF] which are less sensitive to motion artifacts and which reduce the scanning time to 18 seconds to acquire the dorsi-plantarflexion cycle.

Furthermore, post-processing of acquired dynamic images remains a challenging task due to low resolution and motion artifacts. In [START_REF] Sheehan | Kinematic determinants of anterior knee pain in cerebral palsy : a case-control study[END_REF], the displacement of specific points on each bone were quantified by integrating the velocity data using a three-dimensional Fourier tracking method [START_REF] Zhu | Three-dimensional motion tracking with volumetric phase contrast MR velocity imaging[END_REF]. Clarke et al. [START_REF] Clarke | A non-invasive, 3D, dynamic MRI method for measuring muscle moment arms in vivo : demonstration in the human ankle joint and achilles tendon[END_REF] determined in vivo muscle moment arms during ankle movement using a 3D mesh-based registration method. They estimated bone motions by registering the high-resolution joint geometry from MRI scans of the stationary joint to low-resolution geometries from ultrafast MRI scans of the slowly moving joint. However, in addition to the segmentation of the high-resolution static scan, it was needed to manually segment each low-resolution time-frame which was time-consuming and prone to segmentation bias.

In this study, the spatial information of the high-resolution static MRI data and the temporal information of the dynamic data were combined using a log-euclidean polyrigid framework (LEPF). The purpose of this study was twofold : 1) to estimate bone motion, allowing to derive the 3D kinematics of the joint, 2) to compute temporal dense deformation fields of the joint, allowing to reconstruct a high-resolution dynamic MRI sequence from a low-resolution dynamic MRI sequence and one high-resolution MR image.

The use of the LEPF proposed by Arsigny et al. [START_REF] Arsigny | Polyrigid and polyaffine transformations : a novel geometrical tool to deal with non-rigid deformations-application to the registration of histological slices[END_REF][START_REF] Arsigny | A log-euclidean polyaffine framework for locally rigid or affine registration[END_REF][START_REF] Arsigny | A fast and log-euclidean polyaffine framework for locally linear registration[END_REF] has been investigated for registration of structures such as the hip and knee joint in mouse CT images [START_REF] Papademetris | Articulated rigid registration for serial lower-limb mouse imaging[END_REF], human hands from 2D hand radiographs [START_REF] Martín-Fernández | A log-euclidean polyaffine registration for articulated structures in medical images[END_REF], bones in lower-abdomen area [START_REF] Commowick | An efficient locally affine framework for the smooth registration of anatomical structures[END_REF], and also intra-subject mandible 3D datasets [START_REF] Seiler | Capturing the multiscale anatomical shape variability with polyaffine transformation trees[END_REF].

In this work, we adapted this framework to deal with 4-dimensional MRI data by proposing a fast and accurate algorithm to compute the dense deformation fields in regular grids. In fact, we propose to compute the exponential mapping of the LEPF using matrix eigendecomposition instead of the scaling and squaring method used in the literature [START_REF] Arsigny | A fast and log-euclidean polyaffine framework for locally linear registration[END_REF][START_REF] Commowick | An efficient locally affine framework for the smooth registration of anatomical structures[END_REF].

Materials and methods

Subject recruitment

Six typically developing children with a mean age of 12 years and with a mean weight of 35.8 kg participated in this study which was approved by the regional ethics committee.

Children were selected with no contraindications to MRI and with no history of lower limb musculo-skeletal injury or surgery in the past six months. MRI data were acquired in a single visit after parents signed informed consent forms.

Data Acquisition

MRI data have been acquired using a 3T MR scanner (Achieva dStream, Philips Medical Systems, Best, Netherlands). An MRI-compatible orthotic fixture was used to perform dynamic movements inside the scanner(Figure 1). Individual range-of-motion was checked by the cilnician outside the MRI scanner and stops were placed on wire guides for controlling the range-of-motion inside the scanner. Each dynamic scan lasted for 18 seconds and ankle joint was moved (actively or passively) through the range-of-motion for each subject within 18 seconds with a rotation speed of 4•/s to 5•/s (depending on the range). Both passive and active acquisitions were standardized for all subjects as detailed in Annexe A. Dynamic images were acquired in the oblique axis of the ankle motion. This was done by acquiring axial dynamic scans first and then obliquely orienting the scanning plane for images to be acquired in sagittal plane.

The scanning protocol included one high-resolution static 3D scan of the ankle joint with a resolution of 0.26 × 0.26 × 0.8mm (T1-weighted gradient-echo, flip angle 10 • , matrix Guide wires helped the technician to control and operate the passive plantar-dorsi flexion movements of the ankle joint while rest of the limb is fixed using straps at foot, tibia and mid thigh locations. Position of the ankle rest can be adjusted based on the limb length. Knee angle can be adjusted from full extension to 45°flexion. Cushioned ankle rest supports the foot to be imaged.

Data post-processing

Temporal tracking of the ankle bones is a challenging task because of nonlinearly articulated joint motion. Image registration is an iterative process maximizing a similarity measure between two images (source and target) in order to find the optimal geometric transformation that best align them. In this work, we used a locally-linear intensity-based registration method to estimate bone motion, and then we fused these local transformations to compute temporal dense deformation fields.

Motion estimation

The bones of interest were manually segmented in the high-resolution static image (Figure 2) and then automatically propagated throughout the dynamic low-resolution sequence using intensity-based registration.

The proposed approach for motion estimation consisted in two steps : 1) estimating a mapping for each bone between static MRI data and dynamic sequences, 2) estimating the relative bone motions using the dynamic sequence. step1 : The first step was to compute the transformations from the static image to the dynamic images. In order to make this step robust, the multi-resolution registration was carried out in two steps : First, the static image was globally (and rigidly) registered to each low-resolution time frame {D k } k∈1...K . This provided a set of rigid transforms {T S→D k } k∈1...K .

Second, initialized with the set of global rigid transforms {T S→D k } k∈1...K , the temporal position of each bone was refined thanks to a local rigid registration using bone masks. In this step, the static image is registered to each dynamic image {D k } k∈1...K using input weights on each of the bones separately. This provided another set of rigid transforms {T i S→D k } i=1,...,3 k=1,...,K .

To maximize the robustness of the algorithm, an image D k * from the dynamic sequence whose the content is the closest to the static image, was automatically detected. The underlying assumption is that the registration process is much simpler to converge for transformations close to the identity. To this end, the image of the sequence for which the overlap between regions of interest before and after bone-dependent registration was maximum, was selected as the reference image which served as a starting point for the motion estimation within the dynamic sequence. ). This allows to identify the dynamic MR image which is closest to the high resolution image and which will serve as a starting point for the motion estimation within the dynamic sequence. Last, local and rigid registrations are conducted from the reference towards the sequence first and last images.

Computation of temporal dense deformation fields

A dense deformation map of the joint from the static to each time frame was obtained by fusing rigid transforms corresponding to the bones of interest. For each bone, a non-negative weighting function was computed based on a distance to the mask [START_REF] Commowick | An efficient locally affine framework for the registration of anatomical structures[END_REF]. Each weighting function reflects the local influence of the rigid transform of each bone. For each voxel x, the weighting functions were computed as follows :

w i D k (x) = 1 1 + dist 2 (x, M i D k ) (1)
where : dist(x, M i D k ) is the Euclidean distance between the point x and the mask of the bone i propagated onto the dynamic image D k .

Given a point x in the high-resolution static image, the target location of this point in the dynamic image D k can be computed with the following equation [START_REF] Arsigny | Polyrigid and polyaffine transformations : a novel geometrical tool to deal with non-rigid deformations-application to the registration of histological slices[END_REF] :

T S→D k (x) = exp ( N ∑ i=1 wi D k (x) log(T i S→D k ) ) .x (2)
where T S→D k is the dense invertible deformation field from the static image S to the dynamic image D k , wi D k is a normalized weighting function (i.e.

∑ N i=1 wi D k (x) = 1, ∀x ∈ D k ).

Fast computation of dense deformation fields

In this section, we describe an efficient algorithm to compute deformation fields in dense regular grid based on the algebraic properties of the Log-Euclidean polyrigid framework.

In [START_REF] Arsigny | A log-euclidean polyaffine framework for locally rigid or affine registration[END_REF][START_REF] Arsigny | A fast and log-euclidean polyaffine framework for locally linear registration[END_REF], the matrix exponential was recursively approximated with a certain level of accuracy using the scaling and squaring method by taking into account that the matrix exponential is much simpler to compute for matrices close to zero via the Padé approximant. In this work, we computed the exponential mapping of equation 2 exactly using eigen decomposition. Assuming the log-euclidean mean of linear transformations

L(x) = ∑ N i=1 wi D k (x) log(T i S→D k )
to be diagonalizable (i.e. L(x) can be expressed in homogenous coordinates as : L(x) = P.D.P -1 (x), where the columns of P contains the eigenvectors {v k (x)} k∈{0...4} corresponding to the complex eigenvalues of L(x), {λ 1 (x), λ 2 (x), λ 3 (x), 1}), we broadcasted the exponentiation of transformation eigenvalues over all grid points using the following Equation :

e L(x) = P e D P -1 (x) = P                 e λ 1 0 0 0 0 e λ 2 0 0 0 0 e λ 3 0 0 0 0 e                 P -1 (x) (3)

Determination of ankle joint biomechanics 2.4.1. Definition of anatomically based coordinate systems

Anatomical coordinate systems {R i } i=1,...,N were defined on each bone in the highresolution image S following the same protocols used in [START_REF] Sheehan | In vivo talocrural and subtalar kinematics : a non-invasive 3D dynamic MRI study[END_REF] as illustrated in Figure 4, and then mapped into the neutral position D n using the estimated local-rigid transforms {T i S→Dn } i=1,...,N . The neutral dorsi-plantarflexion position was defined at an ankle angle of 90°as recommended by the ISB standards committee [START_REF] Wu | ISB recommendation on definitions of joint coordinate system of various joints for the reporting of human joint motion-part i : ankle, hip, and spine[END_REF]. However, the above defined neutral position cannot be always achieved during dynamic scanning. Therefore, the image in which the foot position was closest to a dorsi-plantarflexion position of zero degrees was selected as the neutral position.

(a) (b) (c)

Quantification of ankle joint kinematics

In this study, the kinematics were defined through the 3D registration of the entire bone volume. More precisely, rotations and translations of each bone i were determined from estimated rigid transforms {T i D k →D k+1 } i=1,...,N . Ranges of motion were estimated using low-resolution image data on multiple dynamic images (15 ankle positions). Kinematics of talocrural joint, subtalar joint, and calcaneal-tibial complex were then derived from

{T i D k →D k+1 } i=1,...,N .
Given the image coordinate system R I = (O, ⃗ x, ⃗ y, ⃗ z), we defined the transformation, M R I →R j , that mapped R I to an anatomical bone-based coordinate system R j = (O j , ⃗ x j , ⃗ y j , ⃗ z j ). This transformation changed a representation from the R j system to the R I system.

Individual bone transformation matrices expressed in R I were then converted into the transformation matrices expressed in the new anatomical coordinate systems R j as follows :

T i j (k) = (M R I →R j × T i D k →D k+1 × M -1 R I →R j ) k=1,.
..,K-1; i=1,...,N ; j=1,...,N where :

-T i j (k) : 4 × 4 rigid transformation matrix of the i th bone relative to the j th local bone coordinate system at time k.

-M R I →R j : the change of basis matrix.

Method evaluation

To evaluate the proposed registration approach for bone motion estimation, the propagated mask of each bone through the dynamic sequence was compared with a manually delineated bone on each time frame. The RMS error was computed as follows :

RM SE(B i k , B i m k ) = √ 1/n c . ∑ nc xc=1 dist 2 (x c , ζ B i k )
where n c is the total number of voxels of ground-truth contours (i.e. contours of B i m k ) and dist(x c , ζ B i k ) is the Euclidean distance between x c and the contour of B i k .

Results

Results on bone motion estimation were reported in Figures 89for the six subjects of the pilot study. In most of the cases, the DICE coefficient was greater than 0. Three-dimensional kinematics of the joint under both active and passive ankle dorsiplantarflexion movement using the current 3D dynamic MRI method were presented in Figures 5 to 7. Kinematics of all ankles (expect A5) were derived from the estimated temporal rigid transforms of each bone of interest, converted into the matrices of the calcaneus relative to the tibia, the talus relative to the tibia, and the calcaneus relative to the talus.

The talus mistracking for the ankle A5 was caused by the considerable missing portion from the bone in intermediate time frames due to the inappropriate initial placement of the flex coils.

The mean rotation about the Inferior-Superior and about the Posterior-Anterior axes (averaged across time for all subjects) was close to 0 degrees for both active and passive motion.

The mean range of passive rotation about the Medial-Lateral axis for talocrural an calcanealtibial joints was between -11 degrees of maximal dorsiflexion and +32 degrees of maximum plantarflexion with full extended knee. This demonstrates excellent correspondance with the mean range of passive rotation reported in [START_REF] Alanen | Ankle joint complex mobility of children 7 to 14 years old[END_REF] and validated for 245 healthy children between 7 and 14 years old.

Discussion

In this study, we have presented a method to determine full 3D-kinematics of the ankle joint from dynamic anatomical MRI data. An intensity-based registration pipeline has been proposed to estimate bone motion through the dynamic sequence. The robustness of this approach for both passive and active motion has been evaluated on in vivo data acquired for a pilot control cohort study. Although based on a successive estimation of transformations, the proposed approach demonstrates advantages in efficiency, reliability and robustness for all subjects.

The proposed algorithm is robust enough to image quality as well as motion artefacts in a way that it can determine kinematics even for unconstrained motion. This was evaluated by changing the delta-time ∆k separating each couple of source/target images during the successive rigid registrations and the algorithm is still robust even for ∆k = 5.

Normative passive kinematics show that the subjects maintained a consistent movement trajectory. Three-dimensional kinematics of the joint under passive ankle movement using the proposed method (Figures 5 to 7) were similar to ankle joint kinematics reported in other studies using velocity data [START_REF] Sheehan | In vivo talocrural and subtalar kinematics : a non-invasive 3D dynamic MRI study[END_REF]. Normative active kinematics show that the subjects were able to perform voluntary plantar-dorsiflexion between the extreme positions with the same temporal regularity as for passive motion. The comparison between passive and active kinematics shows closest temporal averages despite some temporal fluctuation of the standard deviation across subjects.

Although the fact that both our method and other motion methods tend toward the same goal which is bone motion tracking, there are some technical differences making it difficult to perform an objective comparison based on a common error parameter.

In [START_REF] Clarke | A non-invasive, 3D, dynamic MRI method for measuring muscle moment arms in vivo : demonstration in the human ankle joint and achilles tendon[END_REF], the authors were more interested in the Achilles tendon moment arms and the accuracy and precision of bone motion tracking were not reported. Clarke et al. [START_REF] Clarke | A non-invasive, 3D, dynamic MRI method for measuring muscle moment arms in vivo : demonstration in the human ankle joint and achilles tendon[END_REF] have discussed some limitations in their methods which can be particularly useful in adults. They have explicitly mentioned that their protocols require some modifications to be useful in the smaller joints of children. To satisfy the mesh density requirements for using the ICP algorithm, they have recommended the reduction of the slice gap for the dynamic scanning protocol while increasing the number of slices per time frame so that the slices span the joint width.

However, assuming that they have used ultrafast MRI sequences (with 8 sagittal slices/time frame compared to 6 slices in our work, and taking around 2 min to acquire 10 to 20 time frames compared to only 18 sec to acquire 15 time frames in our work), increasing the number of slices will increase the scanning time of dynamic data to more than 2 min which will be uncomfortable especially for the pathological subjects. The authors have also notified that the rotation speed for their study was slow (mean 1 • /s compared to 4.2 • /s in our work) in order to reduce the effects of motion artifacts, so the motion itself is reduced to almost static or quasi-static nature where no real functionality can be assessed. Our choice for FFE sequences for dynamic data acquisition was motivated by the fact that these techniques are less sensitive to motion artifacts [START_REF] Stehning | Fast isotropic volumetric coronary MR angiography using free-breathing 3D radial balanced FFE acquisition[END_REF] and more than 5 times faster than ultrafast MRI.

In the work of Sheehan et al. [START_REF] Sheehan | In vivo talocrural and subtalar kinematics : a non-invasive 3D dynamic MRI study[END_REF], the displacements of three to five points on each bone were calculated through Fourier integration and then converted into three-dimensional orientation angles. As part of the tracking process, they identified a series of vertices, defining regions on each bone of interest in the first time frame (i.e. three to five vertices per bone). The distance between these vertices in the first time frame was considered the absolute distance under the assumption that it should remain fixed throughout the movement because bones are rigid. Thus, for each bone, the difference in distance between vertices in each time frame relative to the absolute distance was defined as the tracking error. A potential limitation of this study [START_REF] Sheehan | In vivo talocrural and subtalar kinematics : a non-invasive 3D dynamic MRI study[END_REF] was the out of plane accuracy of the bone motion because the accuracy of cine-PC MRI is independent of the shape of the bone being tracked [START_REF] Behnam | Assessing the accuracy and precision of musculoskeletal motion tracking using cine-PC MRI on a 3.0 T platform[END_REF]. Another clinical limitation is that subjects with musculoskeletal disease who cannot complete a large number of repeated motions cannot be studied with cine-PC MRI techniques. Although the vertices were well-defined from 3D high-resolution adult data in previous work of Sheehan et al. [START_REF] Sheehan | In vivo talocrural and subtalar kinematics : a non-invasive 3D dynamic MRI study[END_REF],

it was hard to reproduce their exact locations from developing children data for comparison purposes because of the differences in bone size and topology between children and adults.

From another point of view, it was not sufficient to consider the distance between only some vertices as system bias and thus we have considered all the bone segmentations (i.e. all the volumetric information) for method evaluation.

The proposed method remains robust as long as the field-of-view (FOV) is well adjusted (i.e. when the FOV covers either the full or nearly-full anatomy of the ankle bones throughout the entire joint trajectory).

We have also extended the log-euclidean framework to estimate temporal dense deformation fields from multiple rigid transforms. The polyrigid approach provides a way to combine high resolution spatial information with temporal dynamics of joints. The output of this process is a series of high-resolution anatomical images portraying the different phases of the movement cycle. In this context, we proposed to compute the exponential mapping of the LEPF in an efficient and elegant way using matrix diagonalization-based techniques rather than using techniques involving approximation theory. The preservations of the bone shapes after the polyrigid fusion were checked by computing the Jacobian maps associated to the estimated deformation fields (as shown in figure 11).

To conclude, our work is the first effort to track ankle joint motion and deformation from volumetric image data using intensity-based non-rigid image registration rather than using explicit mechanistic models. It therefore achieves its goals without having to perform time-consuming manual segmentations. In immediate future, the proposed approach will be applied to compare the ankle joint dynamics of children with spastic equinus with age-matched healthy children. Since spastic equinus is typically defined as the inability to dorsiflex the foot above plantigrade, with the hindfoot in neutral position and the knee in extended position [START_REF] Kay | Outcome of gastrocnemius recession and tendo-achilles lengthening in ambulatory children with cerebral palsy[END_REF], the foot was constrained to a specific path using the fixture in this work. However, this set-up is problem-specific and can be changed or removed while evaluating the pathomechanics of other joint(s). We will also explore this method to extract fine biomechanical parameters of tendon, and cartilage contact mechanics (e.g. temporal joint space width) of the tibio-talar joint which is the primary joint responsible for plantarflexion and dorsiflexion of the ankle. Furthermore, these techniques can be applied to other joints and related musculoskeletal disorders. This suggests that anatomical dynamic MRI and dedicated image processing techniques can open a new way to study in vivo human joints.

Acknowledgments

The research leading to these results has received funding from Région Bretagne, Chaire

  576×576, FOV 150mm×150mm, TR/TE 7.81/2.75 ms, mean acquisition duration : 424.32 s), and three low-resolution dynamic sequences per child : two passive sequences for repeatability measures and one active sequence, all acquired with knee angle kept at full extension (approximately between 0 • to 10 • ). Each sequence is composed of 15 time frames with a voxel size of 0.57 × 0.57 × 8mm (flip angle 15 • , matrix 352×352, FOV 200mm×200mm, TR/TE 20.61/1.8 ms, acquisition duration : 18.98 s).

Figure 1 :

 1 Figure1: Orthotic fixture specially designed to acquire MRI data on pediatric ankle joint. The entire fixture is made out of MRI-compatible material. Guide wires helped the technician to control and operate the passive plantar-dorsi flexion movements of the ankle joint while rest of the limb is fixed using straps at foot, tibia and mid thigh locations. Position of the ankle rest can be adjusted based on the limb length. Knee angle can be adjusted from full extension to 45°flexion. Cushioned ankle rest supports the foot to be imaged.

Figure 2 :

 2 Figure 2 : Bones of interest : calcaneus (red), talus (green) and tibia (cyan). (a) : Mid-sagittal image from the high-resolution static scan ; (b) : Mid-coronal image from the high-resolution static scan ; (c) : Three dimensional rendering of segmented bones.

Figure 3 :

 3 Figure 3 : Proposed pipeline for ankle motion estimation : Bones of interest are segmented in the highresolution static image. This high-resolution MR image is globally (and rigidly) registered to each MR image of the dynamic sequence (1.), and the position of each bone is refined thanks to a local rigid registration using bone segmentations (2.). This allows to identify the dynamic MR image which is closest to the high resolution image and which will serve as a starting point for the motion estimation within the dynamic sequence. Last, local and rigid registrations are conducted from the reference towards the sequence first and last images.

Figure 4 :

 4 Figure 4 : Anatomical coordinate systems for one subject's calcaneus (a), talus (b), and tibia(c) as per ISB.

8

 8 and the RMSE was smaller than 0.6 indicating accurate propagations of bone masks over the dynamic sequence, for both passive and active motions. From the set of rigid transforms of each bone over the dynamic sequence, it was then possible to compute temporal dense deformation fields to synthesize a high-resolution dynamic MRI sequence. Using the fast computation technique described in Section 2.3.3, a dense deformation field was computed in 3min on an Intel ® Xeon ® Processor E3-1271 v3 3.60 GHz on a 576 × 576 × 90 regular grid and in 15min on a very high dimensional space (576 × 576 × 202 ≃ 67 million deformation vectors) which required a high capacity RAM because computations were performed in the complex domain for the purpose of enhancing the accuracy. Figure10illustrates, for one subject (A6), such high-resolution reconstructed data.

Figure 5 :

 5 Figure 5 : 3D normative kinematics of the calcaneal-tibial complex using the proposed image processing method. The kinematics of the calcaneus are represented relative to the tibial coordinate system defined in neutral position. Average rotation and translation were computed for both passive and active motion for the studied ankles. Standard deviation above and below the average line are shown (hatched shaded area for passive motion and green shaded area for active motion).

Figure 6 :

 6 Figure 6 : 3D normative kinematics of the talocrural (talar-tibial) joint using the proposed image processing method. The kinematics of the talus are represented relative to the tibial coordinate system in neutral position. Average rotation and translation were computed for both passive and active motion for the studied ankles. Standard deviation above and below the average line are shown (hatched shaded area for passive motion and green shaded area for active motion).

Figure 7 :

 7 Figure 7 : 3D normative kinematics of the subtalar joint using the proposed image processing method. The kinematics of the calcaneus are represented relative to the talar coordinate system in neutral position. Average rotation and translation were computed for both passive and active motion for the studied ankles. Standard deviation above and below the average line are shown (hatched shaded area for passive motion and green shaded area for active motion).

Figure 8 :

 8 Figure 8 : 3D DICE coefficients between manual and automatic bone segmentation maps for the six subjects of the pilot study taking account of the whole set of bone trajectories. Results for passive (left column) and active motion (right column) are separated. A dice value of 1 indicates perfect geometric alignment between automatic and ground truth segmentations.

Figure 9 :

 9 Figure 9 : Temporal evolution of the root-mean-square error (RMSE in voxels) between manual and automatic bone segmentations for subjects of the pilot study. Errors are represented using error bars. The average is shown with a dotted line with one standard deviation above and below this line. Results for passive (left column) and active motion (right column) are separated.

Figure 10 :

 10 Figure 10 : High resolution temporal reconstruction for one subject (subject 6 : passive motion). First row : original dynamic images, second row : corresponding high-resolution reconstructed images, downsampled to resemble to the original low-resolution images in order to validate the reconstruction accuracy. Each column corresponds to one time frame. Contours of the first time frame show the reconstruction accuracy (column 1) and the joint motion (columns 2 to 4).

  d'excellence INSERM-IMT Atlantique, Fondation de l'Avenir, Paris, France, and Fondation Motrice, Paris, France.

Figure 11 :

 11 Figure 11 : Jacobian map of the joint deformation fields that maps D 1 to D 2 for subject 4. The Jacobian determinant J(x) at a voxel x measures how the voxel volume changes after registration. It indicates a volume increase if > 1, and a volume decrease if < 1. The Jacobian determinant is equal to 1 inside bone segmentations (no voxel volume changes) as bones only perform linear rigid transformations (rotations and translations).

  The accuracy of the estimations of the successive rigid transforms was assessed by computing the DICE coefficient Dice(B i k , B i m k ) and the RMS error RM SE(B i k , B i m k ) where B i k was the propagated mask of the bone i on dynamic image D k and B i m k was a manually delineated mask of the same bone on the dynamic image D k . A DICE value close to one indicates that the bone masks have been well propagated through the entire dynamic sequence.

Algorithm 1 Bone motion estimation

1 : Input : segmentation of bones of interest in S. 2 : Mapping estimation from S to {D k } :

-Estimate rigid transforms {T S→D k } k=1,...,K -Estimate rigid transform {T i S→D k } k=1,...,K using bone mask M i , initialized with

..,2; i=1,...,N . 4 : Output : compute bone rigid transforms {T i S→D k } k=1,...,K; i=1,...,N by combining steps 2 and 3 through which he/she could hear a metronome. Fixture was then moved in a consistent passive plantar-dorsiflexion by the technician at the beat of the metronome and using the guide wires attached to the fixture. For active movement, no technician was present and children were asked to perform voluntary plantar-dorsiflexion between the extreme positions on the beat of the metronome. Meaning, the active motion was completely voluntary with all the muscles dynamically taking part in the production of the plantar-dorsiflexion movement.

Annexe B. Algorithm

The overall algorithm for bone motion estimation is described in Algorithm 1 and in Figure 3, using the following notations : S : high-resolution static image, {D k } k=1,...,K : set of low-resolution dynamic images, {M i } i=1,...,N : set of mask of bones of interest, T A→B : rigid transform from image A to image B, T i A→B :

rigid transform from image A to image B for the bone i, Dice(A,B) is the DICE coefficient which measures the overlap between segmented regions (Dice(A, B) = 2|A∩B| |A|+|B| ). Rigid registrations were performed using FSL-flirt [START_REF] Jenkinson | Improved optimization for the robust and accurate linear registration and motion correction of brain images[END_REF] with correlation ratio as similarity measure. A voxel-wise weighting based on the segmentation masks was used when estimating bone-dependent motion.

Code is available online at https://github.com/rousseau/dynMRI.