# A Novel Energy-Efficient Cross-Layer Design for Scheduling and Routing in
6TiSCH Networks

Ahlam Hannachi, Wael Jaafar, Salim Bitam, and Nabil Ouazene  Ahlam Hannachi
and Salim Bitam are with the Department of Computer Science, University of
Biskra, Biskra, Algeria (e-mail<EMAIL_ADDRESS>s.bitam@univ-
biskra.dz). Wael Jaafar is with the Department of Software and IT Engineering,
École de Technologie Supérieure (ETS), Montreal, Canada (e-mail:
wael.jaafar@etsmtl.ca). Nabil Ouazene is with the Laboratory of Applications
of Mathematics to Computing and Electronics, University of Batna 2, Batna,
Algeria (e-mail: nabil.ouazene@univ-batna2.dz).

###### Abstract

The 6TiSCH protocol stack was proposed to ensure high-performance
communications in the Industrial Internet of Things (IIoT). However, the lack
of sufficient time slots for nodes outside the 6TiSCH’s Destination Oriented
Directed Acyclic Graph (DODAG) to transmit their Destination Advertisement
Object (DAO) messages and cell reservation requests significantly hinders
their integration into the DODAG. This oversight not only prolongs the
device’s join time but also increases energy consumption during the network
formation phase. Moreover, challenges emerge due to the substantial number of
control packets employed by both the 6TiSCH Scheduling Function (SF) and
routing protocol (RPL), thus draining more energy resources, increasing medium
contention, and decreasing spatial reuse. Furthermore, an SF that overlooks
previously allocated slots when assigning new ones to the same node may
increase jitter, and more complications ensue when it neglects the state of
the TSCH queue, thus leading to packet dropping due to queue saturation.
Additional complexity arises when the RPL disregards the new parent’s schedule
saturation during parent switching, which results in inefficient energy and
time usage. To address these issues, we introduce in this paper novel
mechanisms, strategically situated at the intersection of SF and RPL that are
designed to balance the control packet distribution and adaptively manage
parent switching. Our proposal, implemented within the 6TiSCH simulator,
demonstrates significant improvements across vital performance metrics, such
as node’s joining time, jitter, latency, energy consumption, and amount of
traffic, in comparison to the conventional 6TiSCH benchmark.

###### Index Terms:

Industrial IoT, TSCH, 6TiSCH, SF, RPL network formation, scheduling, routing.

##  I Introduction

Amid current digital advancements, the presence of the Internet permeates
nearly every aspect of our lives, spanning domains from smart cities and
agriculture to industry and medicine. At the heart of this expansive reach is
the Internet of Things (IoT), poised to interlink billions of devices with the
traditional Internet. Many of these devices, while cost-effective, come with
inherent constraints related to energy, computation, storage, and
communication capabilities. The overarching success of IoT hangs on these
resource-limited objects’ ability to operate reliably and sustainably while
maintaining connectivity [1].

In IoT systems, energy efficiency stands at the forefront of technological
challenges, with researchers globally striving to meld enhanced performances
with reduced or equilibrated energy consumption. Specifically, the medium
access control (MAC) layer is pivotal in controlling the energy drain since it
ensures the network’s reliability by regulating access to the shared
communication medium; typically a wireless channel. The IEEE.802.15.4e
standard has been crafted to cater to industrial applications. It introduces
different operational MAC modes, such as Time Slotted Channel Hopping (TSCH),
Deterministic and Synchronous Multi-channel Extension (DSME), and Low Latency
Deterministic Network (LLDN) [2], [3]. Among these modes, TSCH stands out as
the leading one and is included in several standards such as ANSI/ISA100.11a
and WirelessHart IEC62591 [4], [5].

Figure 1: TSCH schedule of nodes A, B, and Y.

The core concept of TSCH lies in the integration of Time Division Multiple
Access (TDMA) with channel hopping techniques. This fusion results in a dual
medium access mechanism based on time and channel. In this context, “time” is
referred to by time slot, and “frequency” is referred to by channel. A time
slot within a specific channel is defined as a cell. For every time slot,
shown in Fig. 1, an IoT node specifies its activity in a given channel,
choosing to transmit to a node (Tx), receive from a node (Rx), or sleep, with
the latter not requiring any specific cell selection [6]. The efficiency of
the TSCH schedule is reflected in the consistent actions of devices. Thus,
each device has to judiciously select the combination of (state, cell,
neighbor), where the state$\in\left\\{Tx,Rx,Sleep\right\\}$. For instance, in
node A’s TSCH schedule in Fig. 1, the cell (slot 1, channel 2) is set for
receiving from node B. Conversely, in node B’s schedule, the same cell is
allocated for transmitting to node A. Moreover, A is set to sleep mode during
time slots 2 and 4, while B enters sleep mode in time slot 2.

As the transition from IPv4 to IPv6 becomes more pressing due to IPv4’s
address limitations and IPv6’s advanced capabilities, compatibility issues
have risen. To tackle them, several bridging technologies have been developed
to cope with the protocols’ distinct structures and behaviors. The latter
became the de facto tools of the protocols’ coexistence, and thus ensuring a
smooth transition from IPv4 towards IPv6. Within this framework, the Internet
Engineering Task Force (IETF) took the initiative of forming the 6TiSCH
working group, aiming to establish protocols that facilitate the
implementation of IPv6 over the TSCH mode. The proposed protocol stack is
depicted in Fig. 2. Specifically, the 6top layer holds the Scheduling Function
(SF) responsible for TSCH cell allocation and traffic demands/performance
metrics monitoring. The conventional SF is the Minimal Scheduling Function
(MSF) [7].

Moreover, the Routing Protocol for Low-Power and Lossy Networks (RPL) layer,
which encompasses the Objective Function (OF), plays a pivotal role in
establishing and sustaining the routing topology. It does so by constructing a
Destination-Oriented Directed Acyclic Graph (DODAG), a structured network
topology that enables optimized routing. Through DODAG, RPL efficiently
propagates routing information, tailoring pathways based on various critical
metrics to meet the specific demands of the network. This mechanism is
fundamental in addressing the complex routing challenges inherent in low-power
and lossy networks, particularly within the Industrial IoT (IIoT) framework.
The latter allows IoT nodes to discover their best potential parents and
establish a route towards the network root. The conventional OF of RPL is the
Objective Function Zero (OF0) [8].

Figure 2: The 6TiSCH protocol stack (Adapted from [9])

Both RPL and 6top utilize control packets to maintain and optimize the IoT
network’s operation, namely DODAG Information Solicitation (DIS), DODAG
Information Object (DIO), and Destination Advertisement Object (DAO) for RPL,
and 6P transactions in 6top. The relation between RPL and 6top is reciprocal.
For instance, while RPL alters the routing topology by changing a node’s
parent, 6top dynamically adjusts the scheduled cells to accommodate this
change, i.e., releasing assigned cells to the old parent and negotiating new
cell allocation with the new parent. RPL determines the next-hop parent for
data transmission based on its topology, while 6top ensures cell availability
for this specific transmission according to the TSCH schedule.

However, this cross-layer interaction between the RPL and 6top layers (shown
by red lines in Fig.2) controls implicitly very critical performance factors.
To the best of our knowledge, very few works have investigated the nature of
this interaction. Consequently, we propose a novel approach that improves the
RPL-6top interaction, while aiming to reduce energy consumption and enhance
the spatial reuse of resources. Our contributions can be summarized as
follows:

  * •

We aim to reduce the number of control packets by merging two types of control
packets and to simplify the cell negotiation and assignment process.
Specifically, upon parent switching, RPL control messages will be used for
cell allocation rather than transmitting separate 6P transaction messages.

  * •

We eliminate instances of parent switching failures, prevalent when the new
parent schedule is at full capacity. This is achieved by embedding a list of
available slots within the DODAG Information Object (DIO) packets.
Consequently, a transition to a parent offering a superior rank will only be
executed if the list of its free slots satisfies the requirement of the new
child.

  * •

We improve the jitter and latency by adopting proximity-based slot selection
rather than the conventional random selection.

  * •

To preempt packet drop events that occur when the TSCH queue reaches its
limit, we propose setting a predefined queue size to solicit new cells from
the parent node.

  * •

Through extensive simulations over the 6TiSCH simulator, we prove the superior
performance of our proposed method in terms of joining time, energy
consumption, traffic load, delivery ratio, and jitter, compared to the MSF
benchmark.

The structure of this paper is outlined as follows: Section II presents a
comprehensive overview of the research context and discusses related works.
Section III delves into the intricate details of the proposed solution.
Through simulations, section IV evaluates the performances of our approach and
compares them to those of the MSF benchmark. Finally, Section V concludes the
paper.

##  II Context and Related Work

In the field of 6TiSCH networks, ongoing research is dedicated to developing
effective, reliable, and energy-efficient communication strategies. This often
involves examining and refining existing standards. In this section, we first
outline the key conventional protocols within the 6TiSCH stack, thus providing
the foundational framework for researchers. Then, we highlight recent related
works.

###  II-A Context: 6TiSCH Background

The 6TiSCH standards lay the groundwork for the operation of IIoT networks. We
present their main components relevant to our research, namely the Minimal
Configuration (MC) and the Minimal Scheduling Function (MSF). The MC details
the procedure for the node joining process [10], while the MSF defines the
scheduling mechanism, which outlines how communication cells are assigned and
released [11].

####  II-A1 Minimal Configuration

The MC defines the sequential steps that each node must undertake before
becoming operative and start transmitting data. Its steps, shown in Fig. 3,
are described as follows:

Figure 3: 6TiSCH network formation (MC).

  * •

Synchronization: Initially, each node is in a “pledge” state, passively
listening on randomly selected channels to anticipate Enhanced Beacons (EBs)
sent by nodes that are already part of the network, including the root node
that initiate the network formation. Upon the reception of a predefined number
(nb) of EB packets, a pledge synchronizes with the IoT network. For instance,
node A in Fig. 3 starts as a “pledge” at time $T$ and synchronizes at time
$T^{A}_{\rm Sync}$, where energy is expanded inefficiently during the interval
$T_{1}=T^{A}_{\rm Sync}-T$.

  * •

DODAG joining: After synchronization, if secure joining is required, the node
sends a “JOIN REQUEST” to the EB source and awaits the “JOIN RESPONSE”.
Otherwise, it directly sends a DODAG Information Solicitation (DIS) packet.
The DIS mechanism operates in two distinct modes:

    * –

Unicast DIS: In this mode, the receipt of a DIS message prompts the targeted
node to respond with a DIO message.

    * –

Broadcast DIS: In this mode, receiving a DIS message leads nodes to reset
their DIO trickle timers, thereby influencing the timing of subsequent DIO
message dissemination.

The DODAG joining is driven by the periodic multicasting of DIO messages by
nodes that have previously integrated the DODAG. DIO messages encapsulate
vital information regarding the DODAG, including its identification (ID),
version, and the node’s rank within the structure. DIO packets are delineated
into:

    * –

Periodic DIO: These messages are needed to establish and sustain the topology
of a DODAG. They provide critical information about the DODAG’s configuration,
including the sender’s rank and set of nodes utilized for parent selection
within the DODAG. The period of the DIO transmission is defined by the Trickle
algorithm, and according to the network’s consistency [12].

    * –

