On the Generalized Delay-Capacity Tradeoff of Mobile Networks with L\'evy Flight Mobility

In the literature, scaling laws for wireless mobile networks have been characterized under various models of node mobility and several assumptions on how communication occurs between nodes. To improve the realism in the analysis of scaling laws, we propose a new analytical framework. The framework is the first to consider a L\'{e}vy flight mobility pattern, which is known to closely mimic human mobility patterns. Also, this is the first work that allows nodes to communicate while being mobile. Under this framework, delays ($\bar{D}$) to obtain various levels of per-node throughput $(\lambda)$ for L\'evy flight are suggested as $\bar{D}(\lambda) = O(\sqrt{\min (n^{1+\alpha} \lambda, n^2)})$, where L\'evy flight is a random walk of a power-law flight distribution with an exponent $\alpha \in (0,2]$. The same framework presents a new tighter tradeoff $\bar{D}(\lambda) = O(\sqrt{\max (1,n\lambda^3)})$ for \textit{i.i.d.} mobility, whose delays are lower than existing results for the same levels of per-node throughput.

I. INTRODUCTION
Since the work in [1] that showed that mobility can be exploited to improve network throughput, there has been a plethora of work on this subject. A major effort in this direction has been in the design of delay tolerant networks (DTNs). However, this benefit in throughput comes at a significant delay cost. The amount of delays required to achieve a level of throughput for various mobility models such as i.i.d. mobility, random waypoint (RWP), random direction (RD), and Brownian motion (BM) have been extensively studied in [2]- [6]. Specifically, the delay required for constant pernode throughput has been shown to grow as Θ(n), which scales as fast as the network size n, for most mobility models including i.i.d. mobility, RWP, RD, and BM [2], [5], [7], [8]. Despite significant advances in the development of delaycapacity scaling laws, there has been considerable skepticism regarding the applicability of the results to real mobile networks because of various simplifying assumptions used in the analysis.
In this paper, we address two issues towards making the delay-capacity tradeoff analysis more realistic: 1) contacts among nodes in the middle of their movements and 2) Lévy mobility patterns of nodes in the network. In the literature, for mathematical simplicity, existing results have assumed that This work was supported in part by NSF grants CNS-1065136, CNS-1012700, CNS-0910868 and CNS-1016216, and ARO MURI Award W911NF-08-1-0238. Our considerations for improving the realism of delay-capacity tradeoff: (a) Conventionally, nodes A and B are assumed to have no contact opportunity, but we consider they can meet each other even during their movements. (b) We also consider that nodes follow Lévy flight mobility. nodes show slotted movements, and they do not communicate with each other while being mobile. Thus, they make contacts with other nodes and transfer data only at the edge of time slots. In other words, as shown in Fig. 1 (a), the opportunity for meeting other nodes during mobility has been ignored, although such opportunities can substantially reduce packet delivery delays. Also, in this work we focus on the Lévy flight model, which is widely accepted to closely mimic the actual movement of humans [9], [10]. The trajectory of this model is illustrated in Fig. 1 (b). To enhance the realism in the analysis of delay-capacity tradeoff, we develop a new analytical framework which takes both of these factors into account by developing a technique that characterizes the distribution of "first meeting time" among nodes conforming to Lèvy flight mobility in a two-dimensional space. It is important to note that the exact distribution of the first meeting time of Lévy flight even in a one-dimensional space has been an open problem even though it has applicability in a diverse set of research problems (e.g., characterization of particle movements and animal movements) in physics and mathematics. It is also informative to note that the distribution of the first meeting time of BM, which can be considered as an extreme case of Lévy flight, is also an open problem as noted in [11], [12].
Lévy flight, the mobility model we focus on in this paper, is a subset of Lévy mobility in which a node moves from position to position in a constant time. Another special case, Lévy walk, in which a node moves from one position to another in time proportional to the distance between the positions.
Except for the notion of the time required for each movement, Lévy flight and Lévy walk are fundamentally the same random walk whose flight length distribution asymptotically follows a power-law f α (z) ≈ 1/z α+1 , where z and α denote the flight length (i.e., moving distance of each slotted movement) and the power-law slope ranging 0 < α < 2, respectively. The heavy-tailed movements of Lévy mobility render the delay characterization extremely challenging. Our framework addresses these challenges using theories from stochastic geometry and probability, and provides a delay-capacity tradeoff for Lévy flight. Also, for a simpler i.i.d. mobility model, we provide a tighter delay-capacity tradeoff compared to existing studies using the same framework. Fig. 2 and Table I summarize the new tradeoffs identified using our analytical framework. The results show that the tradeoff for Lévy flight followsD(λ) = O( min{n 1+α λ, n 2 }) to obtain a per-node throughput of λ = Θ(n −η ) (0 ≤ η ≤ 1/2) as shown in Figs. 3 (a) and (b). These results are well aligned with the critical delay (i.e., minimum delay required to achieve λ = ω(1/ √ n)) suggested in [13]. Our tradeoffs show an important finding that the delay required to obtain constant per-node throughput (i.e., λ = Θ(1)) can be smaller than Θ(n) in mobile networks with mobility models such as Lévy flight with α < 1 and i.i.d. mobility. This is an important observation given that most of the existing studies present the delay required to obtain constant per-node throughput to be Θ(n) for almost all mobility models including the i.i.d. mobility.
Our tradeoff for Lévy flight becomes especially more interesting when we input α values from measurements presented in Table II into the tradeoff. This gives us a hint on how the performance of the network will scale in reality when the network consists of devices mainly carried (or driven) by humans. For α values between 0.53 and 1.81, the delays to obtain λ = Θ(1) are expected to lie between O(n 0.77 ) and O(n). This implies that in reality, a DTN mainly operated by human mobility may indeed experience less than Θ(n) delay in some areas. This observation of smaller delay suggests that mobile networks relying on opportunistic transmissions may have higher practical values in reality given that the delays have been overestimated by mobility and contact models with less realism.  The rest of the paper is organized as follows. We overview a list of related work in Section II and introduce our system models and definitions of performance metrics in Section III.    We then provide the intuition on how our analytical framework evaluates delay-capacity tradeoff using the properties of first meeting time of random walks in Section IV. Based on the understanding in Section IV, we analyze the tradeoffs of Lévy flight and i.i.d. mobility models in Sections V and VI, respectively. After briefly concluding our work in Section VII, we provide full details of proofs used in Sections V and VI through Appendices.

II. RELATED WORK
In [14], it was shown that the per-node throughput of random wireless networks with n static nodes scales as O(1/ √ n). The result was later enhanced to Θ(1/ √ n) by exercising individual power control [15]. Grossglauser and Tse [1] proved that constant per-node throughput is achievable by using mobility when nodes follow ergodic and stationary mobility models. This contradicted the conventional belief that node mobility negatively impacts network capacity due to interruptions in connectivity.
Many follow-up studies [2], [4], [5], [7], [16]- [18] have been devoted to characterize and exploit the delay-capacity tradeoffs. In particular, the delay required to obtain constant per-node throughput has been studied under various mobility models [2]- [6]. The key message is that the delay of 2-hop relaying proposed in [1] is Θ(n) for most mobility models such as i.i.d. mobility, RD, RWP, and BM.
More realistic mobility models, Lévy mobility models known to closely capture human movement patterns, were first analyzed in [13] for a special case of the tradeoff. Using spatio-temporal probability density functions, the critical delay defined by the minimum delay required to achieve larger throughput than Θ(1/ √ n) is identified for Lévy flight as well as Lévy walk.
Existing results on delay-capacity tradeoffs for mobile networks have been built under the assumption that nodes are able to communicate with each other only at the edge of time slots for slotted movements. Also, there has been no framework which can fully understand the delay-capacity tradeoff for Lévy mobility. In this paper, we develop an analytical framework which handles both of these issues and are able to use this framework to characterize the delaycapacity tradeoff.

A. Network Model
We consider a wireless mobile network indexed by n ∈ N {1, 2, . . .}, where, in the nth network, n nodes move on a completely wrapped-around disc D (⊂ R 2 ) whose radius scales as Θ( √ n). 1 Without loss of generality, we set the radius and the center of the disc D as √ n and 0 (0, 0), respectively, i.e., D = {x ∈ R 2 | |x| ≤ √ n}. We assume that the density of the network is fixed to 1 as n increases. 2 We also assume that all nodes are homogeneous in that each node generates data with the same intensity to its own destination. The packet generation process at each node is independent of node mobility. The generated packets are assumed to have no expiration until their delivery and the size of each node's buffer is assumed to be unlimited. Each packet can be delivered by either direct one-hop transmission or over multiple hops using relay nodes.
To model interference in wireless networks, we adopt the following protocol model as in [6], [14]. Let X i (t) denote the location of node i (= 1, . . . , n) at time t (≥ 0). Let L (i,j) (t) |X i (t)−X j (t)| denote the Euclidean distance between nodes i and j at time t. Under the protocol model, nodes transmit packets successfully at a constant rate W bits/sec, if and only if the following is satisfied: for a transmitter i, a receiver j and every other node u = i, j transmitting simultaneously, where ∆ is some positive constant. In addition, the distance between nodes i and j at time t should satisfy L (i,j) (t) ≤ r, where r (> 0) denotes the maximum communication range. We assume the fluid packet model [6], which allows concurrent transmissions of node pairs (with the rate divided by the number of pairs) interfering each other. We denote by Π the class of all feasible scheduling schemes conforming our descriptions.

B. Mobility Model
In this subsection, we mathematically describe the Lévy flight model and the i.i.d. mobility model. At time t = 0, node i chooses its location uniformly on the disc D (i.e., X i (0) ∼ Uniform(D)), which is independent of the others X j (0) for j = i. We assume that time is divided into slots of unit length and is indexed by k ∈ N. At the beginning of the kth slot (i.e., at time t = k −1), node i chooses its next slotted location X i (k) according to the associated mobility model. During the kth slot (i.e., during time t ∈ (k − 1, k]), node i moves from X i (k−1) to X i (k) with a constant velocity. Thus, and X i (k) as follows: (1)

