# Towards more realistic models of genomes in populations: the Markov-
modulated sequentially Markov coalescent

Julien Y. Dutheil  dutheil@@evolbio.mpg.de RG Molecular Systems Evolution, Max
Planck Institute for Evolutionary Biology, August-Thienemann-Str. 2, 24306
Plön, GERMANY

The development of coalescent theory paved the way to statistical inference
from population genetic data. In the genomic era, however, coalescent models
are limited due to the complexity of the underlying ancestral recombination
graph. The sequentially Markov coalescent (SMC) is a heuristic that enables
the modelling of complete genomes under the coalescent framework. While it
empowers the inference of detailed demographic history of a population from as
few as one diploid genome, current implementations of the SMC make unrealistic
assumptions about the homogeneity of the coalescent process along the genome,
ignoring the intrinsic spatial variability of parameters such as the
recombination rate. Here, I review the historical developments of SMC models
and discuss the evidence for parameter heterogeneity. I then survey approaches
to handle this heterogeneity, focusing on a recently developed extension of
the SMC.

##  1\. Modelling the evolution of genomes in populations

When modelling the evolution of large genomic sequences at the population
level, in particular for sexually reproducing species, a key biological
mechanism to account for is meiotic recombination, which shuffles genetic
material at each generation. We first introduce the concept of the ancestral
recombination graph, needed to represent the complete genealogy of a sample
undergoing recombination. We then review the statistical approaches used to
fit models accounting for recombination to population genomics data.

###  1.1. The ancestral recombination graph

The evolution of the set of sequences carried by all individuals forming a
population, generation after generation, can be modelled by a stochastic
process, where each individual leaves a variable number of descendants in the
next generation. As a result, at any position of the sequence, the _genealogy_
of a sample of $n$ individuals can be described by a tree (Figure 1A) [35].
The tips of the tree represent the sampled individuals and the inner nodes
their common ancestors. In the case of sexually reproducing organisms, which
will be the focus of this chapter, the genealogy is not identical for every
position in the sequence. During _sexual reproduction_ , two individuals
contribute part of their sequence to their descendant(s) in the next
generation. The mechanism of _recombination_ is responsible for randomly
sampling the new sequence from the two parental ones (Figure 1B). How often
and where the recombination points occur will be discussed in Section 2. The
consequences of the recombination process can be stated as: (1) the genealogy
of the sequence on the left of a recombination breakpoint potentially differs
from the genealogy of the sequence on the right, (2) the genealogy at two
positions in the sequence are more likely to differ as the distance between
the two points is large, and (3) the genealogy of the complete sequence can no
longer be described by a single tree, but by a collection of such trees and
associated breakpoints. This tree and breakpoints collection can be
represented as a single graph, called the _ancestral recombination graph
(ARG)_ (Figure 1C) [25]. The complexity of the ARG grows with the number of
individuals (which dictates the size of the underlying trees) and the number
of recombination events (which determines how many trees are needed to
represent the history of the full sequences). The ARG represents the complete
history of the sampled individuals, where the trees at each position (referred
to as the “marginal genealogies”)are embedded [49]. It describes the history
of each segment of the sampled sequences, tracing back their ancestors in
potentially distinct individuals. Such segments, which have left descendants
in the sample, are termed _ancestral_. Contained in the ARG is also the
history of some _non-ancestral segments_ , which did not leave a descendant in
the sample, but were once part of a sequence that contained both ancestral and
non-ancestral segments (Figure 1D).

The characteristics of the ARG are determined by the demography of the
population (the history of population size changes), the recombination
landscape (where do the recombination events occur), but also the selective
forces acting on the sequence, as natural selection influences the
distribution of the number of descendants for each individual, based on the
nature of the sequences themselves. While the ARG contains the signature of
the biological processes that shaped the genome sequences, it is unfortunately
not directly accessible. In order to access the embedded information, it is
necessary to model the evolution of sequences in populations.

Figure 1.  Genealogies and recombination: relationships between individuals
and along sequences. A) example genealogy of a sample of five individuals at a
given position in the sequences, under a Kingman coalescent. Tip nodes
represent the samples (1-5) and inner nodes their common ancestors (6-9). B)
illustration of the recombination tree resulting from a process without
coalescence (large number approximation): two chromosomes (solid black and
dotted grey) are paired during sexual reproduction and exchange segments at a
breakpoint $x_{1}$. At the next generation, the descendant sequence recombines
with another sequence (solid grey) at another breakpoint $x_{2}$, etc. The
sequence of any sampled individual is therefore a mosaic of segments with
distinct ancestors separated by a series of breakpoints ($x_{1}$, $x_{2}$,
$x_{3}$). C) a simple ancestral recombination graph (ARG) representing the
genealogy of a sample of size three with one recombination event. The ARG is a
combination of a coalescence tree (as in A) and a recombination tree (as in
B). For ease of interpretation, two mutation events, $\bullet a$ and $\bullet
b$ have been added. The relative coordinate of the recombination event is also
indicated: $x=0.4$, assuming a total sequence length of 1. D) partial graph
showing the different classes of recombination events. Ancestral segments are
depicted as filled rectangles, while non-ancestral segments are shown as
simple lines (subfigure created after Figure 1 in [47]).

###  1.2. The coalescent with recombination as a chronological process

When modelling evolution, the most intuitive approach is to consider the
process chronologically, that is, to model the state of the system generation
after generation. One of the most simple models, the so-called _Wright–Fisher
process_ , considers that the gametes forming one generation are a random
sample of the gametes produced at the previous generation, that is,
reproduction is a purely random process where each individual has the same a
priori chance to contribute to the next generation. In addition, the
population has a finite, constant size. A similar model, termed the _Moran
process_ , considers a slightly different set-up with overlapping generations
[50]. The Wright–Fisher and Moran processes can both accommodate
recombination, modelled by randomly choosing a breakpoint along the genome and
exchanging the parental genetic segments (see the contribution of Baake and
Baake [6] in this volume). In such processes, the fate of a genetic variant is
purely stochastic and governed only by the population size.

When conditioning on a sample of the result of the evolutionary process, a
backwards-in-time modelling is used. Each sequence in the sample represents a
lineage, and the aim of the model is to determine which lineages find a common
ancestor in the past and when. Every time two lineages _coalesce_ into a
common ancestor, the number of lineages to model is reduced by one, until the
last common ancestor of the sample is reached. This process of lineages
merging backwards in time is termed the _coalescent_ [35].

The probability that two lineages merge at a given generation back in time
depends on the population size. When the population size is constant in time,
the number of generations until coalescence follows a geometric distribution
with parameter $\frac{1}{2\cdot N_{e}}$, where $N_{e}$ is the _effective
population size_. Assuming a large $N_{e}$, this discrete-time coalescent is
well approximated by a continuous-time coalescent, where the divergence time
between two sequences follows an exponential distribution with average $2\cdot
N_{e}$ generations. For convenience, time is, therefore, measured in
“coalescence” units equal to $2\cdot N_{e}$, so that the mean divergence time
between two sequences in a sample is equal to 1.

In the _coalescent with recombination_ process, recombination events are
modelled in addition to coalescence events [30]. Backwards-in-time, a lineage
undergoing a recombination event splits in two, the left and right sequences
having distinct ancestors. Since the rates of coalescence and recombination
events, at any time point, depend only on the current lineages, the process is
Markovian in time [63]. This property enabled the development of simulation
procedures and inference methods, allowing the estimation of various
parameters by integrating over the unknown genealogy of a sample (e.g. [2,
14]). Such methods, however, do not scale well with the length of the modelled
sequences, as the number of events in the underlying ARG grows with the
sequence length [20], preventing efficient integration even with Markov chains
Monte-Carlo [75]. These methods are, therefore, restricted to small samples
with relatively few loci.

###  1.3. The coalescent with recombination as a sequential process

Following the initial work by Simonsen and Churchill [63], Wiuf and Hein
extended the two-loci model of coalescence with recombination to multiple loci
[82]. The resulting process models the ARG sequentially along the genome
rather than chronologically. The resulting _sequential coalescent with
recombination_ aims at modelling the genealogy of the sample at position $i$
given the genealogies at previous positions. In fact, genealogies at two
distinct positions in the sequence are not independent: they are identical if
no recombination event occurred between the two positions since the last
common ancestor of the sample and can only differ if at least one
recombination event occurred. Despite this intuitive correlation structure,
the coalescent with recombination is not Markovian along the sequence.
Computing the probability distribution of the marginal genealogy at a given
position proved to be quite challenging because of long-range dependencies,
the genealogy at position $i$ depending not only on the genealogy at position
$i-1$, but on the genealogy at all positions $1$ to $i-1$.

With the goal to simplify the likelihood calculation under the coalescent with
recombination, McVean and Cardin proposed an approximation where certain types
of coalescence events are ignored [49]. An intuitive description of the
simplified process was provided by Marjoram and Wall [47], who recognised five
types of recombination events on the ARG, based on the type of segments in the
parental sequences on both sides of the recombination event (Figure 1D): type
1 events occur in ancestral segments (events at $x_{3}$ and $x_{4}$ in Figure
1D) while types 2-5 occur in non-ancestral segments. Type 2 events occur in
so-called _trapped genetic material_ [57], that is, non-ancestral segments
flanked on both sides by ancestral segments (event at $x_{1}$ in Figure 1D).
Events of types 3, 4, and 5 occur in non-ancestral segments only flanked by
non-ancestral segments on one or both sides (e.g. event at $x_{2}$ in Figure
1D). Such events (3, 4 and 5) do not affect the sample generated by the
corresponding ARG, and therefore do not impact the likelihood of the sample
given the ARG. They can therefore be ignored without introducing any
additional hypotheses, see [82] (types 4 and 5) and [31] (types 3, 4 and 5).
The process of McVean and Cardin, which was further improved by Marjoram and
Wall [47] and Hobolth and Jensen [29] also ignores type 2 recombination
events, that is, recombination events occurring in trapped genetic material
[57]. Doing so implies ignoring potential long-range dependencies between
loci, and the distribution of samples generated by this approximated process
is, therefore, different from that of the standard coalescent with
recombination. The approximated process, however, has the additional property
that the distribution of genealogies at position $i$ is only dependent on the
genealogy at position $i-1$, and is, therefore, Markovian along the sequence.
Such a process is referred to as the _sequentially Markov coalescent (SMC)_
[49, 47]. Importantly, the SMC process generates samples with patterns of
genetic diversity that are very similar to the ones generated by the full
coalescent process [49]. The SMC, in particular, can be seen as a first-order
Markov approximation of the true coalescent with recombination process [81],
and higher order extensions have been introduced [66]. Furthermore, the Markov
property enables very efficient likelihood calculation using dynamic
programming algorithms to integrate over all ARGs. Such methodology comes from
the field of _hidden Markov models (HMM)_ , which we introduce in the next
section.

###  1.4. Coalescent hidden Markov models

Figure 2.  Chronology of sequentially Markov coalescent (SMC) and coalescent
hidden Markov models (CoalHMM). PSMC: pairwise SMC. MSMC: multiple SMC. ILS:
incomplete lineage sorting. I: isolation model. IM: isolation with migration
model.

Because of the SMC approximation, likelihood calculation under a coalescent
with recombination process represents a classical bioinformatic problem where
the probability of an observed state in the sequence depends on an unobserved
state, which is then said to be hidden. In the case studied here, the observed
states are sequence polymorphisms (between 2 or more individual sequences) and
the hidden states are the underlying marginal genealogies. HMMs have been
broadly used in sequence analysis [15]. _Coalescent hidden Markov models
(CoalHMM)_ refer to HMMs where the hidden states are genealogies. It was
introduced by Hobolth et al [27] as a name of the first model developed, which
we introduce later in this section, but was then extended to generally
encompass a full class of models [65] (Figure 2).