DIO as a response: It is a unicast response to a DIS packet sent by a
soliciting node asking for relevant routing information. For instance, node A
in Fig. 3 joins the DODAG at time $T^{A}_{\rm Join}$, where the time interval
$T_{2}=T^{A}_{\rm Join}-T^{A}_{\rm Sync}$ corresponds to the joining
procedure.

  * •

Cell allocation: Once synchronized, every node reserves its first cell,
referred to as the “minimal cell”, at time slot $0$ and channel $0$. Following
this initial reservation, when a node receives a DIO message, it selects a
parent among neighbors that advertised the DIO messages, according to the
RPL’s objective function. Then, the node triggers the process of cell
reservation with the new parent by initiating a 6P transaction and transmits a
DAO packet to the DAG root to indicate its parent, so that the DAG root can
build/update the topology tree. If the parent is the DAG root itself, then the
node sends both the 6P transaction and DAO packet to it. DAO packets are
destined to the DAG root to facilitate the update of the routing tree. They
come in two types:

    * –

Periodic DAO: It is transmitted periodically and used to maintain and refresh
the route information at the DAG root.

    * –

Parent Switching DAO: It is sent immediately after a parent-switching event to
swiftly notify the DAG root of changes in the routing topology.

In summary, nodes within 6TiSCH networks dynamically alter their operational
states throughout their life cycle, influenced by localized parameters and
network conditions. These transitions and state changes, illustrated in Fig.
4, embody the node’s continual attempt to preserve network connectivity,
synchronization, and data transmission integrity within the inherently
challenging communication environment of 6TiSCH networks. Notably, a node
reverts to its initial “pledge” state if it fails to receive any packet within
a predefined temporal window (Dy_sync), indicative of potential communication
disruptions or isolation from the network. Furthermore, if the node
experiences a linkage failure to its parent or while executing parent
switching (Dis_join), it regresses to a “synchronized” state, reflecting a
lack of association with the DODAG, yet maintaining synchronization with the
network’s time source.

Figure 4: Node status cycle.

####  II-A2 Scheduling Function

The scheduling function, SF, serves as a critical mechanism to allocate time
slots and channels within the TSCH mode. Its primary role is to manage the
scheduling of communications across the time-frequency domain. By dynamically
responding to the varying demands and conditions of the network, SF ensures
the judicious use of resources and maintains dependable communications. As
shown in Fig. 5, the SF conducts the cell reservation procedure using 6P
protocol transactions [13, 11], where the term “Requester” denotes a node in
need of reserving new cells from another node, referred to as the “Responder”.
In this process, the “Requester” initiates the cell reservation process. There
are two distinct approaches: the 2-step and 3-step methods. The key difference
lies in which node proposes the list of cells. In the 2-step transaction
illustrated in Fig. 5(a), the “Requester” proposes the list of cells, which
leads to exchanging two control packets, whereas in the 3-step transaction,
presented in Fig. 5(b), the “Responder” proposes the list of cells, thus
exchanging three control packets. Nevertheless, in both approaches, the
proposed cells are locked during the period $T_{\rm Lock}$.

The conventional SF of the 6TiSCH protocol stack is the Minimal Scheduling
Function, which follows a strategy of random cell selection, both in cell
proposal during 6P transactions and in cell reservation among the
proposed/available cells. Nevertheless, the nodes dynamically determine the
number of cells needed based on their traffic requirements and network
conditions. Although this stochastic approach provides adaptability to varying
traffic conditions, it may face challenges such as internal collisions,
particularly in dense networks. This raises concerns about its applicability
in scenarios with stringent Quality-of-Service (QoS) requirements.

(a) 2-Step 6P Transaction.

(b) 3-Step 6P Transaction.

Figure 5: An Example of 6P Transactions.

###  II-B Related Works

The conventional mechanisms, while foundational, have been scrutinized for
potential limitations and areas of enhancement. This section investigates the
existing literature for proposed methodologies, achievements, and gaps,
thereby contextualizing the motivation and novelty for our work.

Aiming to enhnace fundamental aspects of 6TiSCH, in particular the TSCH
configuration, authors of [14] adopted heterogeneous time slot length in MSF
to minimize communication latency. In [15], the authors introduced a heuristic
for parent selection in multi-PHY TSCH, showing that slot bonding can adapt to
different PHY data rates and achieve high Packet Delivery Ratios (PDRs).

From the perspective of network formation in 6TiSCH, the study conducted in
[16] evaluated the network formation procedure, and pinpointed numerous
limitations within its scope. Then, they proffered several recommendations and
established guidelines for parameter adjustments, tailoring them to several
network topologies. Also, authors of [17] addressed performance degradation
during 6TiSCH network formation. They highlighted the challenges faced by new
nodes joining the network, primarily due to the policy of allocating only one
shared cell per slotframe for control packet transmission. They particularly
emphasized the limitations of prioritizing EBs during network formation. Then,
they proposed an improved approach that effectively integrates control packets
into the network formation process.

Typically, research works that enhance SF fall into three categories: 1)
Centralized approaches such as [18, 19], 2) distributed approaches like [20],
and 3) autonomous approaches, e.g., [21]. In [18], the authors proposed a
fully centralized schedule for data gathering, where each node has to send the
allocation request to the DAG root and receive the answer, thus wasting a
considerable amount of energy and time. Alternatively, in [19], a central
entity is responsible for managing the scheduling of transmissions in the
network and uses a recurrent neural network (RNN) based model to predict the
next transmission slot for each node, aiming to reduce the latency of end-to-
end communications. However, a higher control packets overhead is generated.
Moreover, authors of [22] introduced an Adaptive Static Scheduling (ASS)
technique designed to enhance the network’s energy efficiency. Several works
delved into the scheduling strategies and network formation solutions. They
are comprehensively summarized in surveys, e.g., [23, 24]

Other works focused on optimizing the RPL’s objective function. For instance,
authors of [25] introduced a novel parent selection function that accounts for
the number of neighbors and traffic load and aiming to optimize the PDR, while
minimizing frequent parent switching. In [26], a Q-learning based algorithm is
proposed to refine the RPL’s objective function, while [27] proposed an early
parent switching scheme, based on the remaining queue capacity of the parent
node. Finally, the study in [27] introduced a novel parent switching scheme
predicated on the remaining queue capacity of the parent node, which
effectively improved the PDR performance.

##  III Proposed Solution

Our work extensively focus on the intricate interplay between RPL and SF, and
between SF and packet queue management. The core philosophy of our solution is
as follows: 1) Reduce the transmission frequency of control packets upon
parent switching; 2) Minimize the number of failed parent switching events
that result from cell unavailability; 3) Prevent packet dropping events due to
queue saturation; and 4) Improve the jitter and latency performances by slot
selection. Our approach targets saving valuable energy resources and
mitigating unnecessary medium contention, thus enhancing the spatial reuse.

For understanding purposes of our solution, we highlight the following
technical terms usage:

  * •

Reference to “DAO” is explicitly linked to the DAO packet triggered upon
parent switching, unless otherwise specified. This specificity arises from our
focus on mechanisms and enhancements relevant to the parent-switching events,
while the processes involving periodic DAO packets remain unaltered through
our approach.

  * •

Reference to “DIO” may refer to either periodic or direct response to a DIS
packet, dependent on the context of the mentioned processes or mechanisms.

  * •

“DAO-ACK” and “ACK” are two distinct message. DAO-ACK denotes a packet
transmitted by the DAG root during a to acknowledge DAO messages for reliable
route registration/de-registration. In contrast, a packet “ACK”, known also as
a link-layer or TSCH acknowledgment (shown in Fig. 6), is issued upon packet
reception [28].

Figure 6: Packet TSCH Acknowledgment.

###  III-A Regulated Triggering of Parent Switching

In 6TiSCH networks, parent switching might be challenging. Indeed, the
decision to switch to a potentially superior parent is complex, given the
node’s lack of awareness regarding the new parent’s cell availability. This
uncertainty can lead to inefficient energy usage since the node may continue
control packet transmissions without realizing effective cell allocation.
Also, the RPL specifications [29, 8] did not detail the protocol’s behavior
regarding transmission interruption to the former parent during this process.
Consequently, increased latency and potential packet loss may occur. To tackle
this issue, our solution proposes to modify the DIO packet in order to carry a
list of available slots of the sender. This approach ensures that a node
initiates parent switching only when a shared available slot exists, hence
improving the efficiency of the process. Algo. 1 below delineates this
mechanism, particularly focusing on the node’s decision-making process at the
juncture of parent switching.

Data: $parents$ : DIO’s senders

$min\\_nb$: Minimum number of slots to allocate

$max\\_nb$: Maximum number of slots to allocate

Result: Parent switching decision

$nb$ : Number of slots to reserve

// \------ Initialization

1 $nb\\_parentSlot\leftarrow\text{Number of slots allocated to parent}$