Lévy Flight Model.
At the beginning of the kth slot (i.e., at time t = k − 1), node i chooses flight angle and flight length, denoted by θ i (k) (∈ (0, 2π]) and Z i (k) (> 0), respectively. During the kth slot (i.e., during time t ∈ (k − 1, k]), node i moves from X i (k − 1) to the selected direction θ i (k) of the distance Z i (k). Thus, the location X i (k) is determined as where The flight angle θ i (k) and the flight length Z i (k) are independent of each other and also independent of the previous locations X i (t) for the times t ∈ [0, k − 1] before they are generated. Hence, Each flight angle θ i (k) and flight length Z i (k) are independent and identically distributed across node index i and slot index k. Let θ and Z be a generic random variable for θ i (k) and Z i (k), respectively. Then, the flight angle θ is uniformly distributed over (0, 2π], and the flight length Z is generated from a random variable Z ⋆ having the Lévy α-stable distribution [20] by the relation Z = |Z ⋆ |. The probability density function of Z ⋆ is given by where ϕ Z ⋆ (t) E[e itZ ⋆ ] is the characteristic function of Z ⋆ and is given by ϕ Z ⋆ (t) = e −|st| α . Here, |s| > 0 is a scale factor determining the width of the distribution, and α ∈ (0, 2] is a distribution parameter that specifies the shape (i.e., heavytail-ness) of the distribution. The flight length Z for α ∈ (0, 1) has infinite mean and variance, while Z for α ∈ [1, 2) has finite mean but infinite variance. For α = 2, the Lévy α-stable distribution reduces to a Gaussian distribution with a mean of zero and variance of 2s 2 , for which the flight length Z has finite mean and variance. Due to the complex form of the distribution, the Lévy αstable distribution for α ∈ (0, 2) is often treated as a power-law type of asymptotic form: where we use the notation a(z) ∼ b(z) for any two real functions a(z) and b(z) to denote lim z→∞ [a(z)/b(z)] = 1 [20]. The form (5) is known to closely approximate the tail part of the distribution in (4), and a number of papers in mathematics and physics, e.g., [21], [22], analyze Lévy mobility using the form (5). For mathematical tractability, in our analysis we will also use the asymptotic form (5). Specifically, we assume that there exist constants c (> 0) and z th (> 0) such that i.i.d. Mobility Model. At the beginning of the kth slot (i.e., at time t = k − 1), node i chooses X i (k) uniformly on the disc D, which is independent of its previous locations X i (t) for the times t ∈ [0, k − 1] as well as the others X j (k) for j = i. Thus, X i (k) is independent and identically distributed across node index i and slot index k.

C. Contact Model
In our contact model, nodes are allowed to meet while being mobile. Hence, for a time t ⋆ in a domain {t | t ≥ 0}, we say that nodes i and j meet at time t ⋆ (or are in contact at time t ⋆ ) if they satisfy In the widely adopted contact model where nodes are allowed to meet only at the end of their movements (i.e., at slot boundaries), a meeting event can occur for a time k ⋆ in a domain {k | k ∈ {0} ∪ N} satisfying L (i,j) (k ⋆ ) ≤ r. We call this class of contact model slotted contact model throughout the paper.
Mobile nodes are exposed to more contact opportunities in our contact model compared to the slotted contact model.

D. Performance Metrics
The key performance metrics of our interest are per-node throughput and average delay as defined next: Definition 1 (Per-node throughput). Let Λ π:i (t) be the total number of bits received at the destination node i up to time t under a scheduling scheme π ∈ Π. Let λ π be the per-node throughput under π. Then, Definition 2 (Average delay). Let D π:i,v be the time taken for the vth packet generated from the source node i to arrive at its destination node under a scheduling scheme π ∈ Π. Let D π be the average delay under π. Then, In this paper, we focus on analyzing the scaling property of the smallest average delay achieving per-node throughput λ. We call this minimum average delay optimal delay throughout this paper. We focus on the throughput only in the range from Θ(1/ √ n) to Θ(1), since this range corresponds to the case where mobility can be used to improve the per-node throughput.
Definition 3 (Optimal delay). LetD(λ) be the optimal delay to achieve per-node throughput λ. It is then given bȳ

E. Throughput Achieving Scheme
We now consider a schemeπ that can realize per-node throughput λπ scaling from Θ(1/ √ n) to Θ(1). The schemeπ operates as follows: • When a packet is generated from a source node and the destination of the packet is within the communication range of the source node, the packet is transmitted to the destination node immediately. • Otherwise, the source node broadcasts the packet to all neighboring nodes within its communication range. Note that this broadcast is only performed by the source node when the packet is generated. • Any nodes carrying the packet can deliver the packet to the destination node when they are within the communication range of the destination node. • When one of the packets (including the original packet and the duplicated ones) reaches the destination node, all others are not considered for delivery.
By appropriately scaling r as a function of n, the schemeπ can achieve per-node throughput λπ ranging from Θ(1/ √ n) to Θ(1), as shown in the following lemma. Proof: If the network has been running for a long enough time, all nodes become to work as relay nodes and begin to have packets for all other nodes. Therefore, for a network with n/a n disjoint area where a n = Θ(r 2 ), all areas with more than two nodes can always be activated. Let b n denote the probability of having more than two nodes in an area. We then have b n = 1 − 1 − a n n n − n a n n 1 − a n n n−1 .
In addition, the total network throughput becomes b n n/a n and accordingly the per-node throughput is λπ = b n /a n . Without loss of generality, we assume r = n β (0 ≤ β ≤ 1/4). Then, the per-node throughput λπ is given by In the following, we will show that Hence, the per-node throughput λπ under the schemeπ scales as Θ(n −2β ). Note that Similarly, which gives By applying (8) and (9) to (7a), we have (7b). This complete the proof.
LetDπ(β) be the average delay underπ when r = Θ(n β ). Lemma 1 implies that by setting β = − log n ( √ λ), the schemeπ achieves the per-node throughput λπ = λ. Since the schemeπ is of the class Π, the order ofD(λ) can be obtained fromDπ(β) with the use of β = − log n ( √ λ), i.e., IV. PRELIMINARIES In this section, we provide the key intuition to understand how our analytical framework utilizes the properties of first meeting time in the derivation of delay-capacity tradeoffs under the Lévy flight and the i.i.d. mobility models. We then sketch the challenges residing in our framework and briefly describe our approach to address these challenges.

A. Delay Analysis with First Meeting Time
The first meeting time of two nodes moving in a twodimensional space, which is directly connected toDπ, is defined below: Definition 4 (First meeting time). For i = j, the first meeting time of nodes i and j, denoted by T (i,j) , is defined as Since T (i,j) is independent and identically distributed across pair index (i, j), we use T to denote a generic random variable for T (i,j) .
Let D (s,d) be a random variable representing the time taken by a packet generated from a source node s to arrive at a destination node d. Since the packet generation process is independent of node mobility, we consider that each packet is generated at time t = 0 without loss of generality. Then, the packet delay under the schemeπ, denoted by Dπ :(s,d) , can be expressed in terms of the first meeting time as where I(s) {i | L (s,i) (0) ≤ r} denotes a set of node indices that are within the communication range of the node s at time t = 0. Note that s ∈ I(s) by definition. Hence, the following equation represents the schemeπ described in Section III: From Definition 2, the average delayDπ can be obtained bȳ

B. Distribution of the First Meeting Time
In order to evaluate (11), the distribution of the first meeting time P{T (i,d) < τ } is essential. To obtain the distribution, we start from defining the following: let I (i,j) (k) (i = j, k ∈ N) be a random variable indicating the occurrence of a meeting event between nodes i and j during the kth slot (i.e., time t ∈ (k − 1, k]), i.e., For notational simplicity, throughout this paper, we omit (i, j) in I (i,j) (·) and L (i,j) (·), unless there is confusion. We then define a function H(k, l 0 ) for k ∈ N and l 0 ∈ (r, 2 √ n], which denotes the probability that nodes i and j are not in contact during the kth slot, conditioned on the fact that the initial distance between the nodes was l 0 and after that the nodes have not been in contact by time t = k − 1, i.e., Note that l 0 is upper bounded by 2 √ n since the radius of the disc D is set to √ n. We find that the distribution of the first meeting time T can be obtained from the function H(k, l 0 ) as shown in the following lemma.
The identity P{T > 0} = P{L(0) > r} shown in Lemma 2 has the following implications for P{T > 0}: (i) It is determined by the spatial distribution of nodes at time t = 0 (which is assumed to be uniform on the disc D). Hence, it is invariant for both the Lévy flight and the i.i.d. mobility models. (ii) It represents the probability that two arbitrary nodes are out of the communication range at time t = 0. Since P{T > 0} is frequently used throughput this paper, we define it as P o and summarize its implications using the following lemma.
Then, P o is bounded by Proof: See Appendix A.

C. Technical Challenge and Approach
In our framework, characterizing the function H(k, l 0 ) in (12) which appears in the expression for P{T > τ } in (13) is the key to analyze the optimal delay. The major technical challenge arises from tracking meeting events in the middle of a time slot. The meetings over time are heavily correlated irrespective of mobility models. The correlation can be understood as follows: let us consider two consecutive slots, say the kth and the (k + 1)th slots, for ease of explanation. The occurrence of a meeting event during the kth slot (resp. the (k + 1)th slot) is determined by the locations of nodes i and j at the slot boundaries, i.e., at times t = k − 1, k (resp. at times t = k, k + 1). Hence, both I(k) and I(k + 1) depend on the values of X i (k) and X j (k), and accordingly the sequence {I(k)} k∈N is correlated in our contact model. Due to the complexity involved in the correlation, deriving the exact form of H(k, l 0 ) appears to be mathematically intractable. To address this challenge, we take a detour to derive a bound on H(k, l 0 ) using theories from stochastic geometry and probability. The detailed analysis of H(k, l 0 ) for the Lévy flight model and the i.i.d. mobility model is presented in Lemmas 5 and 10, respectively, which allow us to reach the conclusions of this paper.

V. DELAY ANALYSIS FOR THE LÉVY FLIGHT MODEL
In this section, we analyze the optimal delay under the Lévy flight model. We use the following four steps in our analysis: Step 1 (Formulation of the average delay using the first meeting time distribution): As shown in (11) Under the Lévy flight model, T (i,d) for i ∈ I(s) are heavily correlated since the next slotted location X i (k + 1) depends on the current location X i (k). Note that all the nodes i ∈ I(s) are in proximity of the node s, and thus min(T (i,d) ; i ∈ I(s) \ {d}) is not easily tractable. Therefore, we use the following bound to describē Dπ using T .Dπ For the simpler i.i.d. mobility model, we are able to derive a tighter bound on min(T (i,d) ; i ∈ I(s) \ {d}). We present the result in Step 1 of Section VI. Let ⌈T ⌉ denote the smallest integer greater than or equal to T . Then, since , and E[⌈T ⌉] is an upper bound onDπ, as shown in the following lemma.

