Juhoon Kim
email: jkim@net.t-labs.tu-berlin.de

Luigi Iannone

Anja Feldmann

A Deep Dive into the LISP Cache and what ISPs should know about it

Keywords: Locator/ID Separation, Next Generation Internet, Addressing and Routing Architectures, Measurements

Due to scalability issues that the current Internet is facing, the research community has re-discovered the Locator/ID Split paradigm. As the name suggests, this paradigm is based on the idea of separating the identity from the location of end-systems, in order to increase the scalability of the Internet architecture. One of the most successful proposals, currently under discussion at the IETF, is LISP (Locator/ID Separation Protocol). A critical component of LISP, from a performance and resources consumption perspective, as well as from a security point of view, is the LISP Cache. The LISP Cache is meant to temporarily store mappings, i.e., the bindings between identifiers and locations, in order to provide routers with the knowledge of where to forward packets. This paper presents a thorough analysis of such a component, based on real packet-level traces. Furthermore, the implications of policies to increase the level of security of LISP are also analyzed. Our results prove that even a timeout as short as 60 seconds provides high hit ratio and that the impact of using security policies is small.

Introduction

In the last years, there has been an increasing interest in the Locator/ID Split paradigm, which is recognized to be a strong candidate to become a fundamental building block of the Next Generation Internet. Differently from the current architecture, where one single namespace, namely the IP address space, is used for both indentifying and locating end-systems, in the Locator/ID Split paradigm two different namespaces are used: the ID space and the Locator space, to respectively identify and locate end-systems. Such a separation aims at solving the scalability issues that the current Internet is facing [START_REF] Meyer | Report from the IAB Workshop on Routing and Addressing[END_REF], mainly concerning the continuously increasing BGP routing table [START_REF]BGP Routing Table Analysis Report[END_REF], but also concerning addressing, mobility, multi-homing [START_REF] Quoitin | Evaluating the Benefits of the Locator/Identifier Separation[END_REF], and inter-domain traffic engineering [START_REF] Saucez | Interdomain Traffic Engineering in a Locator/Identifier Separation context[END_REF].

The main effort to tackle these issues has started three years ago, when the Routing Research Group (RRG) has been rechartered to explore the possibility to enhance the Internet architecture by introducing some form of Locator/ID Split. Among the numerous proposals introduced since then [START_REF] Li | Recommendation for a Routing Architecture[END_REF], the most successful so far are ILNP (Identifier/Locator Network Protocol [START_REF] Atkinson | ILNP Concept of Operations[END_REF]), which has been adopted by the RRG, and LISP (Locator/ID Separation Protocol [START_REF] Farinacci | Locator/ID Separation Protocol (LISP). IETF -Internet Engineering Task Force[END_REF]), which has been chartered as Working Group in the IETF and has been already deployed in an international test-bed [START_REF]International LISP Infrastructure[END_REF]. On the one hand, all the proposals improve the Internet architecture, in a way or another, by some form of Locator/ID Split. On the other hand, they introduce the necessity to distribute and store bindings between IDs and Locators (i.e., mappings), which is the key critical component that can make the difference. Hence, the importance of evaluating what is the cost (i.e., memory consumption, overhead, ...) of storing and distributing these mappings. In addition, the way they are managed does also have an impact on the robustness of the proposed architecture with respect to security threats.

Focusing on LISP as reference protocol, in this paper we present a thorough analysis of the mapping cache. In a previous work, Iannone et al. analyzed a Netflow traffic trace from a small/medium sized university campus [START_REF] Iannone | On the Cost of Caching Locator/ID Mappings[END_REF]. We build our work on the same methodology, but develop a new emulator in order to go beyond previous results and provide a deeper analysis. For this purpose, we use two [START_REF] Saucez | LISP Security Threats[END_REF]-hour packet-level traces from a large European ISP, taken at different periods of the year. Further, Saucez et al. [START_REF] Saucez | LISP Security Threats[END_REF], in their security analysis of the LISP protocol, highlighted that a number of attacks can be avoided by slightly modifying the policy that LISP uses to encapsulate and decapsulate packets. Hence, we evaluate what are the implications, from a scalability point of view, of running LISP in such a more secure manner, which we call symmetric LISP.

The contributions of this paper are many-fold. We thoroughly analyze the behavior of a LISP Cache for traffic of a large number of DSL customers. Our results show that it is possible to maintain a high hit ratio with relatively small cache sizes, showing that there is no use in large timeouts. Compared to [START_REF] Iannone | On the Cost of Caching Locator/ID Mappings[END_REF], we show how the LISP Cache has good scalability properties. We analyze what and how much would change if instead of running vanilla LISP, i.e., as defined in the specifications, symmetric LISP is used in order to increase security. Our results are important for ISPs in order to identify and quantify the resources needed to deploy LISP, with respect to the level of resiliency that they want to guarantee.

The remainder of this paper is structured as follows. In Section 2 we present the related work. In Section 3 we provide an overview of LISP in order to highlight how the LISP Cache works and the motivation for symmetric LISP. Section 4 describes how we collected and analyzed the traces used to obtain the results presented in Section 5. Finally, we draw our conclusions in Section 6.

Related Work

The concept of Locator/ID Split was already discussed, but never widely explored, in the late 90s ( [START_REF] Saltzer | On the Naming and Binding of Network Destinations[END_REF], [START_REF] Hiden | New Scheme for Internet Routing and Addressing (ENCAPS) for IPNG[END_REF]); rather used to solve specific issues, like cryptographic security with HIP (Host Identity Protocol [START_REF] Moskowitz | Host Identity Protocol (HIP) Architecture[END_REF]) or multi-homing for IPv6 with Shim6 [START_REF] Nordmark | Level 3 Multihoming Shim Protocol for IPv6[END_REF]. The situation has changed after the rechartering of the RRG, with quite a number of proposals being discussed, either based on tunneling (e.g., LISP [START_REF] Farinacci | Locator/ID Separation Protocol (LISP). IETF -Internet Engineering Task Force[END_REF], HAIR [START_REF] Feldmann | HAIR: Hierarchical Architecture for Internet Routing[END_REF]) or on address rewriting (e.g., ILNP [START_REF] Atkinson | ILNP Concept of Operations[END_REF], Six/One [START_REF] Vogt | Six/One: A Solution for Routing and Addressing in IPv6[END_REF]). Despite the plethora of proposals, very few works have tackled the evaluation of such a critical component as the cache.

Kim et al. [START_REF] Kim | Revisiting Route Caching: The World Should Be Flat[END_REF] studied route caching based on Netflow data, showing how caching is feasible also for large ISPs. However, they used a generic approach, not focused on LISP, and a limited cache size. Iannone et al. [START_REF] Iannone | On the Cost of Caching Locator/ID Mappings[END_REF] have performed an initial study, without considering security issues, based on a single Netflow trace of a small/medium sized university campus. We adopt a similar methodology to provide comparable results, which we discuss in Section 6. Nevertheless, because we base our results on two 24-hour packet level traces captured in different periods we are able to provide a deeper analysis, including security aspects.

In the work of Zhang et al. [START_REF] Chen | Evaluating the Performance on ID/Loc Mapping[END_REF], the authors propose a simple LISP Cache model with bounded size, using a Least Recently Used (LRU) policy to replace stale entries. 1 Their evaluation is based on a 24-hour trace of outbound traffic of one link of the CERNET backbone, the China Education and Research Network. In our analysis we take a different approach, assuming that the cache is not limited in size. This is a reasonable assumption since, as we will show later, even for large ISPs there is no need to use large caches and by tuning the expiration timer, caches do not grow so large as to hit memory limits of current routers.

Jakab et al. [START_REF] Jakab | LISP-TREE: A DNS Hierarchy to Support the LISP Mapping System[END_REF] propose a LISP simulator able to emulate the behavior of different mapping systems. However, they focused on the evaluation of the latency of their own mapping system solution, namely LISP-TREE, and compared it to other approaches but neglecting the analysis of the LISP Cache.

LISP Overview

The main idea of the Locator/ID Separation Protocol (LISP [START_REF] Farinacci | Locator/ID Separation Protocol (LISP). IETF -Internet Engineering Task Force[END_REF]) is to split the IP addressing space in two orthogonal spaces, one used to identify the end-hosts, and one used to locate them in the Internet topology. By re-using IP addresses, LISP is incrementally deployable and meant to be used on the border routers of stub domains, like in the scenario presented in Figure 1. Stub domains use internally an IP prefix, called EID-Prefix as for End-system IDentifier, which is part of the ID space. This prefix does not need to be globally routable, since it is used only for routing in the local domain. On the contrary, the core of the Internet, known as Default Free Zone (DFZ), will use globally routable IP addresses that are part of the locator space. In particular, the IP addresses used by stub domains' border routers on their upstream interfaces (toward the provider) are part of such a space and represent the Routing LOCators (RLOCs), since they allow to locate EID in the Internet. The binding between an EID-Prefix and the set of RLOCs that locate it is a mapping. EID-Prefixes are no longer announced in the DFZ, allowing a size reduction of BGP's routing table. Evaluating such kind of benefits is beyond the scope of this paper. Further information can be found in the work of Quoitin et al. [START_REF] Quoitin | Evaluating the Benefits of the Locator/Identifier Separation[END_REF] and Iannone et al. [START_REF] Iannone | Modeling the Economics of Loc/ID Split for the Future Internet[END_REF].

In the context of LISP, end-hosts generate normal IP packets using EIDs as source and destination addresses, where the destination EID is obtained, for instance, through a DNS query. LISP then tunnels those packets from the border router of the source domain to the border router of the destination domain, using the RLOCs as source and destination addresses in the outer header. For tunneling, LISP uses an IP-over-UDP approach, with the addition of a LISP-specific header between the outer UDP header and the inner (original) IP header. 2 In the LISP terminology, the border router in the source domain performing the encapsulation is called Ingress Tunnel Router (ITR), while the one in the destination domain, performing the decapsulation, is called Egress Tunnel Router (ETR). In general, they are just named xTRs. To perform tunneling operations routers use two data stores, namely the LISP Database and the LISP Cache.

The LISP Database stores mappings that bind the local EID-Prefixes (i.e., inside the local domain) to a set of RLOCs belonging to the xTRs deployed in the domain. The purpose of the LISP Database is two-fold. For outgoing packets, if a mapping exists for the source EID it means that the packet has to be LISP encapsulated and the source RLOC is selected from the set of RLOCs associated to the source EID-Prefix. For incoming packets, if a mapping exists for the destination EID, then the packets are decapsulated. The LISP Database is statically configured on each xTR. Its size is directly proportional to the number of the EID-Prefixes and xTRs that are part of the local domain. Due to its static nature and limited size, the LISP Database does not present any scalability issue and is not further analyzed in this paper.