2 $child\\_slots\leftarrow\text{Node's free slots}$

// \------ Parent Selection and Slot Reservation

3 $selected\\_parent\leftarrow parent\mid parent.rank=\min(p.rank\textbf{ For
}p\textbf{ in }parents\textbf{ and }p.free\\_slot\cap
child\\_slots\neq\emptyset)$

4 if  _$selected\\_parent\neq None$_ then

5 if  _Old parent is None_ then

// Node is joining the DODAG

6 $nb\leftarrow$ min_nb

7

8 else

// Node is switching parent

9 $nb\leftarrow$ nb_parentSlot

10 if  _$nb >max\\_nb$_ then

// Ensure selected slots are not exceeding the max limit

11 $nb\leftarrow max\\_nb$

12

13

14 $free\\_slots\leftarrow child\\_slots\cap parent\\_free\\_slots$

15 $slots\leftarrow\textbf{Select}\ \text{nearest}\ nb\ \text{slots}\
\textbf{from}\ free\\_slots$

16 $child\\_slots\leftarrow child\\_slots-slots$

17 Append $slots$ to DAO packet

18 Append $child\\_slots$ to DAO packet

19 Transmit DAO packet

20

21else

22 Pass

// Retain current situation since there is no slot availability

23

Algorithm 1 Parent switching procedure

###  III-B Minimizing 6P Transaction Packets

In the sequenced control packet transmission process of 6TiSCH networks,
depicted in Fig. 3, when a node decides to change its parent, it communicates
this transition by transmitting a DAO packet to the DAG root and triggers the
cell reservation process with the new parent. This entails dispatching a cell
reservation request using a 6P transaction as shown in Fig 5.

Rather than dedicating an entire packet for this purpose, we propose to
judiciously embed the requisite information for cell reservation within the
RPL control packet DAO. Since the next hop for each child packet is inherently
its parent, all packets, including the modified DAO, will go through the
designated parent and hence can extract the embedded information about cell
reservation and define its relevant cells. This approach is detailed as
follows:

  1. 1.

Autonomous channel selection: Transmission channel is autonomously determined
based on the MAC address of the receiver, thereby enabling nodes to distribute
only the slot offset to represent the cell. To avoid interference and benefit
from multi-channeling, we adopt the same approach as MSF. Specifically, the
current transmission/reception channel ($channel\\_{\rm ID}$) is determined
based on the Absolute Slot Number (ASN) and the channel offset ($cell_{\rm
channel}$) of the current cell as follows:

$channel\\_{\rm ID}=\left(ASN+cell_{\rm channel}\right)\ \%\ nb\\_{\rm
channels},$ (1)

where $nb_{\rm channels}$ is the total number of channels and $\%$ is the
remainder of the Euclidean division.

  2. 2.

RPL packet size management: In 6TiSCH networks where the packet size is
intrinsically constrained, judicious management of the RPL packets DIO, DAO,
and DAO ACK, is important. Consequently, we propose the following
modifications:

    1. (a)

Selective data injection: Conventionally, cells are represented as tuples of
slot-offset and channel-offset. Their injection into the RPL packets leads to
increased packet size and processing time. Aiming to reduce packets’ size, we
propose a simplified integer-only representation of the available cells
instead of using space-greedy tuples.

    2. (b)

Efficient Slot ID encoding: To further enhance the packet size, we introduce
an encoding scheme for slot IDs. Typically, 6TiSCH systems limit the slotframe
length to 101 slots. However, we propose here to extend the slotframe length
to 255 slots, such that the slot ID never exceeds this value. By doing so, we
can encode each slot ID into a single byte (8 bits), thus minimizing the slot
ID data size injected into the packet.

    3. (c)

Compact free slots notation: Within our proposed approach, the DAO and DIO
packets have to convey the node’s free slots. To optimize memory usage, we
propose that the node transmits either an ordered list of occupied slots or of
free slots, depending on which is shorter and anticipating it with a flag (1
bit). Within the 6TiSCH context, slot “0” is reserved by nodes upon joining
the network, hence it will never be included in the ordered list of free
slots. This reserved status allows to efficiently utilize slot “0” as the flag
for the occupied/free slots list. This flexible method guarantees less memory
usage.

The aforementioned techniques align with the fundamental requirements of low-
power and lossy networks (LLNs) since they aim to reduce transmission
overhead, save energy, and enhance the network’s throughput. Moreover, they
simplify demands’ processing, which is critical for devices operating within
constrained environments, such as in 6TiSCH networks.

  3. 3.

TSCH ACK-based optimization of slot allocation: For efficient time slot
allocation, we propose a novel method that uses TSCH ACK for immediate
confirmation of slot reservations upon DAO packet reception. Doing so bypasses
the conventionally incurred delays for reservation confirmation. Leveraging
TSCH ACK ensures quicker communication of available time slots, thus reducing
the locking duration of time slots (shown as $T_{\rm Lock}$ in Fig. 5) and
increasing the network’s flexibility for topological changes. Our technique
embeds a list of available slots (identified with integers) directly into the
TSCH ACK messages. For the sake of data size reduction, we limit the maximum
number of time slots allocated at once to 5 resulting in an additional payload
of merely 5 bytes in TSCH ACK. Moreover, to ensure efficient time slot
reservation, we assume that each node maintains an up-to-date variable
containing its ordered available slots list. This variable is dynamically
updated following each slot allocation or de-allocation. Upon receiving a
packet, the joining node compares the list of available slots embedded within
the packet against its variable of available slots. Hence, the decision to
reserve slots is based on the most current network status, which improves the
responsiveness and adaptability of the 6TiSCH network to real-time demands.

###  III-C Reducing DODAG Joining Time

To reduce packet collisions and reduce the nodes’ joining time, we propose
that the DIO sender sets temporary time slots or cells, ensuring it is
actively listening, and it designates one permanent cell to maintain its
ongoing listening state. In scenarios where a DIO is concurrently received by
multiple nodes, thus causing potential DAO collisions, the proposed mechanism
that uses a combination of temporary and permanent cells can diminish the
impact of this issue as follows: 1) The DAO sender temporarily allocates cells
until it receives the TSCH ACK of the DAO packet. This strategy intentionally
sacrifices idle listening to prioritize energy saving, given that the energy
consumed during idle states is notably lower than that needed to manage
collisions and trigger packet retransmissions [30]. 2) In the initial phase of
network formation where nodes actively compete to join the network and parent
switching events are frequent, our method sets a high number of proposed time
slots in the DIO to facilitate rapid integration and reduced collision events.
As the network stabilizes, we strategically reduce the number of proposed
slots to preserve energy.

###  III-D Scheduling Information Injection

In the context of TSCH, transmitting packets from higher-layer protocols is
not an instantaneous process; rather, it entails a meticulous queuing and
scheduling mechanism. As depicted in Fig. 7, incoming packets are first
enqueued within the TSCH queue and await their transmissions. Since we inject
the proposed time slots in DIO and DAO packets, they will be locked during
time $T_{q}=T_{\rm de}-T_{\rm en}$, where $T_{\rm en}$ is the instant at which
the packet entered the TSCH queue, and $T_{\rm de}$ is the time at which it
left the queue. To avoid locking the proposed time slots for a long time
$T_{q}$, we propose to inject them within the packet only before leaving the
queue, i.e., at instant $T_{de}$, hence $T_{q}\approx 0$. By carefully
orchestrating the data injection in the RPL control packets, we optimize time
slot utilization.

Figure 7: Illustration of a packet time in a 6TiSCH queue.

###  III-E Avoiding TSCH Queue Overflow

In 6TiSCH networks, effective synchronization between SF and TSCH queue
management is essential for network operation. With various protocols filling
the queue with packets, timely processing during the SF-designated
transmission cells is crucial. A mismatch between the queue and cell
allocation can lead to packet loss. To address this issue, we propose the
introduction of an early cell reservation technique. Specifically, as the TSCH
queue begins to be filled and before reaching its limit, our system
proactively initiates 6P transactions with the preferred parent for an
additional cell reservation. This proactive approach offers several benefits:
1) Prevent packet drops: Early reservation of additional cells reduces the
likelihood of packet loss due to queue overflow, and 2) Mitigate resource
urgency: This method allows nodes to secure necessary resources in advance,
thus avoiding the last-minute pressure of queue overflow and network
instability.

###  III-F Adaptive Self-Sufficient Cell Reservation

In a typical 6TiSCH network, all DAO packets are destined towards the DAG
root, thus navigating in a multi-hop fashion. In our proposal, essential
information about cell reservation, which consists of the available slots and
proposed slot lists, is embedded within these DAO packets. When a parent node
receives a DAO packet, as it is the principal conduit to the DAG root, it
reserves new cells for its new child where the number of cells should not
exceed the maximum number of cells to allocate at once. Then, it loads the
list of its free slots (instead of the old one that belongs to the child) and
the number of slots to be selected (slots to be allocated to its new child)
into the DAO packet. This packet is then seamlessly forwarded, with each
intermediate node towards the DAG root executing the same procedure. As shown
in Fig. 8, upon reception of a DAO packet from node $3$, node $i$ undertakes
cell selection for reservation and embeds the list of chosen slots into the
TSCH ACK packet. It then forwards the enhanced packet to its parent. To reduce
latency at every hop, we propose that the receiver of a DAO packet selects the
nearest time slot to the actual one used for the reception. One way to do it
is to rely on the concept of the Low-Latency SF (LLSF) protocol [31] where the
slots are selected according to the gap to the left of the actual slot.
However, we propose here to select the nearest slot to the actual one in any
direction. Accordingly, if the node is joining the DODAG, the actual slot will
be in the list of proposed slots, thus maximizing its chance of selection. If
the node has already joined the DODAG, then the actual slot is the slot that
is already allocated to its parent. Consequently, both energy and latency tied
to time slot reservations can be reduced.

Figure 8: Illustration of the DAO packet transmission and TSCH ACK utilization. TABLE I: Simulation parameters Parameters | Value
---|---
Number of nodes | 50, 100
Number of simulation runs | 100
Simulation duration (min) | 30, 270
Network area (km2) | 1
Service or application packet period (sec) | 5, 15, 45, 60
Slotframe length (in number of slots) | 100
Time slot duration (sec) | 0.01
Secure Joining required | True
TSCH maximum payload length (bytes) | 120
TSCH queue size (number of packets) | 10
DIO packet size (bytes) |  |  MSF |  PB
---|---
76 |  120
DAO packet size (bytes) |  |  20 |  75
---|---

##  IV Performance evaluation

###  IV-A Simulation configuration

To comprehensively assess the efficiency and resilience of our proposed
scheme, we employed the 6TiSCH simulator [32]. For the sake of thorough
analysis, generated plots include outlier results, thus reflecting the full
range of observed results. Our simulations spanned $30$ minutes and $270$
minutes ($4.5$ hours) and were repeated 100 times across two distinct network
configurations:

  * •

Conventional 6TiSCH configuration: It consists of the use of the conventional
MSF as the scheduling function and OF0 as the RPL objective function [11, 8].
Within the results, this configuration is denoted as “MSF”.

  * •

Proposed solution-based configuration: This setup leveraged our proposed
solution with the modified MSF and OF0 as the scheduling and objective
functions, respectively, and labeled as “PB” within the results.

For our simulations, we assume a $1\times 1$ km2 area, where $N=50$ or $N=100$
nodes are randomly deployed [32]. Moreover, we assume that the frequency of
data transmissions for the application or service is in the values
$\\{5,15,45,60\\}$ seconds. Also, the slotframe length has been fixed for all
nodes in the network, either equal to $N$ or to $2N$. The remaining simulation
parameters are presented in Table I, while Table II exposes the specific
parameters of the proposed solution.

The selected value for the “TSCH queue empty places threshold for cell
reservation” is justified as follows. Upon each packet enqueuing operation,
nodes assess the remaining capacity within the TSCH queue. When the available
space in the queue reaches a predefined threshold, the node initiates a cell
reservation request using the 6P protocol, where “Number of cells reserved
upon reservation process” (specified in Table II) cells can be requested at
that time, and subsequent requests will be constrained by the value of the
parameter “Cycle interval before next cell request”. This threshold is set to
$2$ empty places in the queue to accommodate the 6P transaction packet and to
reserve capacity for any imminent packet. Such a configuration prevents queue
overflow and enhances the network’s capability to handle data packets,
especially in dense-traffic scenarios.

TABLE II: Parameters of the proposed solution Parameter | Value
---|---
Number of proposed slots per DIO | 7
DIO cells duration (in number of slotframes) | 10
Number of permanent DIO slots | 1
Slot selection ratio within DIO slots | 3
Maximum number of cells to reserve at once | 5
Initial phase duration (min) | 45
TSCH queue empty places threshold for cell reservation | 2
Cycle interval before next cell request | 2
Number of cells reserved upon reservation process | 1
Figure 9: Joining time vs. nbr. of proposed slots per DIO. Figure 10: Consumed
current vs. nbr. of proposed slots per DIO.

To set the parameter “Number of proposed slots per DIO” in Table II to the
best value, we conducted several experiments with different values.
Specifically, in Figs. 9-10, we evaluated the impact of different numbers of
proposed time slots per DIO on the joining time and current consumption
(equivalent to energy consumption).