Lemma 4. The average delayDπ of the schemeπ under the Lévy flight model is bounded bȳ
where T is the generic random variable for the first meeting time , which gives (16). Since the random variable ⌈T ⌉ takes on only nonnegative integer values, the expectation E[⌈T ⌉] can be obtained by where the second equality comes from the property that P{⌈T ⌉ ≥τ } = P{T >τ − 1} for allτ = 1, 2, . . .. Replacinĝ τ − 1 with τ gives the lemma. Step

Lemma 5.
Under the Lévy flight model, the function H(k, l 0 ) in (12) has the following characteristics: representing a flight differential between nodes i and j. 3 Then, geometrically the function H(1, l 0 ) can be viewed as the probability of the flight differential falling into a set S(l 0 ) (⊂ R 2 ) defined as follows. Let D r (u) (⊂ R 2 ) denote a disc of radius r centered at u ∈ R 2 , i.e., D r (u) {x ∈ R 2 | |x − u| ≤ r}. 3 The existence of the generic random variable for Let (v, w) denote a line connecting two points v, w ∈ R 2 . For a fixed l 0 ∈ (r, 2 √ n], define a set S(l 0 ) as An example of S(l 0 ) is shown in Fig. 4. The set S(l 0 ) has a connection with the function H(1, l 0 ) as follows: (iv) For k = 2, 3, . . ., each function H(k, l 0 ) is also bounded above byP . Thus, for all k ∈ N and l 0 ∈ (r, 2 √ n], we have (v) There exist constants c l , c u (> 0), and n th ∈ N such that for all n ≥ n th ,P is bounded above and below bŷ Based on the formula for P{T > τ } in Lemma 2 and (iv) in Lemma 5, we derive a bound on P{T > τ } in terms ofP as follows: for τ = 1, 2, . . ., we have from (iv) in Lemma 5 that Since (20) holds for all l 0 ∈ (r, 2 √ n], by integrating (20) over l 0 ∈ (r, 2 √ n], we have By combining (21) and Lemma 2, we have Since P{T > τ } = P o for τ = 0, the bound in (22) also holds for τ = 0. The above result is summarized in Lemma 6.

Lemma 6. Under the Lévy flight model, the CCDF of the first meeting time T is bounded by
whereP and P o are defined in (18) and (14), respectively.

Lemma 7. The expectation E[⌈T ⌉] under the Lévy flight model is bounded by
Proof: Using Lemma 6, we can give a bound on E[⌈T ⌉] in Lemma 4 as By (v) in Lemma 5, we haveP < 1 for any r > 0. Thus, the expectation E[⌈T ⌉] is bounded by the geometric series which converges to Po 1−P .
Step 3 (Analysis of the delay scaling): In Lemma 1, we have analyzed the order of the per-node throughput λπ of the schemeπ. The results in Lemma 4,(v) in Lemma 5, and Lemma 7 allow us to analyze the order of the average delayDπ, which is shown in Lemma 8.
In addition, under the Lévy flight model, we have a trivial upper bound for all α ∈ (0, 2] as Combining (23) and (24) yields our lemma.
Step 4 (Analysis of the delay-capacity tradeoff): In the last step, we derive the delay-capacity tradeoff under the Lévy flight model. By combining the capacity scaling in Lemma 1 and the delay scaling in Lemma 8, we get the following theorem.

VI. DELAY ANALYSIS FOR THE i.i.d. MOBILITY MODEL
In this section, we provide detailed analytical steps for obtaining the optimal delay under the i.i.d. mobility model. We again follow the four steps analogous to those used for the Lévy flight model.
where T ⋆ v (v = 1, . . . , m) denotes the first meeting time of the node d and the vth node in the set I(s), provided that |I(s)| = m and d / ∈ I(s). Let T 1 , . . . , T m be m independent copies of the generic random variable T . Then, we can derive a bound on T ⋆ v in terms of T v as follows: Here, i v denotes the vth index in the set I(s) and d = denotes "equal in distribution". The last equation comes from the aforementioned nature of the i.i.d. mobility model in which the locations of nodes are reshuffled at every time slot.
We define a functionŪ : Note that discretization of a random variable T v to ⌈T v ⌉ is for mathematical simplicity and it does not affect the result (i.e., order of the optimal delay) of this paper. The function U (m) works as a tight upper bound onDπ as shown in the following lemma.
Step   (12) has the following characteristics: (i) Let ∆X be a generic random variable for X i (t) − X j (t) representing a location differential between nodes i and j. 4 Then, geometrically the function H(1, l 0 ) can be viewed as the probability of the location differential falling into a set where the definitions of (·, ·) and D · (·) can be found in Lemma 5. An example of S ⋆ (l 0 ) is shown in Fig. 5. The set S ⋆ (l 0 ) has a connection with the function H(1, l 0 ) as follows: (v)P is bounded above and below for all n ∈ N bŷ Similarly to Step 2 in Section V, we derive a bound on P{T > τ } in terms ofP as follows: from (iv) in Lemma 10 and Lemma 2, we have Since P{T > τ } = P o for τ = 0, the bound in (37) also holds for τ = 0. The above result is summarized in Lemma 11. whereP and P o are defined in (18) and (14), respectively.
Proof: Combining Lemma 2 and (iv) in Lemma 10 gives Lemma 11. The detailed derivation was described earlier in (37).
Using Lemma 11, we can give a bound on the function U (m) in Lemma 9 as By (v) in Lemma 10, we have (P ) m < 1 for any r > 0. Thus, U (m) is bounded by a convergent geometric series and we summarize the result in Lemma 12.
where γ is a constant in (0, 1) and γr 2 implies the γ fraction of the average number of nodes within the communication range of a source node. In (40), we used the property that U (m) is a nonincreasing function of m. Hence, by Lemma 9 and (40), the average delayDπ of the schemeπ under the i.i.d. mobility model is bounded by: The results in (41), Lemmas 3 and 12, and (v) in Lemma 10 allow us to analyze the order of the average delayDπ, which is shown in Lemma 13. Proof: Here, we provide a sketch of the proof with details given in Appendix C.

Step 4 (Analysis of the delay-capacity tradeoff):
In the last step, we derive the delay-capacity tradeoff under the i.i.d. mobility model. By combining the capacity scaling in Lemma 1 and the delay scaling in Lemma 13, we get the following theorem.

VII. CONCLUDING REMARKS
In this paper, we developed a new analytical framework that substantially improves the realism in delay-capacity analysis by considering (i) Lévy flight mobility, which is known to closely resemble human mobility patterns and (ii) contact opportunities in the middle of movements of nodes. Using our framework, we obtained the first delay-capacity tradeoff for Lévy flight and derived a new tighter tradeoff for i.i.d. Our result is well aligned with the critical delay suggested in [13]. For i.i.d. mobility, our analysis providesD(λ) = O( max(1, nλ 3 )) as shown in Fig. 2 (b). These tradeoffs are especially remarkable in both Lévy flight and i.i.d. mobility for the constant pernode throughput (i.e., λ = Θ(1)) as they demonstrate that the delay can be less than Θ(n), which has been widely accepted for most mobility models. Our future work includes (i) an extension of our framework to analyze the delay-capacity tradeoff under Lévy walk and (ii) another extension to capture correlated movement patterns among nodes.

APPENDIX A PROOF OF LEMMA 3
By the definition of P o in (14), we have Let F Xi(0) (·) denote the CDF of X i (0). Then, by conditioning on the values of X i (0), the probability P{L (i,j) (0) ≤ r} in (46) can be rewritten as where the last equality comes from the independence between X i (0) and X j (0). Note that, since X j (0) ∈ D with probability 1 and X j (0) ∼ Uniform(D), the probability P{X j (0) ∈ D r (x)} in the integral in (47) is given by where Area(S) denotes the area of a set S ⊂ R 2 . An example of D∩D r (x) is shown in Fig. 6. From the figure, it is obvious that Area(D ∩D r (x)) is nonincreasing as x approaches to the boundary of the disc D. Hence, (48) is bounded above by In addition, it is bounded below by where θ 2 cos −1 ( r 2 √ n ) (See Fig. 6). Since √ n ≥ r, we have θ ≥ 2 cos −1 ( r 2r ) = 2π 3 . Hence, the inequality in (50) is further bounded by By substituting (49) and (51) into (47), we have which, combined with (46), gives

APPENDIX B PROOFS OF LEMMAS FOR THE LÉVY FLIGHT MODEL
Here, we give detailed proofs of Lemmas 5 and 8, which are used for analyzing the optimal delay under the Lévy flight model. To prove Lemma 5, we need the following Lemmas 14, 15, and 16.

Lemma 14.
For i = j and k ∈ N, let where V · (k) (representing the kth flight of a node ·) is defined in (3). Then, under the Lévy flight model, ∆V (i,j) (k) has the following properties: and slot index k. Hence, we use ∆V to denote a generic random variable for ∆V (i,j) (k).  (1) and (2)). Since V i (k) is independent of F u (k − 1), it is independent of X u (t). By the same reason, V j (k) is independent of X u (t). Therefore, the difference V i (k)−V j (k) is independent of X u (t).
(ii) Since each of the flight angle θ u (k) and the flight length An example of vectors satisfying the condition in (54): for a given (v i , v j , v), the angle ∠∆V (i,j) is determined by the angle ∠V j . Since Z u (k) is independent and identically distributed across node index u and slot index k, the random variable V u (k) ( (Z u (k) cos θ u (k), Z u (k) sin θ u (k))) is also independent and identically distributed across u and k. Therefore, the difference V i (k) − V j (k) is identically distributed across pair index (i, j) and slot index k. However, it is not necessarily independent across (i, j) while it is independent across k for a fixed (i, j). (iii) To prove (iii), it suffices to show that for any v ≥ 0, where 0 < θ ≤ 2π. In the following, we will prove (52). For simplicity, we omit the slot index k in V · (k) and ∆V (i,j) (k) in the rest of this proof. By conditioning on the values of (|V i |, |V j |), we can rewrite the probability on the left-hand side of (52) as follows: For a fixed v ≥ 0, consider an event An example satisfying (54) is shown in Fig. 7. Under the con- is determined by the angle ∠V j as the figure shows. Since Since the above equality holds for any (v i , v j ) satisfying (54) for a given v, the probability in (53) boils down to the following: This completes the proof.