The LISP Cache temporarily stores the mappings for EID-Prefixes that are not part of the local domain. This is necessary to correctly encapsulate outgoing packets, in particular to select the RLOC to be used as destination address in the outer header. Mappings are stored only for the time that are used to encapsulate packets, otherwise, after a timeout, they are purged from the cache. While critical for the dimensioning of the system, the LISP specification does not provide any value for this timeout, leaving the choice (or responsibility) to the implementers and system administrators. It should be clear that since the entries in the cache can expire, they are also entered in an on-demand fashion. This makes the cache so critical, since its content, size, and efficiency is totally traffic driven. In particular, the first outgoing packet, destined to an EID for which there is no mapping in the cache, triggers a cache-miss. Such an event, in turn, triggers a query message that the ITR sends to the Mapping Distribution System. 3 The latter is a lookup infrastructure designed to retrieve the mapping for the destination EID of the packet that triggered the query. This is the same principle of DNS, which allows to retrieve the IP address(es) of a server from its FQDN (Fully Qualified Domain Name). Insofar, several Mapping Distribution Systems have been proposed, but only the BGP-based ALT [START_REF] Farinacci | LISP Alternative Topology (LISP+ALT). IETF -Internet Engineering Task Force[END_REF] is an official working item of the LISP Working Group. It is out of the scope of this paper to compare different mapping systems. Further information can be found in the work of Jakab et al. [START_REF] Jakab | LISP-TREE: A DNS Hierarchy to Support the LISP Mapping System[END_REF] and Mathy et al. [START_REF] Mathy | LISP-DHT: Towards a DHT to map Identifiers onto Locators[END_REF].

                                        
The attentive reader should have noticed that there is a fundamental asymmetry in LISP when performing encapsulation and decapsulation operations. Indeed, while doing encapsulation the ITR uses the cache to locate what to put in the outer header, on the contrary, the ETR, while decapsulating the packet, does not use the cache. Rather, it just performs the checks described previously for the LISP Database and eventually decapsulates the packet.

This asymmetry leaves the LISP protocol vulnerable against specific attacks on the ETR, as pointed out by the work of Saucez et al. [START_REF] Saucez | LISP Security Threats[END_REF]. The authors describe attacks exploiting data packets as well as attacks that leverage on fields of the LISP header. Attacks are not limited to classic DoS and spoofing, but regard also cache poisoning, i.e., injecting wrong information into the cache, or cache overflow, i.e., saturating the cache and increasing the number of cache-misses. The analysis they performed, led the authors to conclude that to increase the level of security of the whole LISP architecture, the best solution is to use a symmetric model with a drop policy. Put differently, the authors suggest to only allow encapsulation as well as decapsulation if mappings are present in both the LISP Cache and the LISP Database. Otherwise, when the mapping is not in the cache a miss is generated and the packet is dropped. The rationale behind this suggestion is the assumption that the mapping distribution system is secured and trusted, hence, sanity checks can be performed at both ends, when encapsulating and when decapsulating. In particular, for the latter operation it means that if the packet contains information that is not coherent with the content of the cache, it is dropped. The implications of introducing the symmetric model to increase security are two-fold: on the one hand, performing additional checks when performing decapsulation may reduce the performances, however this is a common price for security mechanisms. On the other hand, this means that the size of the cache, its dynamics, and the control traffic overhead is increased.

In the second half of this paper, we will refer to the symmetric model as symmetric LISP and we will assess its impact as compared to vanilla LISP to evaluate what are the trade-offs to increase security in the LISP protocol.

LISP Emulation

To evaluate LISP we implemented pcap2lisp, which is meant to emulate the behavior of LISP's xTRs. The emulator is essentially designed to be fed with pcap-formatted traffic data and mimics as much as possible the LISP architecture; hence we based its implementation on two main modules (cf., Figure 2). The LISP Database, which is a manually configured list of internal network prefixes (EIDs). The LISP Cache, which stores EID-Prefixes and related statistic. Besides these two modules there is a central logic that creates the correct statistics, periodical reports that are written in logs, and that overviews the correct management of the cache timeouts. In addition, we use a local BGP prefixes database, fed with the list of BGP prefixes published by the iPlane Project [START_REF] Madhyastha | iPlane: An Information Plane for Distributed Services. 7th symposium on Operating Systems Design and Implementation[END_REF]. The database is used to group EID-to-RLOCs mappings with the granularity of existing BGP prefixes, because, as for today, there is no sufficient information to predict what will be the granularity of mappings in a LISP-enabled Internet. This BGP granularity follows the methodology proposed in [START_REF] Iannone | On the Cost of Caching Locator/ID Mappings[END_REF], thus making the results comparable. Furthermore, such an approach allows using the BGP database as a black box returning the mappings needed to feed the cache.