These experiments have been conducted over an observation time of $30$
minutes, with an application period of $5$ seconds, and for different network
sizes, i.e., $N_{1}=50$ and $N_{2}=100$ nodes. As shown in Fig. 9, for any
number of nodes, the joining time reduces when the number of proposed slots is
higher. This is expected since more proposed slots mean that a novel node can
rapidly be configured and join the DODAG. However, since the number of
proposed slots cannot be very high, as it may increase the size of packets,
our preference is for the value $7$ since it provides similar results to the
system with the value $13$, in particular for $N_{2}$, i.e., for a large
network. A similar explanation can be provided for Fig. 10 where the value $7$
achieves low current consumption for the large network. Indeed, a very high
value, such as $13$, can result in inefficient idle listening, thus increasing
power consumption. Consequently, we recommend the number of proposed slots per
DIO to be $7$ as the best value since it is the lowest one that achieves good
performances in terms of joining time and consumed power.

(a) $N=50$

(b) $N=100$

Figure 11: Joining time vs. length of service packet period (in sec) (first 30
minutes).

###  IV-B Simulation Results

The initial phase of the network’s life is known as network formation. This
phase is pivotal for establishing the network’s parameters and topology. To
demonstrate the effectiveness of our solution in accelerating the network
formation timeline and minimizing power consumption [30], we conducted
simulations in Figs. 11-12 under the same conditions as the previous figures,
while we varied the value of the parameter “Service packet period”. Also, we
compared the results of our proposal, denoted PB, with those of the benchmark
MSF.

(a) $N=50$

(b) $N=100$

Figure 12: Consumed current vs. length of service packet period (in sec)
(first 30 minutes).

(a) $N=50$

(b) $N=100$

Figure 13: Total nbr. of sent/received packets (different service packet
periods, 4.5 hours).

As shown in Fig 11, our proposed solution PB achieves lower joining time for
any $N$ and service packet period duration. For instance, PB achieves an
average joining time of 13 minutes for a service packet period equal to $15$
and $N=100$, while MSF needs approximately 14 minutes in the same conditions.
Indeed, a node within the PB framework uses proposed slots in the received DIO
to define the cells during which it sends the DAO, unlike MSF which relies on
sending a 6P request to its preferred parent and then waiting for a response
for cell reservation. In PB, the DAO packet serves a dual purpose:
Communicating with the RPL to establish the DODAG structure and participating
in SF. As the service packet period increases, we notice a degradation in the
joining time. This is expected since a node with a low (resp. high) service
packet period will transmit EB packets more (resp. less) often, which
increases (resp. reduces) the likelihood that other nodes synchronize and join
the DODAG.

In 6TiSCH networks, the process of network formation often expends a
significant amount of energy, even before useful data transmission. In Fig.
12, we present the consumed current related to the scenarios of Fig.11. The PB
scheme performs better than the MSF one, for any $N$ and service packet period
duration. The performance gap increases as $N$ becomes higher. This result
proves the robustness of our scheme and its superiority when operating in
dense networks. Finally, we notice that the service packet period has a slight
impact on the total consumed current.

After validating the improvements of our solution on network formation, we
assess the performances of the schemes for long-term operation. We evaluate
the performances in terms of the number of packets transmitted/received, the
end-to-end (E2E) latency, network jitter, and consumed charge. For the next
results, experiments have been conducted for 4.5 hours.

Fig. 13 illustrates the number of transmitted (Tx) and received (Rx) packets
for different service packet periods and $N$, where PB and MSF are compared.
For any $N$ and service packet period, our solution outperforms the MSF
benchmark. In other words, the network that uses our approach realizes more
frequent transmissions, and nodes within it successfully receive more packets.
When the service packet period increases, i.e., a longer service is running,
the number of exchanged packets drops. such a result is expected since during
the observation time of 4.5 hours, a longer service time impacts the frequency
of its execution over the observation period. Moreover, when $N=100$, the
number of exchanged packets almost doubles compared to the case with $N=50$.
Indeed, the availability of more nodes joining the DODAG makes it possible to
transmit/receive more packets, within the limits of available channels and
time slots in 6TiSCH.

(a) $N=50$

(b) $N=100$

Figure 14: E2E latency vs. service packet period (in sec) (4.5 hours).

(a) $N=50$

(b) $N=100$

Figure 15: Median and avg. jitter for transmitted packets during 4.5 hours
(different service packet periods).

Fig. 14 shows the network’s E2E latency for different service packet periods
and $N$. End-to-end latency is defined as the duration measured from the ASN
at which a source node transmits a service/application packet, to the ASN at
which this packet is successfully received and acknowledged by the
destination. As it can be seen, the PB scheme achieves lower E2E latency in
any network condition than the MSF benchmark. This result underscores the
efficiency of our method. Even when the service packet period increases, the
E2E latency follows in a logarithmic trend with values stabilizing at high
values. For instance, for $N=100$, the average E2E latency is almost the same,
about 2.1 sec, for service packet periods $45$ and $60$.

In Fig. 15, we plot the median and average jitter performance for the
transmitted packets during the observation time, for different $N$ and service
packet periods. Jitter refers to the variation in packet transmission latency
over time within a network. It is the inconsistency observed in the arrival
time of packets, primarily caused by varying network congestion levels, route
changes, or other unforeseen factors. To obtain the median and average
results, we generated $100$ scenarios for each packet order point. The median
offers insight into the typical experience of a node without the influence of
outliers, while the average jitter aggregates these variations to provide a
general sense of the network’s jitter over time. According to Fig. 15, MSF
shows an initial peak (median and average jitters), due to the network’s
formation phase, which might disrupt real-time services. In contrast, PB
consistently maintains lower jitter performances than MSF, which indicates a
stable packet delivery time. The median jitter of MSF exhibits high
variability, thus unsuitable for services requiring consistent data rates,
while PB follows a more predictable jitter performance, which makes it fit for
real-time applications. As $N$ increases from $50$ to $100$, and for a given
scheme and service packet period, we notice that the jitter performance is
stable and does not degrade. For instance, for a service packet period of $5$
seconds and PB scheme, the average jitter is at 0.6 seconds for $N=50$ and
$N=100$. However, we notice that for $N=50$, variations in median jitter is
slightly higher than for $N=100$. This might be due to the fixed slotframe
length that might not be suitable for any network size. For a given $N$, the
jitter degrades with the service packet period. Indeed, a longer service time
impacts packet transmissions and thus might cause additional delay and/or
instability. In summary, the jitter is slightly impacted by the size of the
network, while the influence of the service packet delay is more significant
on both the average and median jitters.

(a) $N=50$

(b) $N=100$

Figure 16: CDF of consumed charge.

(a) $N=50$

(b) $N=100$

Figure 17: CDF of nodes’ lifetime.

Fig. 16 presents the Cumulative Distribution Function (CDF) as a function of
the charge consumed (in Coulomb) for the PB (blue) and MSF (purple) solutions,
for different $N$ and service packet periods. For any scenario, the charge
consumption profiles of PB and MSF are almost identical. This result
emphasizes the efficiency of our solution in conserving energy similarly to
MSF, while realizing significant improvements in terms of other performance
metrics, such as joining time, E2E latency, and jitter. We note that for a
given $N$, when the service packet period increases, the CDF’s slope is almost
the same. This means that the amount of charge consumed by a node in the
network is not affected by the service packet period. Also, the CDF between
$N=50$ ad $N=100$ at a fixed service packet period are very similar. Hence,
$N$ has a negligible impact on the consumed charge per node.

In Fig. 17, we evaluate the CDF of the lifetime of the network’s nodes for the
PB (blue) and MSF (purple) solutions, for different $N$ and service packet
periods. For any $N$ and service packet period value, the lifetime of nodes in
the PB-based network is slightly lower than that of the MSF-based one. For
instance, for $N=50$ and service packet period equal to $45$, we see that
$67\%$ of the nodes will have a lifetime of less than 5 years using PB, while
this proportion is only $58\%$ for MSF. When $N$ increases, the gap between PB
and MSF gets lower, which highlights the capability of our method to save
energy in dense networks.

##  V Conclusion

Our research introduces a novel approach to enhance 6TiSCH network
performances by carefully optimizing the interaction between SF and RPL.
Through a multitude of modifications within the exchanged messages in 6TiSCH
and improvements in SF and RPL operations, our proposed solution is capable of
enhancing several performance metrics, including joining time,
transmitted/received packet number, E2E latency, and jitter, while keeping the
energy usage at levels similar to the MSF benchmark and nodes’ lifetime.
Despite the realized improvements, several challenges remain to be solved. For
instance, the fixed duration of time slots prevents the network from achieving
its full potential.

## Acknowledgment

This work was supported by funding from the Algerian Ministry of Higher
Education and Scientific Research.