Lemma 15.
Suppose k ∈ N and l ∈ (r, 2 √ n]. Then, for any sets L(·) ⊂ [0, 2 √ n] satisfying we have under the Lévy flight model the following: The definitions of ∆V and S(l) can be found in Lemma 5.

A. Proof of (i)
To prove (i), we need the following: for any x ∈ [0, 1], The proof of (87) is given at the end of this section. From (19a) in Lemma 5 with r = n β , we have for all n ≥ n th the following: Since n β−1/2 2 ∈ [0, 1] for any β ∈ [0, 1/4] and n ∈ N, we further have from the lower inequality in (87) that Hence, we have Using a similar approach as above, from (19b) in Lemma 5 and the upper inequality in (87), we have for all n ≥ n th the following: Hence, we have lim sup Combining (88) and (89) proves (i).

B. Proof of (ii)
Without loss of generality, we assume P{Z α > z th } = 1. (In this proof, subscript α is added to all random variables to specify the underlying parameter α of the Lévy flight model.) Then, from (6), we have P{Z α > z} = ( zth z ) α for all z ≥ z th , which gives for any 0 < α 1 ≤ α 2 ≤ 2 and z ≥ z th the following: The inequality in (92) shows that for any t 2 > t 1 ≥ 0 having a sufficiently small difference ǫ t 2 − t 1 > 0, we get Note that since P{T α1 > t} = P{T α1 > t, T α1 > t − ǫ} for t ≥ ǫ, we can express P{T α1 > t} in a nested form as Using the nested form continuously, we have Hence, by applying (93) to (94), we have Note that P{T α1 > 0} = P{L α1 (0) > r}. In addition, since X i (0) ∼ Uniform(D) for all i = 1, . . . , n regardless of α, we have P{T α1 > 0} = P{T α2 > 0}. Thus, the right-hand side of (95) boils down to P{T α2 > t}, and consequently P{T α1 > t} ≤ P{T α2 > t} for all t ≥ 0.
where X · (t) denotes the location of a node · at time t. Then, under the i.i.d. mobility model, ∆X (i,j) (t) has the following properties: (i) ∆X (i,j) (k) (k ∈ N) is independent of X u (t) for all u = 1, . . . , n and t ∈ [0, k − 1]. (ii) ∆X (i,j) (t) is identically distributed across pair index (i, j) and time t (≥ 0). Hence, we use ∆X to denote a generic random variable for ∆X (i,j) (t).
(ii) For any i = j and t ≥ 0, X i (t) and X j (t) are independent and identically distributed. Therefore, the difference X i (t)−X j (t) is identically distributed across pair index (i, j) and time t. However, it is not necessarily independent neither across (i, j) nor across t. (iii) To prove (iii), it suffices to show that for any x ≥ 0, where 0 < θ ≤ 2π. By noting that ∠X i (t) ∼ Uniform[0, 2π] for any i = 1, . . . , n and t ≥ 0 and using a similar approach as in the proof of (iii) in Lemma 14, we can prove (iii) in Lemma 17. Due to similarities, we omit the details.
we have under the i.i.d. mobility model the following: The definitions of ∆X and S ⋆ (l) can be found in Lemma  Proof: Using a similar approach as in the proof of Lemma 15, we can prove Lemma 18. The difference is that the key idea of this proof is to use the following equality: for any k ∈ N, l ∈ (r, 2 √ n], θ ∈ (0, 2π], and F (k − 1), we have where the definition of F (k − 1) can be found in (57). Then, similarly to the proof of Lemma 15, using the key equality in (99) we can prove Lemma 18. Due to similarities, we omit the details.
In the following, we show (99). We first consider the event {I(k) = 0}. Since (61) also holds for the i.i.d. mobility model, by the same reason in the proof of Lemma 15, the event {I(k) = 0} occurs if and only if the following event occurs: We next consider the event {I(k) = 0} conditioned by ∆X(k − 1) = le jθ and F (k − 1). Under these conditions, (100) is reduced to the following: An example of S ⋆ (l, θ) is shown in Fig. 11. Hence, the probability on the left-hand side of (99) becomes By (i) in Lemma 17, ∆X(k) is independent of ∆X(k − 1) and F (k − 1), and thus we have In addition, by (ii) in Lemma 17, P{∆X(k) ∈ S ⋆ (l, θ)} = P{∆X ∈ S ⋆ (l, θ)}.
In addition, for any u ∈ D and sufficiently small ǫ (> 0), the area Area(R (x,ǫ) (u)) is calculated as By applying the combined result of (107) and (108) to (106), we have This proves the lemma.

A. Proof of (i)
Similarly to the proof of (i) in Lemma 5, we can prove (i) in Lemma 10 by applying Lemma 18 with k = 1, l = l 0 , and L(0) = (r, 2 √ n]. Due to similarities, we omit the details.
Hence, we have From Lemma 3, we also have P c o = 1 − P o ≤ r 2 n . Hence, the term 2(n − 2)P c o in (131) is bounded above by 2(n − 2)P c o ≤ 2(n − 2) r 2 n ≤ 2r 2 , from which we have Thus, by (132) and (133), the argument of the exponential function in (131)


extfieldsofstudy:
Computer Science
Medicine

text:
BCI Competition IV – Data Set I: Learning Discriminative Patterns for Self-Paced EEG-Based Motor Imagery Detection

Detecting motor imagery activities versus non-control in brain signals is the basis of self-paced brain-computer interfaces (BCIs), but also poses a considerable challenge to signal processing due to the complex and non-stationary characteristics of motor imagery as well as non-control. This paper presents a self-paced BCI based on a robust learning mechanism that extracts and selects spatio-spectral features for differentiating multiple EEG classes. It also employs a non-linear regression and post-processing technique for predicting the time-series of class labels from the spatio-spectral features. The method was validated in the BCI Competition IV on Dataset I where it produced the lowest prediction error of class labels continuously. This report also presents and discusses analysis of the method using the competition data set.