We note as $\\{\boldsymbol{\mathcal{A}}_{i}\\}_{1\leq i\leq L}$ the site-
specific random variable of observed states in a sample of $M$ sequences of
length $L$. Such states (noted ${\\{A_{g}\\}}_{1\leq g\leq S}$) are, in the
general case, a combination of the four nucleotides $A$, $C$, $G$ and $T$ (one
per modelled sequence), with the possibility to additionally account for
missing data (coded as $N$), so that ${\\{A_{g}\\}}_{1\leq g\leq
S}\in{\\{A,C,G,T,N\\}}^{M}$ and $S<5^{M}$ because of symmetry relationships
between trees making some of them unidentifiable. We note as
$\\{x_{i}\\}_{1\leq i\leq L}$ a particular realisation of
$\\{\boldsymbol{\mathcal{A}}_{i}\\}_{1\leq i\leq L}$, that is, the sequence
data. Furthermore, we note as $\\{\boldsymbol{\mathcal{H}}_{i}\\}_{1\leq i\leq
L}$ the site-specific random variable describing the marginal genealogies at
each position in the sequences. In the general case, such genealogies are
rooted trees with $M$ leaves.

In HMM terminology, the probabilities of observing the sequence data $x_{i}$
at a given position $i$ given a realisation of $\boldsymbol{\mathcal{H}}_{i}$,
$\operatorname{Pr}(\boldsymbol{\mathcal{A}}_{i}\mid\boldsymbol{\mathcal{H}}_{i})$,
are called the _emission probabilities_. $\boldsymbol{\mathcal{H}}_{i}$ is a
random variable that has a continuous distribution. To make likelihood
calculations tractable, this distribution is discretized, so that
$\boldsymbol{\mathcal{H}}_{i}$ can take a finite number $n$ of hidden states,
$\\{H_{j}\\}_{1\leq j\leq n}$. Under a discretised distribution of hidden
states, the emission probabilities for each position $i$, hidden state $k$ can
be more explicitly written as:

(1)
$e_{i,k}(x)=\operatorname{Pr}(\boldsymbol{\mathcal{A}}_{i}=x\mid\boldsymbol{\mathcal{H}}_{i}=H_{k}).$

We further introduce the so-called _transition probability_ of a genealogy
$H_{j}$ at position $i-1$ to a genealogy $H_{k}$ at position $i$ as

(2)
$q_{i,j,k}=\operatorname{Pr}(\boldsymbol{\mathcal{H}}_{i}=H_{k}\mid\boldsymbol{\mathcal{H}}_{i-1}=H_{j}).$

Given the set of emission and transition probabilities we can then write the
likelihood of the data by recursion. We define $F_{i,k}$ the joint probability
of the data (observed states) $x_{1},\ldots,x_{i}$ at positions 1 to $i$ and
the ancestral genealogy $H_{k}$ at position $i$ as:

(3)
$F_{i,k}=\operatorname{Pr}(x_{1},\ldots,x_{i},H_{k})=\left\\{\begin{array}[]{ll}f_{k}&\mathrm{if}\,i=0\\\
e_{i,k}(x_{i})\cdot\sum_{j}q_{i,j,k}\cdot F_{i-1,j}&\mathrm{if}\,i>0\\\
\end{array}\right.,$

where $\\{f_{k}\\}_{1\leq k\leq n}$ denotes some initial conditions. These
conditions may be set to $1/n$, or to the stationary distribution of the
Markov chain (providing it exists). Equation (3) is called the _forward
algorithm_ and allows for the computation of the likelihood of the sequences
as

(4) $\mathcal{L}=\sum_{k}F_{L,k}.$

This recursion is an example of dynamic programming, allowing for the
integration over all possible ARGs very efficiently, as it scales in
$\mathcal{O}(L\cdot n^{2})$. The symmetry relationships in the transition
matrix and the relatively low frequency of the observed heterozygous states,
however, allow for further improvement of the run time of the forward
algorithm [26, 76, 59].

The likelihood function depends on a set of parameters $\Theta$, which
includes the demographic model (effective population size and its variation)
and the recombination rate. More complex models can be implemented, for
instance allowing for population structure. These parameters can affect either
the emission probabilities $e_{i,k}(x)$, the transition probabilities
$q_{i,j,k}$, or both. The emission and transition probabilities define the
type of model used. Importantly, most models assume that the process is
homogeneous along the sequence, so that both $e_{i,k}(x)$ and $q_{i,j,k}$ are
actually independent of $i$. In the following, we will review several examples
of coalescent HMM models.

###  1.5. The two-genome case

Figure 3.  Demographic models and hidden states for CoalHMM with two and three
sequences. A: Pairwise SMC. Hidden states correspond to discretised divergence
times between two sequences. Parameters of the model potentially contain a
species divergence time $\tau$ and several epochs of constant effective
population sizes $\theta_{k}$. The CoalHMM model of Mailund et al. [42] uses
only one epoch and $\theta$, the corresponding states $t_{j}$ being,
therefore, drawn from the exponential distribution with mean $2\cdot\theta$,
shifted by $\tau$. The PSMC model of Li and Durbin [39], assumes a skyline
model of multiple epochs, yet with individuals from the same population
($\tau=0$). B: CoalHMM with three sequences and ILS. The hidden states
correspond to four genealogies that differ both in time and order of the
coalescence events. The model assumes constant but distinct ancestral
effective population sizes $\theta_{12}$ and $\theta_{123}$, as well as the
two species divergence times $\tau_{12}$ and $\tau_{123}$.

When the genome sample consists only of two genomes, the marginal genealogies
have a more simple encoding consisting of a single (continuous) number
representing the time to the most recent common ancestor (TMRCA) of the two
sequences. The TMRCA can be further discretised into $n$ hidden states, each
represented by a mean value $(t_{j})_{1\leq j\leq n}$. The transition
probabilities between the hidden states can be calculated under the SMC.
Variants of this model were developed independently by Li and Durbin [39] and
Mailund et al [42]. In the latter, the two genomes come from two distinct
populations that diverged at a time $\tau$ units ago (Figure 3A). The common
ancestral population is assumed to have a constant effective population size
$\theta_{\text{anc}}$. The TMRCA $(t_{j})_{1\leq j\leq n}$ follows an
exponential distribution shifted by an amount of $\tau$. Mailund and
collaborators applied this model to the newly sequenced genomes of two
Orangutan subspecies in order to estimate their ancestral effective population
size and the time of their last genetic exchange. They further extended this
model to allow for a period of gene flow after the initial separation of the
two populations [43].

In the Li and Durbin model, named pairwise sequentially Markov coalescent
(PSMC), the two genomes come from a single population. This approach was
further improved by Schiffels and Durbin [61], who used a more accurate
recombination model. The authors consider a demographic model where the
effective population size is piecewise constant over a given number of epochs
(Figure 3A). The parameters of the model consist of the set of ancestral
sizes, as well as the recombination rate, presumed to be constant along the
sequences. While the epochs of the demographic model and the discretisation
scheme used for the divergence time are distinct aspects, it is convenient to
have some overlap between the two, providing that there are at least as many
hidden states as epochs (otherwise some parameters would become
unidentifiable). Li and Durbin proposed to consider one hidden state per
epoch, so that each segment is represented by one value of $t_{j}$ and one
value of $\theta_{j}$ (Figure 3A). By estimating one $\theta$ per epoch, the
PSMC model allows the reconstruction of a “skyline” plot where population size
varies in time, from present to the distant past (Figure 4). This method was
applied to data from the 1000 Genomes Project [1] in order to infer the
demographic history of distinct populations, which show the signature of the
out-of-Africa bottleneck.

The two approaches of Li and Durbin and Mailund et al. further differ in their
calculation of the emission probabilities. Focusing on the population level
and a relatively short time span of $\sim 1$ Myr, Li and Durbin consider an
_infinite sites model_ where only one mutation per site can happen. They
further consider all types of mutations as equally probable and ignore the
biochemical nature of the underlying nucleotides. This reduces the number of
observed states to three types: homozygous (the two sequences are identical at
a given position), heterozygous (the two sequences differ at a given
position), and unknown (at least one sequence has an unresolved state at that
position). The emission probabilities then take the simple form:

(5) $\displaystyle e_{i,j,\mbox{homozygous}}$ $\displaystyle=$
$\displaystyle\exp(-\theta\cdot t_{j})$ (6) $\displaystyle
e_{i,j,\mbox{heterozygous}}$ $\displaystyle=$ $\displaystyle
1-\exp(-\theta\cdot t_{j})$ (7) $\displaystyle e_{i,j,\mbox{unknown}}$
$\displaystyle=$ $\displaystyle 1,$

where $\theta=4\cdot N_{e}\cdot u$ denote the population mutation rate, and
$u$ the per nucleotide, per generation molecular mutation rate.

Comparing genomes from two distinct (sub)species, therefore potentially
encompassing larger time scales, Mailund et al. used a fully parametrised
substitution model as used in parametric phylogenetic reconstruction methods
[22]. The mutation process is then a continuous time, discrete-state Markov
model with generator $U\\!$, and the emission probabilities are given by
$\exp(U\cdot t_{j})$ for each hidden state $j$. In both models, the emission
probabilities only depend on the observed states and are independent of the
actual position in the sequence, assuming a homogeneous mutation process along
the genome.

Figure 4.  Demographic inference under the PSMC. The MSMC2 software was used
independently on 20 diploid individuals from the 1000 genomes project [1], 10
from the CEU population (Utah residents with European ancestry) and 10 from
the LWK population (Luhya in Webuye, Kenya). MSMC2 was run on data from
chromosome 9 only, with default parameters. The results show that individuals
with European ancestry underwent a stronger bottleneck between 50 and 100 ky
ago, corresponding to the out-of-Africa event.

###  1.6. The three-genome case

In 2007, Hobolth et al. introduced the first CoalHMM model [27], by modelling
the possible genealogical relationships between three species: “$(A,B),C$”,
“$A,(B,C)$” and “$(A,C),B$” (Figure 3B). Considering the two speciation events
that separate first the ancestor of $A$ and $B$ from the ancestor of $C$, and
then the ancestor of $A$ from the ancestor of $B$, the probability that an
individual sequence from $A$ coalesces with a sequence from $B$ within the
$AB$ ancestral species depends on the ancestral population size and the time
between the two speciation events [19]. Backwards in time, if the two
corresponding lineages do not coalesce until the most ancient speciation time,
any of them can coalesce with a sequence from species $C$ before coalescing
with each other. This phenomenon, which results in the marginal genealogy
being incongruent with the phylogeny, is termed _incomplete lineage sorting
(ILS)_. Using coalescent theory, Hobolth et al. derived relationships between
the transition probabilities and used them to infer ancestral effective
population sizes as well as the dates of species divergence, the so-called
speciation times. In this first model the hidden states differ in tree
topology and divergence times. The first hidden state corresponds to the case
where the two lineages A and B coalesce within the ancestral population of A
and B, leading to a genealogy congruent with the phylogeny. The three other
topologies denote cases where the A and B lineages did not coalesce within the
AB ancestor, so that the three lineages A, B and C were already present within
the ancestral population of the three species. These three states correspond
to the cases where A and B, A and C, or B and C coalesce first, respectively
(Figure 3B). The model assumes constant ancestral population sizes and the
divergence times for each topology are reduced to the averages of the
corresponding exponential distributions. The proportion of ILS topologies
directly depends on the time separating the speciation events $\Delta_{T}$ and
the effective size of the ancestral population $\theta_{\text{anc}}$ [19]:

(8)
$\operatorname{Pr}(ILS)=\frac{2}{3}\exp\left(-2\cdot\frac{\theta_{\text{anc}}}{\Delta_{T}}\right),$

allowing the estimations of these parameters from the patterns of topology
variation. The three-species CoalHMM model was applied to genome sequences of
Great Apes: Orangutan [28], Gorilla [60], Bonobo [56], Baboons [58], in order
to infer the patterns of ILS and the ancestral effective population sizes in
this group of species (reviewed in [44]). It was also applied to species of
fungal pathogens where it was used to infer recombination rates [73].

###  1.7. The multiple-genome cases

With the development of sequencing technologies and the increasing sample size
of population genomic datasets, models able to extract the genealogical
information contained in multiple genomes are needed. Building CoalHMM models
with more than two or three sequences poses, however, a computational
challenge because of the underlying combinatorics of marginal genealogies. The
number of possible topologies increases hyper-exponentially with the number
$M$ of sampled sequences, and there is an infinite number of possible
genealogies with a given topology due to the continuous nature of branch
lengths (divergence times). Further approximations are, therefore, required to
scale CoalHMM approaches with larger datasets.

####  1.7.1. Using conditional sampling distributions

In a series of articles [54, 70, 62, 69], Song and collaborators developed an
approach based on the so-called _conditional sampling distribution_ (CSD)
introduced by Li and Stephens [40]. This approach stems from the chain rule of
conditional probabilities, allowing the expression of the likelihood of a
sample of $M$ sequences $S_{1},\ldots,S_{M}$ as a product of conditional
likelihoods:

(9)
$\operatorname{Pr}(S_{1},S_{2},\ldots,S_{M}\mid\Theta)=\operatorname{Pr}(S_{1}\mid
S_{2},\ldots,S_{M},\Theta)\cdot\operatorname{Pr}(S_{2},\ldots,S_{M}\mid\Theta)\\\
=\operatorname{Pr}(S_{1}\mid
S_{2},\ldots,S_{M},\Theta)\cdot\operatorname{Pr}(S_{2}\mid
S_{3},\ldots,S_{M}\mid\Theta)\cdot\operatorname{Pr}(S_{3},\ldots,S_{M}\mid\Theta)\\\
=\operatorname{Pr}(S_{M}\mid\Theta)\prod_{k=1}^{M-1}\operatorname{Pr}(S_{k}\mid
S_{k+1},\ldots,S_{M},\Theta),$

where $\Theta$ denotes the parameter vector. The conditional likelihoods,
however, are approximated, so that the resulting likelihood is a product of
approximate conditionals (PAC) [40], which depends on the order by which the
sequences are treated in the product chain. This is usually accommodated by
permutations and averaging [40], or by a composite likelihood approach such as
the “leave-one-out” strategy [62]:

(10)
$\operatorname{Pr}(S_{1},S_{2},\ldots,S_{M}\mid\Theta)\simeq\prod_{i=1}^{M}\operatorname{Pr}\big{(}S_{i}\mid\left\\{S_{j}\right\\}_{j\neq
i},\Theta\big{)}.$

The CSD are computed under an SMC model, given a piecewise constant
demographic model, as in the PSMC. The model was further extended to allow
more complex demographic scenarios with population structure and migration
[69].

####  1.7.2. Modelling the most recent coalescence events

Schiffels and Durbin [61] developed the _multiple sequentially Markov
coalescent (MSMC)_ , which models only the most recent coalescent event in the
sample. The underlying rationale was that the PSMC is lacking resolution in
the more recent past, due to the very small number of mutations and
coalescences happening in the most recent epochs. Combining multiple samples,
therefore, has the potential to compensate for the lack of information in a
single pair of genomes. The MSMC approach is elegant as it reduces the
combinatorics of the hidden states to one continuous variable (which is
discretised, as in the PSMC): the time of coalescence, together with the index
of the two genomes in the sample that are coalescing, bringing the number of
hidden states to $n={M\choose 2}\cdot K$, where $K$ is the number of discrete
classes used for the distribution of divergence times. The efficiency of the
MSMC approach is, however, paradoxical: by gaining resolution in the present
as the sample size increases, the method progressively loses power as the
number of modelled sequences becomes larger (see [17] for an illustration). In
practice, the authors showed that for a human dataset, the maximum resolution
is obtained for eight haploid genomes [61].

####  1.7.3. Using a composite likelihood

In the MSMC2 approach [45, 78], Schiffels and collaborators proposed to
approximate the likelihood of a sample of $M$ genomes by independently
considering all pairs of genomes. The likelihood of the sample is then
approximated by the product of all pairwise likelihoods, each computed under
the PSMC model. While the pairwise likelihoods are exact under the SMC, the
likelihood of the sample is an example of _composite likelihood_ [36]:

(11)
$\operatorname{Pr}(S_{1},S_{2},\ldots,S_{M}\mid\Theta)\simeq\prod_{i=1}^{M-1}\prod_{j=i+1}^{M}\operatorname{Pr}(S_{i},S_{j}\mid\Theta).$

The likelihood here is an approximation since the divergence times between
pairs of sequences in a genealogy are not independent. The MSMC2 approach is
therefore better described as a “multiple pairwise SMC”. It was shown to
display good resolution in both the past and present time, efficiently making
use of the increasing quantity of signal as the sample size increases.

####  1.7.4. Augmenting the PSMC with site frequency spectra

While approaches like diCal and MSMC2 allow for the efficient modelling of the
evolution of multiple genomes, they are intrinsically limited as the
computational cost become prohibitive for samples of more than one or two
dozen genomes (at least for genomes with a size of the order of that of
humans). Terhorst and colleagues introduced a hybrid approach combining the
power of the SMC, which makes efficient use of linkage patterns, with that of
classical site frequency spectrum (SFS) based approaches [76]. This modelling
framework, termed SMC++, considers a “focus” diploid individual that is
modelled with a PSMC approach. The observed states are then augmented by
taking into account additional genomes to compute an SFS. The emission
probabilities are calculated as the probability of observing the local SFS
given the genealogy at the focus individual, and the authors proposed an
approach to compute such a conditional site frequency spectrum (CSFS). The
resulting SMC++ model can accommodate hundreds of individual genomes. Another
innovation introduced in this approach is the abandonment of the “skyline”
model of piecewise constant effective population size in favour of a spline
model. While divergence times are still discretised, the corresponding times
for each category are derived from a spline curve whose parameters are
estimated. This reduces the number of parameters to estimate and ensures
smoother inferred demographies. The SMC++ approach has been applied to human
data as well as other species, including Drosophila and Zebra finch [76].

##  2\. Heterogeneity of processes along the genome

In all models that we evoked so far, evolutionary processes have been
considered to be homogeneous along sequences. In this section, I review
evidence that these assumptions are at odds with our current knowledge of the
biology of genomes.

###  2.1. Variation of the recombination rate

Recombination rates can vary extensively between species [68], between sexes
[37] and within genomes. At the molecular scale, multiple levels can be
distinguished: recombination rate correlates negatively with chromosome size,
a pattern attributed to the mechanism of meiosis and crossing-over
interference [32]. Given that the rate of crossing-over events is low, this
leads to a higher recombination rate in small chromosomes. Recombination rates
vary also within chromosomes: in Primates, it is generally higher at the start
and end of the chromosomes (the so-called telomeric regions) [72], while in
Drosophila the opposite pattern is observed [13, 10]. In many species
recombination events have an increased chance to occur in particular regions,
called hotspots [55, 52, 74] (but see [77] for a counterexample).

The variation of recombination rate has two types of consequences on the
patterns of sequence diversity. Because the molecular mechanisms of
recombination are tightly linked to DNA repair, recombination itself can be
mutagenic and locally increase sequence variability [38, 4]. Furthermore, in
many species, the repair mechanisms involve gene conversion between homologous
sequences, as one chromosome is used as a template to repair the other one.
However, this mechanism is biased in many species: in the case of heterozygous
positions, the “C” or “G” nucleotides are preferred over “A” and “T”
nucleotides, potentially resulting in large scale variations of GC content
[16] mirroring the variations in recombination rate. The recombination rate
also has indirect effects on genetic diversity: because it breaks down genetic
linkage, recombination counteracts the reduction of diversity at sites linked
to loci under selection, both negative (background selection [11]) and
positive (genetic hitch-hiking [9]). By modulating the local effective
population size, variation of recombination rate along the genome has a strong
impact on the underlying genealogy.

###  2.2. Variation of the mutation rate

Finally, the rate at which mutations occur can vary extensively along the
genome [7]. Mutations can occur via direct modification of the DNA or
indirectly, via errors in the replication or repair mechanisms. Particular
sequence motifs, such as CpG dinucleotides are known to undergo comparatively
higher mutation rates, via the methylation of the cytosine, which is then
deaminated into a thymine, leading to a TpG dinucleotide. In addition to the
potentially mutagenic effect of recombination, which also plays a role in the
repair of DNA damage, the replication machinery itself is error-prone. This
error rate is position dependent: it increases with the replication time,
being lower close to the origins of replication [67, 3, 79]. Under a neutral
scenario, the mutation process is independent of the coalescent process [30],
and, therefore, has no impact on the underlying genealogies. Yet, for
inference models, mutation rate variation acts as a confounding factor, as a
high sequence divergence can be either explained by an ancient coalescence
time or a high mutation rate. In CoalHMM models, the mutation rate will have
an impact on the emission probabilities, that is, the probability of observing
the observed sequence diversity given a genealogy.

##  3\. Existing approaches to account for spatial heterogeneity

In this section, I review the approaches that have been developed to cope with
the heterogeneity of evolutionary processes acting along the genome.

###  3.1. Inferring sequential heterogeneity alone

A large body of work is built on the idea that, if a parameter affects certain
patterns of genetic diversity, it should be possible to use these patterns to
recover the underlying variation of the parameter. The most studied case in
this respect is the recombination rate, via its impact on linkage
disequilibrium. Given an a priori known demographic scenario, it is possible
to compute the likelihood of the data for any given recombination rate, and
use it to estimate the most likely recombination rate value. Due to the
complexity of the underlying model, however, approximations are required to
apply these methods to large genomic datasets. McVean, Awadalla and Fearnhead
[48] introduced the use of a composite likelihood, approximating the full
likelihood by the product of the likelihoods of all pairs of positions within
a minimum distance of each other. This approach is the basis of several
popular methods for recombination rate inference such as LDhat [5] and
LDhelmet [10]. Further developments of these models allowed for the
incorporation of variable population sizes [33, 64]. The underlying
demography, however, has to be estimated independently from the data. Li and
Stephens [40], on the other hand, used the conditional sampling distribution
and the product of approximate conditionals (PAC) to approximate the
likelihood. An application of this method also includes the reconstruction of
haplotypes from genotypic data, a problem known as _phasing_ [71].

###  3.2. Inference using piecewise-homogeneous processes

The most simple approach to infer heterogeneous processes along the genome
while jointly accounting for demography is to use a window-based approach,
consisting of dividing the genome into segments of fixed sizes and estimating
model parameters independently in each resulting window. This strategy was
used by Stukenbrock et al. [73] to use the patterns of ILS and a CoalHMM model
to estimate the recombination rate in 100 kb windows along the genome of the
fungal pathogen Zymoseptoria tritici. In most cases, however, SMC models
require long genome sequences to be able to confidently estimate parameters,
and cannot be run in windows of small sizes, in particular for models at the
population level. Furthermore, window-based approaches raise the issue of the
window size and boundaries to use.

###  3.3. Using sequentially heterogeneous simulation procedures

While the computation of the likelihood of the data under a sequentially
heterogenous process is notoriously difficult, simulating under the
corresponding model can be comparatively easy. Software like the Markovian
coalescent simulator (MaCS) [12], the sequential coalescent with recombination
model (SCRM) [66], fastsimcoal [21] and MSprime [34] allow the simulation of
population genomic data sets under models with variable recombination rate.
Owing to their high computational efficiency, they can be used within an
approximate Bayesian computation (ABC) framework in order to estimate
demographic parameters under realistic recombination models [80]. This
possibility, however, has to date been underexploited, as demographic
inference is so far conducted with data simulated under a homogeneous
recombination landscape (see for instance [41]). While no ABC method has been
developed with the goal to infer the variation of population genomic
parameters along the genome, Gao et al. [24] introduced a machine learning
approach to infer recombination rates. The underlying simulations, however,
are conducted under a model with constant recombination rate.

###  3.4. A posteriori inference of heterogeneous processes

The HMM methodology allows, via the forward algorithm, to compute the
likelihood of the data given a specified demographic model by efficiently
integrating over the unknown underlying ARG. The HMM toolbox further allows
for the computation of the a posteriori probability of each marginal genealogy
for each position [17]:

(12) $\operatorname{Pr}(\boldsymbol{\mathcal{H}}_{i}=H_{j}\mid
x_{1},\ldots,x_{L})=\frac{\operatorname{Pr}(x_{1},\ldots,x_{L},\boldsymbol{\mathcal{H}}_{i}=H_{j})}{\operatorname{Pr}(x_{1},\ldots,x_{L})}.$

The denominator of the ratio is the likelihood of the data, $\mathcal{L}$. In
order to compute the numerator, we need to introduce the so-called _backward
algorithm_ [15]:

(13)
$B_{i,j}=\operatorname{Pr}(x_{i+1},\ldots,x_{L}\mid\boldsymbol{\mathcal{H}}_{i}=H_{j})=\left\\{\begin{array}[]{ll}1&\mathrm{if}\,i=L\\\
\sum_{k}e_{i+1,k}(x_{i+1})\cdot q_{i+1,j,k}\cdot
B_{i+1,k}&\mathrm{if}\,i<L,\\\ \end{array}\right.$

The posterior probability of hidden state $H_{j}$ can therefore be computed as

(14) $\operatorname{Pr}(\boldsymbol{\mathcal{H}}_{i}=H_{j}\mid
x_{1},\ldots,x_{L})=\frac{F_{i,j}\cdot B_{i,j}}{\mathcal{L}}.$

This formula allows for the reconstruction of the most probable marginal
genealogy at each position $i$ by taking the maximum posterior probability

(15)
$\widehat{\boldsymbol{\mathcal{H}}}_{i}=\underset{j}{\operatorname{arg\,max}}\left(\operatorname{Pr}(\boldsymbol{\mathcal{H}}_{i}=H_{j}\mid
x_{1},\ldots,x_{L})\right),$

a procedure called _posterior decoding_. The posterior decoding is performed
after fitting the model parameters by maximizing the likelihood. It is
therefore an example of _empirical Bayesian_ inference [46].

Posterior probabilities of marginal genealogies can also be used to obtain
posterior estimates of biological quantities of interest, accounting for the
uncertainty on the underlying genealogy. The posterior mean estimate
$\widehat{\lambda}_{i}$ at position $i$ of a property $\Lambda(H_{j})$ can be
obtained by

(16)
$\widehat{\lambda}_{i}=\sum_{j}\operatorname{Pr}(\boldsymbol{\mathcal{H}}_{i}=H_{j}\mid
x_{i},\ldots,x_{L})\cdot\Lambda(H_{j}).$

If $\Lambda$ is the coalescence time between two sequences, this formula can
be used to get posterior estimates of sequence divergence along the genome
[53]. More complex examples of functions include, for instance, whether
$\boldsymbol{\mathcal{H}}_{i}$ is distinct from
$\boldsymbol{\mathcal{H}}_{i-1}$, or, in other words, whether a recombination
event occurred between positions $i$ and $i-1$. This allows for the
reconstruction of a recombination map, integrating over the ARG. Such approach
was notably used by Munch et al [51] to reconstruct the recombination map of
the human-chimpanzee ancestor. Posterior estimates are rather robust to the
specified input model and can therefore offer a powerful approach to infer
aspects of the process that are not directly accounted for by the model.
However, because some model properties are intrinsically confounded, such as
local divergence and mutation rate, ignoring spatial heterogeneity might
result in biased inference [8].

##  4\. The integrative sequentially Markov coalescent

In order to account for heterogeneous processes along the genome, we recently
developed the integrative sequentially Markov coalescent (iSMC) [8], an
extension of the SMC. In this framework, parameters of the original SMC vary
along the genome in a Markovian manner, allowing for the modelling of genome
heterogeneity in addition to demographic processes. I illustrate this approach
with results from a recent application of this framework to infer
recombination landscapes and further discuss possible extensions.

###  4.1. The Markov-modulated sequentially Markov coalescent

Whilst the framework can be applied to cases where more than one parameter
varies along the genome, for simplicity, we here consider the case where one
parameter only varies, which we label $\boldsymbol{\mathcal{R}}$,
$\boldsymbol{\mathcal{R}}_{i}$ denoting the values of
$\boldsymbol{\mathcal{R}}$ at position $i$ in the sequences. We assume that
$\boldsymbol{\mathcal{R}}$ follows an a priori known discrete distribution
with $n^{R}$ categories, each with mean value $R_{k}$, with $1\leq k\leq
n^{R}$. In the iSMC framework, the transition and/or emission probabilities
are functions of $\boldsymbol{\mathcal{R}}$ and are, therefore, noted as
$e^{\text{SMC}}_{i,k}(x,R)$ and $q^{\text{SMC}}_{i,j,k}(R)$, respectively. The
key assumption is then to consider that the variation of
$\boldsymbol{\mathcal{R}}$ along the genome can be modelled as a Markov model,
that is, there is a matrix of probabilities $q^{R}$ defined as

(17)
$q^{R}_{i,j,k}=\operatorname{Pr}(\boldsymbol{\mathcal{R}}_{i}=R_{k}\mid\boldsymbol{\mathcal{R}}_{i-1}=R_{j}).$

The forward recursion of the CoalHMM can then be written as

(18)
$F^{\mathrm{SMC}}_{i,j,k}=\operatorname{Pr}(x_{1},\ldots,x_{i},H_{j},R_{k})=\\\
\left\\{\begin{array}[]{ll}f_{j}\cdot
f^{R}_{k}&\mathrm{if}\,i=0\\\\[5.69054pt]
e^{\mathrm{SMC}}_{i,j}(x_{i},R_{k})\cdot\sum_{u}\sum_{v}q^{\mathrm{SMC}}_{i,u,j}(R_{k})\cdot
q^{R}_{i,v,k}\cdot F^{\mathrm{SMC}}_{i-1,u,v}&\mathrm{if}\,i>0,\\\
\end{array}\right.$

where $f^{R}_{k}$ denotes the a priori probability
$\operatorname{Pr}(\boldsymbol{\mathcal{R}}=R_{k})$. Because the SMC now
depends on a parameter that itself follows a Markov process, the resulting
process can be described as a Markov-modulated Markov chain (MMMC). As an MMMC
is itself a Markov process [23], we can rewrite the forward recursion as:

(19)
$F^{\text{iSMC}}_{i,k}=\operatorname{Pr}(x_{1},\ldots,x_{i},H^{\text{iSMC}}_{k})=\left\\{\begin{array}[]{ll}f^{\text{iSMC}}_{k}&\mathrm{if}\,i=0\\\\[5.69054pt]
e^{\text{iSMC}}_{i,k}(x_{i})\cdot\sum_{j}q^{\text{iSMC}}_{i,j,k}\cdot
F^{\text{iSMC}}_{i-1,j}&\mathrm{if}\,i>0.\\\ \end{array}\right.$

In the iSMC hidden Markov model, the hidden states $H^{\text{iSMC}}$ consist
of all possible pairs of genealogies and $\boldsymbol{\mathcal{R}}$:
$(R_{a},H_{b})$, with $1\leq a\leq n^{R}$ and $1\leq b\leq n$. The emission
probabilities $e^{\text{iSMC}}_{i,j}(x)$ now depend on the pair $(R,H)_{j}$,
and the initial probabilities of the hidden states are
$f^{\text{iSMC}}=f^{R}\otimes f$. Similarly, the transition probabilities of
the Markov-modulated SMC (MMSMC) can be written as a function of the
transition probabilities of the two Markov chains:

(20) $q_{i}^{\text{iSMC}}=\left(\begin{array}[]{ccc}q^{R}_{i,1,1}\cdot
q_{i}^{\text{SMC}}(R_{1})&\cdots&q^{R}_{i,1,n^{R}}\cdot
q_{i}^{\text{SMC}}(R_{1})\\\ \vdots&\ddots&\vdots\\\ q^{R}_{i,n^{R},1}\cdot
q_{i}^{\text{SMC}}(R_{(n^{R})})&\cdots&q^{R}_{i,n^{R},n^{R}}\cdot
q_{i}^{\text{SMC}}(R_{(n^{R})})\end{array}\right).$

(As we consider the process modelling the variation of
$\boldsymbol{\mathcal{R}}$ to be itself homogeneous along the genome, we have
$\forall
i_{1},i_{2},\;q^{\text{iSMC}}_{i_{1},j,k}=q^{\text{iSMC}}_{i_{2},j,k}=q^{\text{iSMC}}_{j,k}$.)
The iSMC model can therefore be analysed with standard HMM methodology, just
like the homogeneous SMC. The number of hidden states, however, is now
$n^{R}\cdot n$, meaning that the complexity of the likelihood calculation
becomes $\mathcal{O}\big{(}L\cdot(n\cdot n^{R})^{2}\big{)}$.

The iSMC model adds relatively few extra parameters to the SMC: the transition
probabilities $q^{R}$, which can be reduced to one parameter (see below), and
parameters of the a priori distribution of $\boldsymbol{\mathcal{R}}$. While
parameters of the distribution of $\boldsymbol{\mathcal{R}}$ are generally not
of direct biological interest, the posterior decoding of the HMM allows for
the inference of the underlying landscape of the heterogeneous parameter.
Distinct decoding procedures can be performed in the case of Markov-modulated
HMMs:

  1. (1)

A full decoding, where the most likely pair $(R,H)$ at each position is
reconstructed:

(21)
$\widehat{\left(\boldsymbol{\mathcal{R}},\boldsymbol{\mathcal{H}}\right)}_{i}=\underset{j,k}{\operatorname{arg\,max}}\left(\operatorname{Pr}(\boldsymbol{\mathcal{H}}_{i}=H_{j},\boldsymbol{\mathcal{R}}_{i}=R_{k}\mid
x_{i},\ldots,x_{L})\right),$

  2. (2)

A partial decoding of genealogies, where the most likely genealogy is
inferred, summing over all heterogeneous parameters:

(22)
$\hat{\boldsymbol{\mathcal{H}}_{i}}=\underset{j}{\operatorname{arg\,max}}\Big{(}\sum_{k}\operatorname{Pr}(\boldsymbol{\mathcal{H}}_{i}=H_{j},\boldsymbol{\mathcal{R}}_{i}=R_{k}\mid
x_{i},\ldots,x_{L})\Big{)},$

  3. (3)

A posterior mean estimation of the heterogeneous variable. Setting
$\Lambda(H_{j},R_{j})=R_{j}$ and applying equation 16, we get:

(23)
$\hat{\boldsymbol{\mathcal{R}}_{i}}=\sum_{j}\sum_{k}\operatorname{Pr}(\boldsymbol{\mathcal{H}}_{i}=H_{j},\boldsymbol{\mathcal{R}}_{i}=R_{k}\mid
x_{i},\ldots,x_{L})\cdot R_{k}\\\
=\sum_{k}R_{k}\cdot\sum_{j}\operatorname{Pr}(\boldsymbol{\mathcal{H}}_{i}=H_{j},\boldsymbol{\mathcal{R}}_{i}=R_{k}\mid
x_{i},\ldots,x_{L}).\\\ $

The partial decoding of genealogies enables the reconstruction of the ARG
while accounting for the heterogeneity of the SMC along the genome. The
posterior estimates of the heterogeneous variable allows the reconstruction of
the variation along the genome while accounting for the genealogy and its
uncertainty. In the next section, we apply this framework to model the
variation of the recombination rate along the genome.

###  4.2. A case study: inference of recombination rate variation

Recombination is the best documented heterogeneous process along the genome.
It can be measured experimentally or indirectly using genomic approaches
(reviewed in [5]). In the context of the sequential coalescent and the SMC
approximation, the local recombination rate affects the probability of
transition from one genealogy to another, which increases with higher
recombination rates. To model variable recombination rates in iSMC, we
considered that the local population recombination rate $\rho=4\cdot
N_{e}\cdot r$, where $r$ is the molecular recombination rate in cM / bp per
generation, is the product of a genome average $\rho_{G}$ and a local modifier
$r^{\rho}$. This modifier follows a prior discrete distribution of mean 1 and
with $n^{\rho}$ categories, for instance a discretised Gamma distribution with
shape parameter $\alpha$. We further considered a simple model for the
transition probabilities between the $r^{\rho}$ classes, assuming equal
probabilities of change, $\gamma$. The transition probability matrix
$q^{\rho}$ takes the form:

(24)
$q^{\rho}=\left(\begin{array}[]{cccc}1-\gamma&\gamma/(n^{\rho}-1)&\cdots&\gamma/(n^{\rho}-1)\\\
\gamma/(n^{\rho}-1)&1-\gamma&\ddots&\vdots\\\
\vdots&\ddots&\ddots&\gamma/(n^{\rho}-1)\\\
\gamma/(n^{\rho}-1)&\cdots&\gamma/(n^{\rho}-1)&1-\gamma\\\
\end{array}\right).$

Using the forward equation 19, we can compute the likelihood of the parameters
of the iSMC model. Using optimisation procedures, estimates of the $\alpha$
and $\gamma$ parameters together with the average $\rho_{G}$ and the
demography parameters can be obtained by maximizing the likelihood function.
The likelihood calculation can also be used to perform model comparisons and
test for the heterogeneity of the coalescent process along the genome, for
instance using Akaike’s information criterion (AIC).

Figure 5.  Posterior estimates of recombination rates from a single diploid
genome using iSMC. A 30 Mb region was simulated using the SCRM program [66]
and a variable recombination rate, and then inferred with iSMC, as described
in [8]. Recombination rates were averaged in windows of 200 kb.

A posterior decoding approach (equation 23) can then be used to obtain
estimates of site-specific recombination rates. Simulations under controlled
recombination landscape and demography can be used to assess the accuracy of
the iSMC inference [8]. Figure 5 shows that iSMC recovers the underlying
recombination landscape with good accuracy, despite generally underestimating
high recombination rates. A possible explanation for this is the
discretisation procedure, as the posterior mean estimate is bounded by the
class with the highest mean recombination value. Allowing for more
recombination classes allows for a wider range of values and can potentially
reduce this bias, at the cost of increasing the running time and memory usage.
Because it can recover the recombination landscape from a single diploid only,
the iSMC model was used to generate recombination maps from extinct hominids
from their ancient genome sequences [8].

###  4.3. Extension of iSMC: multiple genomes and multiple heterogeneous
parameters

The Markov-modulated Markov model framework underlying the iSMC approach can
be applied to other SMC models, such as the MSMC [61]. Hidden states of the
resulting CoalHMM are combinations of TMRCA, pairs of genome indices
undergoing the most recent coalescent event, and classes of heterogenous
parameters such as the recombination rate. Extension to multiple genomes can
also be achieved using a composite likelihood approach, as implementated in
MSMC2 (see 1.7.3). The likelihood of the dataset is then approximated by the
product of the likelihoods of each pair of genomes, which are modelled
separately with their own process. In the case of the iSMC approach, this
implies considering that the heterogeneous parameters vary independently along
each pair of genomes; for a model with variable recombination rate, this is
equivalent to estimating a distinct recombination map for each pair of
genomes. This assumption is clearly incorrect for the vast majority of
positions within genomes from the same population. Extensions enforcing a
common map while allowing the coalescent processes to be independent will,
therefore, be instrumental in efficiently scaling the iSMC approach to larger
sample sizes.

The iSMC framework further allows the joint modelling of the variation of
multiple parameters along the genome, such as the mutation and recombination
rates. In such approach, the hidden states are a combination of TMRCA and
classes for each discretised parameter. The addition of any additional
heterogeneous parameter multiplies the complexity of the likelihood
calculation by ${n^{\phi}}^{2}$, where $n^{\phi}$ is the number of discrete
classes considered for the parameter distribution. Besides the increased
complexity, identifiability issues may also arise, since local patterns of
diversity may be equally explained by variation of demography, recombination
rate or mutation rate.

However, when these parameters vary at a scale larger than the variation of
the TMRCA along the genome and when very large genome sequences are analyzed,
increasingly complex models may be successfully fitted.

##  5\. Conclusions

The availability of complete genome data opened the floodgates for the
detailed inference of the demographic history of species. A new generation of
coalescent-based models permits the extraction of demographic signal from the
patterns of genetic linkage along sequences. Such models, however, largely
ignore fundamental aspects of genome biology, that is, that processes such as
recombination and mutation are highly heterogeneous along genomes. Extending
these approaches to account for such heterogeneity not only potentially
improves demographic inference, but also allows to reconstruct the underlying
genomic landscape.

### Acknowledgement

I would like to thank Gustavo V. Barroso for critical reading of this
manuscript and for providing the simulation data plotted in Figure 5. I am
deeply indebted to Ellen Baake, as well as Jeffrey P. Spence and an anonymous
reviewer for their careful reading of the manuscript, for finding several
mistakes and typos, and for their suggestions on how to improve its clarity.

## References

  * [1] 1000 Genomes Project Consortium (and 8 coauthors), A map of human genome variation from population-scale sequencing, Nature 467 (2010), 1061–1073.
  * [2] A.M. Adams and R.R. Hudson, Maximum-likelihood estimation of demographic parameters using the frequency spectrum of unlinked single-nucleotide polymorphisms, Genetics 168 (2004), 1699–1712.
  * [3] N. Agier and G. Fischer, The mutational profile of the yeast genome is shaped by replication, Mol. Biol. Evol. 29 (2012), 905–913.
  * [4] I. Alves, A.A. Houle, J.G. Hussin, and P. Awadalla, The impact of recombination on human mutation load and disease, Philos. Trans. R. Soc. Lond. B, Biol. Sci. 372 (2017), 20160465.
  * [5] A. Auton and G. McVean, Estimating recombination rates from genetic variation in humans, Methods Mol. Biol. 856 (2012), 217–237.
  * [6] E. Baake and M. Baake, Ancestral lines under recombination, this volume.
  * [7] C.F. Baer, M.M. Miyamoto, and D.R. Denver, Mutation rate variation in multicellular eukaryotes: causes and consequences, Nature Reviews Genetics 8 (2007), 619–631.
  * [8] G.V. Barroso, N. Puzovic, and J. Dutheil, Inference of recombination maps from a single pair of genomes and its application to archaic samples, PLoS Genet. 15 (2019), e1008449.
  * [9] N.H. Barton, Genetic hitchhiking, Philos. Trans. R. Soc. Lond. B, Biol. Sci. 355 (2000), 1553–1562.
  * [10] A.H. Chan, P.A. Jenkins, and Y.S. Song, Genome-wide fine-scale recombination rate variation in Drosophila melanogaster, PLoS Genet. 8 (2012), e1003090.
  * [11] B. Charlesworth, M.T. Morgan, and D. Charlesworth. The effect of deleterious mutations on neutral molecular variation, Genetics 134 (1993), 1289–1303.
  * [12] G.K. Chen, P. Marjoram, and J.D. Wall, Fast and flexible simulation of DNA sequence data, Genome Res. 19 (2009), 136–142.
  * [13] J.M. Comeron, R. Ratnappan, and S. Bailin, The many landscapes of recombination in Drosophila melanogaster, PLoS Genet. 8 (2012), e1002905.
  * [14] A.J. Drummond, A. Rambaut, B. Shapiro, and O.G. Pybus, Bayesian coalescent inference of past population dynamics from molecular sequences, Mol. Biol. Evol. 22 (2005), 1185–1192.
  * [15] R. Durbin, S.R. Eddy, A. Krogh, and G. Mitchison, Biological Sequence Analysis: Probabilistic Models of Proteins and Nucleic Acids, Cambridge University Press, Cambridge, 1998.
  * [16] L. Duret and N. Galtier, Biased gene conversion and the evolution of mammalian genomic landscapes, Annu. Rev. Genomics Hum. Genet. 10 (2009), 285–311.
  * [17] J.Y. Dutheil, Hidden Markov models in population genomics, Methods Mol. Biol. 1552 (2017), 149–164.
  * [18] J.Y. Dutheil, G. Ganapathy, A. Hobolth, T. Mailund, M.K. Uyenoyama, and M.H. Schierup, Ancestral population genomics: the coalescent hidden Markov model approach, Genetics 183 (2009), 259–274.
  * [19] J.Y. Dutheil and A. Hobolth,  Ancestral population genomics, Methods Mol. Biol. 856 (2012), 293–313.
  * [20] S.N. Ethier and R.C. Griffiths, On the two-locus sampling distribution, J. Math. Biol. 29 (1990), 131–159.
  * [21] L. Excoffier and M. Foll, fastsimcoal: a continuous-time coalescent simulator of genomic diversity under arbitrarily complex evolutionary scenarios, Bioinformatics 27 (2011), 1332–1334.
  * [22] J. Felsenstein, Inferring Phylogenies, 2nd ed., Sinauer, Sunderland, MA, 2003.
  * [23] N. Galtier and A. Jean-Marie, Markov-modulated Markov chains and the covarion process of molecular evolution, J. Comput. Biol. 11 (2004), 727–733.
  * [24] F. Gao, C. Ming, W. Hu, and H. Li, New software for the fast estimation of population recombination rates (FastEPRR) in the genomic era, G3 6 (2016), 1563–1571.
  * [25] R.C. Griffiths and P. Marjoram, Ancestral inference from samples of DNA sequences with recombination, J. Comput. Biol. 3 (1996), 479–502.
  * [26] K. Harris, S. Sheehan, J.A. Kamm, and Y.S. Song, Decoding coalescent hidden Markov models in linear time. In: Research in Computational Molecular Biology. RECOMB2014, R. Sharan (ed.), Springer, Cham, 2014, pp. 100-114.
  * [27] A. Hobolth, O.F. Christensen, T. Mailund, and M.H. Schierup, Genomic relationships and speciation times of human, chimpanzee, and gorilla inferred from a coalescent hidden Markov model, PLoS Genet. 3 (2007), e7.
  * [28] A. Hobolth, J.Y. Dutheil, J. Hawks, M.H. Schierup, and T. Mailund, Incomplete lineage sorting patterns among human, chimpanzee, and orangutan suggest recent orangutan speciation and widespread selection, Genome Res. 21 (2011), 349–356.
  * [29] A. Hobolth and J.L. Jensen, Markovian approximation to the finite loci coalescent with recombination along multiple sequences, Theor. Popul. Biol. 98 (2014), 48–58.
  * [30] R.R. Hudson, Properties of a neutral allele model with intragenic recombination, Theor. Popul. Biol. 23 (1983), 183–201.
  * [31] R.R. Hudson, Generating samples under a Wright–Fisher neutral model of genetic variation, Bioinformatics 18 (2002), 337–338.
  * [32] D.B. Kaback, Chromosome-size dependent control of meiotic recombination in humans, Nat. Genet. 13 (1996), 20–21.
  * [33] J.A. Kamm, J.P. Spence, J. Chan, and Y.S. Song, Two-locus likelihoods under variable population size and fine-scale recombination rate estimation, Genetics 203 (2016), 1381–1399.
  * [34] J. Kelleher, A.M. Etheridge, and G. McVean, Efficient coalescent simulation and genealogical analysis for large sample sizes, PLoS Comput. Biol. 12 (2016), e1004842.
  * [35] J.F.C. Kingman, The coalescent, Stochastic Processes Appl. 13 (1982), 235–248.
  * [36] F. Larribe and P. Fearnhead, On composite likelihoods in statistical genetics, Statistica Sinica 21 (2011), 43–69.
  * [37] T. Lenormand and J. Dutheil, Recombination difference between sexes: a role for haploid selection, PLoS Biol. 3 (2005), e63.
  * [38] M.J. Lercher and L.D. Hurst, Human SNP variability and mutation rate are higher in regions of high recombination, Trends Genet. 18 (2002), 337–340.
  * [39] H. Li and R. Durbin, Inference of human population history from individual whole-genome sequences, Nature 475 (2011), 493–496.
  * [40] N. Li and M. Stephens, Modeling linkage disequilibrium and identifying recombination hotspots using single-nucleotide polymorphism data, Genetics 165 (2003), 2213–2233.
  * [41] S. Li and M. Jakobsson, Estimating demographic parameters from large-scale population genomic data using Approximate Bayesian Computation, BMC Genet. 13 (2012), 22–37.
  * [42] T. Mailund, J.Y. Dutheil, A. Hobolth, G. Lunter, and M.H. Schierup, Estimating divergence time and ancestral effective population size of Bornean and Sumatran orangutan subspecies using a coalescent hidden Markov model, PLoS Genet. 7 (2011), e1001319.
  * [43] T. Mailund, A.E. Halager, M. Westergaard, J.Y. Dutheil, K. Munch, L.N. Andersen, G. Lunter, K. Prüfer, A. Scally, A. Hobolth, and M.H. Schierup, A new isolation with migration model along complete genomes infers very different divergence processes among closely related great ape species, PLoS Genet. 8 (2012), e1003125.
  * [44] T. Mailund, K. Munch, and M.H. Schierup, Lineage sorting in apes, Annu. Rev. Genet. 48 (2014), 519–535.
  * [45] A.-S. Malaspinas (and 74 coauthors), A genomic history of Aboriginal Australia, Nature 538 (2016), 207–214.
  * [46] J.S. Maritz and T. Lwin, Empirical Bayes Methods, Routledge, London, 2018.
  * [47] P. Marjoram and J.D. Wall, Fast “coalescent” simulation, BMC Genet. 7 (2006), 16.
  * [48] G. McVean, P. Awadalla, and P. Fearnhead, A coalescent-based method for detecting and estimating recombination from gene sequences, Genetics 160 (2002), 1231–1241.
  * [49] G.A.T. McVean and N.J. Cardin, Approximating the coalescent with recombination, Philos. Trans. R. Soc. Lond. B, Biol. Sci. 360 (2005), 1387–1393.
  * [50] P.A.P. Moran, Random processes in genetics, Math. Proc. Camb. Phil. Soc., 54 (1958), 60–71.
  * [51] K. Munch, T. Mailund, J.Y. Dutheil, and M.H Schierup, A fine-scale recombination map of the human-chimpanzee ancestor reveals faster change in humans than in chimpanzees and a strong impact of GC-biased gene conversion, Genome Res. 24 (2014), 467–474.
  * [52] S. Myers, L. Bottolo, C. Freeman, G. McVean, and P. Donnelly, A fine-scale map of recombination rates and hotspots across the human genome, Science 310 (2005), 321–324.
  * [53] P.F. Palamara, J. Terhorst, Y.S. Song, and A.L. Price, High-throughput inference of pairwise coalescence times identifies signals of selection and enriched disease heritability, Nat. Genet. 50 (2018), 1311–1317.
  * [54] J.S. Paul, M. Steinrücken, and Y.S. Song, An accurate sequentially Markov conditional sampling distribution for the coalescent with recombination, Genetics 187 (2011), 1115–1128.
  * [55] T.D. Petes, Meiotic recombination hot spots and cold spots, Nature Reviews Genetics 2 (2001), 360–369.
  * [56] K. Prüfer (and 40 coauthors), The bonobo genome compared with the chimpanzee and human genomes, Nature 486 (2012), 527–531.
  * [57] M.D. Rasmussen, M.J. Hubisz, I. Gronau, and A. Siepel, Genome-wide inference of ancestral recombination graphs, PLoS Genet. 10 (2014), e1004342.
  * [58] J. Rogers (and 41 coauthors), The comparative genomics and complex population history of Papio baboons, Sci. Adv. 5 (2019), eaau6947.
  * [59] A. Sand, M. Kristiansen, C.N.S. Pedersen, and T. Mailund, zipHMMlib: a highly optimised HMM library exploiting repetitions in the input to speed up the forward algorithm, BMC Bioinformatics 14 (2013), 339–348.
  * [60] A. Scally (and 70 coauthors), Insights into hominid evolution from the gorilla genome sequence, Nature 483 (2012), 169–175.
  * [61] S. Schiffels and R. Durbin, Inferring human population size and separation history from multiple genome sequences, Nat. Genet. 46 (2014), 919–925.
  * [62] S. Sheehan, K. Harris, and Y.S. Song, Estimating variable effective population sizes from multiple genomes: a sequentially markov conditional sampling distribution approach, Genetics 194 (2013), 647–662.
  * [63] K.L. Simonsen and G.A. Churchill, A Markov chain model of coalescence with recombination, Theor. Popul. Biol. 52 (1997), 43–59.
  * [64] J.P. Spence and Y.S. Song, Inference and analysis of population-specific fine-scale recombination maps across 26 diverse human populations, Sci. Adv. 5 (2019), eaaw9206.
  * [65] J.P. Spence, M. Steinrücken, J. Terhorst, and Y.S. Song, Inference of population history using coalescent HMMs: review and outlook, Curr. Opin. Genet. Dev. 53 (2018), 70–76.
  * [66] P.R. Staab, S. Zhu, D. Metzler, and G. Lunter, scrm: efficiently simulating long sequences using the approximated coalescent with recombination, Bioinformatics 31 (2015), 1680–1682.
  * [67] J.A. Stamatoyannopoulos, I. Adzhubei, R.E. Thurman, G.V. Kryukov, S.M. Mirkin, and S.R. Sunyaev, Human mutation rate associated with DNA replication timing, Nat. Genet. 41 (2009), 393–395.
  * [68] J. Stapley, P.G.D. Feulner, S.E. Johnston, A.W. Santure, and C.M. Smadja, Variation in recombination frequency and distribution across eukaryotes: patterns and processes, Philos. Trans. R. Soc. Lond. B, Biol. Sci. 372 (2017), 20160455.
  * [69] M. Steinrücken, J. Kamm, J.P. Spence, and Y.S. Song, Inference of complex population histories using whole-genome sequences from multiple populations, Proc. Natl. Acad. Sci. U.S.A. 116 (2019), 17115–17120.
  * [70] M. Steinrücken, J.S. Paul, and Y.S. Song, A sequentially Markov conditional sampling distribution for structured populations with migration and recombination, Theor. Popul. Biol.. 87 (2013), 51–61.
  * [71] M. Stephens and P. Scheet, Accounting for decay of linkage disequilibrium in haplotype inference and missing-data imputation, Am. J. Hum. Genet. 76 (2005), 449–462.
  * [72] L.S. Stevison (and 9 coauthors) Great Ape Genome Project, The time scale of recombination rate evolution in great apes, Mol. Biol. Evol. 33 (2016), 928–945.
  * [73] E.H. Stukenbrock, T. Bataillon, J.Y. Dutheil, T.T. Hansen, R. Li, M. Zala, B.A. McDonald, J. Wang, and M.H. Schierup, The making of a new pathogen: insights from comparative population genomics of the domesticated wheat pathogen Mycosphaerella graminicola and its wild sister species, Genome Res. 21 (2011), 2157–2166.
  * [74] E.H. Stukenbrock and J.Y. Dutheil, Fine-scale recombination maps of fungal plant pathogens reveal dynamic recombination landscapes and intragenic hotspots, Genetics 208 (2018), 1209–1229.
  * [75] M.P.H. Stumpf and G.A.T. McVean, Estimating recombination rates from population-genetic data, Nat. Rev. Genet. 4 (2003), 959–968.
  * [76] J. Terhorst, J.A. Kamm, and Y.S. Song, Robust and scalable inference of population history from hundreds of unphased whole genomes, Nat. Genet. 49 (2017), 303–309.
  * [77] A. Wallberg, S. Glémin, and M.T. Webster, Extreme recombination frequencies shape genome variation and evolution in the honeybee, Apis mellifera, PLoS Genet. 11 (2015), e1005189.
  * [78] K. Wang, I. Mathieson, J. O’Connell, and S. Schiffels, Tracking human population structure through time from whole genome sequences, PLoS Genetics 16 (2020), 1–24
  * [79] C.C. Weber, C.J. Pink, and L.D. Hurst, Late-replicating domains have higher divergence and diversity in Drosophila melanogaster, Mol. Biol. Evol. 29 (2012), 873–882.
  * [80] D. Wegmann, C. Leuenberger, S. Neuenschwander, and L. Excoffier, ABCtoolbox: a versatile toolkit for approximate Bayesian computations, BMC Bioinformatics 11 (2010), 116–123.
  * [81] P.R. Wilton, S. Carmi, and A. Hobolth, The SMC’ is a highly accurate approximation to the ancestral recombination graph, Genetics 200 (2015), 343–355.
  * [82] C. Wiuf and J. Hein, Recombination as a point process along sequences, Theor. Popul. Biol. 55 (1999), 248–259.


\mathcal{F}_k(\rho_{k-1}(k)) = f_{\rho}\underbrace{(\rho_{k-1}(k))}_{2^{k-1}d  -dimensional} + f_{\mathbb{I}_d} \underbrace{\frac{\mathbb{I}_d}{d}}_{d -dimensional} \otimes \underbrace{\text{Block-Tr}_{k-1}(\rho_{k-1}(k))}_{2^{k-1} -dimensional}
\end{equation}

\begin{equation}
\mathcal{R}_k(\rho_{k-1}(k)) = r_{\rho}\underbrace{(\rho_{k-1}(k))}_{2^{k-1}d  -dimensional} + r_{\mathbb{I}_d} \underbrace{\frac{\mathbb{I}_d}{d}}_{d -dimensional} \otimes \underbrace{\text{Block-Tr}_{k-1}(\rho_{k-1}(k))}_{2^{k-1} -dimensional}
\end{equation}

Using <ref> and <ref> in the result above:

\begin{align}
(\mathbb{I}_d \otimes \bra{+}^{\otimes k-1})\left(f_{\rho}(\rho_{k-1}(k)) + f_{\mathbb{I}_d}\frac{\mathbb{I}_d}{d}\otimes\text{Block-Tr}_{k-1}(\rho_{k-1}(k))\right)(\mathbb{I}_d \otimes \ket{+}^{\otimes k-1})\\
+ 2\sqrt{\theta\overline{\theta}} (\mathbb{I}_d \otimes \bra{+}^{\otimes k-1})\left(r_{\rho}(\rho_{k-1}(k)) + r_{\mathbb{I}_d}\frac{\mathbb{I}_d}{d}\otimes\text{Block-Tr}_{k-1}(\rho_{k-1}(k))\right)(\mathbb{I}_d \otimes \ket{+}^{\otimes k-1})
\end{aligned}
\end{bmatrix}\\
&\text{rearranging the terms, we get:}\nonumber\\
(f_{\rho} + 2\sqrt{\theta\overline{\theta}}r_{\rho})(\mathbb{I}_d \otimes \bra{+}^{\otimes k-1})(\rho_{k-1}(k))(\mathbb{I}_d \otimes \ket{+}^{\otimes k-1})\\
+  (f_{\mathbb{I}_d}+ 2\sqrt{\theta\overline{\theta}}r_{\mathbb{I}_d})(\mathbb{I}_d \otimes \bra{+}^{\otimes k-1})(\frac{\mathbb{I}_d}{d}\otimes\text{Block-Tr}_{k-1}(\rho_{k-1}(k))(\mathbb{I}_d \otimes \ket{+}^{\otimes k-1})
\end{aligned}\right\}\\
&\text{substituting } (\mathbb{I}_d \otimes \bra{+}^{\otimes k-1})(\rho_{k-1}(k))(\mathbb{I}_d \otimes \ket{+}^{\otimes k-1}) = M_{k-1}(\rho_{k-1}(k))\nonumber\\
&= \frac{1}{2}\left\{
(f_{\rho} + 2\sqrt{\theta\overline{\theta}}r_{\rho})M_{k-1}(\rho_{k-1}(k))
+  (f_{\mathbb{I}_d}+ 2\sqrt{\theta\overline{\theta}}r_{\mathbb{I}_d})(\mathbb{I}_d \otimes \bra{+}^{\otimes k-1})(\frac{\mathbb{I}_d}{d}\otimes\text{Block-Tr}_{k-1}(\rho_{k-1}(k))(\mathbb{I}_d \otimes \ket{+}^{\otimes k-1})\right\}
\end{align}

Thus, we get a recursive recipe to get the density operator after k iterations:

\begin{equation}
M_{k}(\rho_{k}(k)) = \frac{1}{2}\left\{
(f_{\rho} + 2\sqrt{\theta\overline{\theta}}r_{\rho})M_{k-1}(\rho_{k-1}(k))
+  (f_{\mathbb{I}_d}+ 2\sqrt{\theta\overline{\theta}}r_{\mathbb{I}_d})(\mathbb{I}_d \otimes \bra{+}^{\otimes k-1}\otimes\text{Block-Tr}_{k-1}(\rho_{k-1}(k))\otimes \ket{+}^{\otimes k-1}\right\}
\end{equation}

After that, we can analyze the density operator obtained for success probability after k Grover iterations.


11institutetext: Department of Electrical Engineering, Universidad de Chile
11email<EMAIL_ADDRESS>
22institutetext: Advanced Mining Technology Center, Universidad de Chile
22email<EMAIL_ADDRESS>
Authors’ Instructions

# Near Real-Time Object Recognition for Pepper based on Deep Neural Networks
Running on a Backpack

Esteban Reyes 11 Cristopher Gómez 11 Esteban Norambuena 11 Javier Ruiz-del-
Solar 11 2 2

###### Abstract

The main goal of the paper is to provide Pepper with a near real-time object
recognition system based on deep neural networks. The proposed system is based
on YOLO (_You Only Look Once_), a deep neural network that is able to detect
and recognize objects robustly and at a high speed. In addition, considering
that YOLO cannot be run in the Pepper’s internal computer in near real-time,
we propose to use a Backpack for Pepper, which holds a Jetson TK1 card and a
battery. By using this card, Pepper is able to robustly detect and recognize
objects in images of 320x320 pixels at about 5 frames per second.

###### Keywords:

Pepper robot, YOLO, Jetson TK1, ROS

##  1 Introduction

Environments where service and social robots operate/live are highly dynamic,
in part because of the people living in those environments constantly interact
with each other and carry out daily life activities. Therefore, service and
social robots require perception systems that are highly robust, but at the
same time are able to operate in real-time111In the service robot domain, we
understand as real-time tasks that have a reaction span that looks natural to
people, i.e. $\sim$5 Hz.. Object detection and recognition are some of the
vision tasks that require at least near real-time operation.

The commercial Pepper robot [1] is a social robot used to research on human-
robot interaction in real-world human environments (e.g. it is the official
platform of the RoboCup@Home Standard Platform League), but its operation is
constrained to the fact that it lacks the computational power necessary to run
state-of-the-art vision algorithms.

On the other hand, the uprising of Deep Neural Networks (DNNs) has lead
researchers to use them in the development of models that can quickly
recognize objects and persons in a robust manner. However, they are very
expensive in terms of computational power, so they cannot be run directly on
typical service robots. In particular, most of state-of-the-art DNN models
cannot run on Pepper’s internal computer, a situation that stands as the
problem we address in this paper. Therefore, the main goal of this paper is to
propose a solution to this problem that considers two components: First, the
selection of YOLO (_You only look Once_) [2], a DNN with real-time detection
and recognition of objects and persons, as the most suitable DNN to this task;
and second, the development and implementation of an add-on for Pepper, a
backpack that permit the attachment of a single board computer onto Pepper,
particularly a Nvidia Jetson TK1, which can run YOLO at about 5 FPS when
processing images of 320x320 pixels. The use of an external enhancing device,
attached to the robot without modifying its structure, comes as an inspiration
from similar projects developed for the NAO robot [3], where a fully
replicable backpack was built.

In this work we provide details on how to reproduce the backpack for Pepper,
and on how to install and implement YOLO on it, making the backpack CAD model,
hardware, software specifications and an installation guide available for
replication.

The remainder of this paper is structured as follows. In Section 2, we present
an overview of the YOLO network. Subsequently, in Section 3 we present the
process done to adapt and use YOLO on the Pepper’s internal computer and in
the Jetson TK1 card. Then, in Section 4 we address the mechanical design of
the Pepper Backpack. Finally, in Section 5 we draw some conclusions of this
work.

##  2 Robust and Fast Object Recognition using YOLO

_You Only Look Once_ (YOLO) [2], is a computer vision system capable of
detecting a wide variety of objects in a single image, with an accuracy
similar to RetinaNet [4], but with a superior speed of inference when compared
to others state of the art systems such as SSD [5], R-FCN [6] and FPN FRCN
[7]. Its speed makes it one of the best-suited systems for real-time object
detection needed in systems such as service robots.

One of the main features of YOLO is to treat the detection of objects as a
regression problem, where the model is trained to identify bounding boxes
(BB), along with probabilities of certainty, in areas that might contain an
object. Unlike other systems, where _proposals_ or _regions of interest_ are
generated explicitly as object candidate windows, on which a class inference
is executed. This difference gives YOLO advantages in terms of speed, having
to process each image only once to perform multiple detections, instead of
processing _proposals_ individually. Also, the system processes the complete
image and not just a region of it, which makes its inferences contemplate the
global context of the image, making it less likely to detect background
content as an object, which translates into a lower number of false positives
in comparison to other systems.

Another distinctive feature of this model is its end-to-end training process,
meaning that it has a unique _pipeline_ that is trained jointly. Unlike other
systems that have different components and need to be trained separately, such
as Faster-RCNN [8].

Since the first release of YOLO, there have been 3 major versions of the
algorithm, each one aiming to improve accuracy with respect to previous
versions. The first version, YOLOv1 [2], achieved 63.4% mean average precision
(mAP) over PASCAL VOC 2007 [9], with an inference speed of 45 FPS.
Introduction of a fully-convolutional model [5], multi-scale training [10],
batch normalization [11], BB dimensions priors [10], among other techniques,
raised YOLOv2 [10], which gets an mAP of 48.1% on COCO [12] dataset and 78.6%
mAP on PASCAL VOC 2007, while working at 40 FPS. The latest version of YOLO,
YOLOv3 [13], includes a larger model with 75 convolutional layers that use
residual blocks [14], prediction of BB across 3 different scales by using a
procedure similar to feature pyramid networks [7], among other improvements
that result in an mAP of 57.9% on COCO, at 20 FPS on the same TitanX GPU that
all models where tested. Differences among YOLO versions clearly depicts that
performance can be sacrificed to gain processing speed, allowing to choose the
best-suited version for a given application. Moreover, there are reduced
versions of YOLOv1 and YOLOv2 which are even faster. A review of most versions
of YOLO is shown in Table 1.

Table 1: Performance of different YOLO versions. Inference speed results (FPS) where obtained when running on a Titan X GPU. Fastest YOLO version is tiny-YOLO v2 at 207 FPS, and most accurate version is YOLOv3 with a 57.9% mAP on COCO dataset. Model | Input size | Train set | Test set | mAP | FPS
---|---|---|---|---|---
YOLOv1 | 448x448 | VOC 2007+2012 | VOC 2007 | 63.4% | 45
Fast YOLOv1 | 448x448 | VOC 2007+2012 | VOC 2007 | 52.7% | 155
YOLOv2 | 416x416 | VOC 2007+2012 | VOC 2007 | 76.8% | 67
tiny-YOLOv2 | 416x416 | VOC 2007+2012 | VOC 2007 | 57.1% | 207
YOLOv2 | 608x608 | COCO | COCO | 48.1% | 40
YOLOv3 | 608x608 | COCO | COCO | 57.9% | 20

##  3 Adapting YOLO to be used with Pepper Robots

In spite of YOLO’s exceptional performance and speed, it is necessary to take
into account that the reported speeds were measured using a platform with a
powerful GPU. When running YOLO on a CPU-only system, such as the Pepper’s
internal computer, the speed of operation decreases considerably. To deal with
this situation is that we base our work on the fastest versions of the model
(tiny-YOLO), which has only 15 convolutional layers.

###  3.1 YOLO for Pepper

To run YOLO on the Pepper’s on-board computer, the Darknet framework [15] must
be first installed. Darknet is written in C, a low-level language, so it is
easy to port to different platforms. Thus, compilation for Pepper’s computer
is straightforward.

To easily integrate YOLO with other modules of Pepper, the darknet_ros package
from ETH Autonomous System Lab [16] is used. This package operates with an
older version of Darknet, however, it is compatible with the most recent
versions of Darknet. The YOLO ROS package implements all the tools needed to
feed the system with a standard ROS stream of images through a defined topic.
Moreover, the information of the detected objects is also published through a
ROS topic.

To compile the YOLO ROS package for the Pepper’s computer, mild modifications
are needed to source and compilation files. A guide to the process can be
found in the _Supplementary Material_ Section.

###  3.2 YOLO for Jetson TK1

As an alternative strategy to the Pepper’s on-board processing of tiny-YOLO,
an external processing unit is used, where the only task that directly
involves Pepper is to publish on a ROS topic the images from his camera at a
rate that depends on the resolution of the images. The maximum resolution that
enables a rate of $\sim 30$Hz is 640x480. Jetson TK1 development card is used
for the external processing. We select this unit because it has an integrated
2GB Nvidia GPU, which allows the use of parallel computing platform CUDA [17],
accelerating the calculation of certain Darknet operations written exclusively
for parallel processing in CUDA.

As a first approach, the original version of Darknet is compiled, it is
important to note that its compilation with GPU is immediate, but not the use
of the cuDNN library [18] of efficient computing for deep learning. This
incompatibility raises because of the CUDA version available for the Jetson
TK1 is CUDA 6.5, and the library cuDNN from Darknet uses the CUDA 7.0 version,
thus processing speeds will be lower.

Afterwards, YOLO ROS is compiled on the Jetson card, for which ROS Indigo was
previously installed by following the ROS Jetson guide [19]. The YOLO ROS
compilation without GPU is straightforward, but to enable its use some
modifications need to take place. The changes made to YOLO ROS focuses on the
file CMakeList.txt, the scripts in C++, the configuration of the CUDA paths
and extensions of some ROS Indigo files. Due to the extensive and tedious
modifications made to the code, they will be omitted, but these and the entire
installation process in the Jetson TK1 are reviewed in detail in the
_Supplementary Material_ Section.

Finally, a _ROS Network_ is configured to connect Pepper to the Jetson through
an Ethernet cable. Configuring Pepper as MASTER and Jetson as a HOST. This
enables the access of Jetson TK1 card to the topic where images from Pepper’s
camera are published, then process them with tiny-YOLO to finally publish
detection information on another topic that Pepper can subscribe to and use as
desired. A diagram that depicts the connections made is shown in Fig. 1.

Figure 1: Pepper to Jetson TK1 communication diagram. Through the ROS/NaoQi
driver, Pepper publishes images from his camera. The Darknet ROS node runs the
YOLO network that feeds from the images. The detected objects and their
classes are published in a custom ROS message, which can be easily accessed by
software modules running in Pepper’s internal computer.

###  3.3 Results

The training process of tiny-YOLO considers input images that vary in
resolution. In the process, the images are resized every ten training
iterations, to resolutions equal to $320+32n$, where $n\in\\{0,1,...,9\\}$ is
chosen randomly. By doing this, images from 320x320 pixels up to 608x608
pixels are considered in the training process. This multi-scale property
allows to choose the input image resolution at runtime; so the input image
size is used as a system parameter to manage the trade-off between inference
speed and detection accuracy. The original tiny-YOLO configuration, that
considers 416x416 pixel input images, is first tested on a notebook with a
quadcore i5-4210U 1.70GHz CPU, and the inferences speed reached roughly 1 FPS.
Then, input images are reduced to 160x160 pixels, allowing detection of big
objects at $\sim$6 FPS, while small ones were not detected. Afterwards, we
compared tiny-YOLO and tiny-YOLO ROS in the notebook, Pepper, and Jetson at a
fixed input size of 160x160, these results are summarized in Table 2.

Table 2: Comparison between different tiny-YOLO systems and platforms. The inference speed and image size are evaluated. Model | Computer |  | Input image size
---
| Inference
---
speed [FPS]
tiny-YOLO |  | Notebook CPU
---
quadcore i5-4210U 1.70GHz
160x160 |  $\sim$6
tiny-YOLO ROS |  | Notebook CPU
---
quadcore i5-4210U 1.70GHz
160x160 |  $\sim$3
tiny-YOLO | Pepper | 160x160 |  $\sim$0.6
tiny-YOLO ROS | Pepper | 160x160 |  $\sim$0.3
tiny-YOLO ROS |  | ROS Network
---
Pepper-Jetson
160x160 |  $\sim$16

To consider the object detection method as a near real-time system, it must
have an inference speed of at least 5 FPS. Table 2 makes clear that the
Pepper’s on-board computer cannot be used to run the tiny-YOLO ROS system, so
incorporation of the Jetson TK1 unit is needed. In Table 2 is shown that tiny-
YOLO ROS running on Jetson TK1 jointly with Pepper reaches 16 FPS, which gives
a margin to enlarge input image size and improve detections accuracy at cost
of decreasing inference speed. Experiments on how inference speed decreases
when making input image size bigger, are shown in Table 4, where mAP over VOC
2007 test set is calculated to depict that a greater image size means better
detection, specifically of small objects. The highest resolution that could be
tested was 384x384 because a greater one leaves Jetson TK1 without
computational resources.

Table 3: Speed of tiny-YOLO ROS for Pepper-Jetson TK1 at different input image size. Inference speed decreases as input size increases, but larger images get higher mAP on VOC 2007 test set. | Input image
---
size
| Inference
---
speed [FPS]
mAP
160x160 |  $\sim$15.5 | 24.60%
224x224 |  $\sim$5.9 | 37.10%
288x288 |  $\sim$5.8 | 44.72%
320x320 |  $\sim$4.8 | 47.69%
352x352 |  $\sim$4.5 | 50.32%
384x384 |  $\sim$3.6 | 52.62%

Table 4: Class-wise mean average precision (mAP) on VOC 2007 test set for tiny-YOLO ROS in Pepper-Jetson TK1, with an input image size of 320x320. aero | bike | bird | boat | bottle
---|---|---|---|---
47 | 60.1 | 39.5 | 30.6 | 15.8
bus | car | cat | chair | cow
61.2 | 57.8 | 65.6 | 23.9 | 43.1
table | dog | horse | m bike | person
49 | 59.2 | 66.3 | 63.5 | 52.3
plant | sheep | sofa | train | tv
21.9 | 42.2 | 48.1 | 59.8 | 47.1

Table 4 clearly show that tiny-YOLO ROS for Pepper-Jetson TK1 achieves a near
real-time inference speed of $\sim$4.8 FPS at an input image size of 320x320
while getting 47.69% mAP on VOC 2007 test set. Using this same data, class-
wise mAP is calculated on Table 4. This gives a detailed insight on the
performance of the system at 320x320 resolution, where it is clear that
correct detection of small objects like bottles or potted plants is a
challenge for the system.

After evaluating the system on classic data sets, we jumped to try it in an
indoor environment, a place closer to the reality of what a service robot has
to deal with. We analyzed different sizes of the input images, in order to
select the one that allows obtaining an acceptable detection rate as well as
processing speed. Fig. 2(a) shows how for an image size of 160x160, the YOLO
based detection system is unable to detect a bottle and it also throws false
positives by detecting a non-existing car. On the other hand, when using a
resolution of 384x384 pixels (Fig. 2(b)) the YOLO based detection system is
able to correctly detect the bottle and also better fits the BB for a person.

((a))

((b))

Figure 2: Detection examples of tiny-YOLO ROS for Pepper-Jetson TK1, at
different resolution inputs. Color BB with labels shows detected objects. (a)
Input image size of 160x160 pixels, where the model is unable to detect small
objects like the bottle and mistakenly detect a car. (b) Input image size of
384x384 accurately detects multiple small to medium size objects.

When testing the detection system with input image of 320x320 pixels, it
achieved near real-time inference speed, and it showed robust object detection
of people and infrequent false positives. These results are shown in Fig 3,
where multiple sights of the indoor environment where shown to the system.
Numerous object detection can be seen in Fig. 3(a) where 2 people and a chair
are correctly detected. On the other hand, false positives that momentarily
happen are depicted in Fig. 3(b), by detection of a non-existing car. Small
object detection difficulties can be seen in Fig. 3(c) and 3(d), where the
bottle is not detected at first, but its recognized after moving it around,
straighten it up and putting it closer to the camera of the Pepper robot.

Therefore, we conclude that is very relevant to select an appropriate image
size, and that an image size of 320x320 pixels allows obtaining a good
tradeoff between detection accuracy and processing speed.

((a))

((b))

((c))

((d))

Figure 3: Detection examples of tiny-YOLO ROS for Pepper-Jetson TK1, at
$\sim$4.8 FPS and input image size of 320x320. This resolution allows
detection of people and medium to large size objects, while struggles with
small ones (a) Correct detection of multiple objects. (b) Example of an
incorrect car detection, false positives are infrequent. (c) Correct detection
of a person and a chair, but the bottle is ignored. (d) Correct detection of a
small object when moving bottle from (c) closer to the camera.

##  4 Pepper Backpack

Considering that using an external processing card (Jetson TK1) is the best
alternative for implementing a YOLO based object detection system, we decided
to design and built a backpack for Pepper, which can hold the card as well as
its battery.

###  4.1 Mechanical Design

The most important restriction to design the Pepper Backpack was to not modify
the robot structure. For this reason, the attachment of the backpack to the
robot should be non-invasive, therefore, an attachment using a suction cup is
proposed.

The main enclosure of the Jetson computer corresponds to one acrylic plate in
the base and a 3D printed case with a battery compartment mounted on the front
of the board. A custom joint was 3D printed to connect the enclosure of the
Jetson and the suction cup structure. Whole backpack renders can be seen in
Fig. 4, while perspectives of the backpack attached to Pepper robot through a
5.8 $mm$ of diameter commercial suction cup can be seen in Fig. 5.

The CAD models, list of components and materials required to build the
backpack can be found in the Supplementary Material section.

((a))

((b))

((c))

Figure 4: 3D Renders of the Pepper Backpack. The main enclosure (white) is 3D
printed. A commercial suction cup (black) is used to attach the backpack to
Pepper. (a) Front view. (b) Lateral View. (c) Back view.

((a))

((b))

Figure 5: Pepper robot with the Jetson TK1 backpack attached to its lower
back. (a) Perspective view. (b) Lateral close-up tho backpack.

###  4.2 Hardware

####  4.2.1 Nvidia Jetson TK1

The Nvidia Jetson TK1 main processor is a Nvidia Tegra K1 which is a
CPU+GPU+ISP single chip. The existence of a GPU in the Jetson allows it to run
state-of-the-art Deep Learning algorithms.

The Nvidia Jetson TK1 comes with a power supply that can provide 12 $volts$
and 5 $ampere$ maximum. These requirements are fulfilled by a 3-Cell LiPo
battery. The battery is connected to the Jetson through a standard 2.1x5.5$mm$
barrel jack. The typical power-draw of the Jetson does not surpass 10 $watts$
[20]. It is important to note that Pepper does not provide a power output,
thus using it as a source of power for the Jetson card is not possible without
structural modifications of Pepper.

To connect the Jetson computer to Pepper, and enable data stream between
platforms, the Ethernet Gigabit port of both is used.

##  5 Conclusions and Future Work

In this work, we studied the use of YOLO on the Pepper’s computer, by testing
different versions of the model and variation of their parameters. We
concluded that it is not possible to surpass the 1 FPS minimum limit in order
to achieve near real-time object detection. As an alternative, we introduced
external image processing, choosing the Jetson TK1 computer as the device to
run tiny-YOLO at an input image resolution of 320x320, which demonstrated to
be the best-suited model to reach high-speed processing of $\sim$4.8 FPS. By
using a smaller input image size we gained speed, at cost of performance, this
was reflected on an mAP of 47.69% over VOC 2007 test set, and a low capacity
of detection for small objects on real-world indoor environments.

To enable high-speed communication between Pepper and the Jetson TK1 computer
we use an Ethernet Gigabit connection. More important is the fact that we
directly attached the board onto Pepper through a custom made backpack, which
does not affect the movement of the robot.

As a future work task, we propose usage of another external computer besides
Jetson TK1, because it may not be the most suitable platform to best exploit
inference speed of tiny-YOLO. This statement lies in the fact that the Jetson
TK1 is a 32-bit system that supports up to CUDA 6.5 version, which does not
allow usage of the deep neural networks dedicated library cuDNN. Platforms
such as Jetson TX1 or Jetson TX2 which allow higher versions of CUDA and thus
cuDNN usage may outperform results presented in this work.

Finally, all the necessary files to replicate the project will be publicly
available.

## Supplementary Material

All the necessary to replicate the project resides in two GitHub repositories.
One repository provides the CAD models and list of components to build the
Pepper Backpack: https://github.com/uchile-robotics/pepper-
backpackhttps://github.com/uchile-robotics/pepper-backpack. The other
repository corresponds to a fork of darknet_ros with instructions to run tiny-
YOLO ROS on the Pepper-Jetson TK1: https://github.com/uchile-robotics-
forks/darknet_roshttps://github.com/uchile-robotics-forks/darknet_ros

## Acknowledgements

This research was partially funded by FONDECYT Project 1161500.

## References

  * [1] SoftBank Robotics.  Pepper www.ald.softbankrobotics.com/en/robots/pepper.  Last accessed, 20, 2017.
  * [2] Joseph Redmon, Santosh Divvala, Ross Girshick, and Ali Farhadi.  You only look once: Unified, real-time object detection.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 779–788, 2016.
  * [3] Matías Mattamala, Gonzalo Olave, Clayder González, Nicolás Hasbún, and Javier Ruiz-del Solar.  The nao backpack: An open-hardware add-on for fast software development with the nao robot.  arXiv preprint arXiv:1706.06696, 2017.
  * [4] Tsung-Yi Lin, Priya Goyal, Ross Girshick, Kaiming He, and Piotr Dollár.  Focal loss for dense object detection.  arXiv preprint arXiv:1708.02002, 2017.
  * [5] Wei Liu, Dragomir Anguelov, Dumitru Erhan, Christian Szegedy, Scott Reed, Cheng-Yang Fu, and Alexander C Berg.  Ssd: Single shot multibox detector.  In European conference on computer vision, pages 21–37. Springer, 2016.
  * [6] Jifeng Dai, Yi Li, Kaiming He, and Jian Sun.  R-fcn: Object detection via region-based fully convolutional networks.  In Advances in neural information processing systems, pages 379–387, 2016.
  * [7] Tsung-Yi Lin, Piotr Dollár, Ross Girshick, Kaiming He, Bharath Hariharan, and Serge Belongie.  Feature pyramid networks for object detection.  In CVPR, volume 1, page 4, 2017.
  * [8] Shaoqing Ren, Kaiming He, Ross Girshick, and Jian Sun.  Faster r-cnn: Towards real-time object detection with region proposal networks.  In Advances in neural information processing systems, pages 91–99, 2015.
  * [9] M. Everingham, L. Van Gool, C. K. I. Williams, J. Winn, and A. Zisserman.  The PASCAL Visual Object Classes Challenge 2007 (VOC2007) Results.  http://www.pascal-network.org/challenges/VOC/voc2007/workshop/index.html.
  * [10] Joseph Redmon and Ali Farhadi.  Yolo9000: better, faster, stronger.  arXiv preprint, 2017.
  * [11] Sergey Ioffe and Christian Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  arXiv preprint arXiv:1502.03167, 2015.
  * [12] Tsung-Yi Lin, Michael Maire, Serge Belongie, James Hays, Pietro Perona, Deva Ramanan, Piotr Dollár, and C Lawrence Zitnick.  Microsoft coco: Common objects in context.  In European conference on computer vision, pages 740–755. Springer, 2014.
  * [13] Joseph Redmon and Ali Farhadi.  Yolov3: An incremental improvement.  arXiv, 2018.
  * [14] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.
  * [15] Joseph Redmon.  Darknet: Open source neural networks in c.  http://pjreddie.com/darknet, 2016, 2013.
  * [16] ETH Autonomous System Lab.  Darknet ros.  https://github.com/leggedrobotics/darknet_ros.  Accessed: 2018-03-30.
  * [17] Jason Sanders and Edward Kandrot.  CUDA by example: an introduction to general-purpose GPU programming.  Addison-Wesley Professional, 2010.
  * [18] Sharan Chetlur, Cliff Woolley, Philippe Vandermersch, Jonathan Cohen, John Tran, Bryan Catanzaro, and Evan Shelhamer.  cudnn: Efficient primitives for deep learning.  arXiv preprint arXiv:1410.0759, 2014.
  * [19] ROS wiki.  Jetson ros.  http://wiki.ros.org/NvidiaJetsonTK1.  Accessed: 2018-03-30.
  * [20] elinux.org.  Jetson tk1 powerdraw.  https://elinux.org/Jetson/Computer_Vision_Performance#Power_draw_during_computer_vision_tasks.  Accessed: 2018-03-30.