## References

  * [1] A. Baiyere, H. Topi, V. Venkatesh, J. Wyatt, and B. Donnellan, “The internet of things (iot) - a research agenda for information systems,” _Communications of the Association for Information Systems_ , vol. 47, 05 2020\.
  * [2] D. De Guglielmo, S. Brienza, and G. Anastasi, “IEEE 802.15.4e: A survey,” _Computer Communications_ , vol. 88, pp. 1–24, 2016.
  * [3] D. Guglielmo, G. Anastasi, and A. Seghetti, “From ieee 802.15.4 to ieee 802.15.4e: A step towards the internet of things,” _Advances in Intelligent Systems and Computing_ , vol. 260, pp. 135–152, 03 2014.
  * [4] J. Adriano, E. C. Rosário, and J. Rodrigues, “Wireless sensor networks in industry 4.0: Wirelesshart and isa100.11a,” in _Proc. IEEE Int. Conf. Indus. Appl. (INDUSCON)_ , 11 2018, pp. 924–929.
  * [5] N. Choudhury, R. Matam, M. Mukherjee, and J. Lloret, “A performance-to-cost analysis of ieee 802.15.4 mac with 802.15.4e mac modes,” _IEEE Access_ , vol. 8, pp. 41 936–41 950, 2020.
  * [6] “Ieee standard for low-rate wireless networks,” _IEEE Std 802.15.4-2015 (Revision of IEEE Std 802.15.4-2011)_ , pp. 1–709, 2016.
  * [7] T. Chang, M. Vučinić, X. Vilajosana, S. Duquennoy, and D. R. Dujovne, “6TiSCH minimal scheduling function (MSF),” RFC 9033, May 2021. [Online]. Available: https://www.rfc-editor.org/info/rfc9033
  * [8] P. Thubert, “Objective Function Zero for the Routing Protocol for Low-Power and Lossy Networks (RPL),” RFC 6552, Mar. 2012. [Online]. Available: https://www.rfc-editor.org/info/rfc6552
  * [9] A. Tabouche, B. Djamaa, and M. R. Senouci, “Traffic-Aware Reliable Scheduling in TSCH Networks for Industry 4.0: A Systematic Mapping Review,” _IEEE Communications Surveys & Tutorials_, pp. 1–1, 2023. [Online]. Available: https://ieeexplore.ieee.org/document/10208136/
  * [10] X. Vilajosana, K. Pister, and T. Watteyne, “Minimal IPv6 over the TSCH Mode of IEEE 802.15.4e (6TiSCH) Configuration,” RFC 8180, May 2017. [Online]. Available: https://www.rfc-editor.org/info/rfc8180
  * [11] T. Chang, M. Vučinić, X. Vilajosana, S. Duquennoy, and D. R. Dujovne, “6TiSCH Minimal Scheduling Function (MSF),” Internet Engineering Task Force, Internet-Draft draft-ietf-6tisch-msf-10, 2019, work in Progress. [Online]. Available: https://datatracker.ietf.org/doc/html/draft-ietf-6tisch-msf-10
  * [12] P. Levis, T. H. Clausen, O. Gnawali, J. Hui, and J. Ko, “The Trickle Algorithm,” RFC 6206, Mar. 2011. [Online]. Available: https://www.rfc-editor.org/info/rfc6206
  * [13] Q. Wang, X. Vilajosana, and T. Watteyne, “6TiSCH Operation Sublayer (6top) Protocol (6P),” RFC 8480, Nov. 2018. [Online]. Available: https://www.rfc-editor.org/info/rfc8480
  * [14] M. Rady, Q. Lampin, D. Barthel, and T. Watteyne, “6DYN: 6TiSCH with heterogeneous slot durations,” _Sensors_ , vol. 21, no. 5, p. 1611, Feb. 2021\.
  * [15] G. Daneels, D. Van Leemput, C. Delgado, E. De Poorter, S. Latré, and J. Famaey, “Parent and PHY selection in slot bonding IEEE 802.15.4e TSCH networks,” _Sensors_ , vol. 21, p. 5150, Jul. 2021.
  * [16] D. Fanucchi, B. Staehle, and R. Knorr, “Network formation for industrial IoT: Evaluation, limits and recommendations,” in _Proc. IEEE Int. Conf. Emerg. Technol. Factory Automation (ETFA)_ , vol. 1, 2018, pp. 227–234.
  * [17] A. Kalita and M. Khatua, “Opportunistic transmission of control packets for faster formation of 6TiSCH network,” _ACM Trans. IoT_ , vol. 2, no. 1, pp. 1–29, Feb. 2021.
  * [18] Y. Tanaka, P. Minet, M. Vučinić, X. Vilajosana, and T. Watteyne, “Ysf: A 6tisch scheduling function minimizing latency of data gathering in iiot,” _IEEE Internet of Things Journal_ , vol. 9, no. 11, pp. 8607–8615, 2022.
  * [19] G. Daneels, B. Spinnewyn, S. Latré, and J. Famaey, “ReSF: recurrent low-latency scheduling in IEEE 802.15.4e TSCH networks,” _Ad Hoc Netw._ , vol. 69, pp. 100–114, 2018.
  * [20] W. Yang, Y. Cao, H. Wang, Z. Zhang, and C. Wu, “EDSF: efficient distributed scheduling function for IETF 6TiSCH-based industrial wireless networks,” _Mob. Netw. & Appl._, Jun. 2022.
  * [21] S. Jeong, J. Paek, H.-S. Kim, and S. Bahk, “TESLA: Traffic-aware elastic slotframe adjustment in TSCH networks,” _IEEE Access_ , vol. 7, pp. 130 468–130 483, Sept. 2019.
  * [22] X. Fafoutis, A. Elsts, G. C. Oikonomou, R. J. Piechocki, and I. Craddock, “Adaptive static scheduling in IEEE 802.15.4 TSCH networks,” in _Proc. IEEE World Forum IoT (WF-IoT)_ , 2018, pp. 263–268.
  * [23] A. R. Urke, O. Kure, and K. Ovsthus, “A survey of 802.15.4 TSCH schedulers for a standardized industrial internet of things,” _Sensors_ , vol. 22, no. 1, 2022.
  * [24] S. Kharb and A. Singhrova, “A survey on network formation and scheduling algorithms for time slotted channel hopping in industrial networks,” _J. Netw. & Computer Appl._, vol. 126, pp. 59–87, 2019.
  * [25] B. Aydın, S. Görmüş, H. Aydın, and S. Külcü, “A new routing objective function for IETF 6TiSCH protocol,” in _Proc. Sig. Process. Commun. Appl. Conf. (SIU)_ , 2022, pp. 1–4.
  * [26] T. Bekar, S. Görmüş, B. Aydın, and H. Aydın, “Q-learning algorithm inspired objective function optimization for IETF 6TiSCH networks,” in _Proc. Int. Conf. Smart Appl. Commun. Network. (SmartNets)_ , 2023, pp. 1–6.
  * [27] A. Kalita, A. Hazra, and M. Gurusamy, “Efficient schemes for improved performance in 6TiSCH networks,” in _Proc. IEEE Conf. Computer Commun. Wrkshps. (INFOCOM WKSHPS)_ , 2023, pp. 1–6.
  * [28] T. Watteyne, M. R. Palattella, and L. A. Grieco, “Using IEEE 802.15.4e Time-Slotted Channel Hopping (TSCH) in the Internet of Things (IoT): Problem Statement,” RFC 7554, May 2015. [Online]. Available: https://www.rfc-editor.org/info/rfc7554
  * [29] R. Alexander, A. Brandt, J. Vasseur, J. Hui, K. Pister, P. Thubert, P. Levis, R. Struik, R. Kelsey, and T. Winter, “RPL: IPv6 Routing Protocol for Low-Power and Lossy Networks,” RFC 6550, Mar. 2012. [Online]. Available: https://www.rfc-editor.org/info/rfc6550
  * [30] X. Vilajosana, Q. Wang, F. Chraim, T. Watteyne, T. Chang, and K. Pister, “A realistic energy consumption model for tsch networks,” _Sensors Journal, IEEE_ , vol. 14, pp. 482–489, 02 2014.
  * [31] T. Chang, T. Watteyne, Q. Wang, and X. Vilajosana, “LLSF: Low latency scheduling function for 6TiSCH networks,” in _Proc. Int. Conf. Dist. Comput. Sensor Syst. (DCOSS)_ , 2016, pp. 93–95.
  * [32] E. Municio, G. Daneels, M. Vučinić, S. Latré, J. Famaey, Y. Tanaka, K. Brun, K. Muraoka, X. Vilajosana, and T. Watteyne, “Simulating 6TiSCH networks,” _Trans. Emerg. Telecommun. Technol._ , vol. 30, no. 3, p. e3494, 2019.


# SHIELD: Secure Haplotype Imputation Employing Local Differential Privacy

Marc Harary<EMAIL_ADDRESS>Yale University

###### Abstract

We introduce Secure Haplotype Imputation Employing Local Differential privacy
(SHIELD), a program for accurately estimating the genotype of target samples
at markers that are not directly assayed by array-based genotyping platforms
while preserving the privacy of donors to public reference panels. At the core
of SHIELD is the Li-Stephens model of genetic recombination, according to
which genomic information is comprised of mosaics of ancestral haplotype
fragments that coalesce via a Markov random field. We use the standard
forward-backward algorithm for inferring the ancestral haplotypes of target
genomes—and hence the most likely genotype at unobserved sites—using a
reference panel of template haplotypes whose privacy is guaranteed by the
randomized response technique from differential privacy.

##  1 Introduction

In the context of biomedical analyses of large patient cohorts, whole-genome
sequencing still remains prohibitively expensive for existing high-throughput
technology. On the other hand, array-based genotyping platforms provide a more
efficient method of collecting data for large-scale studies of human disease,
albeit at the expense of the statistical power of genome-wide association
(GWA) studies that intend to fine-map causal variants or facilitate meta-
analyses [1, 2, 3, 4, 5].

One solution is genotype imputation, a preliminary stage in many GWA studies
that consists of inferring the genotype for a given target genome at loci that
have not been directly assayed, essentially expanding the dimensionality of
the original dataset [6, 5, 7, 8, 9, 10, 2]. Employing a reference panel of
donated haplotypes sequenced via higher-quality technology and at a far denser
set of variants, imputation algorithms like MaCH [7], Minimac [8], BEAGLE [9],
PLINK [10], fastPHASE [11], and IMPUTE [2] have been demonstrated to reliably
augment both the coverage and statistical power of GWA analyses and hence
become an essential component of many clinical studies [6].

Further to this end, public databases like the UK biobank (UKB) [12], All of
Us research program [13], Haplotype Reference Consortium [14], and 1,000
Genomes Project (1KG) [15] have been made available to facilitate genomic
research in part by offering standardized and readily accessible reference
panels [16]. In cases where running imputation algorithms using large
reference panels is impractical on local hardware or the direct access to the
biobank data is prohibited, public web services like the Michigan Impute
Server [17] are often established to answer queries to clients submitting
target haplotypes for imputation.

Unfortunately, as part of a growing literature on privacy concerns in genomic
research, it has also been documented that coordinated attacks on the part of
cryptographic adversaries are capable of compromising the privacy of research
subjects that donate to public reference panels [18, 19, 20, 21]. For example,
attackers have been able to exploit ancestral data [22] or other personally
identifying information [23] to reconstruct reference genomes. An urgent
challenge is therefore to develop a suite of imputation algorithms that can
simultaneously facilitate high-utility, statistically reliable GWA studies
while protecting the privacy of contributors to reference haplotype panels
[24, 18, 25].

One solution is the technique of differential privacy, which has rapidly
become the “gold-standard” for statistical queries by being able to provide
both robust privacy guarantees for participants in studies and meaningful
results for researchers in commercial and scientific settings [26, 27, 28]. At
the crux of the technique is a rigorous mathematical formalization of privacy
that quantifies the extent to which adding pseudorandom noise to the results
of computations can protect the anonymity of members of a database [29].

The following work introduces Secure Haplotype Imputation Employing Local
Differential privacy (SHIELD), a program that employs the Li-Stephens model of
genetic recombination [30, 5] to impute missing haplotype variants in target
genomes while incorporating differential privacy techniques to protect
reference panel donors. Specifically, SHIELD proceeds in two stages: (i)
initial input perturbation to guarantee local differential privacy [31] via
randomized response [32, 33] and (ii) fitting a hidden Markov model [34] to
each subsequent client query via the forward-backward algorithm [35]. In an
experiment that closely simulates a real-world use case for haplotype
imputation, we show that SHIELD is able to obtain state-of-the-art imputation
accuracy while providing mathematically formalized privacy guarantees.

##  2 Results

###  2.1 Overview

The setting for which SHIELD is intended consists of a client user uploading
target genomes to a public imputation server [18]. In the standard imputation
workflow, contributors to a biobank upload their sequenced genomic data to a
central, publicly available server, where the data are then collated to create
a haplotype reference panel to pass as an argument to an imputation algorithm
[12, 15, 14]. Subsequently, client researchers may then upload target genomes
as part of a clinical study to the server, where the targets are imputed using
the private haplotype reference panel and, most often, an algorithm based in
hidden Markov models [34, 7, 8, 9, 10, 2] and the forward-backward algorithm
[35]. At no point in the workflow is the haplotype reference panel directly
visible to client researchers submitting jobs to the server. However, while
the privacy of the contributors to the reference panel may appear guaranteed,
it has been demonstrated that adversarial attacks employing carefully
coordinated queries to the server can divulge the sequences of reference
haplotypes [18].