INTRODUCTION
Self-paced brain-computer interfaces (BCIs) have received increasing attention in recent years in the BCI community (Mason and Birch, 2000;Blankertz et al., 2002;Millan and Mourino, 2003;Scherer et al., 2008;Zhang et al., 2008). Conventionally, BCIs often require the users to follow a specific computer-generated cue before performing a specific mental control task. By contrast, self-paced BCIs allow the users to perform the control at will at anytime, by detecting specific brain signals associated with the mental control at each and every time point (Blankertz et al., 2002(Blankertz et al., , 2007(Blankertz et al., , 2008a. Providing the user with continuous control is also important for efficient control and user-training applications: first, it means that the user can perceive the system's response in a continuous and real-time manner so as to plan the mental task to activate desired BCI actions; second, it enables real-time feedback training in which the users (ALS patients; Kübler et al., 2005), for example) can learn to regulate brain waves so as to improve the BCI performance.
Various EEG modalities have been demonstrated for self-paced BCIs, such as P300 , motor imagery (Townsend et al., 2004), or finger movement related signals (Mason and Birch, 2000). In this work we focus on motor imagery (MI) which is the mental rehearsal of a motor act without any real motor output. It provides an important means for the BCIs that directly compensate lost motor functions for physically disabled (Pfurtscheller et al., 1997). It has also been shown that naïve subjects can operate a MI-BCI .
Numerous signal processing and pattern recognition techniques have been developed for classification of two or multiple MI classes (e.g., imaginary movements of left hand, right hand, tongue, or foot). For example, the common spatial pattern method (CSP; Pfurtscheller et al., 1997) and its various extensions (Lemm et al., 2005;Blankertz et al., 2008b;Wu et al., 2008) are widely used for extracting discriminative spatial (or joint spatio-spectral) patterns that contrast the power features of spatial patterns in different MI classes. For tackling multi-class problems, an information theoretic feature extraction method (Grosse-Wentrup and Buss, 2008) and other extensions of CSP (Dornhege et al., , 2004 have been proposed. And various classifiers have also been studied for MI classification . In self-paced MI-BCIs, the system not only needs to differentiate between specific MI activities, but also it has to detect them against a not-so-well-controlled class called non-control which is the aggregate of all user states other than the MI activities. For example in Millan and Mourino (2003), a local neural classifier was used to reject non MI signals. In Townsend et al. (2004), CSP features were combined with a linear discriminant analysis method to produce a scalar feature, which determined the userstate (e.g., left hand MI or NC) via thresholds. Special processing of the output incorporating a dwell and refractory period was also studied. More recently, it was reported in Scherer et al. (2008) that able-bodied subjects were able to navigate through a virtual environment using a self-paced MI-BCI with three bi-polar EEG channels only. It worked by combining two classifiers: one for discrimination between MI tasks and the others for detecting specific motor activities in the brain. It was also reported that the second classifier was sensitive to the non-stationarity of EEG, and this may indicate how challenging it is to differentiate MMI EEG against NC.
In this paper we present and study a MI detection method that won the first place in the BCI Competition IV, Dataset I. In particular, we develop a robust machine learning technique for spatio-spectral feature extraction and selection to differentiate multi-class MI and NC. We also employ a non-linear regression machine to predict the class labels at each time from the EEG features. And we propose a non-linear regression method www.frontiersin.org to post-process the time-series of predicted class labels so as to improve prediction accuracy. We conduct an offline analysis using 5-fold cross-validation on the calibration data. The method yields a mean-square-error (MSE) for class label prediction in the range from 0.20 to 0.29 for the subjects. The method then produces an average MSE of 0.38 on the evaluation data.
The rest of the paper is organized as follows. Section 2 briefly describes the EEG data. Section 3 overviews the system and presents detailed descriptions of the essential components. We study the system using the EEG data in Section 4, and present the conclusion in Section 5.

THE EEG DATA
The proposed method was evaluated using the BCI Competition IV Dataset I (Blankertz et al., 2007), which was recorded from 4 human subjects performing motor imagery tasks. A few computergenerated artificial data were also present in the dataset, though they are not considered both in this study and in the competition. The EEG recordings included a total of fifth-nine channels that mostly were distributed over or around sensorimotor areas. Each subject contributed in two sessions: a calibration session and an evaluation session (Scherer et al., 2008).
In the calibration session, each subject chose to perform two classes of motor imagery tasks from left hand, right hand, or foot imaginary movements. At the beginning of each task, a visual cue was displayed in a computer screen to the subject who then started to perform a motor imagery task accordingly for 4 s. Each subject performed a total of 200 motor imagery tasks that were balanced between the two classes. The motor imagery tasks were interleaved with breaks of 4-s long. The evaluation session followed a different protocol. The subjects followed the soft voice commands from an instructor to perform motor imagery tasks of varying time length between 1.5 and 8 s. Consecutive tasks were also interleaved with a varying time length interval also between 1.5 and 8 s.
The four human subjects are referred to as "a," "b," "f," and "g." For the sake of computational efficiency, we use the 100-Hz version of the data in this work.

Overview
As illustrated in Figure 1, the system consists of three processing components, namely, m-class feature extraction based on the filter-bank common spatial pattern (FBCSP) technique (Ang et al., 2008), information theoretic feature selection and non-linear regression for sample-based prediction, and post-processing. The input is a stream of multi-channel EEG waveforms, while the output is a sequence of scalar values in the range of [−1 1]. And the output is expected to approximate the ordinal class labels at each time point: −1 for the first motor imagery class (MI-1), 0 for NC, and 1 for the second motor imagery class (MI-2). Please note that the accuracy of the system is measured by the mean-square-error between the outputs and the true class labels.
The first two processing components work on individual windowed EEG data that we refer to as EEG epochs. Basically, the epochs are created by a shifting short-time window of a specific length (we will discuss the window length later) running through the data stream. The ending point of the window is the current sample whose class label is to be predicted. The third component, i.e., post-processing, looks into not only the current short-time window, but also the history of EEG data so as to improve the prediction. The following subsections will describe the three components successively.

m-class FBCSP for feature selection
The primary phenomenon of MI EEG is event-related desynchronization (ERD) or event-related synchronization(ERS; Pfurtscheller et al., 1997;Müller-Gerking et al., 1999), the attenuation, or increase of the rhythmic activity over the sensorimotor cortex generally in the μ (8-14 Hz) and β (14-30 Hz) rhythms. The ERD/ERS can be induced by both imagined movements in healthy people or intended movements in paralyzed patients (Dornhege et al., 2004;Kübler et al., 2005;Grosse-Wentrup and Buss, 2008).It is noteworthy that another neurological phenomenon called Bereitschafts potential is also associated with MI EEG but non-oscillatory . In this works we consider ERD/ERS features only.
Feature extraction of ERD/ERS is, however, a challenging task due to its poor, low signal to noise ratio. Therefore, spatial filtering in conjunction with frequency selection (via processing in either temporal domain or spectral domain) in multi-channel EEG has been highly successful for increasing the signal to noise ratio (Ramoser et al., 2000;Lemm et al., 2005;Dornhege et al., 2006;Blankertz et al., 2008b;Zhang et al., 2011). This technique of mclass FBCSP is developed from Ang et al. (2008), but account for the NC class.
Suppose there is an EEG epoch given in form of a matrix where n c is the number of channels, and n t is the number of time samples. The precise frequency band which is responsive to MI activities can vary from one subject to another. From the viewpoint of learning, multiple possible frequency bands have to be examined. Therefore, we process the EEG epoch using an array of FIGURE 1 | Online processing system. It consists of three processing steps to map continuous EEG data into a final output in the range of [−1 1]. Here FBCSP stands for filter-bank common spatial pattern (see Section 2), and TDNN for time-delay neural network (see Section 4).

Frontiers in Neuroscience | Neuroprosthetics
filter-banks, each of them is a particular band-pass filter and all together they cover a continuous frequency range, e.g., from 4 to 32 Hz. In this work, a total of 8 zero-phase filters based on Chebyshev Type II filters are built with central frequencies from 8 to 32 Hz at a constant interval in the logarithm domain. Consequently, the center frequencies are respectively 8, 9.75, 11.89, 14.49, 17.67, 21.53, 26.25, and 32 Hz. The filters have a uniform Q factor (bandwidth-to-center frequency) of 0.33 as well as an order of 4.
Let us useX n to denote a processed EEG epoch by an arbitrary filter-bank, say the n-th one. A set of spatial filters then apply to the matrix to extract the spatio-spectral characteristics of MI activities in that frequency band. Each spatial filter transformsX n into a time-series of the same length by where nm is the index of the spatial filter, t :1 ≤ t ≤ n t is the index of time sample, andx n (t ) is the column vector of matrix X n at time t.
The principle of the spatial filtering is to maximize the contrast between two classes in y in terms of Rayleigh coefficient where R 1 and R 2 are the covariance matrices of all epochs ofx n in the two classes respectively. The maximization of the Rayleigh coefficient is achieved by solving the following generalized eigenvalue problem (Ramoser et al., 2000) R 1 w = γ R 2 w. (4) The system selects the maximal 2 and the minimal 2 eigenvalues and the corresponding eigenvectors for w.
It should be noted that the feature used for BCI are the short-time power of the filtered signal y where l defines the length of the short-time window. Different from the conventional 2-class BCIs, the present system has to account for MI classes in addition to NC. And it takes a pair-wise approach. The first pair compares MI-1 vs. MI-2. The second and the third pairs compare MI-1 vs NC and MI-2 vs NC respectively. The construction of the spatial filters for MI-1 vs. MI-2 is performed with the original CSP method, while the spatial filters for the other two pairs are obtained in a different way as below.
To contrast MI-1 vs NC or MI-2 vs NC, we need to account for large within-class variations in NC, where the brain activity is not as well controlled as in motor imagery. Therefore, we consider that NC epochs can be further categorized into sub-states or modes. And each sub-state may exhibit different spatio-spectral characteristics, and form a particular cluster in EEG observations. Therefore, we employ a clustering approach to identify the sub-states of NC EEG epochs. First, we reduce the dimensionality of EEG by principal component analysis (PCA), and we use the energy vector u from the top few PCA components to represent each NC epoch: (the filter-bank index n is dropped for simplicity) wherex denotes the mean value vector of x in an epoch, and Q is the matrix of top eigenvectors for the maximal eigenvalues. Then we cluster the feature vectors of u from all the training data using the standard k-means algorithm into n k clusters, by minimizing the total intra-cluster variance.
where k denotes the index of the k-th cluster of NC epochs, and S k denotes the set of all feature vector indices belonging to the cluster andū k is their mean feature vector. We will examine different sub-state number n k in the experiment later.
To contrast a motor imagery class, say MI-1, against NC, we have to compare each and every sub-state in NC against MI-1. Again, we use the same CSP method that now maximizes the Rayleigh coefficient between a NC sub-state and MI-1.
With all the CSP filters constructed for the three class-pairs in every filter-bank, we will run them and aggregate the outcomes to form a joint, raw feature vector for each input EEG epoch. It can be seen that the size of the raw feature vector is partially determined by the number of clusters for learning NC. Then we consider how to further select a robust feature vector from the raw feature vector and map it to the desired output in the next subsection.

Information theoretic feature selection and regression
We denote the raw feature vector variable by A, and the selected feature vector variable by A η . And consider that generally the class label as a discrete random variable C with value from 1 to N c . We will use a η and c to represent a particular selected feature vector and its class label.
For modeling the dependency between A η and C, the mutual information is an important quantity that measures the mutual dependence of the two variables according to information theory (Papoulis, 1984). Mathematically, it is given by where H (A η ) denotes the entropy of the random feature vector, and H (A η | C) is the conditional entropy From the viewpoint of information theory, the optimum feature set is the one that carries the most mutual information about the www.frontiersin.org class label. This is exactly the idea of the maximum mutual information (MMI) criterion, which has been established as the basis for discriminative learning procedures in various machine learning techniques. Therefore, our system seeks the optimal feature vector A η which maximizes the mutual information.
The objective above involves joint probability density functions (PDFs) that need to be estimated from a given training data set. Please note that, to simplify the descriptions, we will omit the symbol η in the expressions unless otherwise specified. Using kernel density estimation, the PDF of a is given by where a i is a training sample of the feature vector, and ϕ is a smoothing kernel that takes a Gaussian form here.
where ψ is the covariance matrix that is assumed diagonal and estimated from training data according to the normal optimal smoothing strategy (Bowman and Azzalini, 1997). We can then adopt a method proposed in Viola and Wells (1997) to approximate the entropy H (A) with a given set of samples.
Combining the above equations, the entropy of the random vector A can be approximated by The within-class entropy H (A | c) can be similarly estimated with the training samples from the class c only.
With the above developments, we can compute the mutual information estimate for any subset of features. And we will select the subset which yields the largest mutual information estimate. Now we consider the mapping from the features to the desired outputs of class labels as a regression problem. The features are linearly normalized to the range [−1 1] using their empirical upper and lower bounds. To account for possibly non-linearity in the mapping, we employ a generalized regression neural network with non-linear hidden neurons. Briefly, the network consists of three layers of neurons: the second (hidden) layer contains radial basis function neurons, while the third layer has linear neurons with normalized inputs (Wasserman, 1993). This mapping can be cast as a general regression neural network (GRNN). Here in this work we use the SD of all the training set features for the spread parameter that basically defines the kernel width for the radial basis functions.

Post-processing and optimization
The above procedure addresses the problem of predicting class labels from individual EEG segments. On the other hand, postprocessing of the predicted label series allows us to explore information in the dynamics of brain activities and corresponding EEG observations. Particularly, we introduce a time-delay neural network (TDNN; Clouse et al., 1997). Here the network has three layers: the input layer acts as a 4-s buffer of epoch-based label predictions from training data; the output is the sequence of desired true label sequence in the same time frame; and the hidden layers consists of a few neurons with radial basis transfer functions. Again, this network can be cast as a GRNN.
There are two parameters to be optimized during system calibration, including the time interval for motor imagery segment extraction for training; the number of sub-states for NC segments clustering (hereafter n NC ). It is worthwhile to note that the training data usually are still trial-based for calibration, and we may take advantage of the timing information to extract most effective time intervals related to motor imagery EEG. Particularly, the system will examine the following time intervals: four 2-s-long intervals starting at 0.5, 1, 1.5, and 2 s from the trial-start cue, and three 2.5-s-long intervals starting at 0.5, 1, 1.5 s from the cue. For the submission to competition, we used a 5-fold cross-validation to optimize the two parameters in a subject-dependent manner.
Since the true class labels are from the limited number set {−1, 0, 1}, the output of the system is finally cut to the range of [−1 1]. This, however, did not produce any significant effects in the performance measure in our tests below.

RESULTS
First we would like to present the results of cross-validation before post-processing. For the sake of computational efficiency, the prediction of class label is performed every 0.1 s. Figure 2 plots the mean-square-error (MSE) of regression with respect to the selection of time interval T itv for training the feature extractor (see Section 1). Note that for each subject and each time interval, the MSE presented in the graph is the lowest MSE with the optimal T itv . Interestingly, it can be seen that all of the lower MSEs were achieved by 2.5-s-long intervals instead of 2-slong ones. Overall, the interval of 1.5-4.0 s gave rise to the smallest regression error. Nonetheless, it seems that the regression error is not very sensitive to the selection of time interval, as the variation of the error on each subject's data is small across different time intervals. Figure 3 plots MSE with respect to n k : the number of NC substates (modes, or clusters) in calibration. Note that for each subject and each n k , the MSE presented in the graph comes from the time interval with the lowest MSE. And, there is a slightly decrease in the error on two subjects ("a" and "f ") by 2-mode NC learning. Now let us assess the post-processing technique. With the best time interval and NC sub-state number obtained from the study Frontiers in Neuroscience | Neuroprosthetics  above, the MSEs before and after the post-processing are compared in Table 1. It can be seen that on every subject, the technique effectively reduced MSE by 0.02 or 0.03, which is approximately 10% of the original MSE.
Finally, we would like to give the evaluation results, which we submitted to the competition, on the independent data sets (i.e., the evaluation sets in the competition) in Table 2. www.frontiersin.org

DISCUSSION
In summary, we have presented a computational method for motor imagery detection. It can extract and learn effective spatio-spectral features to discriminate between three EEG classes including noncontrol and 2 motor imagery classes. It is expected that the use of dwell and refractory periods (c.f. Townsend et al., 2004) may further improve the performance. Furthermore, comparing the performance on the calibration data and that on the evaluation data, it seems that the performance degrades significantly. And this is due to the effect of session-to-session transfer that often cause considerable changes in the EEG characteristics, especially in the competition data where the calibration and the evaluation sessions employed different protocols. Finally, since this method uses an autonomous learning framework such that it does not rely on ad hoc tuning, it can serve as a favorable baseline for future research in motor imagery detection.


Revisiting Timed Specification Theories: A Linear-Time Perspective

We consider the setting of component-based design for real-time systems with critical timing constraints. Based on our earlier work, we propose a compositional specification theory for timed automata with I/O distinction, which supports substitutive refinement. Our theory provides the operations of parallel composition for composing components at run-time, logical conjunction/disjunction for independent development, and quotient for incremental synthesis. The key novelty of our timed theory lies in a weakest congruence preserving safety as well as bounded liveness properties. We show that the congruence can be characterised by two linear-time semantics, timed-traces and timed-strategies, the latter of which is derived from a game-based interpretation of timed interaction.

Introduction
Component-based design methodologies can be encapsulated in the form of compositional specification theories, which allow the mixing of specifications and implementations, admit substitutive refinement to facilitate reuse, and provide a rich collection of operators. Previously [1], we developed a linear-time specification theory for reasoning about untimed components that interact by synchronisation of input and output (I/O) actions, inspired by interface automata [2]. Models can be specified operationally by means of transition systems augmented by an inconsistency predicate on states, or declaratively using traces. The theory admits non-determinism, a refinement preorder based on traces, and the operations of parallel composition, conjunction and quotient. The refinement is strictly weaker than alternating simulation and is actually the weakest precongruence preserving inconsistent states. This implies that our refinement is substitutive, meaning component A refines component B iff A can replace B in any environmental context without introducing additional errors.
In this paper we target component-based development for real-time systems with critical timing constraints. We formulate a timed extension of the lineartime specification theory of [1], by allowing for both operational descriptions of components, as well as declarative specifications based on traces. Our operational models are based on a variant of timed automata with I/O distinction (although we do not insist on input-enabledness, cf [3]), augmented by two special states: ⊥ for safety and bounded-liveness errors, and for timestop. Tracebased declarative specifications are shown to be a suitable semantic domain for the operational models. In addition to timed-trace semantics, we present timedstrategy semantics, which coincides with the former but relates our work closer to the timed-game frameworks used by [4] and [5]. The substitutive refinement of our framework gives rise to the weakest congruence preserving ⊥, and is shown to coincide across all our formalisms.
Amongst notable works in the literature, we briefly mention a theory of timed interfaces [5] and a theory of timed specifications [4]. Timed interface theory contributes a framework based on timed games to formalise notions such as interfaces and compatibility, and also provides a parallel composition operator. However, the work cannot be considered a specification theory as it does not deal with the notion of refinement for component substitution or the operations of conjunction, disjunction and quotient. In this respect, [4] provides a complete theory; however, the refinement is a timed version of the alternating simulation originally defined for interface automata [2]. Consequently, it is too strong for determining when a component can be safely substituted with another (cf the example in Figure 3).
Outline. In Section 2 we introduce timed I/O automata, their semantic mapping to timed I/O transition systems, and supply the operational definitions for the operations of parallel composition, conjunction, disjunction and quotient. In Section 3 we use the timed-game framework to introduce timed-strategy semantics, which we relate to the operational framework. Similarly in Section 4, we present timed-trace semantics and relate these to the operational definitions. Section 5 discusses related work, and finally Section 6 concludes.

Formal Framework
In this section we introduce timed I/O automata, timed I/O transition systems and a semantic mapping from the former to the latter. Timed I/O automata are compact representations of timed I/O transition systems. Our theory will be developed using timed I/O transition systems, which are endowed with a richer repertoire of semantic machinery.

Timed I/O Automata
Clock constraints. Given a set X of real-valued clock variables, a clock constraint over X , cc : CC (X ), is a boolean combination of atomic constraints of the form x d and x − y d where x , y ∈ X , ∈ {≤, <, =, >, ≥}, and d ∈ N.
A clock valuation over X is a map t that assigns to each clock variable x in X a real value from R ≥0 . We say t satisfies cc, written t ∈ cc, if cc evaluates to true under valuation t. t + d denotes the valuation derived from t by increasing the assigned value on each clock variable by d ∈ R ≥0 time units. t[rs → 0] denotes the valuation obtained from t by resetting the clock variables in rs to 0. Sometimes we use 0 for the clock valuation that maps all clock variables to 0. -C ⊆ X is a finite set of clock variables -A (= I ∪ O) is a finite alphabet, where I and O are disjoint sets of input actions and output actions respectively -L is a finite set of locations l 0 ∈ L is the initial location -AT ⊆ L × CC (C ) × A × 2 C × L is a set of action transitions -Inv : L → CC (C ) and coInv : L → CC (C ) assign invariants and coinvariants to states, each of which is a downward-closed clock constraint.
We use l , l , l i to range over L and use l g,a,rs − −−− → l as a shorthand for (l , g, a, rs, l ) ∈ AT . g : CC (C ) is the enabling guard of the transition, a ∈ A the action, and rs the subset of clock variables to be reset.
Our TIOAs are similar to existing variants of timed automata with input/output distinction, except for the introduction of co-invariants and noninsistence on input-enabledness. While invariants specify the bounds beyond which time may not progress, co-invariants specify the bounds beyond which the system will time-out and enter error states. Our TIOAs can be used to describe both the assumptions made by the component on the inputs, together with the guarantees provided by the component on the outputs. Such assumptions and guarantees can be time constrained: guards on output transitions express safety timing guarantees, while guards on input transitions express safety timing assumptions; invariants (urgency) express liveness timing guarantees on outputs while co-invariants (time-out) express liveness timing assumptions on inputs.
When components interact together, we check whether the guarantees they provide meet the assumptions they make on each other. If not, there are two types of errors: -An input arrives in a state and at a time when it is not expected (i.e. not satisfying the guards on the input transitions). This is a safety error. -An input does not arrive in a state within a time bound (specified by a co-invariant) as expected. This is a bounded-liveness error.
Example. Figure 1 depicts TIOAs representing a job scheduler together with a printer controller. The invariant at location A of the scheduler forces a boundedliveness guarantee on outputs in that location. As time must be allowed to progress beyond t = 100, the start action must be fired within the range 0 ≤ t ≤ 100. After start has been fired, the clock x is reset to 0 and the scheduler waits (possibly indefinitely) for the job to finish. If the job does finish, the scheduler is only willing for this to take place between 5 ≤ t ≤ 8 after the job started (safety assumption), otherwise an unexpected input error will be thrown. The controller waits for the job to start, after which it will wait exactly 1 time unit before issuing print (forced by the invariant y ≤ 1 on state 2 and the guard y = 1). The controller now requires the printer to indicate the job is printed within 10 time units of being sent to the printer, otherwise a time-out error on inputs will occur (co-invariant y ≤ 10 in state 3 as liveness assumption). After the job has finished printing, the controller must indicate to the scheduler that the job has finished within 5 time units.

Timed Actions and Words
In this section we introduce some notation relating to timed actions and timed words that will be of use to us in later sections.
Timed actions. For a set of input actions I and a set of output actions O, define tA = I O R >0 to be the set of timed actions, tI = I R >0 to be the set of timed inputs, and tO = O R >0 to be the set of timed outputs. We use symbols like α, β, etc. to range over tA.
Timed words. A timed word (ranged over by w , w , w i etc.) is a finite mixed sequence of positive real numbers (R >0 ) and visible actions such that no two numbers are adjacent to one another. For instance, 0.33, a, 1.41, b, c, 3.1415 is a timed word denoting the observation that action a occurs at 0.33 time units, then another 1.41 time units lapse before the simultaneous occurrence of b and c, which is followed by 3.1415 time units of no event occurrence. The empty word is denoted by .
Operations on timed words. We use last(w ) to denote the last element in the sequence w , and l (w ) to indicate the length, which is obtained as the sum of all the reals in w . Concatenation of timed words w and w is obtained by appending w onto the end of w and coalescing adjacent reals (summing them). For instance, a, 1.41 0.33, b, 3.1415 = a, (1.41 + 0.33), b, 3.1415 = a, 1.74, b, 3.1415 . Prefix/extension are defined as usual by concatenation, and we use ≤ for the prefix partial order. We write w tA 0 for the projection of w onto timed alphabet tA 0 , which is defined by removing from w all actions not inside tA 0 and coalescing adjacent reals.

Semantics as Timed I/O Transition Systems
The semantics of TIOAs are given as timed I/O transition systems, which are a special class of infinite labelled transition systems. The states of the TIOTS for a TIOA capture the configurations of the automaton, i.e. its location and clock valuation. Therefore, each state of the TIOTS is a pair drawn from L × R C , which we refer to as the set of plain states, denoted P . In addition, we introduce two special states ⊥ and , which are required for the semantic mapping of disabled inputs/outputs, invariants and co-invariants.
⊥ is called the inconsistent state, representing safety and bounded-liveness errors.
is the so-called timestop state, representing the magic moment from which no error can occur. 1 An intuitive way to understand and ⊥ is from an input/output game perspective. The component controls output and delay while the environment controls input. ⊥ is the losing state for the environment. So a disabled input at a state p is equated to an input transition from p to ⊥. is the losing state for the component. So a disabled output/delay at p is equated to an output/delay transition from p to . Thus we can have two semantics-preserving transformations on TIOTSs.
The ⊥-completion of a TIOTS P, denoted P ⊥ , adds an a-labelled transition from p to ⊥ for every p ∈ P P and a ∈ I s.t. a is not enabled at p. ⊥completion will make a TIOTS input-receptive, i.e. input-enabled at all states. The -completion of a TIOTS P, denoted P , adds an α-labelled transition from p to for every p ∈ P P and α ∈ tO s.t. α is not enabled at p.
Furthermore, for technical convenience (e.g. ease of defining time additivity), the definition of TIOTSs requires that 1) is a quiescent state, i.e. a state in which the set of outgoing transitions are all self-loops, one for each d ∈ R >0 , and 2) ⊥ is a chaotic state, i.e. a state in which the set of outgoing transitions are all self-loops, one for each α ∈ tA. The set of all possible states is denoted S = P {⊥, }. We use p, p , p i to range over P while s, s , s i range over S .
The transition relation → of the TIOTS is derived from the execution semantics of the TIOA.
Note that our semantics tries to minimise the use of transitions leading to /⊥ states. Thus there are no delay transitions leading to . This creates implicit timestops, which we capture using the concept of semi-timestop (i.e. semi-). We say a plain state p is a semi-iff 1) all output transitions enabled in p or any of its time-passing successors lead to the state, and 2) there exists Thus a semi-is a state in which it is impossible for the component to avoid the timestop without suitable inputs from the environment.
− → s for some s. In the sequel of this paper we only consider TIOTSs that are time-additive.
We say a TIOTS is deterministic iff there is no ambiguous transition in the Given a TIOTS P, a timed word can be derived from a finite execution of P by extracting the labels in each transition and coalescing adjacent reals. The timed words derived from such executions are called traces of P. We use tt, tt , tt i to range over the set of traces and use s 0 tt = ⇒ s to denote a finite execution that produces trace tt and leads to s.