We fed pcap2lispwith two sets of anonymized traffic data collected within a large European ISP for 24 hours in April 2009 (APR09) and in August 2009 (AUG09) a vantage point covering more than 20,000 DSL lines and already validated in other studies ( [START_REF] Ager | Revisiting cacheability in times of user generated content[END_REF][START_REF] Kim | Today's usenet usage: Characterizing NNTP traffic[END_REF]).

Results

Similarly to the work in [START_REF] Iannone | On the Cost of Caching Locator/ID Mappings[END_REF], we used three different cache timeout values, respectively 60 seconds, 180 seconds (i.e., three minutes), and 1,800 seconds (i.e., 30 minutes). The reason why we choose 60 seconds, instead of 300 minutes like in [START_REF] Iannone | On the Cost of Caching Locator/ID Mappings[END_REF], is because that work already proved that the 300 minutes timeout value is inefficient considering the hit ratio vs. the cache size.

It is useful starting by identifying the working set, in this case represented by the number of observed BGP prefixes in our measurement environment. Figure 3 depicts the total number of contacted prefixes per minute, as well as the breakdown of incoming, outgoing, and bi-directional traffic. What we can identify in Figure 3 is that the majority of the observed prefixes (i.e., 70.1% in average in both traces) are bi-directional. One thing to keep in mind is that in the case of vanilla LISP, incoming packets do not play any role in the LISP cache, as explained in Section 3. On the contrary, when using symmetric LISP, incoming prefixes have an impact on the cache since the ETR does need a mapping. Hence, the 11.4% (average value in both traces) of incoming prefixes are a key differentiation between the two versions of LISP.

Vanilla LISP

A key benefit of LISP is the reduction of the routing table size in the DFZ [START_REF] Quoitin | Evaluating the Benefits of the Locator/Identifier Separation[END_REF], however, it exists a trade-off between the reduction of the routing table size and the LISP Cache size. The analysis of the caching cost is important for estimating the actual benefits of LISP. Figure 4 illustrates the number of entries and the size of the cache measured per one-minute time bin. Note that there are two measurement scales for the y-axis in the plot. While the left scale indicates the number of cache entries, the right scale indicates the size of the cache expressed in MBytes. The drop of the cache size at the end of the plot is due to the fact that we run the emulation until all cache entries are expired due to the timeout.

We can observe that the number of entries and the size of the cache are pretty low for the 60 and 180 seconds timeout cases (respectively slightly more than 10,000 and 20,000) and their curves show a spiky behavior. The 1,800 seconds timeout is much higher but with smoother changes, however, it takes about 30 minutes to reach the stable working set (almost 60,000 entries).

In the plot, and differently from [START_REF] Iannone | On the Cost of Caching Locator/ID Mappings[END_REF], we calculate the size of the cache using the size of the data structure of a real LISP implementation, namely Open-LISP [START_REF]The OpenLISP Project[END_REF]. In OpenLISP each entry of the cache is a radix node containing all the information of the mapping and pointing to a chained list of RLOCs nodes. Thus, the size C size of the cache is given by:

C size = N E × (Radix s + N R × RLOC s ).
Where N E and N R represent respectively the number of entries in the cache and the number of RLOCs per entry. Radix s is the size of a radix node (56 bytes in OpenLISP), while RLOC s is the size of the RLOC node (48 bytes in OpenLISP). For the results in Figure 4 we assume two RLOCs per mapping.

Figure 5 depicts traffic volume and the traffic overhead of LISP expressed in MBytes/sec. This plot is based on 60 seconds timeout, because, as we will show next, the traffic overhead is in inverse proportion to the timeout value. Therefore, the estimation of the traffic overhead reported in this analysis shows the maximum values. The line plotted in positive values indicates outgoing traffic, whereas the line plotted in negative values indicates incoming traffic. The shading over and under the traffic line indicates the traffic overhead. The overall overhead O can be evaluated by adding the overhead due to the LISP encapsulation, i.e., the number of packets N P multiplied the size of the prepended header E H and the volume of traffic generated to request and receive the missing mappings. Assuming the one single request/reply exchange is performed for each miss, the latter is given by the number of misses N M multiplied by the size of the Map-Request message M REQ , for outgoing traffic, or the size of the Map-Reply, for incoming traffic, which is given by a fixed size M REP plus the size of each RLOCs record R times the number of RLOCs N R . Note that the Map-Request and Map-Reply are the standard messages defined by LISP to request and receive a mapping from the mapping distribution system, independently from the specific instance of this latter. Putting everything together, for outgoing traffic we have:

O out = N Pout × E H + (N M × M REQ )
; while for the incoming traffic we obtain:

O in = N P in × E H + (N M × (M REP + N R × R)).
From the LISP specifications [START_REF] Farinacci | Locator/ID Separation Protocol (LISP). IETF -Internet Engineering Task Force[END_REF] the size of E H is 36 bytes, the size of M REQ is (without any trailing data) 24 bytes, the size of M REP (without any trailing data) is 28 bytes, to which we need to add 12 bytes for each RLOC record R. With these numbers is possible to plot the overhead in Figure 5, which is between 3.6% and 5.2% for APR09 and 3.6% and 4.6% for AUG09.

As the previous analysis shows, the selection of an appropriate timeout value is of prime importance for the efficiency of the cache. Thus, we further analyze the implications of different timeout values. We first evaluate the efficiency of cache by investigating the ratio between cache misses and cache hits. depicts the number of hits and misses in different timeout values. Even with the y-axis in logarithmic scale it is not possible to distinguish the lines of the cache hits, since they are overlapping each other. Hence, we need to look at the difference in the number of cache misses. Despite the identifiable difference, it is hard to avoid the claim that the benefit from the longer timeout value is insignificant compared to its cost in terms of size. Indeed, the hit ratio is always higher than 99% for every timeout value, while the size can be more than 5 times bigger when we compare the 60 seconds case with the 1,800 seconds case.

To confirm that a large cache is not useful, in Figure 7 we plot the CDF of the traffic volume forwarded by each entry, showing that the vast majority of cache entries carry less than 1 MBytes of traffic. Furthermore, Figure 8 shows how at least 50% of the entries have a lifetime only slightly longer than the timeout value. Both CDFs are the consequence of the well-known Internet phenomena that a small number of prefixes are responsible for the majority of Internet traffic. Coupled with the evidence mentioned above, we draw the inference that the minimum timeout value (60 seconds) is the most cost-beneficial.

Symmetric LISP

With respect to the security threats discussed in Section 3, we evaluate the extra cache size and the extra traffic overhead when the symmetric model is used. Since the extension of the mapping mechanism to incoming traffic is the main idea of the symmetric model, an increase in the number of entries and the size of the cache is expected and confirmed by the plot in Figure 9. For a clear comparison, the cache size shown in Figure 4 is presented again in gray color (lighter gray in black and white print). From Table 1, we see that there is a 6.8% to 9.9% increase in the cache size, when the symmetric model is used. Interestingly, the higher the timeout value, the lower the increase; we argue that it is due to the fact that longer timeouts increase entries' reuse probability. Now, we expand the analysis to include the extra traffic overhead. Recall that each cache-miss, even for incoming traffic, triggers at least two messages (Map-Request and Map-Reply). Figure 10 compares the number of cache-miss under vanilla LISP and symmetric LISP. The most interesting result shown in the plot is the fact that the number of miss caused by incoming packets does not correspond to the increase in the number of cache-miss. Indeed, we observe that 2.8% to 13.6% of cache-miss are additionally produced under symmetric LISP, but the percentage of cache-miss due to incoming packets is much higher, while the number of cache-miss due to outgoing packets shrinks.

When translated into traffic overhead, the increase of miss generates a small increase in bandwidth consumption. Nevertheless, after applying the same computation like for vanilla LISP, the measured increase in bandwidth is always lower than 0.5%, hence not critical if not negligible.

Conclusion

The scalability issues of the current Internet have triggered an important amount of research on the Locator/ID Split paradigm. LISP, in particular, has gained momentum and seems to be a strong candidate to be widely adopted. The critical component in the LISP architecture, from a scalability and security point of view, is the LISP Cache. By using two 24-hour traces of a large European ISP, we thoroughly analyzed this component. In the present work we adopted the same methodology as in [START_REF] Iannone | On the Cost of Caching Locator/ID Mappings[END_REF], in this way, the two works nicely complement and validate each other. An interesting comparison with the previous work is about the size of the cache vs. the number of users. In our case the size of the cache is almost doubled, however, in our traces there is more than three times the number of users than in the traces analyzed in [START_REF] Iannone | On the Cost of Caching Locator/ID Mappings[END_REF]. This suggests that the average size of the LISP cache does not grow linearly with the number of end-systems. This is an important result, since it confirms that the LISP cache has good scalability properties. Future work will focus on further analyzing this aspect.

We went further in the analysis, mainly, but not only, in two points: the analysis of very short cache timeout and the analysis of the symmetric LISP model, which improves security. When the timeout value is as small as 60 seconds, the measurements prove that efficiency of the cache is still very high, with more than 99% hit-ratio, while the size reduces almost by half compared to the 180 seconds timeout case, and it is almost five times smaller compared to than the 1,800 seconds timeout case.

Another main contribution of the present work concerns the symmetric LISP model. The presented analysis shows that the increase in the size of the cache and the generated overhead is the order of 13% (for the 60 second timeout case). This looks like a very low cost compared to the security benefits, which have a very high value for ISPs and vendors.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Example of LISP deployment and communication.

Fig. 3 .

 3 Fig. 3. Report of the number of correspondent prefixes per minute.

Fig. 4 .

 4 Fig. 4. Number of entries and cache size (assuming two RLOCs) per minute.

Fig. 5 .

 5 Fig. 5. Traffic volume and overhead due to LISP (60 seconds timeout value).

Fig. 6 .

 6 Fig. 6. Number of hits and misses per minute.

Fig. 7 .

 7 Fig. 7. Cumulative Distribution Function of the traffic volume per cache entry.

Fig. 8 .

 8 Fig. 8. Cumulative Distribution Function of entries' lifetime.

Fig. 9 .Fig. 10 .

 910 Fig. 9. Extra overhead caused by the security enhancement.

Table 1 .

 1 Cache size (in Mbytes) increase due to security enhancement (two RLOCs)

			APR09			AUG09
	Timeout	Vanilla	Symmetric	Difference	Vanilla	Symmetric	Difference
	60 sec	3.47	3.76	+8.36%	3.04	3.34	+9.87%
	180 sec	5.32	5.74	+7.89%	4.85	5.29	+9.07%
	1800 sec	11.29	12.06	+6.82%	10.99	11.76	+7.01%

In the rest of the paper we will use the term LISP Cache and cache interchangeably.

The LISP header contains information about RLOCs' reachability and traffic engineering. Details on this header can be found in the protocol specification ([START_REF] Farinacci | Locator/ID Separation Protocol (LISP). IETF -Internet Engineering Task Force[END_REF],[START_REF] Iannone | LISP Map Versioning[END_REF]).

[START_REF]The OpenLISP Project[END_REF] Note that, in case of cache-miss, the packet that triggered it cannot be encasulated, since there is no mapping available. The LISP specifications do not explicitily describe what to do with the packet, however, it is out of the scope of the present work to evaluate this particular aspect.

Acknowledgments

This work is partially funded by the Trilogy Project (http://www.trilogy-project.org), a research initiative (ICT-216372) of the European Community under its Seventh Framework Program.


Paulo Salvador
email: salvador@ua.pt

António Nogueira
email: nogueira@ua.pt

Markov Modulated Bi-variate Gaussian Processes for Mobility Modeling and Location Prediction

Keywords: mobility modeling, location prediction, Markov model, modulated bi-variate Gaussian process

A general-purpose and useful mobility model must be able to describe complex movement dynamics, correlate movement dynamics with the nodes geographic position and be sufficiently generic to map the characteristics of the movement dynamics to general geographic regions. Moreover, it should also be possible to infer the mobility model parameters from empirical data and predict the location of any node based on known positions. Having the ability to model movement and predict future positioning of mobile nodes in complex environments can be very important to several operational and management tasks. Existing mobility modeling approaches are based on simple and limited models, specifically designed for particular application scenarios or requiring the complete knowledge of the mobility environment. These features make them unusable in complex scenarios with no (or partial) knowledge of the environment. This paper presents a discrete time Markov Modulated Bi-variate Gaussian Process that is able to characterize the position and mobility of any mobile node, assuming that the position within a generic sub-region can be described by a bi-variate Gaussian distribution and the transition between sub-regions can be described by an underlying (homogeneous) Markov chain. With this approach, it is possible to describe the mobile node movement within and between a set of geographic regions determined by the model itself and, due to the Markovian nature of the model, it is also possible to capture complex dynamics and calculate the future probabilistic position of a mobile node. The proposed approach can be applied to scenarios where the possible pathways are unknown or too complex to consider in a real model that must make a prediction in a very short time. The results obtained by applying the proposed model to real and publicly available data demonstrate the accuracy and utility of this approach: the model was able to efficiently describe the movement patterns of mobile nodes and predict their future position. Besides, the model has also revealed higher performances when compared to other modeling approaches.

Introduction

Being able to describe the movement and predict future positioning of mobile nodes in complex environments can be very important to different operational and management tasks: evaluation of different system design alternatives and implementation costs, simulation and study of new routing protocols for Mobile and Delay Tolerant Networks, knowledge of the nodes distribution for performance evaluation and message delivery optimization, prediction of the nodes preferred locations for connectivity issues, etc. Other research topics, like mobile computing, can also greatly benefit from the ability to track and predict the location of mobile devices: an accurate location predictor can significantly improve the performance or reliability of many applications in pervasive computing.

There are two different strands of the mobile modeling problem: the ability to describe or track the movement of mobile nodes, which is mainly used for simulation purposes, and the ability to predict the mobile nodes position, which can be useful to several management tasks. An efficient mobility model must be able to describe complex movement dynamics, correlate the movement dynamics with the nodes position and be sufficiently generic to map movement's dynamics characteristics to general geographic regions where the road/path topology is unknown or imprecise. The parameters of the mobility model should be easily inferred from real mobility data and the model should provide an almost realtime prediction of the nodes position. This paper proposes a discrete time Markov Modulated Bi-variate Gaussian Process to characterize the position and mobility of any mobile node, assuming that the position can be described by a bi-variate Gaussian distribution and the transition between regions can be described by an underlying (homogeneous) Markov chain. This model is able to describe the mobile node movement within and between a set of geographic regions determined by the model itself. Due to the Markov nature of the model, it is possible to determine the future probabilistic position of a mobile node.

The results obtained by applying this model to real mobility data prove that this approach can accurately predict the future position of mobile nodes presenting, at the same time, higher performance than other prediction models that were selected for comparison.

The paper is organized as follows: section 2 presents some of the most relevant works related to the proposed model; section 3 presents the mathematical background on the multidimensional discrete time Markov Modulated Bi-variate Gaussian Process (dMMBVGP) and the description of the proposed model; section 4 presents the algorithmic details of the inference procedure; section 5 presents the results obtained in a proof of concept and a discussion of those results; finally, section 6 presents the main conclusions.

Related work

Multiple proposals for modeling movement of mobile nodes have already been published [START_REF] Camp | A survey of mobility models for ad hoc network research[END_REF]. The proposed models can be mapped into four categories: (i) purely random; (ii) with temporal dependency, when the movement of a mobile node is likely to be affected by its movement history; (iii) with spatial dependency between subjects, when mobile nodes tend to travel in a correlated manner, and (iv) with geographical dependency/restriction, when the movement of nodes is bounded by streets, freeways or obstacles.

In purely random mobility models, mobile nodes move randomly and freely without restrictions, that is, the destination, speed and direction are all chosen randomly and independently of other nodes. The Random Waypoint model [START_REF] Broch | A performance comparison of multi-hop wireless ad hoc network routing protocols[END_REF] and two of its variants, named Random Walk model [START_REF] Camp | A survey of mobility models for ad hoc network research[END_REF] and Random Direction model [START_REF] Royer | An analysis of the optimum node density for ad hoc mobile networks[END_REF], are the most important examples of this type of models. Obviously, these models fail to represent some mobility characteristics that are likely to exist in mobile ad hoc networks: the temporal and spatial dependency of velocity and the geographic restrictions of movement. Since, in general, the velocities of a single node at different time slots are "correlated", other models were proposed to capture this temporal dependency behavior: the Gauss-Markov mobility model [START_REF] Liang | Predictive distance-based mobility management for PCS networks[END_REF], where the velocity of a mobile node is assumed to be correlated over time and modeled as a Gauss-Markov stochastic process and the Smooth Random Mobility Model [START_REF] Bettstetter | Mobility modeling in wireless networks: Categorization, smooth movement and border effects[END_REF], where the speed and direction of the node movement is incrementally and smoothly changed in order to generate more realistic movement behaviors. In [START_REF] Hsu | Modeling spatial and temporal dependencies of user mobility in wireless mobile networks[END_REF], the authors proposed a time-variant community mobility model that realistically captures spatial and temporal correlations. Synthetic traces that match the characteristics of different mobility traces, including wireless LAN traces, vehicular mobility traces and human encounter traces, were generated, demonstrating the flexibility of the proposed model. Several Markovian approaches have also been proposed to address the mobility problem: reference [START_REF] Hyytia | A markovian waypoint mobilitymodel with application to hotspot modeling[END_REF] proposed the Markovian Waypoint model, where mobile nodes move along a straight line segment from one waypoint to another at a certain speed. In this approach, waypoints are obtained from a discrete time Markov process having a subset of IR 2 as the state space. Other mobility models, such as [START_REF] Bettstetter | Mobility modeling in wireless networks: Categorization, smooth movement and border effects[END_REF], are able to produce more smoothly turning trajectories.

In most real life applications, a node's movement is subject to the environment. Therefore, nodes may move in a pseudo-random way on predefined pathways in the environment field. Mobility models with geographical restrictions have been proposed to describe this kind of environments: in the Pathway Mobility Model [START_REF] Tian | Graph-based mobility model for mobile ad hoc network simulation[END_REF] the node movement is restricted to pathways in a map, while in the Obstacle Mobility Model [START_REF] Johansson | Scenariobased performance analysis of routing protocols for mobile ad-hoc networks[END_REF] the effects of obstacles on node mobility and radio propagation are taken into account. All these classes of models are unable to constrain/define the movement dynamics as a function of the mobile node positioning.

There has also been a significant interest in mobility models specially designed for partitioned networks, where end-to-end connectivity rarely or never exists. In [START_REF] Bittner | The area graph-based mobility model and its impact on data dissemination[END_REF], authors proposed the Area Graph Based Mobility Model that considers the inhomogeneity of spatial node distribution by using a directed graph with areas of different densities; in [START_REF] Hsu | Weighted waypoint mobility model and its impact on ad hoc networks[END_REF] similar approaches, called Weighted Waypoint Mobility and Clustered Mobility Models, respectively, were proposed to capture node preferences in choosing possible destinations; the Community Based Mobility Model, proposed in [START_REF] Spyropoulos | Performance analysis of mobilityassisted routing[END_REF], was designed to recreate the characteristics of mobility in sparse networks. All these models exhibit a skewed distribution of location visits, which could result in forming islands of connectivity at popular locations and being unable to stabilize. The Heterogeneous Random Walk model, proposed in [START_REF] Piorkowski | A Parsimonious Model of Mobile Partitioned Networks with Clustering[END_REF], is a very parsimonious model where clusters emerge because of low node speed in popular areas. Musolesi et al. [START_REF] Musolesi | Designing mobility models based on social network theory[END_REF] presented a mobility model that incorporates sociological relationships between mobile nodes, resulting in a quite complex model. Recently, some studies have been devoted to model the realistic social dimension of human mobility. A human mobility model based on heterogeneous centrality and overlapping community structure in social networks was proposed in [START_REF] Yang | Using social network theory for modeling human mobility[END_REF], being able to match characteristics observed from real human motion traces, especially heterogeneous human mobility popularity. This type of models requires pre-established division of the network geographic environment and are not able to adapt themselves to unknown generic scenarios.

Several approaches have been proposed to predict the future location of mobile entities. In [START_REF] Liu | Mobility modeling, location tracking and trajectory prediction in wireless networks[END_REF], authors proposed a method recording a sequence of cell identifiers to predict future location of users in wireless ATM networks, while [START_REF] Akyildiz | The predictive user mobility profile framework for wireless multimedia networks[END_REF] use logical information for location prediction. Other approaches use road topology information to predict user location [START_REF] Soh | Dynamic bandwidth reservation in cellular networks using road topology based mobility predictions[END_REF]. In [START_REF] Lee | Modeling steady-state and transient behaviors of user mobility: formulation, analysis, and application[END_REF], a semi-Markov model is used to describe the state transitions between stationary behavior and movement in a WLAN cell, where all model parameters can be obtained from movement monitoring. Most of the existing methods are restricted to a specific purpose, leading to straightforward implementation but limited effectiveness in generic network management tasks. Designing a robust flexible mobility and location model that can be used in generic environments is the motivation for this paper.

Model and Positioning Prediction

Markov Modulated Bi-variate Gaussian Processes Model

The proposed discrete time Markov Modulated Bi-variate Gaussian Process (dMMBVGP) model characterizes position and mobility of a subject based on the following assumptions: (i) within a restricted generic geographic region the position may be described by a bi-variate Gaussian distribution, (ii) the geographical space can be divided in N distinct regions (non-determined a priori ) and (iii) the transition between regions can be described by an underlying (homogeneous) Markov chain, where each state maps the movement characteristics of a specific region. The dMMBVGP can then be described as a random process (Z) with a bi-variate Gaussian distribution, which determines a position in a two-dimensional environment, whose parameters are a function of the state (J) of the modulator Markov chain (Z, J) with N states. More precisely, the (homogeneous) Markov chain (Z, J) = {(Z k , J k ), k = 0, 1, . . .} with state space IR 2 × U , with U = {1, 2, . . . , N }, is a dMMBVGP if and only if for k = 0, 1, . . .,

P (Z k+1 = z, J k+1 = n|J k = m) = p mn Γ n (z) (1)
where p mn represents the probability of a transition from state m to state n of the underlying Markov chain in time interval [k, k + 1], and Γ n (z) =

1 2π|Σn| 1/2 e -1 2 (z-Λn) T Σ -1
n (z-Λn) is the bi-variate Gaussian distribution in region n centered in Λ n and having covariance matrix Σ n .

The dMMBVGP may be also regarded as a Markov random walk [START_REF] Pacheco | Markov-modulated processes & semiregenerative phenomena[END_REF] where the two-dimensional positions increments (∆) in each instant k have a distribution given by the difference of two bi-variate Gaussian distributions whose parameters are a function of the state of the modulator Markov chain

P (Z k+1 =z + ∆, J k+1 = n|Z k = z, J k = m) = = p mn 1 2π |Σ n + Σ m | 1/2 e -1 2 (∆-Λn+Λm) T (Σn+Σm) -1 (∆-Λn+Λm)
Whenever ( 1) holds, we say that (Z, J) is a dMMBVGP with a set of modulating states with size N and parameter matrices P, S, and L. The matrix P is the transition probability matrix of the modulating Markov chain J,

P =     p 11 p 12 . . . p 1N p 21 p 22 . . . p 2N . . . . . . . . . . . . p N 1 p N 2 . . . p N N     (2)
matrix L defines the centers of each region,

L = Λ 1 Λ 2 . . . Λ N , Λ n = λ (x) n λ (y) n (3)
and matrix S contains the covariance (sub-)matrices of each region:

S = Σ 1 Σ 2 . . . Σ N , Σ n = (σ (x) n ) 2 ρ n σ (x) n σ (y) n ρ n σ (x) n σ (y) n (σ (y) n ) 2 . ( 4
)
Moreover, we denote by Π = [π 1 , π 2 , . . . , π N ] the stationary distribution of the underlying Markov chain. Specifying z = x y T , where x and y represent longitude and latitude, respectively, and defining the general (not weighted) bi-variate Gaussian distribution within a region n by

Γ n (x, y) = 1 2πσ (x) n σ (y) n 1 -ρ 2 n e Qn (5)
with

Qn = -1 2(1-ρ 2 n )   (x-λ (x) n ) 2 (σ (x) n ) 2 + (y-λ (y) n ) 2 (σ (y) n ) 2 - 2ρn (x-λ (x) n )(y-λ (x) n ) σ (x) n σ (y) n  
, the stationary probability density map of position (x, y) can be determined by a weighted sum of N bi-variate Gaussian distributions:

f N (x, y) = N n=1 π n Γ n (x, y) (6)
obtaining a position probability density map.

Positioning prediction

Defining q k = {q k,1 , q k,2 , . . . , q k,N }, k = 0, 1, . . . , K, where q k,n is the probability that, at time interval k (within a total of K observations), the mobile node is in region n (i.e. the probability of the underlying Markov chain of the mobility model is in state n at time instant k) and based on equation ( 6) we can define the predicted position probability map at time instant k as:

fk (x, y) = N n=1 q k,n Γ n (x, y) (7)
with

q k = q 0 P k (8)
where q 0 represents the probability vector of the last known position (x 0 , y 0 ) belonging to the different regions, and can be inferred as

q 0,n = π n Γ n (x 0 , y 0 ) f (x 0 , y 0 ) (9)

dMMBVGP Inference Procedure

Hereafter, we will refer to the empirical position probability density map as f e (x, y) and n-1 (x, y) = f e (x, y)-f n (x, y) will represent the difference function between the empirical and fitted position probability maps with n bi-variate Gaussian distributions (f n ).

The inference procedure can be divided in two steps: (i) inference of matrices L and S and vector Π by adjusting f e (x, y) by a weighted sum of bi-variate Gaussian distributions and (ii) inference of matrix P by analyzing the transitions between the underlying Markov chain inferred in step (i).

The first step of the inference procedure relies on the minimization process defined by min |f n (x, y) -f e (x, y)| where f n (x, y), as defined in [START_REF] Hsu | Weighted waypoint mobility model and its impact on ad hoc networks[END_REF], is a weighted sum of n bi-variate Gaussian distributions and n is determined as part of the inference process and does not need to be defined a priori.

The iterative fitting procedure starts with n = 1, finds the global maximum of 0 (x, y) = f e (x, y) and defines that point as the center of the first mobility region. Each mobility region (n) will have associated a bi-variate Gaussian distribution with parameters λ

(x) n , λ (y)
n , σ

(x) n , σ (y)
n , ρ n and weight π n . The n-th mobility region center (i.e. the bi-variate Gaussian distribution mean) is determined by:

(λ (x) n , λ (y) n ) = -1 n-1 (max( n-1 (x, y))) , n ≥ 1. ( 10
)
The inference of the remaining parameters of the bi-variate Gaussian distributions are estimated by analyzing the bi-dimensional local decay behavior of f e (x, y) in the vicinity of the local maximum points along an orthogonal plane. Therefore, based on the bi-variate Gaussian distribution function [START_REF] Camp | A survey of mobility models for ad hoc network research[END_REF], it is possible to determine the ratio between the partial derivative of π n Γ n (x, y) with respect to the variable x in the plane defined by y = λ (y)

n and π n Γ n (x, λ (y)
n ) as

πn∂Γn(x,λ (y) n ) ∂x π n Γ n (x, λ (y) n ) = - (x -λ (x) n ) (1 -ρ 2 n )(σ (x) n ) 2 ⇔ (σ (x) n ) 2 = γ (x) n 1 -ρ 2 n ( 11
)
and the ratio between the partial derivative of π n Γ n (x, y) with respect to the variable y in the plane defined by x = λ

(x)
n and π n Γ n (λ

(x) n , y) as πn∂Γn(λ (x) n ,y) ∂y π n Γ n (λ (x) n , y) = - (y -λ (y) n ) (1 -ρ 2 n )(σ (y) n ) 2 ⇔ (σ (y) n ) 2 = γ (y) n 1 -ρ 2 n ( 12
)
where

γ (x) n = - (x -λ (x) n )Γ n (x, λ (y) n ) ∂Γn(x,λ (y) n ) ∂x and γ (y) n = - (y -λ (y) n )Γ n (λ (x) n , y) ∂Γn(λ (x)
n ,y) ∂y

. [START_REF] Musolesi | Designing mobility models based on social network theory[END_REF] Based on definitions ( 5) and ( 13) it is possible to redefine the n-th weighted bi-variate Gaussian distribution function

Γ n (x, y) = Γ n (λ (x) n , λ (y) n )e Q n ( 14
)
with Γ n (λ

(x) n , λ (y) n ) = π n √ 1-ρ 2 n 2π √ γ (x) n γ (y) n
and

Q n = -1 2 (x-λ (x) n ) 2 γ (x) n + (y-λ (y) n ) 2 γ (y) n - 2ρn(x-λ (x) n )(y-λ (x) n ) √ γ (x) n γ (y) n .
¿From equation [START_REF] Pacheco | Markov-modulated processes & semiregenerative phenomena[END_REF], it is possible to define the n-th bi-variate Gaussian distribution correlation:

ρ n = γ (x) n γ (y) n (x -λ (x) n )(y -λ (y) n ) r n (x, y) (15)
where

r n (x, y) = ln Γn(x,y) Γn(λ (x) n ,λ (y) n ) + (x-λ (x) n ) 2 2γ (x) n + (y-λ (y) n ) 2 2γ (y) n
, and the n-th bi-variate Gaussian distribution weight (i.e. the stationary probability of the underlying Markov chain n-th state) is given by:

π n = 2πΓ n (λ (x) n , λ (y) n ) γ (x) n γ (y) n 1 -ρ 2 n ( 16
)
It is assumed that within a square region with width 2∆ and center at the n-th local maximum of f e (x, y) (equation ( 10)), the position is modeled exclusively by the n-th weighted bi-variate Gaussian distribution, i.e.,

Γ n (x, y) ≈ n-1 (x, y), x ∈ [λ (x) n -∆ (x) , λ (x) n +∆ (x) ]∧y ∈ [λ (y) n -∆ (y) , λ (y)
n +∆ (y) ]. Therefore, it is possible to find estimators for parameters γ (x) n and γ (y) n based on the average value of equations ( 13) calculated in the vicinities of the f e (x, y) local maximums:

γ(x) n = - 1 2∆ (x) λ (x) n +∆ (x) λ (x) n -∆ (x) (x -λ (x) n ) n-1 (x, λ (y) n ) ∂ n-1(x,λ (y) n ) ∂x dx ( 17
)
and

γ(y) n = - 1 2∆ (y) λ (y) n +∆ (y) λ (y) n -∆ (y) (y -λ (y) n ) n-1 (λ (x)
n , y)

∂ n-1(λ (x)
n ,y) ∂y dy.

(

) 18
The estimator for the n-th bi-variate Gaussian distribution correlation is defined as the average value of ( 15) calculated over a square region in the vicinity of the f e (x, y) local maximums:

ρn = γ (x) n γ(y) n 4∆ (x) ∆ (y) λ (x) n +∆ (x) λ (x) n -∆ (x) λ (y) n +∆ (y) λ (y) n -∆ (y) rn (x, y) (x -λ (x) n )(y -λ (y) n ) dydx (19)
where rn (x, y) = ln n-1 (x,y)

n-1(λ (x) n ,λ (y) n
) + (x-λ (x) n ) 2 2γ (x) n + (y-λ (y) n ) 2 2γ (y) n .
The weight of the n-th bi-variate Gaussian distribution can be calculated using the estimators obtained with equations ( 17), ( 18) and ( 19):

πn = 2π n-1 (λ (x) n , λ (y) n ) γ(x) n γ(y) n 1 -ρ2 n (20)
The n-th bi-variate Gaussian distribution variances can be obtained based on equations ( 11), ( 12), ( 17), ( 18) and ( 19):

(σ (x) n ) 2 = γ(x) n (1 -ρ2 n ) and (σ (y) n ) 2 = γ(y) n (1 -ρ2 n ) (21)
The final iteration step is to re-calculate the difference function and update n = n + 1:

n (x, y) = n-1 (x, y) -πn Γ n (x, y). (22)
The iterative process finishes after a predefined number of iterations or when min |f n (x, y) -f e (x, y)| reaches a predefined fitting error; if that goal is not achieved, the value n is increased by one and the iterative process continues recomputing equations [START_REF] Lee | Modeling steady-state and transient behaviors of user mobility: formulation, analysis, and application[END_REF] to [START_REF] Tian | Graph-based mobility model for mobile ad hoc network simulation[END_REF]. At the end of the iteration process it is possible to determine the model matrices L and S.

The final step of the inference procedure is to infer matrix P, i.e. the transition probabilities between the regions/states defined in the first step. The task is achieved by probabilistically mapping each known position (x k , y k ), k = 0, 1, . . . , K to one region and then average the probabilistic transitions between states.

The position in time interval k (x k , y k ) will be probabilistically mapped to one of the N regions according to a probability vector q k = {q k,1 , q k,2 , . . . , q k,N }, k = 0, 1, . . . , K defined by:

q k,n = π n Γ n (x k , y k ) f (x k , y k ) (23)
Since the underlying Markov chain at instant k is in state m with probability P (J k = m), it is possible to define a partial estimator for p mn at all k time instants as (taking equations 1 and 2 into account):

pmn,k = P (J k = m)p mn = P (J k = m)P (J k+1 = n|J k = m) = = P (J k = m)P (J k+1 = n) = q k,m q k+1,n (24)
Therefore, the overall estimator for p mn can then be inferred by averaging over k the partial estimators defined in (24),

pmn = 1 K -1 K-1 k=0 pmn,k , m, n = 1, . . . , N (25)
Rewriting (25) in matrix form, it is possible to obtain matrix P:

P = 1 K -1 K-1 k=0 q T k × q k+1 (26)

Results

As proof of concept, we applied the proposed model to mobility data of taxicabs in San Francisco, USA. Data was provided by Exploratorium, the San Francisco museum of science, art and human perception, through the cabspotting project (http://cabspotting.org). A pre-processed data set is freely available [START_REF] Piorkowski | CRAW-DAD data set epfl/mobility (v. 2009-02-24[END_REF] and contains GPS coordinates, sampled at approximately 1 per minute, of 500 taxis collected over 30 days of May 2008 in the San Francisco Bay Area. Fig. 1 depicts a sample trajectory of one taxi in the San Francisco Bay area. In this paper, in order to infer and test the proposed model we have only considered the movements within a rectangular region in the San Francisco downtown, which is depicted in Fig. 1 within the white rectangle, with longitude limits between -122.447 o and -122.388 o and latitude limits between 37.77 o and 37.81 o . The data subset was randomly divided in two sub-sets (of 250 cars each), one to infer the model parameters and the other to test the model. Applying the proposed mobility model to the test data sub-set, a model with 64 states (sub-regions) was obtained. The model inference time interval was less than 30 minutes with an octave/Matlab implementation in a generic PC. Fig. 2 shows the position probability distribution map for a taxi in San Francisco downtown and Fig. 3 shows the stationary position probability distribution map of the dMMBVGP inferred model. It is possible to observe that the dMMBVGP model was able to incorporate and model the complexity of the cars movement/location in San Francisco Downtown. With a more detailed analysis, it is possible to observe that the inference algorithm centered the sub-regions in the more visited/crowded San Francisco's landmarks.

In order to evaluate the location prediction performance of the model, we used the test sub-set to predict the node future position based on a known position, using different time-lags (1 to 20 timeslots, approximately 1 to 20 minutes). The performance metric was the (average) distance (in meters) between the predicted and real future node position. After inferring the model parameters, position prediction using the dMMBVGP, equation [START_REF] Hsu | Modeling spatial and temporal dependencies of user mobility in wireless mobile networks[END_REF], is performed in hundreds of microseconds (average value ∼ 921µs), using also an octave/Matlab implementation in a generic PC. Note that the prediction time interval depends on the number of future time slots for which the prediction is made.

As a complementary evaluation of the proposed model performance, we also evaluate other four mobility and location prediction models in the same scenario:

-Complete random model, where node location follows a uniform distribution over the complete geographic region under study. -Grid random model, where node location is probabilistically determined in two phases; firstly, the location sub-region is determined according to the empirical location distribution map inferred over a uniform grid and, secondly, the exact position within the chosen sub-region is determined according to a uniform distribution. -Random walk model [START_REF] Camp | A survey of mobility models for ad hoc network research[END_REF] with reflection at the region borders [START_REF] Papoulis | Probability, random variables, and stochastic processes[END_REF], where the node direction of movement is uniformly distributed in [0, 2π] and the speed has a Gaussian distribution with mean and standard deviation inferred from empirical data. -Markov grid model, which is a sub-region based Markov family model [START_REF] Song | Evaluating location predictors with extensive wi-fi mobility data[END_REF] where location is determined by an underlying Markov chain that determines the transition between sub-regions. We choose to define the sub-regions as an uniform grid. Within a sub-region, node location has an uniform distribution.

-
The complete random model has no geographic dependency and has no information of the node movement history. The model location prediction does not depend on the current node position. This model is used as the worst case reference value. The grid random model prediction does not depend also on the Being able to model movement and predict future positioning of mobile nodes in complex environments can be very important to several operational and management tasks, like evaluation of different system design alternatives, simulation and study of new routing protocols, message delivery and nodes connectivity optimization or mobile computing. Since existing mobility modeling approaches are based on simple and limited models that were specifically designed for particular application scenarios or require the complete knowledge of the mobility environment, this paper presented a discrete time Markov Modulated Bi-variate Gaussian Process that is able to characterize the position and mobility of any mobile node by assuming that the position within a generic sub-region can be described by a bi-variate Gaussian distribution and the transition between subregions can be described by an underlying (homogeneous) Markov chain. This approach can be applied to scenarios where the possible pathways are unknown or too complex to consider in a real model that must make a prediction in a very short time. The results obtained by applying the proposed model to real and publicly available mobility data of San Francisco taxi-cabs demonstrate the accuracy and utility of this approach: the model was able to efficiently describe the movement patterns of the mobile nodes and predict their future position. Besides, the model has also revealed higher performances when compared to other modeling approaches.

Fig. 1 .

 1 Fig. 1. San Francisco map with one sample taxi movement (base map source: GoogleMaps).

  Real position probability map, San Francisco downtown only.

	probability
	37.809						-122.389
	37.803						-122.399
	37.794						-122.411
	latitude	37.786	37.778		-122.435	-122.423	longitude
			37.770	122.447
	37.770 Fig. 2. -122.447 37.778 37.786 37.794 37.803 37.809 latitude probability	-122.435	-122.423	-122.411 longitude -122.399	-122.389
	Fig. 3. Fitted dMMBVGP position probability map, San Francisco downtown only.

current node position, however it incorporates some geographic dependency inferred from empirical historic data. A 8x8 grid (64 equal sub-regions) was used to infer the empirical location distribution map. The random walk model has no geographic dependency but the prediction process utilizes the current node position and extrapolates the future position based on a known constant speed and direction. The model speed was modeled by a Gaussian distribution with mean 302.1m/minute (≈ 18.1Km/h or ≈ 11.2mph) and standard deviation of 346.7m/minute. For the inference of the Markov grid model, a 8x8 grid (64 equal sub-regions) was used in order to obtain a model with the same level of complexity as the dMMBVGP. Note that the Markov grid model differs from the dMMBGP model in the way regions are defined and position is modeled within each region; the time dependability (based on an underlying Markov chain) is modeled and inferred in a similar way in both models. Figure 4 depicts the models prediction performance. It is possible to observe that the dMMBGP outperforms all the other models and it is possible to obtain relatively small prediction errors, in a region so large and complex as the one under study, even in long range predictions (20 timeslots ∼ 20 minutes). As expected, the worst performing models are the ones that do not receive any real-time input at the time of the prediction, and therefore the prediction error is independent of the time-lag. The random walk model is able to surpass the others in very short time-lag predictions (≤ 2 minutes); however the model is not able to capture the complex dynamics of the taxis movements and the prediction error greatly increases when the prediction range is larger than 1 slot and rapidly approximates the error of a completely random location choice. The Markov grid model, due to its nature, is able to capture the complexity and dynamism of the movements; however, the definition of generic regions restricts the model performance. Therefore, the dMMBVGP model is able to model the complexity of the movements due to its Markovian nature, while its adaptable definition of sub-regions allows it to maintain the model complexity at low levels and achieve very good location prediction results.


Nicholas C Valler
email: nvaller@cs.ucr.edu

B Aditya Prakash

Hanghang Tong
email: htong@us.ibm.com

Michalis Faloutsos
email: michalis@cs.ucr.edu

Christos Faloutsos
email: christos@cs.cmu.edu

Epidemic Spread in Mobile Ad Hoc Networks: Determining the Tipping Point

Short-range, point-to-point communications for mobile users enjoy increasing popularity, particularly with the rise in Bluetooth-equipped mobile devices. Unfortunately, virus writers have begun exploiting lax security in many mobile devices and subsequently developed malware exploiting proximity-based propagation mechanisms (e.g. Cabir or CommWarrior). So, if given an ad-hoc network of such mobile users, will a proximityspreading virus survive or die out; that is, can we determine the "tipping point" between survival and die out? What effect does the average user velocity have on such spread? We answer the initial questions and more. Our contributions in this paper are: (a) we present a framework for analyzing epidemic spreading processes on mobile ad hoc networks, (b) using our framework, we are the first to derive the epidemic threshold for any mobility model under the SIS model, and (c) we show that the node velocity in mobility models does not affect the epidemic threshold. Additionally, we introduce a periodic mobility model and provide evaluation via our framework. We validate our theoretical predictions using a combination of simulated and synthetic mobility data, showing ultimately, our predictions accurately estimate the epidemic threshold of such systems.

Introduction

The prevalence and increased functionality of mobile phones present an unique opportunity for malicious software writers. According to popular reports, 45.4 million U.S. This material is based upon work supported by the Army Research Laboratory under Cooperative Agreement No. W911NF-09-2-0053, the National Science Foundation under Grants No. IIS-1017415, CNS-0721736 and CNS-0721889, NETS-0721889 and NECO-0832069 and a Sprint gift. Any opinions, findings, and conclusions or recommendations in this material are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the Army Research Laboratory, the U.S. Government, the National Science Foundation, or other funding parties. The U.S. Government is authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation here on. citizens own smartphones in 2010 [START_REF] Gonsalves | Android Phones Steal Market Share[END_REF], roughly 15% of the total population. These devices are becoming more and more essential for the daily lives of end-users, especially as smartphones offer more and more capabilities. At the same time, mobile devices are equipped with device-to-device (a.k.a. point-to-point) communication technologies (such as Bluetooth), where the communication does not use the phone's service provider's infrastructure. There are estimated over 920 million Bluetooth-equipped devices worldwide shipped in 2008, making Bluetooth the most common point-to-point communication protocol in today's smartphones. To date, at least two smartphone worms have been found in the wild, Cabir and CommWarrior, which spread using Bluetooth.

Given a system of mobile agents, such as smartphones, what can we say about the propagation of a virus4 within the system? A key question is to identify the tipping point, known as the system's epidemic threshold, or take-off point, below which a virus is guaranteed to "die out." For the epidemic models, we focus on one of the most popular models, the flu-like one susceptible-infected-susceptible -SIS (see Section 5 where we handle other models). There, agents maintain no immunity, and become susceptible, immediately after they heal. Our key contributions are as follows:

1. Framework and Formula: We present a framework for estimating the epidemic threshold on any, arbitrary mobile ad hoc network model. To the best of our knowledge, ours is the first theoretical study for the general mobility case. The idea is to derive a sequence of adjacency matrices, and then compute the first eigenvalue of the so-called system matrix (see Theorem 1). There, we show that the epidemic threshold depends only on this first eigenvalue, and nothing else. 2. Closed Formulas: We show how to use our framework to derive simple, approximate (but accurate) formulas for several, popular special cases (Random Walk model, Levy flight model). 3. Insensitivity to Velocity: To our surprise, our results show that the epidemic threshold does not depend on the node velocity (v > 0). Our experiments confirm the accuracy of our approximations, as well as our 'insensitivity' observation.

Jumping ahead, Figure 1 showcases the accuracy of our results (Lemma 3) of the epidemic threshold for the SIS (=flu-like) model on the so-called 'Levy Flight' mobility model. See section 4 for more details -but the point is that our prediction for the take-off point (= epidemic threshold, indicated with a black arrow) is exactly where all curves take off.

We have two additional contributions: through extensive simulation experiments, we show that similar insensitivity results hold for other popular mobility models like Levy flight, Random Waypoint etc.; and, moreover, we introduce the periodic mobility model, which is very popular in biological virus epidemiology [START_REF] Anderson | Infectious Diseases of Humans[END_REF][START_REF] Barrett | EpiSimdemics: an Efficient Algorithm for Simulating the Spread of Infectious Disease over Large Realistic Social Networks[END_REF], and show how to use our framework to estimate its epidemic threshold.

The rest of this paper has the typical organization: background (Section 2, proposed framework and theorems (Section 3), experiments (Section 4), additional observations (Section 5), and finally, related work and conclusions. (Sections 6 and 7). Fig. 1. Accuracy of our results (Lemma 3), for the 'Levy-Flight' model. Take-off plot, plotting the max number of infections vs. strength of the virus. Notice that our predicted take-off point (black arrow) agrees with the simulations, for several node velocities.

Preliminaries

In this section, we present a general background on proximity-based epidemic spreading models and formulate our problem statement.

Epidemic Model: SIS (flu-like)

The SIS epidemic model resembles a flu-like virus, where nodes have no immunity. Healthy ('S' = susceptible) nodes become sick ('I' = infected) stochastically from their infected neighbors with a probability β. Alternatively, a sick node becomes healthy (and open to re-infection), with a probability δ. These two parameters are also referred to as the birth rate (β) and death rate (δ) of the virus.

The tipping point τ , or epidemic threshold, of an SIS epidemic model is the condition under which an infection will die out exponentially quickly irrespective of initial infection, as opposed to spreading out, causing and epidemic (technically, a pandemic). For a survey on SIS and numerous other epidemic models, see Hethcote [START_REF] Hethcote | The mathematics of infectious diseases[END_REF], or [START_REF] Chakrabarti | Epidemic thresholds in real networks[END_REF][START_REF] Ganesh | The effect of network topology in spread of epidemics[END_REF].

Problem Formulation

Using the background discussed above, we now formulate our problem statement. See Table 1 for definitions of various symbols. In this paper, we consider an epidemic on a mobile network, which provides an underlying contact structure for the virus to use as it propagates. By doing so, at any point in time, the system is non-homogenous, as nodes may only transmit the virus to its neighbors. Given:

1. Mobile ad hoc network mobility models, described below.

2. The SIS model parameters, i.e., the virus birth and death probabilities β and δ.

Find:

The epidemic threshold τ or tipping point for the system such that for τ < 1 an infection will die out quickly, irrespective of initial conditions. Our problem naturally leads to other issues like the effect of node velocity in models on the threshold, giving approximations in specific cases etc. We elaborate on them in the upcoming sections.

Framework

In this section, we detail our framework for analyzing mobility models and then move on to specific approximations and questions arising out of the framework. We will present extensive simulations demonstrating the results later in Section 4. Also please see Section 4.2 for a description of the mobility models. Note that node-to-node contacts at a particular time can be represented by an adjacency matrix A. We next provide a general theorem expressing the epidemic threshold for mobility models.

Epidemic Thresholds on Mobility Models

Theorem 1 (Mobility model threshold). If a mobility model can be represented as a sequence of connectivity graphs L = {A1, A2, . . . , AT }, one adjacency matrix At for each time step t ∈ {1..T }, then the epidemic threshold is:

τ = λS (1)
where λ S is the first eigenvalue of matrix S and S = i Si and ∀i ∈ {1..T } Si = (1 -δ)I + βAi (I is the standard N × N identity matrix).

Proof. If the mobility model can be represented as a sequence of graphs, then the epidemic threshold depends on the first eigenvalue of the system matrix [START_REF] Prakash | Virus Propagation on Time-Varying Networks: Theory and Immunization Algorithms[END_REF]. Hence,

τ = λ i ((1-δ)I+βA i ) .
We can now give a simpler closed-form approximation for the threshold in Equation 1 in the following lemma: Lemma 1 (EAAM Approximation for Threshold). Under the same conditions as in Theorem 1, the following is an approximation for the epidemic threshold:

τ ≈ β δ × λ Aavg ( 2
)
where Aavg = i Ai/T is the average adjacency matrix.

Proof. Note that,

S = i ((1 -δ)I + βAi) = (1 -δ) T I + β i Ai + O(β 2 ) + O(β * δ) + O(δ 2 ) ≈ (1 -T δ)I + T βAavg (3)
where we neglected second or lower order terms involving β and δ. Hence, we find that B = (1 -T δ)I + T βAavg is a first order approximation for the S = i Si matrix.

Hence from Theorem 1 we want λB < 1 which implies Equation 2.

We will refer to the above approximation as the 'Eigenvalue of the Average Adjacency Matrix' (EAAM) approximation.

Specific Approximations

Lemma 2 (Random-Walk Threshold). In the random-walk mobility model and under the SIS model, the following is an approximate epidemic threshold:

τ ≈ β/δ × πR 2 × N/A ( 4
)
where R is the radius of influence of each node.

Proof. Under the random-walk model, at every point of time, each node is at a random (x, y) position, uniformly distributed on the field of interest. Each node has a radius of possible connections (like the BlueTooth radius) R. Consequently each node has d = πR 2 × N/A neighbors on average (ignoring boundary effects). The connectivity graph at each time step is roughly a random graph with average degree d. Hence it has first eigenvalue λ1 = d on average. Hence this is approximately equivalent to having a static graph under the SIS model where the epidemic threshold [START_REF] Chakrabarti | Epidemic thresholds in real networks[END_REF] is τ = β/δ × λ1.

We now obtain the lemma after obvious substitutions.

In fact, we can go further and generalize this to any mobility model where the geographic steady state distribution is uniform.

Lemma 3 (Uniformly-Distributed Steady State (USS) Threshold). For any mobility model where the geographic distribution of nodes at the steady state is uniform over the area of interest and under the SIS model, the following is an approximate epidemic threshold:

τ ≈ β/δ × πR 2 × N/A ( 5
)
where R is the radius of influence of each node.

Proof. The proof for Lemma 2 goes through even here precisely because of the geographically uniformly distributed nature of the steady state. Each node has roughly the same number of connections and hence the adjacency graph is approximately a homogenous graph with constant first eigenvalue. The result follows as before.

We will refer to the above approximation as the 'Uniform Steady-State' (USS) approximation. Mobility models like Levy Flight and Random-Walk are examples of models with a geographically uniform distribution of the nodes at the steady state. Lemma 3 allows us to quickly estimate the threshold for these and many other models.

Insensitivity to Node Velocity

As there is no factor depending on the node velocity in Lemma 3, we conclude the following surprising implication: Corollary 1 (Node velocity and threshold). The node velocity (v > 0) does not affect the epidemic threshold in mobility models where the steady state has a geographically uniform steady state distribution like Random Walk, Levy Flight etc.

We conjecture that the velocity does not affect the threshold even for non-geographically uniformly distributed steady state mobility models like Random Waypoint. We provide empirical results supporting this claim later in Section 4.

Conjecture 1 (Effect of velocity). The node velocity (v > 0) does not affect the epidemic threshold in the Random Waypoint mobility model. Does Velocity have an impact at all? The above discussion raises the point whether the node velocity has any effect at all on the dynamics of the epidemic spreading. We expect that the velocity of motion does have an effect, when we are above threshold. Furthermore, simulations resulted in a non-obvious observation. The velocity had an impact on the steady-state number of infected agents in the system. We elaborate more on these issues in Section 5.

Simulation Methodology and Results

Experimental Setup

To facilitate the simulation, we wrote a custom Python2.6 simulation program using the NumPy/SciPy python libraries. All simulations were conducted on a 4 core Intel(R) Xeon(R) CPU operating at 2.53 GHz and 72 GB of memory running CentOS-5.5 (Linux kernel 2.6). We varied the number of agents (nodes) N between 250 and 1500 within a simulation field of area A = 40, 000m 2 (200m by 200m). Thus, node density ρ, commonly defined a N/A, was between 0.125 and .125 nodes per m 2 . All nodes had a transmission range of 5.0 meters. We did not account for signal attenuation, reflection nor other wireless phenomena. Prior to the beginning of the simulation, nodes were distributed on the simulation field in a uniform fashion. Simulations were generally run for a period of 100s with time intervals of ∆T = 0.1 seconds.

We studied three mobility models common to mobile ad hoc networking: Random Walk, Levy Flight and Random Waypoint. In the following sections, we provide detail on each model as well as simulation results. The position Pi,t of each node in the system at time t is a function of mobility model and previous position and time step ∆T , such that Pi,t+1 = M (Pi,t, ∆T ), where M is the mobility model.

The purpose of our simulations was to determine the role of the mobility model in the spread of malware in a point-to-point contact network loosely describing Bluetooth communication technology.

Mobility Models

Random Walk. The Random Walk (RW) mobility model (also referred to as Brownian Motion) was originally formulated to describe the seemingly random motion of particles. Numerous variations exists, here we describe our implementation.

Each node i in the system is parameterized by speed (Vi) and angle (θi). Both Vi and θi are drawn uniformly from systemwide predefined ranges, [vmin, vmax] and [0, 2π), respectively. Clearly, such a system is memoryless. The model we employ varies from the simple RW model by introducing a flight time for each node, Ti. Flight time is drawn uniformly from a range [τmin, τmax]. The spatial distribution of the RW mobility model is uniform over the simulation field. According to our framework, we predict the RW mobility model will follow Lemma 2.

Levy Flight. Levy Flight mobility models have recently attracted attention due to their statistical similarities with human mobility [START_REF] Rhee | On the Levy-Walk Nature of Human Mobility[END_REF]. At the beginning of each flight, each node selects an angle uniformly from within (0, 2 * π], a flight time drawn from some distribution, a flight length and a pause time. Flight length and pause time are drawn from Levy distributions p(t) ∝ |t| -(1+α) and ψ(t) ∝ t -(1+β) , where time t > 0, respectively. When α = 2 and β = 2 the result is a special case of the Levy distribution resulting in a Gaussian distribution. As with the Random Walk, the spatial distribution of the Levy Flight mobility model is uniform over the simulation field.

Random Waypoint. The Random Waypoint (RWP) mobility model is often cited as the de facto mobility model in ad hoc networks. As originally proposed by Johnson et. al [START_REF] Johnson | Dynamic Source Routing in Ad Hoc Wireless Networks[END_REF], the RWP mobility model each node i is described by three parameters: current location (Pi), speed (Vi), waypoint (Wi), and pause time (ρ). In general, the RWP mobility model operates as follows: Initially, a node is stationary. After a pause time ρ, the node selects a waypoint uniformly from the simulation field, then, travels along the shortest path to its waypoint P at a velocity Vi drawn uniformly from (vmin, vmax]. Upon arrival at their waypoint, each node will pause for a time t = ρ. After the pause period is done, each node will repeat the process. The spatial distribution of the RWP mobility model is bell-shaped [START_REF] Bettstetter | The Node Distribution of the Random Waypoint Mobility Model for Wireless Ad Hoc Networks[END_REF].

Summary of Results

Accuracy of Approximations, Random Walk. We present a series of simulation studies of the Random Walk mobility model in Figure 2. In these studies, we varied the birth (β) and death (α) parameters of the SIS infection model. We refer to the resulting plots as "take-off plots," which show the maximum number of infected agents seen in our simulation against our approximated epidemic threshold. For each plot, we labeled the estimated take-off point according to the specific threshold approximation. Figure 2(A), we approximated the epidemic threshold using the first eigenvalue of the average adjacency matrix, Aavg (the EAAM threshold approximation of Lemma 1). We indicate the predicted threshold value at τ = 1. As expected, no epidemic was present at values of τ < 1. At values of τ > 1, we see explosive growth in the max number of infected agents.

In Figure 2(B), we plot the the USS approximation of the epidemic threshold (Lemma 3), i.e. β/δ ×πR 2 ×N/A. As in the EAAM threshold, USS behaves as expected (i.e. no epidemic below τ = 1). In fact, in each of these figures, we note that no infected agents (aside from patient zero) are present for epidemic threshold values below 1. Compared to EAAM, this threshold value takes off at values closer to τ = 1.

Accuracy of Approximations, Levy Flight. As with RW, the spatial distribution of nodes following Levy Flight mobility model is uniformly distributed on the simulation field. Thus, we expect Levy flight to perform similar to RW. Figure 3 presents a take-off plots for Levy Flight simulations. Again, we note that no infected agents exist below either threshold approximations. Furthermore, we find that USS performs better than EAAM.

Accuracy of Approximations, Random Waypoint. The next series of simulations were conducted using the popular Random Waypoint mobility model. We also selected this model specifically because it does not result in a uniform spatial distribution of nodes, therefore does not fall under USS. The results of these simulations are presented in Figure 4. Surprisingly, both threshold approximations perform well against the RWP mobility model. As the earlier mobility models exemplified, USS performs better than EAAM. Insensitivity to Velocity. In order to illustrate Corollary 1 and validate Conjecture 1, we performed a series of simulation in which we varied the node velocity, the results of which are illustrated in Figures 2,3,4.

As expected, for both the RW and Levy Flight mobility models, the take-off points were not greatly affected by increasing the nodes velocity. Furthermore, Figure 4 shows that the take-off point for the RWP mobility model was not affected by velocity, affirming Conjecture 1.

Discussion

We elaborate here on the effect of node velocity on the dynamics of epidemic spreading. We also introduce the periodic mobility model and present an analysis of it via our framework. In addition, we touch upon other epidemic models as well.

More on impact of node velocity

As discussed previously in Section 3, node velocity does not seem to effect the threshold in many models. We now ask whether the velocity affects the epidemic at all? For an "above threshold" system, two more parameters are of interest: (a) steadystate maximum, and (b) warm-up period. The steady-state maximum is the maximum number of infected agents in the system till steady state, whereas the warm-up period is the time necessary to reach steady state. We expect that the velocity of motion does have an effect, when we are above threshold. Clearly, speed will effect the speed of propagation of the virus and thus the warm-up period. Higher velocity means better mixing of agents, and thus faster convergence to the steady state. This observation is also demonstrated through simulations. velocity varied between a fixed 1 m/s fixed and an uniform selected 15 -20 m/s, as indicated in the legends. We performed a longer simulation in order for the systems to settle in a steady state. Less intuitive is that velocity appears to affect the number of infected agents at steady-state. For example, in Figure 5(A), the line corresponding to 1 m/s appears to reach a steady-state of approximately 65 -70% infected agents, whereas, at 15 -20 m/s, the steady-state is roughly 90%. The steady-state for velocities between these two extremes lay in-between. We suspect the degree of mixing, influenced by node velocity, is the root cause of the above observation.

The Periodic Mobility Model

As we indicated in Section 3, our framework predicts the epidemic threshold of mobility models that can be represented as a series of adjacency matrices. The Periodic Mobility Model is a special case of such a series, where a set of k adjacency matrices

{A1, A2, • • • , A k } are repeated periodically.
This is a typical model used in biological virus studies [START_REF] Barrett | EpiSimdemics: an Efficient Algorithm for Simulating the Spread of Infectious Disease over Large Realistic Social Networks[END_REF] to model general movements of a population. As an example, let A1 be an adjacency matrix of people during the day (say, at the office). Let A2 be an adjacency matrix representing contacts/interactions during the evening (say, at home). So the series formed by repeating {A1, A2} represents the daily, repeated interactions of our population. The periodic model offers a realistic, yet general model of mobility, capturing general patterns rather than specific movements of the system. Lemma 4 (Periodic Model Threshold). Under the periodic mobility model with k alternating behaviors repeating periodically, the epidemic threshold is given by τ = λS where S = k i=1 Si and, as before, Si = (1 -δ)I + βAi. Proof. Omitted for brevity, similar to Theorem 1.

As an example, Figure 6 shows the take-off plot of a periodic mobility models where N = 500 nodes, k = 2 adjacency matrices. As predicted by Lemma 4, the max number of infected agents over the simulation period takes off at λS = 1.

Other Epidemic Models

Given recent results on epidemic thresholds on static networks [START_REF] Prakash | Got the Flu (or Mumps)? Check the Eigenvalue![END_REF], we believe that our results will carry through for many other epidemic models as well e.g. SIR (mumpslike), SIRS, SEIR, MSEIR etc [START_REF] Hethcote | The mathematics of infectious diseases[END_REF] which capture differences between the way various diseases spread.

Conjecture 2 (Other Epidemic models). Our results for all the mobility models discussed in this paper hold for all the epidemic models covered in [START_REF] Prakash | Got the Flu (or Mumps)? Check the Eigenvalue![END_REF] as well.

Related Work

Here, we review the related work. It is worth pointing out that while most of existing studies about epidemic spread on mobile networks focus on (1) some particular types of network structures, and/or (2) one specific mobility model; our framework is very general and it applies to arbitrary network structure, and all the three popular mobility model. General Epidemic Modeling. Bailey provides the canonical text on epidemic modeling [START_REF] Bailey | The Mathematical Theory of Infectious Diseases and its Applications[END_REF]. A more recent survey is provided by Hethcote in [START_REF] Hethcote | The mathematics of infectious diseases[END_REF]. Kephart and White [START_REF] Kephart | Directed-Graph Epidemiological Models of Computer Viruses[END_REF][START_REF] Kephart | Measuring and Modeling Computer Virus Prevalence[END_REF] were among the first to propose epidemiology-based models to analyze the spread of computer viruses. The model they suggest provides a good approximation of virus propagation in networks where contact among individuals is essentially homogeneous. Recent discoveries suggest real networks (including social networks [START_REF] Domingos | Mining the Network Value of Customers[END_REF], router and AS networks [START_REF] Faloutsos | On power-law relationships of the Internet topology[END_REF], and Gnutella overlay graphs [START_REF] Ripeanu | Mapping The Gnutella Network: Properties of Large-Scale Peer-to-Peer Systems And Implications For System Design[END_REF]) follow a power-law structure instead, prompting a re-evaluation of the homogeneity assumption common in the works above.

Epidemics on Static Networks. Observation suggests that real networks are not homogeneous, rather, overwhelming evidence suggests real networks follow a power law structure instead. By introducing an underlying structure for a disease to spread, such Predicted Take-off Point as a static network, removes the original homogeneous assumption pioneered by those reference above. Newman [START_REF] Newman | Threshold Effects for Two Pathogens Spreading on a Network[END_REF] studied the epidemic thresholds for multiple competing viruses on special, random graphs. Pastor-Satorras and Vespignani studied viral propagation for such power-law networks [START_REF] Moreno | Epidemic Outbreaks in Complex Heterogeneous Networks[END_REF][START_REF] Pastor-Satorras | Epidemic Dynamics and Endemic States in Complex Networks[END_REF][START_REF] Pastor-Satorras | Epidemic Dynamics in Finite Size Scale-Free Networks[END_REF]. They developed an analytic model for the Barabási-Albert (BA) power-law topology [START_REF] Barabasi | Emergence of Scaling in Random Networks[END_REF]. However, their derivation depends on some assumptions which does not hold for many real networks [START_REF] Kumar | Trawling the Web for Emerging Cyber-Communities[END_REF][START_REF] Faloutsos | On power-law relationships of the Internet topology[END_REF]. Pastor-Satorras et al. [START_REF] Pastor-Satorras | Epidemic Dynamics in Finite Size Scale-Free Networks[END_REF] also proposed an epidemic threshold condition, but this uses the "mean-field" approach, where all graphs with a given degree distribution are considered equal. There is no particular reason why all such graphs should behave similarly in terms of viral propagation. Chakrabarti et. al. [START_REF] Chakrabarti | Epidemic thresholds in real networks[END_REF] observe that epidemic threshold of an arbitrary graph can be captured in a single parameter, the first eigenvalue of the adjacency matrix λ1,A. Their observation was rigorously confirmed in [START_REF] Chakrabarti | Information Survival Threshold in Sensor and P2P Networks[END_REF] and independently by [START_REF] Van Mieghem | Virus Spread in Networks[END_REF]. We again leverage the above observations to formulate our solution in Sect. 3

Epidemics on Mobile Networks. Prompted by the emergence of mobile devices, such as Bluetooth-equipped smartphones, researchers introduced mobility to epidemic spread. Mickens et. al. were among the first to examine device-to-device spreading of malicious software in mobile ad hoc networks [START_REF] Mickens | Analytical models for epidemics in mobile networks[END_REF][START_REF] Mickens | Modeling Epidemic Spreading in Mobile Environments[END_REF]. In their work, they present a queue-based technique for the RWP model to overcome the limitations of the earlier homogeneous models of Kephart and White. In a similar work, Yan et. al. extend the observations of Mickens et. al. by examining additional mobility models and their effect on epidemic spreading of a SIS virus [START_REF] Yan | Bluetooth Worm Propagation: Mobility Pattern Matters! In[END_REF]. Their work is unique in that it models virus propagation, in detail, a Cabir-like bluetooth worm, including the Bluetooth stack and unique worm properties.

Mobility Models. The mobility models used in Sect. 4 are fairly common, with significant literature devoted to the subject. For an overview on mobility models, we refer our readers to the following surveys [START_REF] Camp | A Survey of Mobility Models for Ad Hoc Network Research[END_REF]. The RWP mobility model has been extensively used, despite well know flaws. For a discussion of the merits of RWP, refer to [START_REF] Yoon | Random Waypoint Considered Harmful[END_REF]. The Levy mobility model was first described in [START_REF] Shlesinger | Random Walks with Infinite Spatial and Temporal Moments[END_REF], yet has been used extensively to model human and animal movements [START_REF] Rhee | On the Levy-Walk Nature of Human Mobility[END_REF][START_REF] Rhee | Human Mobility Patterns and their Impact on Delay Tolerant Networks[END_REF].

To conclude, recent malware in the wild, using device-to-device virus propagation schemes, prompted our study of the epidemic threshold in mobile ad hoc networks. Our contributions in this paper are:

1. Framework: We present a framework for the determining the epidemic threshold (for the SIS model) on any mobility model which can be converted into a series of adjacency matrices and give a formula for it (Theorem 1). 2. Closed Formulas: We also give a closed-form approximation for the SIS epidemic threshold on general mobility models (Lemma 1). 3. Insensitivity to Velocity: We analyze the impact of velocity in popular mobility models like Random walk, Levy Flight, Random waypoint etc. and find that it unexpectedly does not affect the threshold (Lemmas 2, 3 and Conjecture 1).

In addition, we introduced the "periodic mobility model," popular in other fields like epidemiology [START_REF] Anderson | Infectious Diseases of Humans[END_REF][START_REF] Barrett | EpiSimdemics: an Efficient Algorithm for Simulating the Spread of Infectious Disease over Large Realistic Social Networks[END_REF], to the networking community and solved it using our framework (Lemma 4). Finally we presented extensive simulations to demonstrate our analysis and results. Future work may concentrate on providing theoretical analysis on the effect of velocities on the steady state behavior of an epidemic on various mobility models.

Fig. 2 .

 2 Fig.2. Accuracy of Framework Approximations, Random Walk Take-off plots for EAAM and USS. Notice that our predictions (black arrows) are accurate.

Fig. 3 .

 3 Fig.3. Accuracy of Framework Approximations, Levy Flight Take-off plots for EAAM and USS. Again our predictions (black arrows) prove accurate.

Figure 5 (Fig. 4 .

 54 Fig. 4. Accuracy of Framework Approximations. Random Waypoint Take-off plots for EAAM and USS.

Fig. 5 .

 5 Fig. 5. Number of Infected Agents vs. Time on the Random Walk and Random Waypoint. Node velocities were varied as indicated in the legends. Steady-State number of infected agents increases with node velocity, while warm-up period shrinks (best viewed in color).

Fig. 6 .

 6 Fig. 6. Take-off Plot for System Matrix Eigenvalue of a Periodic Mobility Model. A1 and A2 had N = 500 nodes spread uniformly across a 200 × 200 simulation field. Notice the accurate our prediction (black arrow).

Table 1 .

 1 Terminology

	Symbol	Definition and Description
	General Terms:
	A, B, . . .	matrices (bold upper case)
	A(i, j)	element at the i th row and j th column of A
	A(i, :)	i th row of matrix A
	A(:, j)	j th column of matrix A
	I	standard n × n identity matrix
	a, b, . . .	column vectors
	I, J , . . .	sets (calligraphic)
	λ B	first eigenvalue (in absolute value) of a
		matrix B
	Mobility Terms:
	M	mobility model
	Pi,t	position of node i at time t
	N	number of nodes
	A	simulation area
	ρ	node density (N/A)
	∆T	Time step
	T	number of different alternating behaviors
	A1, . . . , A T	T corresponding size (n × n) symmetric
		alternating adjacency matrices
	Epidemic Terms:
	β	virus transmission probability in the SIS model
	δ	virus death probability in the SIS model
	τ	epidemic threshold
	Acronyms and Terms:
	USS	Uniform Steady-State Approximation
	EAAM	Eigenvalue of Average Adjacency Matrix
		Approximation
	Take-off Plot	Max number of infected agents vs.
		Epidemic Threshold τ

We focus on proximity propagation, which in mobile networks is affected by mobility. However, virus for smartphones can also propagate through email, mms, or direct access to the web. We do not consider these cases, sine they are not directly affected by mobility