To this end, SHIELD modifies the imputation workflow by leveraging local [31]
differential privacy [26, 27, 28, 36, 32]. Haplotype data can be represented
as a bitstring in which a 1 at the $i$th position in the sequence indicates
that the haplotype possesses the minor allele at the $i$th site and a 0 the
major allele [8]. Prior to submission to the central imputation server,
pseudorandom noise is added to the two bitstrings denoting each individual’s
pair of haplotypes via randomized response, a technique from differential
privacy that simply consists of flipping a random subset of the bits from 0 to
1 and vice versa [32, 33]. The likelihood that a given bit in the haplotype
bitstring is flipped varies as a function of a parameter $\varepsilon$—called
the privacy budget [29]—such that lower values of $\varepsilon$ entail a
higher probability that any bit is flipped and therefore a higher degree of
privacy. The tradeoff, however, is that lower privacy budgets incur a greater
expense to imputation accuracy, rendering it a hyperparameter that the
database curator must carefully adjust to strike an acceptable balance between
donor privacy and client utility. Once all perturbed haplotypes are collected
at the central server, imputation is subsequently performed using the modified
haplotypes as a reference panel.

Privacy is guaranteed by the fact that no contributor’s data will, on average,
be unmodified when input to the imputation algorithm invoked by client
researchers. In this way, no adversary could be certain that the results that
they obtain from an attack accurately reflect the true reference panel. These
privacy guarantees are also local; even if an adversary were to access the
reference panel directly rather than through coordinated queries, the data
obtained would again not perfectly reflect any individual’s true genome [36].

Figure 1: Overview of the SHIELD pipeline, with the key algorithms in orange.
Noise is added once to the reference data (purple) via Perturb, then collated
and stored on the server to guarantee local DP (modified bits in bold). The
client (green) then calls Impute on the server with the target haplotype
(missing sites denoted $\varnothing$) and the reference panel as arguments.

###  2.2 State-of-the-art imputation accuracy

Figure 2: A. The reference haplotype matrix corresponding to the first 128
SNPs on chromosome 20 and 200 haplotypes in 1KG. Empty squares represent the
presence of the major allele, yellow of the minor. B. The same haplotype
matrix perturbed by SHIELD.

To evaluate SHIELD’s performance on a realistic simulation of an imputation
query, we performed an ablation study on the 1KG Phase 3 [15] dataset. We
withheld 100 genomes (equivalent to 200 haplotypes) from the reference panel
to impute via the remaining 2,404 samples. The first 10,000 single-nucleotide
polymorphisms (SNPs) were extracted from 1KG; the remaining were discarded to
render run times more tractable. To simulate an array-based assay of the 200
target haplotypes, we ablated all sites except those included in the Illumina
Human1M-Duo v3.0 DNA Analysis BeadChip manifest, the intersection of which
with the first 10,000 sites in the 1KG data consisted of a total of 253 sites
for an a priori coverage of 2.53%.

To quantify accuracy, we summed the imputed dosages for each pair of
haplotypes to compute a final genotype dosage for each sample, then computed
the coefficient of determination ($R^{2}$) between the genotype dosages and
the ground-truth exome data. Because sites vary massively by minor allele
frequency (MAF), the loci were divided into three bins corresponding to MAFs
of $(0\%,0.5\%)$, $\left[0.5\%,5\%\right)$, and $\left[5\%,50\%\right]$.
Respectively, these bins contained 5,943, 2,157, and 1,900 variants in the
reference set. Accuracy was assessed, by bin, both to compare the performance
of SHIELD to that of Minimac3 [8] and to characterize the effect of the
privacy budget on our method’s accuracy.

Figure 3: A. Comparison between the accuracy by MAF of imputed dosages for
targets withheld from 1KG for both SHIELD (non-differentially private) and
Minimac3. B. SHIELD’s accuracy by MAF versus privacy budget.

Our analyses show nearly identical performance between SHIELD and Minimac3
when no input perturbation is applied, with the former obtaining scores of
0.571, 0.784, and 0.902, respectively, on the three bins enumerated above and
the latter scores of 0.584, 0.787, and 0.901 (Figure 2). SHIELD’s accuracy was
reevaluated at various values of our privacy budget along the interval
$\left[0.01,10\right]$, reflecting the typical range of values that
$\varepsilon$ is assigned in many differentially private algorithms [26].
Expectedly, accuracy exhibits a negative association with $\varepsilon$. At an
upper bound of $\varepsilon=10$, SHIELD performs nearly identically to
Minimac3 (0.564, 0.784, 0.901; Figure 3), while performance degrades
significantly at $\varepsilon=0.01$ (0.014, 0.038, 0.218; Figure 3).

###  2.3 Impact on Markov parameters

As noted above, the parameters for the Markov random field [34] modeling
genomic recombination [30], namely the mutation and recombination rates, were
computed on the unperturbed data by Minimac3 [8]. The rationale was that the
noise added to the reference panel mimicked the behavior of extremely rapid
genomic recombination, causing Minimac3’s expectation-maximization procedure
to dramatically overestimate the recombination rates (5.93 $\times 10^{-3}$
vs. 4.84$\times 10^{-4}$) and, conversely, to underestimate the mutation rates
(Figure 3B). These atypical rates exerted a decidedly negative impact on
imputation accuracy, with performance decreasing by 35.5%, 16.1%, and 5.46%
for each of the three bins, respectively, when the rates were computed on the
reference panel perturbed at $\varepsilon=5.0$. In sum, it is clearly superior
to estimate population parameters a priori, although, notably, doing so on the
reference panel itself is not differentially private and may leak information.

###  2.4 Impact on compression rates

An additional feature of haplotype imputation introduced by Minimac3 was the
M3VCF format for genomic data, which both substantially decreases total file
size over the traditional VCF format and enables the state-space reduction
technique that further improves imputation runtime [8]. The key insight
enabling the format is the observation that, due to identity-by-descent [5],
most haplotypes share identical $k$-mers of genomic material at intervals of
contiguous loci despite being unique overall. In other words, given an
arbitrary interval along the genome, the number of unique $k$-mers
collectively exhibited by the reference panel is almost always smaller than
the total number of reference haplotypes per se. Therefore, it is possible to
implement a compression scheme in which the genome is partitioned into
intervals and only the unique $k$-mer strings are retained, substantially
compressing the original reference panel [8].

An unfortunate consequence of local differential privacy via randomized
response is that, on average, random noise will destroy the exact equality
between haplotypes substrings. From the perspective of a compression algorithm
attempting to identify the set of unique $k$-mers along a given interval, an
apparently larger number of unique fragments will exist, rendering M3VCF-style
compression will less efficient. As an illustration, we partitioned the
genomic data into mutually exclusive, exhaustive blocks of uniform size
ranging from 2 to 500. We then computed the data compression ratio when M3VCF-
style state-space reduction was applied at each block size by dividing the
total $5.008\times 10^{8}$ bits in the uncompressed panel by the number of
bits following compression and plotted the ratio against block size (Figure
3C). Input perturbation resulted in compression rates up to an order of
magnitude smaller.

Figure 4: A. SHIELD’s imputed dosage accuracy $\left(R^{2}\right)$ by MAF
using parameters derived from both the original and perturbed reference
panels. B. The mean recombination and error rates using the original and
perturbed parameters. C. M3VCF-like compression ratio versus haplotype block
size on the original and perturbed parameters.

##  3 Discussion

In this work, we develop Secure Haplotype Imputation Employing Local
Differential privacy (SHIELD), a program for performing genomic imputation
with strong privacy guarantees for reference haplotypes via the randomized
response technique [33]. Analysis shows that SHIELD is able to obtain state-
of-the-art accuracy in realistic experimental settings at typical privacy
budgets.

We note that the strong performance of SHIELD parallels the effectiveness of
RAPPOR [37], a differentially private algorithm for mining strings in
commercial contexts that is also based on randomized response. Unlike SHIELD,
however, RAPPOR is not intended for data that is inherently binary; rather,
arbitrary alphanumeric strings are hashed onto Bloom filters [38] that are
subsequently perturbed. The fact that haplotype data intrinsically consist of
bitstrings makes randomized response particularly convenient in a genomic
context.

But despite the strong performance exhibited in the experiments above, it
should be acknowledged that the privacy guarantees made by our program are
limited to individual variants. In other words, for a given privacy budget
$\varepsilon$ [26, 27, 28], SHIELD can provably ensure protection for each
sample’s genotype at any one site, but not across the entire genome per se.
Certain adversarial attacks are therefore still feasible with SHIELD even
though accurate reconstruction of reference haplotypes is not [19, 20, 21, 22,
23]. Whole-genome privacy would instead require the division $\varepsilon$
across each site (see [27] for a discussion on composition in differential
privacy), which is prohibitively difficult for datasets containing tens of
thousands of variants. On the other hand, such divisions may be possible if a
fairly limited segment of the genome is to be imputed. Future research into
genomic privacy may investigate these scenarios or alternative differentially
private mechanisms.

A second limitation of our program is its dependence on accurate a priori
estimates of population parameters [30, 5, 8], which are non-trivial to
compute while still enforcing local differential privacy. Subsequent work may
inquire into the feasibility of computing population parameters a posteriori
by performing some manner of statistical correction.

Nevertheless, the capacity for basic differentially private mechanisms to
easily provide meaningful results is highly promising for the prospect of
privacy in practical genomic research.

##  4 Methods

The SHIELD algorithm consists of two subroutines, Perturb and Impute, that are
described below. The former is called once on a reference panel $\mathbf{X}$
to produce a locally [31] differentially private [26, 27, 28, 29] reference
panel $\mathbf{\tilde{X}}$ that is stored on the imputation server, whereas
the latter is then called by the client for each subsequent query haplotype
$\mathbf{z}$ using $\mathbf{\tilde{X}}$ as the reference panel.

###  4.1 Differential Privacy and Randomized Response

We derive the privacy guarantees of SHIELD from the notion of differential
privacy [26, 27, 28, 29]. Preliminarily, we develop the notion of neighboring
datasets. Given a universe of datasets $\mathcal{X}$, we say that two datasets
$x,y\in\mathcal{X}$ are neighbors if and only if they differ by at most one
individual sample. We will also call a randomized algorithm
$\mathcal{M}:\mathcal{X}\to\mathcal{F}$, where $\mathcal{F}$ is an arbitrary
probability space, a mechanism. We then say that a mechanism
$\mathcal{M}:\mathcal{X}\to\mathcal{F}$ satisfies
$\left(\epsilon,\delta\right)$-differential privacy if and only if for all
$\mathcal{S}\subseteq\mathcal{F}$ and for all $x,y\in\mathcal{X}$ such that
$x$ are $y$ are neighboring, we have

$P\left(\mathcal{M}\left(x\right)\in\mathcal{S}\right)\leq\exp\left(\varepsilon\right)P\left(\mathcal{M}\left(y\right)\in\mathcal{S}\right)+\delta.$
(1)

Among the most common techniques in differential privacy, randomized response
[32, 33] satisfies $\epsilon$-differential privacy for binary attributes. The
randomized response scheme on a binary attribute $X$ is a mechanism
$\mathcal{M}_{rr}:\left\\{0,1\right\\}\to\left\\{0,1\right\\}$ is
characterized by a $2\times 2$ distortion matrix

$\mathbf{P}=\begin{pmatrix}p_{00}&p_{01}\\\ p_{10}&p_{11}\end{pmatrix},$ (2)

where
$p_{uv}=P\left(\mathcal{M}_{rr}(x_{i})=u|x_{i}=v\right)\quad\left(u,v\right)\in\left\\{0,1\right\\}$.
It can be shown [32] that the highest-utility value for $\mathbf{P}$ is