Operational Specification Theory
In this section we develop a compositional specification theory for TIOTSs based on the operations of parallel composition , conjunction ∧, disjunction ∨ and quotient %. The operators are defined via transition rules that are a variant on synchronised product.
Parallel composition yields a TIOTS that represents the combined effect of its operands interacting with one another. The remaining operations must be explained with respect to a refinement relation, which corresponds to safesubstitutivity in our theory. A TIOTS is a refinement of another if it will work in any environment that the original worked in without introducing safety or bounded-liveness errors. Conjunction yields the coarsest TIOTS that is a refinement of its operands, while disjunction yields the finest TIOTS that is refined by both of its operands. The operators are thus equivalent to the join and meet operations on TIOTSs 2 . Quotient is the adjoint of parallel composition, meaning that P 0 %P 1 is the coarsest TIOTS such that (P 0 %P 1 ) P 1 is a refinement of P 0 . Let 1} be two TIOTSs that are both ⊥ and -completed, satisfying (wlog) S 0 ∩ S 1 = {⊥, }. The composition of P 0 and P 1 under the operation ⊗ ∈ { , ∧, ∨, %}, written P 0 ⊗ P 1 , is only defined when certain composability restrictions are imposed on the alphabets of the TIOTSs. P 0 P 1 is only defined when the output sets of P 0 and P 1 are disjoint, because an output should be controlled by at most one component. Conjunction and disjunction are defined only when the TIOTSs have identical alphabets (i.e. O 0 = O 1 and I 0 = I 1 ). This restriction can be relaxed at the expense of more cumbersome notation, which is why we focus on the simpler case in this paper.
For the quotient, we require that the alphabet of P 0 dominates that of P 1 (i.e.
, in addition to P 1 being a deterministic TIOTS. As quotient is a synthesis operator, it is difficult to give a definition using just statelocal transition rules, since quotient needs global information of the transition systems. This is why we insist on P 1 being deterministic 3 .
-→ is the smallest relation containing → 0 ∪ → 1 , and satisfying the rules: We adopt the notation of s 0 ⊗ s 1 for states, where the associated interpretation is supplied in Table 1. Furthermore, given two plain states p i = (l i , t i ) for i ∈ {0, 1}, we define p 0 × p 1 = ((l 0 , l 1 ), t 0 t 1 ). preorder. Similarly for disjunction. Quotient is the adjoint of parallel composition. If the second component state does not refine the first, the quotient will try to rescue the refinement by producing (so that its composition with the second will refine the first). If the second component state does refine the first, the quotient will produce the least refined value so that its composition with the second will not break the refinement.
An environment for a TIOTS P is any TIOTS Q such that the alphabet of Q is complementary to that of P, meaning I P = O Q and O P = I Q . Refinement in our framework corresponds to contextual substitutability, in which the context is an arbitrary environment.
Definition 5. Let P imp and P spec be TIOTSs with identical alphabets. P imp refines P spec , denoted P spec P imp , iff for all environments Q, P spec Q is ⊥-free implies P imp Q is ⊥-free. We say P imp and P spec are substitutively equivalent, i.e. P spec P imp , iff P imp P spec and P spec P imp .

It is obvious that
induces an equivalance on TIOTSs and no equivalence that preserves the ⊥ state can be weaker than . In the sequel we will give two concrete characterisations of and show that is also a congruence w.r.t. the parallel composition, conjunction, disjunction and quotient operators.
The operational definition of quotient requires that P 1 is determinised, which can be accomplished by a modified subset construction procedure on (P ⊥ 1 ) . If the current state subset S 0 contains ⊥, it reduces S 0 to ⊥; if ⊥ / ∈ S 0 = { }, it reduces S 0 by removing any potential in S 0 . As expected, the determinisation of P, denoted P D , is substitutively equivalent to P. Proposition 1. Any TIOTS is substitutively equivalent to a deterministic TIOTS.
Equipped with determinisation, quotient is a fully defined operator on any pair of TIOTSs. Furthermore, we can give an alternative (although substitutively equivalent) formulation of quotient as the derived operator (P ¬ 0 P 1 ) ¬ , where ¬ is a mirroring operation that first determinises its argument, then interchanges the input and output sets, as well as the and ⊥ states.
Example. Figure 2 shows the parallel composition of the job scheduler with the printer controller. In the transition from B 4 to A1, the guard combines the effects of the constraints on the clocks x and y. As finish is an output of the controller, it can be fired at a time when the scheduler is not expecting it, meaning that a safety error will occur. This is indicated by the transition to ⊥ when the guard constraint 5 ≤ x ≤ 8 is not satisfied.

Timed I/O Game
Our specification theory can be understood from a game theoretical point of view. It is an input-output game between a component and an environment that   -Given two TIOTSs P and Q with identical alphabets, we say P is a partial unfolding [6] of Q if there exists a function f from S P to S Q s.t. 1) f maps to , ⊥ to ⊥, and plain states to plain states, 2) f (s 0 P ) = s 0 Q , and 3) p -We say a simple path L is a run of P if L is a partial unfolding of P.