$\mathbf{P}=\begin{pmatrix}\frac{e^{\varepsilon}}{1+e^{\varepsilon}}&\frac{1}{1+e^{\varepsilon}}\\\
\frac{1}{1+e^{\varepsilon}}&\frac{e^{\varepsilon}}{1+e^{\varepsilon}}\end{pmatrix}.$
(3)

Fixing the number of samples in our reference panel $n$ and the number of
sites $m$, we denote the universe of possible reference panels
$\mathcal{X}=\left\\{0,1\right\\}^{m\times n}$. Because haplotypes are vector-
valued, applying the notion of neighboring datasets is non-trivial. For our
purposes, we will say that two reference panels
$\mathbf{X},\mathbf{X}^{\prime}\in\mathcal{X}$ are neighboring if and only if
their Hamming distance is less than or equal to $1$. In other words, we
consider $\mathbf{X}$ and $\mathbf{X^{\prime}}$ neighbors if and only if
$X_{i,j}\neq X^{\prime}_{i,j}$ for a single marker $i$ and a single individual
$j$ as opposed to a whole-genome interpretation of neighboring datasets in
which $\mathbf{X}$ and $\mathbf{X}^{\prime}$ may differ by an entire row.

It then follows that by applying the randomized response mechanism
$\mathcal{M}_{rr}$ to each entry in a reference panel matrix $\mathbf{X}$, we
may store a perturbed copy $\mathbf{\tilde{X}}$ of the original reference
panel that satisfies entry-wise $\varepsilon$-differential privacy. The
perturbation step of SHIELD then consists of the procedure Perturb. We note
that we use the symbol $\overset{r}{\leftarrow}$ to denote a pseudorandom
sample and $\text{Bern}\left(\vartheta\right)$ to denote a Bernoulli
distribution with parameter $\vartheta$.

Algorithm 1 Applies randomized response mechanism to reference panel.

1:procedure Perturb($\mathbf{X},\varepsilon$)

2: $\mathbf{\tilde{X}}\leftarrow$ empty matrix

3: for $i=1,2,\ldots,n$ do

4: for $j=1,2,\ldots,m$ do

5:
$c\overset{r}{\leftarrow}\text{Bern}\left(\frac{e^{\varepsilon}}{1+e^{\varepsilon}}\right)$

6: if $c=1$ then

7: $\tilde{X}_{i,j}\leftarrow X_{i,j}$

8: else

9: $\tilde{X}_{i,j}\leftarrow\lnot X_{i,j}$

10: end if

11: end for

12: end for

13: return $\mathbf{\tilde{X}}$

14:end procedure

A convenient property of differential privacy is post-processing [26]. If
$\mathcal{M}:\mathcal{X}\to\mathcal{F}$ is an
$\left(\varepsilon,\delta\right)$-differentially private randomized algorithm
and $f:\mathcal{F}\to\mathcal{F}^{\prime}$ is an arbitrary mapping, then
$f\circ\mathcal{M}:\mathcal{X}\to\mathcal{F}^{\prime}$ is
$\left(\varepsilon,\delta\right)$-differentially private. We set
$\mathcal{M}=\mathcal{M}_{rr}$ and define $f$ such that
$f\left(\mathbf{\tilde{X}}\right)=\textsc{Impute}\left(\mathbf{z},\mathbf{\tilde{X}},\boldsymbol{\mu},\boldsymbol{\rho}\right)$
for some fixed values $\mathbf{z}$, $\boldsymbol{\mu}$, and
$\boldsymbol{\rho}$ (see below on the meaning of these parameters). Then by
post-processing, it follows that each call to Impute on the perturbed
reference panel $\mathbf{\tilde{X}}$ will satisfy $\varepsilon$-differential
privacy. In other words, once $\mathbf{\tilde{X}}$ has been collected on the
imputation server and perturbed so as to satisfy local differential privacy,
an unlimited number of queries are able to be made by an algorithmic adversary
without divulging any one haplotype’s value at any one site with a high degree
of certainty.

###  4.2 HMM-based genotype imputation

We will also use the following notation:

  * •

0, 1, and $\varnothing$: the minor allele, major allele, and constant denoting
an unobserved site to be imputed;

  * •

$n$ and $m$: the number of reference samples and reference markers;

  * •

$\left[n\right]$: the set of reference haplotypes, represented as the index
set;

  * •

$\mathbf{X}=\left[\mathbf{x}_{1},\mathbf{x}_{2},\ldots,\mathbf{x}_{m}\right]^{\intercal}\in\left\\{0,1\right\\}^{m\times
n}$: the reference panel haplotype sequences, equivalent to a real (and more,
specifically, binary) matrix;

  * •

$\left(z_{k}\right)_{k=1}^{m}\in\left\\{0,1,\varnothing\right\\}^{m}$: the
sequence corresponding to the observed target haplotype that, because it may
include the missing site letter, is not, strictly speaking, a real vector;

  * •

$\left(\hat{z}_{k}\right)_{k=1}^{m}\equiv\mathbf{\hat{z}}\in\left[0,1\right]^{m}$:
the sequence of imputed haplotype dosages, equivalent to a real vector;

  * •

$\mathbf{y}\in\left[n\right]^{m}$: the site-wise identities of the reference
haplotypes from which $\mathbf{z}$ is descended;

  * •

$\boldsymbol{\rho}\in\left[0,1\right]^{m}$: the recombination rates [30, 5]
such that $\rho_{i}=P(y_{i+1}=j_{2}|y_{i}=j_{1})\quad j_{2}\neq j_{1}$,
meaning that $\boldsymbol{\rho}$ is equivalent to a real vector (we simply let
$\rho_{m}=0$ as a dummy value);

  * •

$\boldsymbol{\mu}\in\left[0,1\right]^{m}$: the mutation rates [30, 5] such
$\mu_{i}=P(z_{i}\neq X_{i,j}|y_{i}=j)$, meaning that $\boldsymbol{\mu}$ is
equivalent a real vector;

  * •

$\mathbf{M}=\left[\mathbf{m}_{1},\mathbf{m}_{2},\ldots,\mathbf{m}_{m}\right]^{\intercal}\in\left[0,1\right]^{m\times
n}$: the emission probabilities in matrix form such that $\mu_{i}=P(z_{i}\neq
X_{i,j}|y_{i}=j)$ such that

$M_{i,j}=\begin{cases}1-\mu_{i}\quad&\text{if}\quad z_{i}=X_{i,j}\\\
\mu_{i}\quad&\text{if}\quad z_{i}=\varnothing\quad\text{or}\quad z_{i}\neq
X_{i,j}\end{cases};$ (4)

  * •

$\boldsymbol{\Gamma}=\left[\boldsymbol{\gamma}_{1},\boldsymbol{\gamma}_{2},\ldots,\boldsymbol{\gamma}_{m}\right]^{\intercal}\in\left\\{0,1\right\\}^{m\times
n}$: the posterior probabilities for haplotype identity for all sites in
matrix form $\Gamma_{i,j}=P\left(y_{i}=j|\left(z_{k}\right)_{k=1}^{m}\right)$;

  * •

$\mathbf{A}=\left[\boldsymbol{\alpha}_{1},\boldsymbol{\alpha}_{2},\ldots,\boldsymbol{\alpha}_{m}\right]^{\intercal}\in\left\\{0,1\right\\}^{m\times
n}$: the forward probabilities [35] for all sites in matrix form such that
$A_{i,j}=P(y_{i}=j|\left(z_{k}\right)_{k=1}^{i})$;

  * •

$\mathbf{B}=\left[\boldsymbol{\beta}_{1},\boldsymbol{\beta}_{2},\ldots,\boldsymbol{\beta}_{m}\right]^{\intercal}\in\left\\{0,1\right\\}^{m\times
n}$: the backward probabilities [35] for all sites in matrix form such that
$B_{i,j}=P\left(y_{i}=j,\left(z_{k}\right)_{k=i+1}^{m}\right)$;

By the Law of the Unconscious Statistician, we compute the expected value
(i.e., dosage) for each $i$th site via the inner product of the posterior
probabilities $\boldsymbol{\gamma}_{i}$ over the reference haplotype space
$y_{i}$ and the values $\boldsymbol{x}_{i}$ that the reference haplotypes have
at the $i$th site:

$\text{E}\left[\hat{z}_{i}|\left(z_{k}\right)_{k=1}^{m};\mathbf{X},\boldsymbol{\mu},\boldsymbol{\rho}\right]=\mathbf{x}_{i}^{\intercal}\boldsymbol{\gamma}_{i}.$
(5)

As is critical to the forward-backward algorithm [35], we note that

$\boldsymbol{\gamma}_{i}\propto\boldsymbol{\alpha}_{i}\circ\boldsymbol{\beta}_{i},$
(6)

where $\circ$ denotes the Hadamard product. However, $\boldsymbol{\gamma}_{i}$
must sum to unity, meaning that to compute the exact posterior distribution we
normalize the two messages via

$\boldsymbol{\gamma}_{i}=\boldsymbol{\alpha}_{i}\circ\boldsymbol{\beta}_{i}\left(\boldsymbol{\alpha}_{i}^{\intercal}\boldsymbol{\beta}_{i}\right)^{-1}.$
(7)

For all sites in matrix form, this is equivalent to

$\text{E}\left[\mathbf{\hat{z}}|\left(z_{k}\right)_{k=1}^{m};\mathbf{X},\boldsymbol{\mu},\boldsymbol{\rho}\right]=\left(\mathbf{X}\circ\boldsymbol{\Gamma}\right)\mathbf{1}$
(8)

where

$\boldsymbol{\Gamma}=\left(\mathbf{A}\circ\mathbf{B}\right)\left(\text{diag}\left(\left(\mathbf{A}\circ\mathbf{B}\right)^{\intercal}\boldsymbol{1}\right)\right)^{-1}.$
(9)

As to computing the forward and backward messages, their recurrence relation
from the Li-Stephens model [30, 5] can be shown to simplify, respectively, to

$\begin{cases}\boldsymbol{\alpha}_{1}=\mathbf{1}\\\
\boldsymbol{\alpha}_{i+1}=\mathbf{m}_{i+1}\circ\left(\frac{\rho_{i+1}}{n}\sum_{j=1}^{n}A_{i,j}+\left(1-\rho_{i}\right)\boldsymbol{\alpha}_{i}\right)\end{cases}$
(10)

and

$\begin{cases}\boldsymbol{\beta}_{m}=\mathbf{1}\\\
\boldsymbol{\beta}_{i}=\mathbf{m}_{i}\circ\left(\frac{\rho_{i}}{n}\sum_{j=1}^{n}B_{i+1,j}+\left(1-\rho_{i+1}\right)\boldsymbol{\beta}_{i+1}\right)\end{cases}.$
(11)

This gives rise to the final implementation of Impute, which computes the
forward and backward messages successively via dynamic programming [39].

Algorithm 2 Uses forward-algorithm to impute dosages according to Li-Stephens
model.

1:procedure Impute($\mathbf{z},\mathbf{X},\boldsymbol{\mu},\boldsymbol{\rho}$)

2: $\mathbf{A},\mathbf{B},\mathbf{M}\leftarrow$ empty matrices

3: for $i=1,2,\ldots,m$ do

4: for $j=1,2,\ldots,n$ do

5: if $z_{i}=X_{i,j}$ then

6: $M_{i,j}\leftarrow 1-\mu_{i}$

7: else

8: $M_{i,j}\leftarrow\mu_{i}$

9: end if

10: end for

11: end for

12: $\boldsymbol{\alpha}_{1}\leftarrow\boldsymbol{1}$

13: for $i=1,2,\ldots,m$ do

14:
$\boldsymbol{\alpha}_{i+1}=\mathbf{m}_{i+1}\circ\left(\frac{\rho_{i+1}}{n}\sum_{j=1}^{n}A_{i,j}+\left(1-\rho_{i}\right)\boldsymbol{\alpha}_{i}\right)$

15: end for

16: $\boldsymbol{\beta}_{m}\leftarrow\boldsymbol{1}$

17: for $i=m-1,m-2,\ldots,1$ do

18:
$\beta_{i}=\mathbf{m}_{i}\circ\left(\frac{\rho_{i}}{n}\sum_{j=1}^{n}B_{i+1,j}+\left(1-\rho_{i+1}\right)\boldsymbol{\beta}_{i+1}\right)$

19: end for

20: $\boldsymbol{\Gamma}\leftarrow\mathbf{A}\circ\mathbf{B}$

21:
$\boldsymbol{\Gamma}\leftarrow\boldsymbol{\Gamma}\left(\text{diag}\left(\boldsymbol{\Gamma}^{\intercal}\boldsymbol{1}\right)\right)^{-1}$

22:
$\mathbf{\hat{z}}\leftarrow\left(\mathbf{X}\circ\boldsymbol{\Gamma}\right)\mathbf{1}$

23: return $\mathbf{\hat{z}}$

24:end procedure

## Acknowledgements

We thank Hoon Cho at the Broad Institute of MIT and Harvard for insightful
discussion.

## References

  * [1] I. H. Consortium _et al._ , “A second generation human haplotype map of over 3.1 million snps,” _Nature_ , vol. 449, no. 7164, p. 851, 2007.
  * [2] J. Marchini, B. Howie, S. Myers, G. McVean, and P. Donnelly, “A new multipoint method for genome-wide association studies by imputation of genotypes,” _Nature genetics_ , vol. 39, no. 7, pp. 906–913, 2007.
  * [3] B. Servin and M. Stephens, “Imputation-based analysis of association studies: candidate regions and quantitative traits,” _PLoS genetics_ , vol. 3, no. 7, p. e114, 2007.
  * [4] D. Taliun, D. N. Harris, M. D. Kessler, J. Carlson, Z. A. Szpiech, R. Torres, S. A. G. Taliun, A. Corvelo, S. M. Gogarten, H. M. Kang _et al._ , “Sequencing of 53,831 diverse genomes from the nhlbi topmed program,” _Nature_ , vol. 590, no. 7845, pp. 290–299, 2021.
  * [5] Y. Li, C. Willer, S. Sanna, and G. Abecasis, “Genotype imputation,” _Annual review of genomics and human genetics_ , vol. 10, pp. 387–406, 2009\.
  * [6] J. Marchini and B. Howie, “Genotype imputation for genome-wide association studies,” _Nature Reviews Genetics_ , vol. 11, no. 7, pp. 499–511, 2010\.
  * [7] Y. Li, C. J. Willer, J. Ding, P. Scheet, and G. R. Abecasis, “Mach: using sequence and genotype data to estimate haplotypes and unobserved genotypes,” _Genetic epidemiology_ , vol. 34, no. 8, pp. 816–834, 2010.
  * [8] S. Das, L. Forer, S. Schönherr, C. Sidore, A. E. Locke, A. Kwong, S. I. Vrieze, E. Y. Chew, S. Levy, M. McGue _et al._ , “Next-generation genotype imputation service and methods,” _Nature genetics_ , vol. 48, no. 10, pp. 1284–1287, 2016.
  * [9] D. L. Ayres, A. Darling, D. J. Zwickl, P. Beerli, M. T. Holder, P. O. Lewis, J. P. Huelsenbeck, F. Ronquist, D. L. Swofford, M. P. Cummings _et al._ , “Beagle: an application programming interface and high-performance computing library for statistical phylogenetics,” _Systematic biology_ , vol. 61, no. 1, pp. 170–173, 2012.
  * [10] S. Purcell, B. Neale, K. Todd-Brown, L. Thomas, M. A. Ferreira, D. Bender, J. Maller, P. Sklar, P. I. De Bakker, M. J. Daly _et al._ , “Plink: a tool set for whole-genome association and population-based linkage analyses,” _The American journal of human genetics_ , vol. 81, no. 3, pp. 559–575, 2007.
  * [11] P. Scheet and M. Stephens, “A fast and flexible statistical model for large-scale population genotype data: applications to inferring missing genotypes and haplotypic phase,” _The American Journal of Human Genetics_ , vol. 78, no. 4, pp. 629–644, 2006.
  * [12] C. Sudlow, J. Gallacher, N. Allen, V. Beral, P. Burton, J. Danesh, P. Downey, P. Elliott, J. Green, M. Landray _et al._ , “Uk biobank: an open access resource for identifying the causes of a wide range of complex diseases of middle and old age,” _PLoS medicine_ , vol. 12, no. 3, p. e1001779, 2015\.
  * [13] A. of Us Research Program Investigators, “The “all of us” research program,” _New England Journal of Medicine_ , vol. 381, no. 7, pp. 668–676, 2019.
  * [14] “A reference panel of 64,976 haplotypes for genotype imputation,” _Nature genetics_ , vol. 48, no. 10, pp. 1279–1283, 2016.
  * [15] N. Siva, “1000 genomes project,” _Nature biotechnology_ , vol. 26, no. 3, pp. 256–257, 2008.
  * [16] S. Das, G. R. Abecasis, and B. L. Browning, “Genotype imputation from large reference panels,” _Annual review of genomics and human genetics_ , vol. 19, pp. 73–96, 2018.
  * [17] P.-R. Loh, P. Danecek, P. F. Palamara, C. Fuchsberger, Y. A Reshef, H. K Finucane, S. Schoenherr, L. Forer, S. McCarthy, G. R. Abecasis _et al._ , “Reference-based phasing using the haplotype reference consortium panel,” _Nature genetics_ , vol. 48, no. 11, pp. 1443–1448, 2016\.
  * [18] N. Dokmai, C. Kockan, K. Zhu, X. Wang, S. C. Sahinalp, and H. Cho, “Privacy-preserving genotype imputation in a trusted execution environment,” _Cell systems_ , vol. 12, no. 10, pp. 983–993, 2021.
  * [19] L. Bonomi, Y. Huang, and L. Ohno-Machado, “Privacy challenges and research opportunities for genomic data sharing,” _Nature genetics_ , vol. 52, no. 7, pp. 646–654, 2020.
  * [20] S. Sankararaman, G. Obozinski, M. I. Jordan, and E. Halperin, “Genomic privacy and limits of individual detection in a pool,” _Nature genetics_ , vol. 41, no. 9, pp. 965–967, 2009.
  * [21] N. Homer, S. Szelinger, M. Redman, D. Duggan, W. Tembe, J. Muehling, J. V. Pearson, D. A. Stephan, S. F. Nelson, and D. W. Craig, “Resolving individuals contributing trace amounts of dna to highly complex mixtures using high-density snp genotyping microarrays,” _PLoS genetics_ , vol. 4, no. 8, p. e1000167, 2008.
  * [22] Y. Erlich, T. Shor, I. Pe’er, and S. Carmi, “Identity inference of genomic data using long-range familial searches,” _Science_ , vol. 362, no. 6415, pp. 690–694, 2018.
  * [23] M. Gymrek, A. L. McGuire, D. Golan, E. Halperin, and Y. Erlich, “Identifying personal genomes by surname inference,” _Science_ , vol. 339, no. 6117, pp. 321–324, 2013.
  * [24] H. Cho, D. J. Wu, and B. Berger, “Secure genome-wide association analysis using multiparty computation,” _Nature biotechnology_ , vol. 36, no. 6, pp. 547–551, 2018.
  * [25] M. Naveed, E. Ayday, E. W. Clayton, J. Fellay, C. A. Gunter, J.-P. Hubaux, B. A. Malin, and X. Wang, “Privacy in the genomic era,” _ACM Computing Surveys (CSUR)_ , vol. 48, no. 1, pp. 1–44, 2015.
  * [26] C. Dwork, “Differential privacy,” in _International colloquium on automata, languages, and programming_. Springer, 2006, pp. 1–12.
  * [27] ——, “Differential privacy: A survey of results,” in _International conference on theory and applications of models of computation_. Springer, 2008, pp. 1–19.
  * [28] F. K. Dankar and K. El Emam, “Practicing differential privacy in health care: A review.” _Trans. Data Priv._ , vol. 6, no. 1, pp. 35–67, 2013.
  * [29] C. Dwork, A. Roth _et al._ , “The algorithmic foundations of differential privacy,” _Foundations and Trends® in Theoretical Computer Science_ , vol. 9, no. 3–4, pp. 211–407, 2014.
  * [30] N. Li and M. Stephens, “Modeling linkage disequilibrium and identifying recombination hotspots using single-nucleotide polymorphism data,” _Genetics_ , vol. 165, no. 4, pp. 2213–2233, 2003.
  * [31] M. Yang, L. Lyu, J. Zhao, T. Zhu, and K.-Y. Lam, “Local differential privacy and its applications: A comprehensive survey,” _arXiv preprint arXiv:2008.03686_ , 2020.
  * [32] Y. Wang, X. Wu, and D. Hu, “Using randomized response for differential privacy preserving data collection.” in _EDBT/ICDT Workshops_ , vol. 1558, 2016, pp. 0090–6778.
  * [33] S. L. Warner, “Randomized response: A survey technique for eliminating evasive answer bias,” _Journal of the American Statistical Association_ , vol. 60, no. 309, pp. 63–69, 1965.
  * [34] L. Rabiner and B. Juang, “An introduction to hidden markov models,” _ieee assp magazine_ , vol. 3, no. 1, pp. 4–16, 1986.
  * [35] L. E. Baum _et al._ , “An inequality and associated maximization technique in statistical estimation for probabilistic functions of markov processes,” _Inequalities_ , vol. 3, no. 1, pp. 1–8, 1972.
  * [36] G. Cormode, S. Jha, T. Kulkarni, N. Li, D. Srivastava, and T. Wang, “Privacy at scale: Local differential privacy in practice,” in _Proceedings of the 2018 International Conference on Management of Data_ , 2018, pp. 1655–1658.
  * [37] Ú. Erlingsson, V. Pihur, and A. Korolova, “Rappor: Randomized aggregatable privacy-preserving ordinal response,” in _Proceedings of the 2014 ACM SIGSAC conference on computer and communications security_ , 2014, pp. 1054–1067.
  * [38] B. H. Bloom, “Space/time trade-offs in hash coding with allowable errors,” _Communications of the ACM_ , vol. 13, no. 7, pp. 422–426, 1970.
  * [39] R. Bellman, “Dynamic programming,” _Science_ , vol. 153, no. 3731, pp. 34–37, 1966.