Strategies.
A strategy G is a deterministic tree TIOTS s.t. each plain state in G is ready to accept all possible inputs by the environment, but allows a single move (delay or output) by the component, i.e. eb G (p) = I mv G (p) s.t. mv G (p) = {a} for some a ∈ O or {} ⊂ mv G (p) ⊆ R >0 , where eb G (p) denotes the set of enabled timed actions in state p of LTS G, and mv G (p) denotes the unique component move allowed by G at p.
A TIOTS P contains a strategy G if G is a partial unfolding of (P ⊥ ) . The set of strategies 4 contained in P is denoted stg(P). Since it makes little sense to distinguish strategies that are isomorphic, we will freely use strategies to refer to their isomorphism classes and write G = G to mean G and G are isomorphic.
Let us give some examples in Figure 3. For the sake of simplicity we use two untimed transition systems P and Q, which have identical alphabets I = {e, f } and O = {a, b, c}, to illustrate the idea of strategies. The transition systems use solid lines while strategies use dotted lines. Plain states are unmarked while the and ⊥ states are marked by and ⊥ resp. 5 We show four strategies of P and  two strategies of Q on the right hand side of P and Q resp. in Figure 3. (They are not the complete sets of strategies for P and Q.) Note that the strategies 3 and 4 own their existence to the completion.
Comparing strategies. When the game is played, the component tries to avoid reaching while the environment tries to avoid reaching ⊥. Different strategies in stg(P) vary in their effectiveness to achieve the objective. Such effectiveness can be compared if two strategies closely resemble each other: we say G and G are affine if s 0 G tt = ⇒ p and s 0 G tt = ⇒ p implies mv G (p) = mv G (p ). Intuitively, it means G and G propose the same move at the 'same' states. For instance, the strategies 1, 3 and A in Figure 3 are pairwise affine and so are the strategies 2, 4 and B .
Given two affine strategies G and G , we say G is more aggressive than G , denoted G G , if 1) s 0 = ⇒ . Intuitively, it means G can reach ⊥ faster but slower than G . forms a partial order over stg(P), or more generally, over any set of strategies with identical alphabets. For instance, strategy A is more aggressive than 1 and 3, while strategy B is more aggressive than 2 and 4.
When the game is played, the component P prefers to use the maximally aggressive strategies in stg(P) 6 . Thus two components that differ only in nonmaximally aggressive strategies should be equated. We define the strategy semantics of component P to be [P] s = {G | ∃ G ∈ stg(P) : G G }, i.e. the upward-closure of stg(P) w.r.t. .
Game rules. When a component strategy G is played against an environment strategy G , at each game state (i.e. a product state p G × p G ) G and G each propose a move (i.e. mv G (p G ) and mv G (p G )). If one of them is a delay and the other is an action, the action will prevail. If both propose delay moves (i.e. mv G (p G ), mv G (p G ) ⊆ R ≥0 ), the smaller one (w.r.t. set containment) will prevail.
Since a delay move proposed at a strategy state is the maximal set of possible delays enabled at that state, the next move proposed at the new state after firing the set must be an action move (due to time additivity). Thus a play cannot have two consecutive delay moves.
If, however, both propose action moves, there will be a tie, which will be resolved by tossing the coin. For uniformity's sake, the coin can be treated as a special component. A strategy of the coin is a function h from tA * to {0, 1}. We denote the set of all possible coin strategies as H .
A play of the game can be formalised as a composition of three strategies, one each from the component, environment and coin, denoted G P h G Q . At a current game state p P × p Q , if the prevailing action is α and we have p P α − → s P and p Q α − → s Q , then the next game state is s P s Q . The play will stop when it reaches either or ⊥. The composition will produce a simple path L that is a run of P Q. Since P Q gives rise to a closed system (i.e. the input alphabet is empty), a run of P Q is a strategy of P Q. This is crucial since it reveals that strategy composition of P and Q is closely related to their parallel composition: stg(P Q) = {G P h G Q | G P ∈ stg(P), G Q ∈ stg(Q) and h ∈ H }.
Parallel composition. Strategy composition, like component (parallel) composition, can be generalised to any pair of components P and Q with composable alphabets. That is, O P ∩ O Q = {}. For such P and Q, G P h G Q gives rise to a tree rather than simple path TIOTS. That is, at each game state p P ×p Q , besides firing the prevailing α ∈ tO P ∪tO Q , we need also to fire 1) all the synchronised inputs, i.e. e ∈ I P ∩ I Q , and reach the new game state s P s Q (assuming p P e − → s P and p Q e − → s Q ) and 2) all the independent inputs, i.e. e ∈ (I P ∪ I Q ) \ (A P ∩ A Q ), and reach the new game state s P × p Q or p P × s Q . It is easy to verify that G P h G Q is a strategy of P Q.
Conjunction/disjunction. Besides strategy composition, strategy conjunction (&) and strategy disjunction (+) are also definable. They are binary operators defined only on pairs of affine strategies. We define G&G = G ∧ G and G + G = G ∨ G . Note that if G and G are not affine, G ∧ G and G ∨ G do not necessarily produce a strategy. For instance the disjunction of the strategies 1 and 2 in Figure 3 will produce a transition system that stops to output after the a transition.
Refinement. Strategy semantics induce an equivalence on TIOTSs. That is, P and Q are strategy equivalent iff [P] s = [Q] s . However, strategy equivalence is too fine for the purpose of substitutive refinement (cf Definition 5). For instance, transition systems P and Q in Figure 3 are substitutively equivalent, but are not strategy equivalent, because 1, 2, 3 and 4 are strategies of Q (due to upwardclosure w.r.t. ), but A and B are not strategies of P.
However, we demonstrate that substitutive equivalence is reducible to strategy equivalence providing we perform disjunction closure on strategies.
Lemma 1. Given a pair of affine component strategies G 0 and G 1 , G 0 h G and G 1 h G are ⊥-free for some environment strategy G and h ∈ H iff G 0 + G 1 h G is ⊥-free.
We say Π + is a disjunction closure of Π iff it is the least superset of Π s.t. G + G ∈ Π + for all pairs of affine strategies G, G ∈ Π + . It is easy to see the disjunction closure operation preserves the upward-closedness of strategy sets. Lemma 2. For -composable TIOTSs P and Q, Thus conjunction and disjunction are the join and meet operations and quotient produces the coarsest TIOTS s.t. (P 0 %P 1 ) P 1 is a refinement of P 0 . Lemma 6. For any TIOTS P, is a congruence w.r.t. , ∨, ∧ and % subject to composability.
Summary. Strategy semantics has given us a weakest ⊥-preserving congruence (i.e. [P] + s ) for timed specification theories based on operators for (parallel) composition, conjunction, disjunction and quotient. Strategy semantics captures nicely the game-theoretical nature as well as the operational intuition of the specification theories. However, in a more declarative manner, the equivalence can also be characterised by timed traces, as we see in the next section.

Declarative Specification Theory
In this section, we develop a compositional specification theory based on timed traces. We introduce the concept of a timed-trace structure, which is an abstract representation for a timed component. The timed-trace structure contains essential information about the component, for checking whether it can be substituted with another in a safety and liveness preserving manner.
Given any TIOTS P = I , O, S , s 0 , → , we can extract three sets of traces from (P ⊥ ) : TP (plain traces) is a set of timed traces leading to plain states, TE (error traces) a set of timed traces leading to ⊥ and TM (magic traces) a set of timed traces leading to . The three sets contain sufficient but not necessary information for our substitutive refinement, which is designed to preserve ⊥ rather than . For instance, adding any trace tt ∈ TE to TP should not change the semantics of the component; similarly it is true for removing any trace tt ∈ TP from TM . Based on a slight abstraction of the three sets we can define a triple-trace structure as the semantics of P. Obviously, TE is extension-closed. TT is non-empty and prefix-closed. TR is prefix-closed and fully branching 7 w.r.t. TT (i.e. tt α ∈ TT for all tt ∈ TR and α ∈ tA). TT \ TR is time-extension closed (i.e. tt ∈ X ⇒ tt d ∈ X ) and any pair of traces from TT \ TR that are related by extension are related by time-extension.
From hereon let P 0 and P 1 be two TIOTSs with triple trace structures The substitutive refinement relation in Section 2.4 can equally be characterised by means of trace containment. Consequently, T T (P 0 ) can be regarded as providing an alternative encoding of the set [P 0 ] + s of strategies.
We are now ready to define the timed-trace structure semantics for the operators of our specification theory. Intuitively, the timed-trace semantics mimic the synchronised product of the operational definitions in Section 2.4. An important fact utilised in formulating these operations on traces is that for any trace tt ∈ tA * and TIOTS P, either tt is a trace of P or there is some prefix tt 0 of tt s.t. tt 0 is an error or magic trace of P.
Parallel composition. The idea behind parallel composition is that the projection of any trace in the composition onto the alphabet of one of the components should be a trace of that component.
The above says tt is an error trace if the projection of tt on one component is an error trace while the projection of tt on the other component is not a magic trace. tt is a realisable trace if tt is either an error trace or a plain trace. tt is a plain trace if the projection of tt on both components are plain traces. Finally, tt is a magic trace if its projection on one component is a magic trace, while the projection of all strict prefixes of tt on the other component is a plain trace.
Disjunction. From any composite state in the disjunction of two components, the composition should only be willing to accept inputs that are accepted by both components, but should accept the union of outputs. After witnessing an output enabled by only one of the components, the disjunction should behave like that component. Because of the way that ⊥ and work in Table 1, this loosely corresponds to taking the union of the traces from the respective components.
Essentially, tt is a magic trace if it is a magic trace on one component while one of its prefixes is a magic trace on the other component. The realisable and error traces are simply the union of the corresponding traces on P 0 and P 1 .
Conjunction. Similarly to disjunction, from any composite state in the conjunction of two components, the composition should only be willing to accept outputs that are accepted by both components, and should accept the union of inputs, until a stage when one of the component's input assumptions has been violated, after which it should behave like the other component. Because of the way that both ⊥ and work in Table 1, this essentially corresponds to taking the intersection of the traces from the respective components.
A trace tt is a magic trace if it is a magic trace on one of the components, and all strict prefixes of the trace are realisable by the other component. The realisable and error traces are simply the intersection of the corresponding traces on P 0 and P 1 .
Quotient. Quotient ensures its composition with the second component is a refinement of the first. Given the synchronised running of P 0 and P 1 , if P 0 is in a more refined state than P 1 , the quotient will try to rescue the refinement by taking as its state (so that its composition with P 1 's state will refine P 0 's). If P 0 is in a less or equally refined state than P 1 's, the quotient will take the worst possible state without breaking the refinement.
The above says tt is an error trace if either 1) tt is an error trace in P 0 , but the projection of any strict prefix of tt on P 1 is not an error trace, or 2) the projection of tt on P 1 is a magic trace, but no strict prefix of tt is a magic trace in P 0 . tt is a magic trace if either 1) tt is a magic trace in P 0 , but the projection of any prefix of tt is not a magic trace in P 1 , or 2) the projection of tt on P 1 is an error trace, but no prefix of tt is an error trace in P 0 .
Mirroring of triple trace structures is straightforward: T T (P 0 ) ¬ = (O 0 , I 0 , TT 0 , TT 0 \ TE 0 , TT 0 \ TR 0 ). This is because dealing with traces means we have implicit determinism, so we can skip the determinisation step. Consequently, quotient can also be defined as the derived operator (T T (P 0 ) ¬ T T (P 1 )) ¬ .

Comparison with Related Works
Based on linear-time, our timed theory owes much to the pioneering work of trace theories in asynchronous circuit verification, such as Dill's trace theory [7]. Our mirror operator is essentially a timed extension of the mirror operator from asynchronous circuit verification. The definition of quotient based on mirroring (for the untimed case) was first presented by Verhoeff as his Factorisation Theorem [8].
Our work is also deeply influenced by the work of [5] on timed games, with some modifications. Firstly, a TIOTS is regarded as a set of component strategies, rather than a timed game graph. We adopt most of the game rules in [5], except that, due to our requirement that proposed delay moves are maximal delays allowed by a strategy, a play cannot have consecutive delay moves. This enables us to avoid the complexity of time-blocking strategies and blame assignment, but does not ensure non-Zenoness 8 . Secondly, we do not use timestop/semi-timestop to model time errors (i.e. bounded-liveness errors). Rather, we introduce the explicit inconsistent state ⊥ to model both time and immediate (i.e. safety) errors.
Timestop is used to model the magic state, which can simplify the definition of parallel, conjunction and quotient and enables us to avoid the complexity of having two transition relations and well-formedness of timed interfaces.
Last but not least, our work is related to [4], as both devise a complete timed specification theory. The major differences lie in the use of timed alternating simulation as refinement in [4], while ours is linear-time. An advantage of our work is that refinement is the weakest congruence preserving inconsistency, while beneficial in [4] is the algorithmic efficiency of branching-time simulation checking. Moreover, [4] has fully implemented the timed-game algorithms.

Conclusions
We have formulated a rich compositional specification theory for components with real-time constraints based on a linear-time notion of substitutive refinement. The operators of hiding and renaming can also be defined, according to our past experiences [13]. We believe that our theory can be reformulated as a timed extension of Dill's trace theory [7]. Future work will include an investigation of realisability and assume-guarantee reasoning.