Inference of demographic histories in structured populations with an application to human evolution.

Résumé

L'inférence de l'histoire démographique des espèces représente l'un des plus grands défis de la génétique des populations. Cette histoire est souvent représentée comme une histoire de changements de taille de populations, sous l'hypothèse que la structure des populations peut être ignorée. Alternativement, lorsqu'on suppose qu'il y a une structure, elle est souvent définie a priori comme un arbre de populations et n'est pas inférée. Cette thèse vise à proposer des méthodes et des outils pour reconstruire l'histoire démographique de populations qui ont été structurées en un nombre inconnu de sous-populations pendant de longues périodes de temps.

Nous présentons deux approches de l'inférence démographique en présence de structure. La première est basée sur l'IICR (Inverse Instantaneous Coalescence Rate), fonction du modèle démographique et du schéma d'échantillonnage, qui peut être estimée pour un seul individu diploïde en utilisant la méthode du PSMC (pairwise sequentially Markovian coalescent). La méthode proposée ajuste les courbes IICR observées avec les courbes IICR exactes obtenues dans le cadre de modèles en 'îles symétriques stationnaires par morceaux, et infère le nombre d'îles, leur taille commune et l'importance du flux génétique à différentes périodes de temps. Notre application à un ensemble de cinq PSMCs humains a fourni des histoires démographiques qui sont en accord avec des études précédentes suggérant une structure humaine ancienne. La seconde approche se base sur des échantillons génétiques multiples, et est centrée sur l'utilisation du SFS (Site Frequency Spectrum) comme statistique pour l'inférence démographique. Nous nous concentrons sur l'efficacité du calcul du SFS exact attendu dans le cadre d'un modèle général de structure de population, et nous montrons que, pour des tailles d'échantillon allant jusqu'à 26 haploïdes, il est possible d'obtenir une bonne précision numérique et de bonnes performances pour les modèles d'îles symétriques en exploitant la propriété de parcimonie de la matrice de transition pour le processus de Markov associé.

Introduction

Population genetics has seen a major transformation in the last few decades with the wide availability of full-scale genomic sequences. Global initiatives like the 1000 Human Genomes Project [START_REF] Siva | 1000 genomes project[END_REF] has made it possible to have a highly detailed insight into the nature of genetic diversity. One of the great challenges facing the field today is analyzing this rapidly increasing mountain of information in a sensible and efficient manner [START_REF] Beaumont | Recent developments in genetic data analysis: what can they tell us about human demographic history?[END_REF][START_REF] Johri | The impact of purifying and background selection on the inference of population history: problems and prospects[END_REF]. The development of mathematical models for the evolution of populations plays a central role in this effort, along with the associated statistical and numerical analysis methods. This thesis aims to contribute methods and tools for reconstructing the demographic history of populations that have long been geographically structured into an unknown number of sub-populations.

Reconstructing demographic histories in the context of population genetics is often achieved by fitting the parameters of a model according to the observed sequenced data. Some parameters of interest are past population sizes [START_REF] Beaumont | Detecting population expansion and decline using microsatellites[END_REF][START_REF] Li | Inference of human population history from individual whole-genome sequences[END_REF][START_REF] Schiffels | Inferring human population size and separation history from multiple genome sequences[END_REF][START_REF] Liu | Exploring population size changes using SNP frequency spectra[END_REF][START_REF] Boitard | Inferring population size history from large samples of genome-wide molecular data-an approximate bayesian computation approach[END_REF] or the dates of major demographic events such as population splits or admixture events and their proportions [START_REF] Gutenkunst | Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data[END_REF][START_REF] De Barros Damgaard | The first horse herders and the impact of early bronze age steppe expansions into asia[END_REF]. One classical example of such parameter inference is Watterson's θ W estimator [START_REF] Watterson | On the number of segregating sites in genetical models without recombination[END_REF], which measures the population mutation rate by counting the number of places in the genome where a sample of individuals exhibit differences. A more recent example can be found in Li and Durbin's PSMC method [START_REF] Li | Inference of human population history from individual whole-genome sequences[END_REF], which estimates the ancestral size history of a population using the full sequence of a single diploid individual, under the operating assumption that such population has remained panmictic over time, i.e., isolated and without internal mating structures.

This assumption of panmixia has been widely prevalent in many studies over the years [START_REF] Li | Inference of human population history from individual whole-genome sequences[END_REF][START_REF] Liu | Exploring population size changes using SNP frequency spectra[END_REF]. Indeed, contending with the presence of structure and migration has the potential to greatly complicate any demographic model. Some studies of human evolution have introduced structure in the form of simplified tree models, where the number of populations is fixed a priori, i.e., the topology is not inferred (other aspects, such as divergence times, are inferred). The branches of these trees often represent large continental regions that are themselves assumed to have remained panmictic over long periods of time [START_REF] Gutenkunst | Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data[END_REF][START_REF] Prado-Martinez | Great ape genetic diversity and population history[END_REF][START_REF] Noskova | GADMA: Genetic algorithm for inferring demographic history of multiple populations from allele frequency spectrum data[END_REF]. These panmictic and tree models are useful approximations, and they have proven their utility in building stories of human expansions and population splits. However, if the inferred parameter values or conclusions may fundamentally change depending on the underlying model assumptions, then the resulting stories can be questioned [START_REF] Wakeley | Nonequilibrium migration in human history[END_REF][START_REF] Mazet | On the importance of being structured: instantaneous coalescence rates and human evolution-lessons for ancestral population size inference?[END_REF][START_REF] Scerri | Beyond multiregional and simple out-of-Africa models of human evolution[END_REF].

The present thesis focuses exclusively on models of structured populations. One of our core objectives in this context is to develop methods of demographic inference where the number of populations can be inferred rather than assumed. The other objective is to understand the extent to which the observed genetic diversity within populations can be attributed to the effects of structure and changes in connectivity rather than to changes in population size. The main theoretical background supporting this approach is the structured coalescent [START_REF] Herbots | Stochastic models in population genetics: genealogy and genetic differentiation in structured populations[END_REF] which provides a Markovian framework for the study of the genealogy of a sample of lineages that are dynamically migrating between a set of islands or demes. In Chapter 1 we present a summary of a few key concepts in population genetics such as the Wright-Fisher model and the associated coalescent process and a few of its extensions, including variable population size, deme-based population structure, and the infinite-sites mutation model. We also introduce the Inverse Instantaneous Coalescence Rate (IICR) function and the site frequency spectrum (SFS), both of which are proposed as summary statistics for inferential approaches.

Chapter 2 describes a method for inferring the number of populations and the changes in connectivity in a piecewise-stationary n-island model. The n-island is the simplest model of demographic structure, where all the islands are assumed of the same size, and connectivity is fully symmetrical in the sense that all islands exchange migrants with all other islands at the same rate. The piecewise-stationarity refers to the fact that this migration rate may change in a discrete manner throughout time. The dates of these events are also inferred. This model has some limitations as it ignores spatial distances and other complexities of real species, but the choice for this first attempt was guided by simplicity and computational feasibility. The method in question uses the IICR of a single diploid sample as input data in order to fit the model parameters. First described in [START_REF] Mazet | On the importance of being structured: instantaneous coalescence rates and human evolution-lessons for ancestral population size inference?[END_REF], the IICR is closely linked with the probability distribution of the coalescence times of one pair of lineages. It is an attractive candidate for demographic inference in our context since it can be always defined independently of the assumed demographic model, and previous simulations [START_REF] Chikhi | The IICR (inverse instantaneous coalescence rate) as a summary of genomic diversity: insights into demographic inference and model choice[END_REF]Rodríguez et al., 2018) have shown that it is sensitive to population structure, sampling patterns and fluctuations of migration rates. The IICR can be estimated from real data using methods such as the PSMC.

The process of fitting the data was accomplished using a general purpose meta-heuristic to explore the parameter space and minimize the distance between a candidate IICR and the target (observed) one. The chapter also details the approach and results of our validation efforts, which consisted in simulating IICR curves under different n-island models and comparing the inferred parameters with the known ones. We then applied the method to human genomic data using five published PSMC curves [START_REF] Prado-Martinez | Great ape genetic diversity and population history[END_REF] and compared the inferred histories between individuals and with previously inferred scenarios by Rodríguez et al. (2018) and [START_REF] Noskova | GADMA: Genetic algorithm for inferring demographic history of multiple populations from allele frequency spectrum data[END_REF]. This study was published in [START_REF] Arredondo | Inferring number of populations and changes in connectivity under the n-island model[END_REF].

Chapter 3 focuses on establishing a framework for demographic inference of structured populations using multi-sample summary statistics. Having multiple samples provides better information resolution in the not-so-distant past of the sample, since in typical ancestral trees, most of the total branch length is distributed close to the samples near the bottom of the tree (leafs). The chapter introduces two candidate multi-sample summary statistics: the IICR k and the SFS. The IICR k is a natural generalization of the IICR that is formulated using the random distribution of the waiting time until the first coalescence in a sample of size k, and can be estimated for real data using methods such as the MSMC [START_REF] Schiffels | Inferring human population size and separation history from multiple genome sequences[END_REF]. On the other hand, the SFS is a histogram of allele frequencies that has been widely applied as an efficient way for summarizing the genetic diversity of a sample [START_REF] Wakeley | Estimating ancestral population parameters[END_REF][START_REF] Griffiths | The age of a mutation in a general coalescent tree[END_REF][START_REF] Nielsen | Estimation of population parameters and recombination rates from single nucleotide polymorphisms[END_REF][START_REF] Gutenkunst | Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data[END_REF][START_REF] Excoffier | Robust demographic inference from genomic and snp data[END_REF].

We show in this chapter that the Markovian framework of the structured coalescent can be generalized to capture the necessary information for computing these functions within arbitrary numerical accuracy. Similar approaches have been avoided in the past due to the computational challenges of handling the state spaces of the associated Markov processes, which grow in size almost exponentially as a function of the number of samples. However, we show that by not tracking the historical origin of the lineages (mainly used to compute the joint SFS of the populations), using appropriate model specializations that take advantage of inherent symmetries, and exploiting of the sparsity patterns that appear in the rate matrix under certain state orderings, it is possible to very quickly compute the expected SFS of samples of size up to k = 26 in models of symmetrical n-island structure. The approach is general in the sense that it allows for many potential types of structured scenarios, however we exclusively showcase the n-island as a model specialization for the implementation and results, once again due to its mathematical simplicity.

All along this work we have strived to not only establish novel methods for approaching the questions of demographic inference in the context of population structure, but also to provide implementations of these methods that are both intuitively usable and highly performant. The codebases for these projects are available in github.com/arredondos/snif and github.com/arredondos/sisifs. Many questions remain unanswered. In the final chapter on conclusions and future work, we attempt to paint a broad map of where the results of this work can fit in a larger framework of demographic inference for structured populations by indicating some of the ways in which the presented methods can be generalized.

Introduction (version française)

La génétique des populations a connu une transformation majeure au cours des dernières décennies, grâce à une large disponibilité de séquences génomiques complètes. Des initiatives mondiales comme le projet 1000 génomes humains [START_REF] Siva | 1000 genomes project[END_REF] ont permis d'avoir un aperçu très détaillé de la nature de la diversité génétique. L'un des grands défis auxquels le domaine est confronté aujourd'hui est l'analyse de cette masse d'information en rapide augmentation, d'une manière pertinente et efficace [START_REF] Beaumont | Recent developments in genetic data analysis: what can they tell us about human demographic history?[END_REF][START_REF] Johri | The impact of purifying and background selection on the inference of population history: problems and prospects[END_REF]. Le développement de modèles mathématiques pour l'évolution des populations ainsi que les méthodes d'analyses statistique et numérique associées jouent un rôle central. Cette thèse a pour but d' apporter des méthodes et des outils de reconstruction de l'histoire démographique de populations qui sont depuis longtemps structurées géographiquement en un nombre inconnu de sous-populations.

La reconstruction de l'histoire démographique dans le contexte de la génétique des populations est souvent réalisée en ajustant les paramètres d'un modèle en fonction des données séquentielles observées. Certains paramètres d'intérêt sont les tailles de population passées [START_REF] Beaumont | Detecting population expansion and decline using microsatellites[END_REF][START_REF] Li | Inference of human population history from individual whole-genome sequences[END_REF][START_REF] Schiffels | Inferring human population size and separation history from multiple genome sequences[END_REF][START_REF] Liu | Exploring population size changes using SNP frequency spectra[END_REF][START_REF] Boitard | Inferring population size history from large samples of genome-wide molecular data-an approximate bayesian computation approach[END_REF] ou les dates d'événements démographiques majeurs tels que les divisions de population ou les événements de mélange et leurs proportions [START_REF] Gutenkunst | Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data[END_REF][START_REF] De Barros Damgaard | The first horse herders and the impact of early bronze age steppe expansions into asia[END_REF]. Un exemple classique de cette inférence de paramètres est l'estimateur θ W de Watterson [START_REF] Watterson | On the number of segregating sites in genetical models without recombination[END_REF], qui mesure le taux de mutation de la population en comptant le nombre d'endroits du génome où un échantillon d'individus présente des différences. Un exemple plus récent se trouve dans la méthode PSMC de Li et Durbin [START_REF] Li | Inference of human population history from individual whole-genome sequences[END_REF], qui estime l'histoire de la taille ancestrale d'une population en utilisant la séquence complète d'un seul individu diploïde, sous l'hypothèse que cette population est restée panmictique au fil du temps, c'est-à-dire isolée et sans structures internes de reproduction.

Cette hypothèse de panmixie a été largement adoptée dans de nombreuses études au fil des ans [START_REF] Li | Inference of human population history from individual whole-genome sequences[END_REF][START_REF] Liu | Exploring population size changes using SNP frequency spectra[END_REF]. En effet, faire face à la présence de la structure et de la migration complique considérablement tout modèle démographique. Certaines études de l'évolution humaine ont introduit une structure sous la forme de modèles d'arbres simplifiés, où le nombre de populations est fixé a priori, c'est-à-dire que la topologie n'est pas inférée (d'autres aspects, comme les temps de divergence, le sont). Les branches de ces arbres représentent souvent de grandes régions continentales qui sont elles-mêmes supposées être restées panmictiques sur de longues périodes de temps [START_REF] Gutenkunst | Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data[END_REF][START_REF] Prado-Martinez | Great ape genetic diversity and population history[END_REF][START_REF] Noskova | GADMA: Genetic algorithm for inferring demographic history of multiple populations from allele frequency spectrum data[END_REF]. Ces modèles panmictiques et arborescents sont des approximations utiles, et ils ont prouvé leur utilité pour construire des histoires d'expansions humaines et de scissions de populations. Cependant, si les valeurs des paramètres déduites ou les conclusions peuvent changer fondamentalement en fonction des hypothèses du modèle sous-jacent, alors les histoires qui en résultent peuvent être remises en question [START_REF] Wakeley | Nonequilibrium migration in human history[END_REF][START_REF] Mazet | On the importance of being structured: instantaneous coalescence rates and human evolution-lessons for ancestral population size inference?[END_REF][START_REF] Scerri | Beyond multiregional and simple out-of-Africa models of human evolution[END_REF].

La présente thèse se concentre exclusivement sur les modèles de populations structurées. Un de nos objectifs principaux dans ce contexte est de développer des méthodes d'inférence démographique où le nombre de populations peut être déduit plutôt que supposé. L'autre objectif est de comprendre dans quelle mesure la diversité génétique observée au sein des populations peut être attribuée aux effets de la structure et des changements de connectivité plutôt qu'aux changements de taille de la population. Le principal fondement théorique de cette approche est le coalescent structuré [START_REF] Herbots | Stochastic models in population genetics: genealogy and genetic differentiation in structured populations[END_REF] qui fournit un cadre markovien pour l'étude de la généalogie d'un échantillon de lignées qui migrent dynamiquement entre un ensemble d'îles ou de dèmes. Dans le chapitre chap :prelim, nous présentons un résumé de quelques concepts clés en génétique des populations tels que le modèle de Wright-Fisher et le processus de coalescence associé, ainsi que quelques-unes de ses extensions, notamment la taille variable de la population, la structure de la population basée sur les dèmes et le modèle de mutation à sites infinis. Nous présentons également la fonction inverse du taux de coalescence instantané (IICR) et le spectre de fréquence des sites (SFS), qui sont tous deux proposés comme statistiques pour les approches inférentielles.

Le chapitre 2 décrit une méthode permettant de déduire le nombre de populations et les changements de connectivité dans un modèle n-île stationnaire par morceaux. Le modèle n-île est le modèle le plus simple de structure démographique, dans lequel toutes les îles sont supposées avoir la même taille, et la connectivité est totalement symétrique dans le sens où toutes les îles échangent des migrants avec toutes les autres îles au même rythme. La stationnarité par morceaux fait référence au fait que ce taux de migration peut changer de manière discrète au cours du temps. Les dates de ces événements sont également déduites. Ce modèle présente certaines limites car il ignore les distances spatiales et d'autres complexités des espèces réelles, mais le choix pour cette première tentative a été guidé par la simplicité et la faisabilité informatique. La méthode en question utilise l'IICR d'un seul échantillon diploïde comme données d'entrée afin d'ajuster les paramètres du modèle. Décrit pour la première fois dans [START_REF] Mazet | On the importance of being structured: instantaneous coalescence rates and human evolution-lessons for ancestral population size inference?[END_REF], l'IICR est étroitement lié à la distribution de probabilité des temps de coalescence d'une paire de lignées. Il s'agit d'un candidat intéressant pour l'inférence démographique dans notre contexte puisqu'il peut toujours être défini indépendamment du modèle démographique supposé, et des simulations antérieures [START_REF] Chikhi | The IICR (inverse instantaneous coalescence rate) as a summary of genomic diversity: insights into demographic inference and model choice[END_REF]Rodríguez et al., 2018) ont montré qu'il est sensible à la structure de la population, aux modèles d'échantillonnage et aux fluctuations des taux de migration. L'IICR peut être estimé à partir de données réelles en utilisant des méthodes telles que le PSMC.

Le processus d'ajustement des données a été réalisé à l'aide d'une méta-heuristique d'usage général pour explorer l'espace des paramètres et minimiser la distance entre un IICR candidat et l'IICR cible (observé). Ce chapitre détaille également l'approche et les résultats de nos efforts de validation, qui ont consisté à simuler des courbes d'IICR sous différents modèles n-îles et à comparer les paramètres déduits avec les paramètres connus. Nous avons ensuite appliqué la méthode à des données génomiques humaines en utilisant cinq courbes PSMC publiées [START_REF] Prado-Martinez | Great ape genetic diversity and population history[END_REF] et comparé les histoires inférées entre les individus et avec les scénarios inférés précédemment par Rodríguez et al. (2018[START_REF] Noskova | GADMA: Genetic algorithm for inferring demographic history of multiple populations from allele frequency spectrum data[END_REF]. Cette étude a été publiée dans [START_REF] Arredondo | Inferring number of populations and changes in connectivity under the n-island model[END_REF].

Le chapitre 3 établit un cadre pour l'inférence démographique de populations structurées à l'aide de statistiques multi-échantillons. Le fait d'avoir plusieurs échantillons permet une meilleure résolution des informations dans un passé pas si lointain de l'échantillon, puisque dans les arbres ancestraux typiques, la majeure partie de la longueur totale des branches est localisée à proximité des échantillons situés au bas de l'arbre (feuilles). Ce chapitre présente deux statistiques multi-échantillons candidates : l'IICR k et le SFS. L'IICR k est une généralisation naturelle de l'IICR formulé en utilisant la distribution aléatoire du temps d'attente jusqu'à la première coalescence dans un échantillon de taille k, et peut être estimé pour des données réelles en utilisant des méthodes telles que le MSMC [START_REF] Schiffels | Inferring human population size and separation history from multiple genome sequences[END_REF]. D'autre part, le SFS est un histogramme des fréquences alléliques qui a été largement appliqué car c'est un moyen efficace de résumer la diversité génétique d'un échantillon [START_REF] Wakeley | Estimating ancestral population parameters[END_REF][START_REF] Griffiths | The age of a mutation in a general coalescent tree[END_REF][START_REF] Nielsen | Estimation of population parameters and recombination rates from single nucleotide polymorphisms[END_REF][START_REF] Gutenkunst | Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data[END_REF][START_REF] Excoffier | Robust demographic inference from genomic and snp data[END_REF].

Nous montrons dans ce chapitre que le cadre markovien du coalescent structuré peut être généralisé pour capturer les informations nécessaires au calcul de ces fonctions avec une précision numérique arbitraire. Des approches similaires ont été évitées dans le passé en raison des difficultés de calcul liées à la manipulation des espaces d'état des processus de Markov associés, dont la taille croît de manière presque exponentielle en fonction du nombre d'échantillons. Cependant, nous montrons qu'en ne suivant pas l'origine historique des lignées (principalement utilisée pour calculer le SFS conjoint des populations), en utilisant des spécialisations de modèle appropriées qui tirent profit des symétries inhérentes, et en exploitant les propriétés de parcimonie qui apparaissent dans la matrice de taux sous certains ordres d'état, il est possible de calculer très rapidement le SFS attendu d'échantillons de taille allant jusqu'à k = 26 dans des modèles de structure symétrique à n-îles. L'approche est générale dans le sens où elle permet de nombreux types de scénarios structurés. Cependant nous présentons exclusivement le modèle en îles pour la mise en oeuvre et les résultats, une fois encore en raison de sa simplicité mathématique.

Tout au long de ce travail, nous nous sommes efforcés non seulement d'établir de nouvelles méthodes pour aborder les questions d'inférence démographique dans le contexte de la structure de la population, mais aussi de fournir des implémentations de ces méthodes qui soient à la fois intuitivement utilisables et hautement performantes. Les programmes sont disponibles dans github.com/arredondos/snif et github.com/arredondos/sisifs.

De nombreuses questions restent sans réponse. Dans le dernier chapitre, consacré aux conclusions et aux travaux futurs, nous tentons de dresser une carte générale de la place des résultats de ce travail dans un cadre plus large d'inférence démographique pour les populations structurées, en indiquant certaines des façons dont les méthodes présentées peuvent être généralisées.

List of publications

The present thesis is based on the following publications and manuscripts:

xkcd.com

Chapter 1.

Preliminaries

In this chapter we briefly summarize some basic elements of population genetics in order to make the document more self-contained, and also to establish a few of the naming conventions and assumptions that will be referenced in later chapters. For the sake of brevity, the material here is discussed quite tersely, and only the theoretical concepts that are most closely related to the themes of the following chapters are covered. For a more complete and through exposition of population genetics theory, see [START_REF] Hein | Gene genealogies, variation and evolution: a primer in coalescent theory[END_REF]; [START_REF] Tavaré | Part I: Ancestral inference in population genetics[END_REF]; [START_REF] Gillespie | Population genetics: a concise guide[END_REF].

The Wright-Fisher model

The Wright-Fisher model in its most simple version describes the evolution of a population of 2N haploid individuals (often interpreted as genes) with minimal assumptions. The generations are discrete and non-overlapping, and generation t + 1 is created by making each individual a copy of one of the 2N from generation t, chosen randomly with probability 1/(2N).

The model in this form makes a number of strong assumptions about the population, which evolves in a sequence of discrete time steps by replacing the entire generation with its descendants in every step. Each descendant is an exact copy if its predecessor, so the model does not have a concept of sexual reproduction or gene recombination. Additionally, since each predecessor is selected randomly with uniform probability from the population, there is no concept of demographic structure, fitness or selection in the model.

Typically, the populations we are modeling consist of the genes or alleles of a genome, instead of individual organisms. In the following sections, we focus on a backwards-in-time view of the Wright-Fisher model, known as the coalescent, and relax some of the assumptions made here in favor of having more useful models.

With this it follows that the probability for a sampled pair of lineages to find a common ancestor in exactly g generations in the past is:

1 - 1 2N g-1 1 2N , (1.2.1)
since they choose different ancestors in the first g -1 generations, and then choose the same ancestor in the next one. In this case we say that the pair of genes coalesces in g generations. We denote the time until coalescence of two lineages by T N 2 , with a probability distribution P (T N 2 = g) as in (1.2.1), which corresponds to a geometric distribution with parameter 1/(2N) and an expected value of 2N.

Generalizing for a sample of size k lineages, we have that the probability of there being no coalescences in the previous generation, again counting favorable and total outcomes, is:

P (T N k > 1) = (2N)(2N -1) • • • (2N -k + 1) (2N) k = 1 - 1 2N 1 - 2 2N • • • 1 - k -1 2N = 1 k-1 + -1 2N 1 k-2 + -2 2N 1 k-2 + • • • + -k + 1 2N 1 k-2 + O 1 N 2 = 1 - 1 2N k-1 i=1 i + O 1 N 2 = 1 - k(k -1) 4N + O 1 N 2 , (1.2.2)
where O 1/N 2 captures all the terms that are divided by N 2 or a higher power. If we introduce the additional assumption that k ≪ N, these terms are negligible and can be discarded. In terms of the model, this is equivalent to assuming that when the size of a random sample is very small compared to the size of the population, it is extremely unlikely that more than two lineages will have a common ancestor in the previous generation, and thus this possibility can be excluded.

With this approximation we have that the probability of observing the first coalescence in a sample of size k in exactly g generations (becoming then a sample of size k -1) is given by:

P (T N k = g) ≈ 1 - k(k -1) 4N g-1 k(k -1) 4N . (1.2.3)
Consequently, the variable T N k can be approximated by a geometric distribution with parameter k(k -1)/(4N) = k 2 /(2N). Since the ancestors of the genes are chosen with uniform probability, the times T N k , T N k-1 , . . . , T N 2 are independent.

The continuous-time approximation

In the previous sections we measure time in discrete units of time (generations), and the corresponding T N k random variables reflect this. However, for large populations sizes and long-running processes (large N and g), the coalescent can be well approximated by a continuous-time process where time is scaled in proportion to the population size. Using this transformation, the coalescence times become independent of the population size. For instance, if we take t = g/(2N), then the probability that a pair of lineages will note coalesce before g generations is:

P (T N 2 > g) = 1 - 1 2N g ⇔ P T N 2 2N > t = 1 - 1 2N 2Nt ≈ e -t ⇔ P T N 2 2N t ≈ 1 -e -t ,
therefore the variable T 2 := T N 2 /(2N) can be approximated by an exponential distribution of parameter 1 for large values of N. Using a similar transformation we can see that the scaled times T k := T N k /(2N) can be approximated by exponential distributions of parameter k 2 :

P (T N k > g) = 1 - k 2 1 2N g ⇔ P (T k t) ≈ 1 -e -( k 2 ) t .
In order to reduce the notation complexity, the approximation T k ≈ exp k 2 will be implicitly assumed in the following equations. With this we can compute the expected value of several quantities related to the genealogy of the samples. For instance, the expected time until the first coalescence among k samples is:

E(T k ) = 1 k 2 = 2 k(k -1)
.

Other quantities of interest are the expected time until the most recent common ancestor of the samples, given by the total height of the tree:

T MRCA = T k + T k-1 + • • • + T 2 : E(T MRCA ) = k i=2 E(T i ) = 2 k i=2 1 i (i -1) = 2 1 - 1 k ; (1.2.4)
and the expected total branch length of the tree, which can be computed by weighing the times T i according to the number of live lineages:

E(BL k ) = k i=2 i E(T i ) = 2 k i=2 1 i -1 = 2 k-1 i=1 1 i . (1.2.5)
We can see from equations (1.2.4) and (1.2.5) that taking large samples has quickly diminishing returns, since the expected tree height is bounded by 2, and the expected total branch length grows on the order of log(k).

An important benefit of looking at the population evolution through the lens of the coalescent is that simulating gene genealogies becomes computationally feasible, even for very long-running processes and very large population sizes. Indeed, given a sample of k genes in a Wright-Fisher population, we can simulate the exponential time T k , and then choose a random pair from the k 2 possible pairs to form a coalescence. Then replace k with k -1 and repeat until the MRCA of the sample is reached. This approach, with many other added capabilities for more general models, is notably used by the ubiquitous simulation software ms [START_REF] Hudson | ms a program for generating samples under neutral models[END_REF].

Demographic models

In this section we discuss some generalizations to the Wright-Fisher model and the associated coalescent process that allow for more flexibility in the population evolution, and thus more realistic and useful models featuring fluctuating population sizes and population structure.

Population size change

In the previous sections we have been assuming that the evolution of the population is constrained to having a constant size throughout time. we now relax this assumption by introducing a population size function. Returning to the discrete-time framework, we denote by N(g) the number of haploid individuals or sequences in the population g generations before the present. Therefore, N(0) = 2N would be the current population size, N(1) the size of the generation prior, and so on. We note that in this model, the probability for two genes sampled in the present to have a common ancestor in the previous generation is

1 N(1)
.

Continuing this argument g generations in the past we have that the probability that a pair of genes have not yet coalesced is given by:

P (T N 2 > g) = g i=1 1 - 1 N(i ) . (1.3.1)
We can consider a continuous-time version of N(g) by rounding up the coalescent times 2Nt to an integer number of generations. With this we can define N c (t) := N(⌈2Nt⌉) as well as the continuous-time relative size function λ(t) = N c (t)/N c (0). We now derive the distribution of the time T 2 for this size changing population. Recall that for x ∈ (0; 1) we have the inequality

x -ln(1 -x) x 1 -x ,
which after substituting x for 1/N(i ) and adding for i = 1, . . . , g reduces to

g i=1 1 N(i ) - g i=1 ln 1 - 1 N(i ) g i=1 1 N(i ) -1 . (1.3.2)
We note that with the change of variable τ = i /2N we can have the transformation

g i=1 1 N(i ) = g i=1 1 N c (i /2N) = g 1 di N c (i /2N) = g/2N 0 dτ λ(τ ) .
Additionally, we know from (1.3.1) that

g i=1 ln 1 - 1 N(i ) = ln g i=1 1 - 1 N(i ) = ln P (T N 2 > g) = ln P (T 2 > g/2N)
therefore (1.3.2) implies that for large N:

-ln P (T 2 > g/2N) ≈ g/2N 0 dτ λ(τ ) ⇔ P (T 2 > t) ≈ exp - t 0 dτ λ(τ ) . (1.3.3)
Equation (1.3.3) indicates that there is a very close relation between the relative size change function λ(t) of a varying size Wright-Fisher population and the distribution of the time until the coalescence of a pair of genes sampled from such a population. We return to this idea in §1.4.3.

A similar result can be had for the distribution of T k , i.e., the time until the first coalescence in a sample of size k. Indeed, by counting favorable and total cases, we have:

P (T N k 1) = N(1)(N(1) -1)(N(1) -2) • • • (N(1) -(k -1)) N(1) k = k-1 i=1 1 - i N(1)
, which, after following a series of derivations similar to those of (1.2.2), conduces to:

P (T N k 1) ≈ 1 - k 2 1 N(1) ⇒ P (T N k g) ≈ g i=1 1 - k 2 1 N(i ) ≈ 1 - k 2 g i=1 1 N(i ) = 1 - k 2 g/2N 0 dτ λ(τ ) ,
from which we get:

P (T k > t) ≈ exp - k 2 t 0 dτ λ(τ ) .
(1.3.4)

Population structure

Another one of the strong assumptions in the basic Wright-Fisher model relates to the fact that parents are chosen at random with equal probability. In a mating population, this would translate to the assumption that every mating pair is equally as likely to produce offspring. In reality, physical distances often impose an important barrier for finding mates, and geographically close pairs are more likely to produce offspring than those that are much further apart. We can model this idea in the Wright-Fisher scheme using the finite Herbots's model of the structured coalescent [START_REF] Herbots | Stochastic models in population genetics: genealogy and genetic differentiation in structured populations[END_REF]. We assume the existence of n populations or demes that behave as haploid Wright-Fisher models of size N i = 2s i N genes each, where s i is the relative deme size and N is large as usual. Migration occurs between demes as in each generation a proportion q ij of lineages migrates from deme i to deme j. We denote by m i the proportion of the population of deme i that was received from other demes in any given generation, such that m i = i =j q ji s j /s i . Herbots showed that measuring time in units of 2N generations and making N go to infinity in such a way that the number of migrants stays bounded, the model converges to a continuous-time Markov process. In this transformation, q ij goes to zero in such a way that the product 2Nq ji s j /s i converges and has limit M ij /2. Thus we can express the transition rates in the rate matrix Q of the Markov process as functions of n, s i and M ij .

There are several different types of demographic models. They differ by the configuration of their deme sizes and the patterns of connectivity. Some notable examples (see Figure 1.2):

• The symmetrical islands models, or n-island model for short [START_REF] Wright | Evolution in Mendelian populations[END_REF]): all islands have the same size, and all the migration rates M ij between any pair of islands i and j are equal, so we use the notation M = (n -1)M ij to denote the migration rate received by any given island.

• Continent-island models: one island (the continent) is significantly larger in size than the rest, which are of the same size. Additionally, migration only occurs from the continent to the islands with rate M 1 . In some cases, migration back from the islands into the continent is modeled with another rate M 2 .

• Stepping-stone models: the islands are spatially arranged, and migration can only occur between neighboring islands, with symmetrical rate M. Many different spatial arrangements are possible, but popular ones are rectangular and linear arrangements.

It is important to note that in general structured scenarios, no closed-form expressions can be obtained for quantities of interest like E(T 2 ) or the site frequency spectrum (SFS), for instance. One recurring theme of this dissertation is exploring such properties in special cases of structured scenarios.

Demographic inference

Demographic inference is one of the central topics of this work. We are interested in extracting information from observed data in the form of estimated model parameters. There are several different coalescent models attempting to describe various aspects of the genealogy and using different sets of assumptions, and there is often a multitude of parameters defining such models. Regardless of these variations, mutations are the source of observed genetic diversity, and any model that attempts to explain real-world data must account for the presence of mutations.

A model of mutation

Here we explore a generalization of the Wright-Fisher model where the assumption that each descendant is exactly identical to its parent is lifted. We introduce genetic types by imposing a mutation process on top of the evolution process, where each time a lineage or gene copies its parent from the previous generation, a mutation is introduced with probability u, and thus the offspring becomes a different genetic type. Alternatively, the new lineage is an exact copy of its ancestor with with probability 1u.

Under the lens of the coalescent we may be interested in knowing, for instance, the average number of generations we need to trace back in the genealogical history of a single lineage in order to find a mutation. In each step backwards in time, the probability that the lineage experienced a mutation has a Bernoulli distribution of parameter u, therefore the probability that the lineage has mutated after exactly g backward generations is given by

P (T N M = g) = u(1 -u) g-1 ,
where T N M is the time until first mutation of a sampled lineage, with a geometric distribution of parameter u. Measuring time in units of 2N generations, we can once again obtain a convenient continuous-time approximation for large N:

P (T N M g) = 1 -(1 -u) g , ⇔ P (T M t) ≈ 1 -e -tθ/2 , (1.4.1)
where T M denotes time in units of 2N generations and θ := 4Nu is the scaled mutation rate. This parameter can be interpreted as the expected number of mutations separating a sample of two lineages, since the expected time until their coalescence is 2N generations, and thus 2Nu mutations are expected on each branch. Equation 1.4.1 indicates that T M can be approximated by an exponential variable of parameter θ/2 for large values of N. A consequence of this is that the number of mutations in a branch of length t in the genealogical tree is Poisson distributed, with intensity tθ/2. Therefore, if we denote this number of mutations M t , we have:

P (M t = k) = (tθ) k k!2 k e -tθ/2 . (1.4.2)
Knowing the distribution of M t allows us to simulate the mutations efficiently. Indeed, given a simulated genealogical tree, we can randomly place M t mutations on each branch of length t to obtain a realization of the coalescent process with mutation.

The SFS

We can model the process in which mutations occur in genomic sequences in a simple way by introducing the assumption that the probability for a mutation to occur in the same position twice is null. This assumption is justified in the presence of very long sequences with low mutation rates, and it is known as the infinite sites model of mutation, in which mutations always happen in new locations. Under this assumption, we can assign binary states to the positions along the genome. The two states, which can be coded as 1 and 0, can be interpreted as the position containing a mutation or not, respectively. We can thus represent a genomic sequence as a string of zeroes and ones. If two genes are observed to have different codes in a given position, we say that the position is a segregating or polymorphic site. We could distinguish a mutated position from a non-mutated one, for instance, by comparing with a sampled sequence from an outgroup species.

A convenient way of summarizing the genetic diversity of a sample under these assumptions is to simply inspect the values of the sampled sequences in the segregating sites, given that the sequences match exactly in all other sites. Consider the following example, where k = 6 sequences were sampled, and 8 segregating sites were observed among them.

The rows of Table 1 contain the sequence encoding for each of the segregating positions (columns, labeled SNP j for single-nucleotide polymorphism number j). For instance, we can see that 4 of the 8 polymorphisms were singletons, meaning that a mutation happened in exactly one of the sampled sequences. We can similarly track the number of sites that exhibited exactly i mutations, and compile this information in a vector ξ = (ξ 1 , . . . , ξ k-1 ), where ξ i is the number of segregating sites at which there are exactly i copies of the mutant

SNP 1 SNP 2 SNP 3 SNP 4 SNP 5 SNP 6 SNP 7 SNP 8 Sample 1 0 1 0 0 0 0 1 0 Sample 2 1 0 1 0 0 0 1 0 Sample 3 0 1 1 0 0 1 0 0 Sample 4 0 0 0 0 1 0 1 1 Sample 5 0 0 1 0 0 0 1 0 Sample 6 0 0 0 1 0 1 1 0 Total 1 2 3 1 1 2 5 1
Table 1.1.: An example for how to compute the observed site frequency spectrum (SFS) in a sample of size k = 6 genes. The columns show the segregating sites (SNPs), which are the positions in the genome where a mutation is observed for at least one of the sampled sequences. With increased sequence length, the table grows horizontally, but the SFS is always of length k -1 because this is the maximum possible number of observed mutations (a column full of 1s is not possible since this would correspond to a mutation taking place before the most recent common ancestor of the sample, and as such it would not visible in the coalescence tree). type. In the example of Table 1, ξ = (4, 2, 1, 0, 1). Under the assumptions of a Wright-Fisher model with constant population size and scaled mutation rate θ, it can be shown that:

E(ξ i ) = θ i , i = 1, . . . , k -1. (1.4.3)
See [START_REF] Tajima | Statistical method for testing the neutral mutation hypothesis by dna polymorphism[END_REF]; [START_REF] Fu | Statistical properties of segregating sites[END_REF]; [START_REF] Griffiths | The age of a mutation in a general coalescent tree[END_REF]; [START_REF] Hudson | A new proof of the expected frequency spectrum under the standard neutral model[END_REF] for various proofs of (1.4.3). The vector ξ is known as the Site Frequency Spectrum, and it has been the basis for many estimators and test statistics for analyzing genomic data. Notable examples are the θ W [START_REF] Watterson | On the number of segregating sites in genetical models without recombination[END_REF] and θ π [START_REF] Tajima | Evolutionary relationship of dna sequences in finite populations[END_REF] estimators of θ; Tajima's D [START_REF] Tajima | Statistical method for testing the neutral mutation hypothesis by dna polymorphism[END_REF]; Fay and Wu's H [START_REF] Fay | Hitchhiking under positive darwinian selection[END_REF]; and the fixation index F ST [START_REF] Wright | Isolation by distance[END_REF].

The IICR

We show in §1.3.1 that the distribution of T 2 times in a size-changing Wright-Fisher population is closely related to the relative size function λ(t). Indeed, in (1.3.3) we can solve for λ(t) to obtain:

λ(t) = P (T 2 > t) f T 2 , (1.4.4)
where f T 2 is the probability density function of the T 2 random variable. There are methods that can estimate the right-hand side of (1.4.4) from real data, most notably the PSMC method from [START_REF] Li | Inference of human population history from individual whole-genome sequences[END_REF] estimates T 2 by analyzing the local density of the number of segregating sites along the genome of a single diploid individual. This method has been used to estimate the function λ(t) in several applications, but it is notable that the quantity P (T 2 > t)/f T 2 can always be defined, regardless of the underlying evolutionary model or its assumptions, therefore it is useful to study it independently of the assumptions of a size-changing Wright-Fisher population. This was first done by [START_REF] Mazet | On the importance of being structured: instantaneous coalescence rates and human evolution-lessons for ancestral population size inference?[END_REF], which referred to it as the Inverse Instantaneous Coalescence Rate, or IICR for short.

It is important to stress that the IICR is not tied to any particular model, panmictic or otherwise. It is defined using the distribution of the coalescence times of a sample of size two as in (1.4.4). It can be approximated to arbitrary numerical precision under many demographic models (Rodríguez et al., 2018); it can also be computed empirically by simulating a sample of coalescent times [START_REF] Chikhi | The IICR (inverse instantaneous coalescence rate) as a summary of genomic diversity: insights into demographic inference and model choice[END_REF]; and it can be read from full sequence genomic data using the appropriate methods, like the previously cited PSMC.

xkcd.com

Chapter 2.

Demographic inference using the IICR

To use the IICR as a summary of genomic information we first assume that an IICR curve can be obtained, which we will use as the target for demographic inference. With simulated data (sequences or T 2 values) this target curve can be obtained under any pre-defined coalescent model that could be expressed with a simulation tool (e.g., the ms program [START_REF] Hudson | ms a program for generating samples under neutral models[END_REF]). With real genome-wide sequence data, the curve can be estimated with the PSMC method of [START_REF] Li | Inference of human population history from individual whole-genome sequences[END_REF]. We then try to identify a piecewise stationary n-island model that generates an IICR that is identical or similar to the target IICR (or PSMC curve). The similarity between the two IICR curves is quantified with a distance metric defined below. We use a genetic algorithm to explore the parameter space (number of populations, migration rate within a time component, and timing of these changes assuming a fixed number of components for each independent analysis) and minimize that distance. We compute the IICR under the non-stationary structured coalescent (NSSC) of Rodríguez et al. (2018).

The piecewise stationary n-island coalescent

The parameter space

We first define the parameter space, as this directly determines the family of demographic histories that we can explore and infer from. The piecewise stationarity refers to the fact that, although migration rate is constant and identical between any pair of islands, this rate may be different between consecutive time periods (components), and there is a fixed number γ that represents the number of demographic events. To say that there are γ 0 changes of gene flow thus means that there are c = γ + 1 components or periods of constant gene flow. Likewise, the deme size, which is the same for all islands, may in theory change through time in the general model presented in Rodríguez et al. (2018). In the present study we focus on models with constant population size but we present a more general model where deme sizes can change between components. In this more general case, the parameter space includes the number of islands n, the times t i for the demographic events, and the values of both the migration rate M i and the local deme size s i at each new demographic period. Note that n is inferred but it does not change through time. We thus assume no extinction, no population split and no creation of new populations.

Given a fixed integer γ of demographic events to consider (γ 0) and a collection of bounds B in the form of:

B = [n min ; n max ], [t 1min ; t 1max ] . . . [t γ min ; t γ max ], [M 0min ; M 0max ] . . . [M γ min ; M γ max ], [s 0min ; s 0max ] . . . [s γ min ; s γ max ] ,
(2.1.1)

we define the parameter space Φ γ,B as:

Φ γ,B = ϕ = (n, t 1 . . . t γ , M 0 . . . M γ , s 0 . . . s γ ) ∈ N × R 3γ+3 , s.t. ∀i : 2 n min n n max ; 0 < t i min t i t i max ; 0 < M i min M i M i max ; 0 < s i min s i s i max . (2.1.2)
We define bounds for each variable because we use a constrained optimization algorithm, for which all parameters must be bounded (see §2.2). Also, since we focus on the case where there are no deme size changes, we enforce this by using B, as making s i min = s i max = 1 for all 0 i γ effectively fixes all deme sizes to 1, and reduces the number of parameters to 2γ + 2.

Computing and scaling the IICR

Given any demographic scenario from Φ γ,B , the associated coalescent process is an instance of the NSSC of Rodríguez et al. (2018). Our main object of interest regarding these scenarios is the IICR. In the supplementary materials we provide a brief overview of how to perform its computation for any given ϕ ∈ Φ γ,B based on the work of Rodríguez et al. (2018).

The computation of the IICR uses functions that receive the time t in units of 2N generations, and return values in units of N generations per coalescence, so these IICR functions are dimensionless in the sense that they operate in a relative frame of reference.

In order to compare the IICR with PSMC inferences, we need to re-scale both the time and the IICR values by a reference deme of size N which specifies how many haploid genes correspond to a local deme size of 1 as follows:

sIICR(g) = N IICR(g/2N);
where sIICR(g) refers to the scaled IICR, and IICR(t) to the unscaled (dimensionless) one. Note that we use g for generations as the variable name for sIICR to further stress the difference. The parameter space for the sIICR can be thought of as a simple one-dimensional addition to Φ γ,B :

Φγ,B = φ = (N, ϕ) ∈ R × Φ γ,B such that 0 < N min N N max .
Considering the c components which constitute the model, we note that during component i , taking 0 i γ = c -1, the underlying coalescent process X t is being governed by an deme (initial state 1) has distribution F same (t) = P (t) (1,3) and density f same (t) = D(t) (1,3) , and for the case where we sample in two different demes (initial state 2), then T 2,diff. would have its distribution given by F diff. (t) = P (t) (2,3) and its density by

f diff. (t) = D(t) (2,3) .
The factor matrices in (2.1.3) can be computed in several ways in the general case (see [START_REF] Herbots | Stochastic models in population genetics: genealogy and genetic differentiation in structured populations[END_REF] or [START_REF] Hobolth | Phase-type distributions in population genetics[END_REF]), but considering this particular instance of size 3 × 3, they may also be computed directly given an arbitrary ∆t and rate matrix Q: (2.1.4) where:

e ∆t Q =   F 11 F 12 1 -F 11 -F 12 F 21 F 22 1 -F 21 -F 22 0 0 1   ,
F 11 = (δ + α -2Ms)exp 2 + (δ -α + 2Ms)exp 1 2δ , F 12 = (n -1)F 21 , F 21 = Ms(exp 1 -exp 2 ) δ , F 22 = (δ + α -2Ms)exp 1 + (δ -α + 2Ms)exp 2 2δ , exp 1 = exp ∆t(δ -α) 2s(n -1) , exp 2 = exp ∆t(-δ -α) 2s(n -1) , δ = α 2 -4Ms(n -1), α = Mns + n -1.
(2.1.5) With this we can efficiently compute either IICR functions:

IICR same (t) = 1 -F same (t) f same (t) , IICR diff. (t) = 1 -F diff. (t) f diff. (t) . (2.1.6)
In contrast with the parameter space for the unscaled IICR (Φ γ,B ) in which there is a one-to-one correspondence between a parameter tuple and the corresponding IICR curve, there are only 3γ + 3 independent degrees of freedom in Φγ,B , even though there are 3γ + 4 parameters. This notion is formalized by the following lemma.

Lemma 1. Given any φ = (N, n, t 1 . . . t γ , M 0 . . . M γ , s 0 . . . s γ ) ∈ Φγ,B , then the parameter tuple:

φ0 = N C , n, Ct 1 . . . Ct γ , M 0 C • • • M γ C , Cs 0 . . . Cs γ is such that: sIICR φ(g) = sIICR φ0 (g) ∀g > 0,
where C is any rescaling factor for which the coordinates of φ0 are within the bounds B.

Proof. Let us denote by π n,M,s (t) the factors e tQ that appear in the definition of P (t) (equations 2.1.3 and 2.1.4). It is easy to verify that for any C > 0, we have:

π n,M,s (t) = π n, M C ,Cs (Ct). (2.1.7)
Indeed, from the expressions in (2.1.5) we can see that the parameter M always appears in the factor Ms, and the parameter t always appears in the factor t/s, which are invariant under the transformation (M, s, t) → ( M C , Cs, Ct).

Next, given any:

ϕ = (n, t 1 . . . t γ , M 0 . . . M γ , s 0 . . . s γ ), ϕ 0 = (n, Ct 1 . . . Ct γ , M 0 C • • • M γ C , Cs 0 . . . Cs γ ),
and any t > 0, we can write P (t) from (2.1.3) as:

P ϕ (t) = i k=1 π n,M k-1 ,s k-1 (t k -t k-1 ) π n,M i ,s i (t -t i ) = i k=1 π n, M k-1 C ,Cs k-1 (Ct k -Ct k-1 ) π n, M i C ,Cs i (Ct -Ct i ) = P ϕ 0 (Ct). (2.1.8)
where i is the largest index such that t i t and subsequently Ct i Ct. We denote now by

F ϕ (t) any of F same (t) = P ϕ (t) (1,3) or F diff. (t) = P ϕ (t) (1,2) . From (2.1.8) we have F ϕ (t) = F ϕ 0 (Ct).
In order to introduce scaling, we consider an arbitrary effective size N and perform the corresponding change of variable t = g/2N. Thus, by having φ = (N, ϕ) and φ0 = (N/C, ϕ 0 ), we get:

F ϕ g 2N = F ϕ 0 Cg 2N ⇒ 1 2N f ϕ g 2N = C 2N f ϕ 0 Cg 2N ⇒ N 1 -F ϕ (g/2N) f ϕ (g/2N) = N C 1 -F ϕ 0 (Cg/2N) f ϕ 0 (Cg/2N) ⇒ sIICR φ(g) = sIICR φ0 (g)
The implication of Lemma 1 is that when trying to infer all the parameters of φ simultaneously, the only parameter for which we may get an absolute estimate is n, as the rest of them can only be distinguished up to an unknown re-scaling factor C. Note that this un-identifiability issue is different from the one identified in [START_REF] Mazet | On the importance of being structured: instantaneous coalescence rates and human evolution-lessons for ancestral population size inference?[END_REF] regarding the inability to discriminate between panmictic and non-panmictic demographies with a single IICR. However, we stress here that in practice this is not necessarily an issue because it suffices to fix one of the model parameters (for instance, s 0 = 1) to be able to uniquely map any sIICR curve to its parameters. In the case of constant size this is even less of an issue, since all deme sizes are fixed to s i = 1 and thus no further considerations are necessary.

Optimization framework: search algorithm and optimality criteria

The search algorithm explores the parameter space and uses an optimality criterion to select the structured scenario that best explains a given target IICR curve, either scaled or unscaled.

We also assume that the underlying coalescence times for these target IICRs have cumulative distribution F 0 and density f 0 . Given a target IICR 0 and a parameter space Φ γ,B , we want to find a parameter tuple ϕ in Φ γ,B such that the exact IICR curve corresponding to the model defined by ϕ (denoted by IICR ϕ ) approximates IICR 0 . We thus want to find the minimal distance:

min ϕ∈Φ γ,B d IICR 0 , IICR ϕ .
(2.2.1)

Regarding the distance d, a straightforward definition would be:

d IICR 0 , IICR ϕ = ∞ 0 IICR 0 (t) -IICR ϕ (t) w (t) dt, (2.2.2)
where w (t) is a weight function that should take into account the natural distribution of the information in an IICR. One reasonable solution for w is to take a quantity proportional to the density f 0 of the coalescence times because it ensures that the integral in (2.2.2) is finite, and also because it assigns more weight to the temporal periods where the target IICR 0 is expected to be more accurate and reliable since more coalescences are likely to have happened.

We thus consider the family of weight functions:

w (t) = f ω 0 (t) f ω 0 1 , (2.2.3)
where • 1 is the L 1 -norm and ω > 0 is a weight-shifting parameter, with the purpose of dampening (if ω < 1) or exaggerating (if ω > 1) the effect of the weight f 0 . Unless otherwise noted, we take ω = 1, which corresponds in practice to giving more weight to recent periods of the IICR in direct proportion to the density f 0 in an n-island model.

In practice, we need to consider that all we know about IICR 0 is a stepwise discretization over a bounded interval of time, so a numerical approximation of the distance (2.2.2) is required. This includes approximating the density f 0 of the underlying T 2 distribution. Given a division of time into I intervals [τ j-1 ; τ j ) for 1 j I, where τ 0 = 0 and τ I < ∞, we can consider a discrete representation of IICR 0 in the form of a collection of I values {y j } such that:

IICR 0 (t) = y j , ∀t ∈ [τ j-1 ; τ j ), 1 j I. (2.2.4)
We can use this form to compute a numerical approximation for the integral in (2.2.2). For instance, a first degree approximation would yield:

d IICR 0 , IICR ϕ = I j=1 y j -IICR ϕ (τ j ) w (τ j-1 ) (τ j -τ j-1 ).
As for the values of w (τ j ), notice that from (1.4.4) we have the identity:

f 0 (t) = exp t 0 -dτ IICR 0 (τ ) IICR 0 (t),
which, using the representation (2.2.4), can be discretized into:

f 0 (0) = 1/y 0 , f 0 (τ j ) = exp j k=1 τ k-1 -τ k y k y j , 1 j I.
We then have that for any given ω, w (τ j ) can be expressed as:

w (τ j ) = f ω 0 (τ j ) I k=1 f ω 0 (τ k ).
An alternative option for the definition of d in (2.2.1) could be one that takes into account the ultimate visual nature of the curve fitting task. Assuming that the points {τ j } are log-distributed and that they will be used for visualization purposes in a horizontally log-scaled plot like Figure 2.11, then the definition:

d(IICR 0 , IICR ϕ ) = I j=1 y j -IICR ϕ (τ j ) ,
(2.2.5) captures the perceived visual difference between the plots of the two curves. We distinguish distance (2.2.2) from (2.2.5) by denoting them d ω and d visual respectively. We keep both definitions because we found that the weighted family of distances generally performs better than the visual distance under certain validation tests, but also that the d visual distance can be used to choose the optimal weight parameter in d ω (see Figure 2.24).

Regarding the optimization problem (2.2.1) itself, we use the Differential Evolution method [START_REF] Storn | Differential evolution-a simple and efficient heuristic for global optimization over continuous spaces[END_REF]. As a genetic meta-heuristics, this algorithm maintains and evolves (using mutation and recombination parameters) a population of solutions iteratively. As a global optimization algorithm, it features mechanisms for escaping local optima of the parameter space. In §2.2.1 we explore the potential effects on the inference results of tuning some of the parameters provided by this algorithms implementation. For our validations, the method runs by using multiple steps of optimization which we refer to as rounds. In addition, we stress that the method should be used multiple times on real data sets. We set a maximum number of allowed rounds, as well as a tolerance ε for the distance which controls the minimum number of rounds.

Comparison of optimization parameters

Here we explore the effect of various parameters of the optimization algorithm on the speed of convergence of the inference process.

The most important parameters that affect the search process and convergence criteria are:

• strategy: it can be one of 12 possible values (see Table 1) and controls how each new generation of solutions are computed from the previous one. Default is 'best1bin'.

• maxiter: maximum number of iterations the algorithm will perform before forcing convergence. Default is 5000.

• popsize: number of simultaneous solutions during any given generation. Default is 15.

• tol: relative tolerance for search convergence. The convergence criteria is met when the standard deviation of the solutions within a generation is smaller than tol times the average energy (in our case, distance) within that generation. Default is 10 -2 .

• mutation: per-generation mutation rate of the solutions. The default behaviour is to draw a random value from a uniform distribution in [0.5, 1] each generation.

• recombination: per-generation recombination rate of the solutions. Default is 0.7.

We selected 10 random simulated scenarios with unscaled IICR and c = 4 components from the set of exact IICR validations that had not converged before 500 rounds (these would correspond to scenarios off the diagonals in Figure A.4). For each one of them, and for each of the 12 possible values for the strategy parameter, we attempted another 100 rounds. Table 1 shows the number of rounds it took for each of the 10 scenarios to converge in each case (a value of 100 means that convergence was not reached). As can be seen, strategy 'best2exp' was best with a combined number of 424 rounds for the 10 scenarios.

Strategy #1 #2 #3 #4 #5 #6 #7 #8 #9 #
Next, using this optimal strategy parameter, we tried one alternative value for the rest of the optimization parameters at a time, again allowing a maximum of 100 rounds. The alternative values were as follows:

• maxiter was changed from 5000 to 10000.

• popsize was changed from 15 to 50.

• mutation was changed from random sample in [0.5, 1] to random sample in [0.5, 1.7].

• recombination was changed from 0.7 to 0.9. The results, shown in

Validation framework and results

We applied the SNIF inferential method to target IICRs generated under piecewise stationary n-island models of increasing complexity (i.e., number of components) and with known parameter values (N, n, t i , M i ) and then compared the inferred parameter values to those actually used.

In what follows we present various ways of generating random demographic scenarios and then computing appropriate IICR curves from them for use in validation.

Φγ,B

Continuous sampling

Discrete sampling

φ (N,n,t i ,M i ) Exact IICR {x i } ∼ T2 110010011• • • 001101101• • • T-sim IICR Seq-sim IICR SNIF φ0 (N 0 ,n 0 ,t 0 i ,M 0 i ) Compare! Figure 2
.2: Flowchart of the validation procedures. Starting from a parameter space Φγ,B we use one of two sampling methods ( §2.3.1) to generate a demographic history φ defined (for the scaled case) by the parameters (N, n, t i , M i ). We then compute the IICR of that demographic history using one of three methods ( §2.3.2) to obtain the target IICR. After that, we run the inference algorithm on this target IICR curve (using wider bounds than those in B) to obtain an estimated (or inferred) demographic history φ0 = (N 0 , n 0 , t 0 i , M 0 i ), which we then compare to the known φ in order to assess the accuracy of the inference methodology ( §2.3.4 and §2.3.5).

Sampling the parameter space

Given a parameter space Φ γ,B (we only discuss the unscaled case here for brevity, but the same principles apply to a scaled parameter space), we sample demographic scenarios from which we compute the corresponding IICRs. We used two sampling strategies which we call continuous and discrete sampling.

Continuous sampling. Assuming that we want to realise L independent tests, this sampling strategy consists in using uniform or log-uniform distributions for each of the 3γ + 3 random variables:

n ∼ U{n min , n min + 1, . . . , n max }, t 1 ∼ LU 10 [t 1min ; t 1max ], . . . , t γ ∼ LU 10 [t γ min ; t γ max ], M 0 ∼ U[M 0min ; M 0max ], . . . , M γ ∼ U[M γ min ; M γ max ], s 0 ∼ U[s 0min ; s 0max ], . . . , s γ ∼ U[s γ min ; M γ max ],
(2.3.1)

where U denotes a uniform distribution (discrete in the case of n and continuous for the rest) and LU 10 denotes a log-uniform distribution of base 10. This distribution is used for sampling the times of changes in a logarithmic space in order to take into account the natural distribution of information in an IICR. After sorting the times t i , we can define the L sampled scenarios by constructing, for 1 j L, the tuple n j , t j 1 . . . t j γ , M j 0 . . . M j γ , s j 0 . . . s j γ . This sampling strategy makes it very unlikely to sample exactly the same parameter values twice or to sample exactly the same M i values in two consecutive components. However, it sometimes produces demographic scenarios in which consecutive t i and/or M i values may be close to each other, and thus difficult to distinguish. This makes it thus harder on our inferential framework compared to cases where we would chose contrasted scenarios with clearly separated events with major changes in migration rates. In other words, our inferential method was sometimes inferring parameters in the case of extremely difficult scenarios as we show below.

In §2.3.4 we show the results obtained using this sampling method with L = 400 scenarios. The bounds for sampling and inferring are shown in Table 2.3. In particular, we note that in practice we disallow deme size changes by fixing the bounds of the s i to 1, which consequently reduces the parameter space to just 2γ + 2 parameters.

Discrete sampling.

Here we sampled L = 100 independent scenarios from the same parameter space, but using the following set of predefined values:

n ∈ n = {2, 5, 10, 15, 20}, t i ∈ t = {1/10, 1/2, 1, 2, 5, 10, 20, 50} ∀i , M i ∈ M = {1/10, 1/5, 1/2, 1, 2, 5, 10, 20, 50} ∀i ,

s i ∈ s = {1} ∀i ;
The inference process was, however, done within the continuous space. For instance, under this validation scheme (see §2.3.5) we only simulated data with 2, 5, 10, 15, 20 islands but the inference process always allowed n to take any value between 2 and 50. The choice of the L independent simulated data sets was done using the following procedure. We first considered the following cartesian product of dimension 2γ + 2:

K = n × t γ × M γ+1 × {1}.
and then uniformly drew L tuples from K without replacement. We then sorted the sampled event times obtaining thus a set of L demographic scenarios. We drew randomly (without replacement) from the set K, rejecting scenarios with identical M i values in two consecutive components, until we reached L accepted scenarios.

The three types of target IICRs

We explored three different types of target IICRs (see Figure 2.3) that could be obtained given a scenario ϕ ∈ Φ γ,B . All IICRs were discretized so as to be comparable to PSMC plots (see equation 2.2.4).

Validating SNIF on PSMC plots across the parameter space described above would be extremely time consuming as it would require simulating genomes and then running the PSMC method (or other related methods) on these genomes before applying our approach. We thus only ran the PSMC method in the case of the scenarios inferred for the human data so as to integrate the uncertainty due to the PSMC inferential process. The issue of uncertainty is crucial but our aim is not to test the PSMC or other inferential methods. To clarify this we explain below the different types of IICR that could be computed given a scenario ϕ ∈ Φ γ,B .

Exact IICR We can compute the IICR for any n-island model at any time value t, but to generate input data we need a discretization as in (2.2.4), so considering that we take a log-distributed sample of size I in the interval [t min , t max ], we end up with a suitable IICR 0 . Note that even though this IICR has been discretized, its values are exact within machine precision, so it is still an artificial product compared to real data.

For the validations using the exact IICR in §2.3.4 we chose for the distance tolerance between a target and an inferred IICR a value of ε = 10 -10 for the unscaled IICRs and an equivalent value of ε = 10 -7 for the scaled IICR (since the simulated N was always 1000). It should be noted that this value of ε is quite small even for double-precision floating-point arithmetic, and thus is only a reasonable choice for validation using exact IICRs (i.e., those where the distance could theoretically be zero).

T-sim IICR

The T-sim IICR is obtained by simulating a finite collection of T 2 realizations using ms and then building an empirical IICR as in [START_REF] Mazet | On the importance of being structured: instantaneous coalescence rates and human evolution-lessons for ancestral population size inference?[END_REF], using the Kaplan-Meier estimator [START_REF] Kaplan | Nonparametric estimation from incomplete observations[END_REF], with log-distributed times. We stress that ms scales time in units of 4N generations whereas our models use a scale of 2N generations (see the example in Figure 2.3), so this must be kept in mind when writing ms commands.

Seq-sim IICR

We simulate genomic sequences with ms and then apply the PSMC method for obtaining the IICR 0 to be used by the inference method. Since simulating genomes and performing PSMC analyses is significantly more time consuming than the other two methods, we limited ourselves to validating the Seq-sim IICR for the human PSMC based scenarios that we obtain after performing the demographic inference described in §2.4. The results of this step are shown in §2.4.2 of the Supplementary Materials.

Quantifying the inference error

In order to quantify the inference error incurred during the continuous-sampling validation phase, we measured the normalised root-mean-square deviation (nRMSD) between the simulated and inferred parameter values for each parameter of the demographic models (c = 1 to c = 6 components). These values can be seen in the lower-right corner of every sub-panel in figures A.1 to A.6 and A.8 to A.13. They are also summarised in Figure 2.4. We note that the number of islands n and the effective size N (in scaled scenarios) is very well inferred regardless of the number of components c. On the other hand, the inference accuracy of the connectivity rates M i does get gradually worse when increasing the number of parameters (see section Validation results using exact target IICRs in the main text). In an effort to quantify the component misidentification phenomenon, we computed the correlation between 100 randomly sampled pairwise values of simulated and estimated parameters (in this case, the simulated M 2 parameter in unscaled scenarios of c = 5 components with the inferred M j parameters of the same scenarios). The key insight underlying this test is that when a parameter is badly estimated, it may be due to the fact that the method estimated the value from the component that is either just before or just after it. Conversely, the value of M from other components should be non-correlated, since values are taken at random within the range of allowed values. Now, we do not know if the method used the component just before or just after, as this may change from simulation to simulation. To solve this problem we computed the correlation (r 2 ) by using either both values or just the value that was closest. We did exactly the same for the values taken from other components that are not expected to be correlated. We find indeed that when no correlation is expected there is no correlation, but when we take the best value the correlation increases, and it increases much more when it is a neighbouring value. Additionally, we found that this effect is significantly amplified (r 2 increased from 0.37 to 0.82) when we exclude from our sample the tests were the inferred rates were within 10% of the simulated values, further indicating that this effect is present mostly when there is a large mismatch between the simulated and inferred M values. These results are displayed in Figure 2 In 'random global' we draw M j from components 0, 1, 3, and 4. In 'random far' we draw M j from components 0 and 4. In 'best far' we draw M j from components 0 and 4 again, but always keep the best fitting value of the two. Finally, in 'best near' we draw M j from components 1 and 3 and always keep the best fitting value. 'Filtered at 10%' indicates that the scenarios where the inferred value was within 10% of the simulated value were excluded.

0 0.1 0.2 Normalized RMSD 0 0.5 1 M 0 M 1 M 2 M 3 M 4 M 5 0 1 2 3 1 2 3 4 5 6 0 0.1 0.2 c n N 1 2 3 4 5 6 0 0.5 1 c 2 3 4 5 6 0 1 2 3 c t 1 t 2 t 3 t 4 t 5

Validation results using exact target IICRs

A first set of figures (Figures 2.6, A.1 to A.6, and A.8 to A.13) represents the simulated and inferred parameter values on the horizontal and vertical axes, respectively, using the continuous sampling strategy. As explained in §2.3.1, the range of possible values in the inference process was always wider than the range used for the simulated values (see Table 2.3 for the exact values). We quantified the inference error for each parameter by computing the Normalized Root-Mean-Square Deviation (nRMSD). This value is displayed in the lower-right corner of each panel of the previously mentioned figures, and summarized for all parameters in Figure 2.4. For reference, we also highlight the y = x line, indicating what would be a perfect inference, and the region corresponding to 10% of relative error around this line (50% of relative error in the case of the t i parameters). A summary of how many tests fall within this margin of error (and others) is displayed in figures A.7 and A.14. Altogether, these figures always show the data points near the y = x, hence suggesting that the inferred parameter was identical or very close to the simulated parameters. This is particularly obvious for all the parameters corresponding to scenarios with up to four components, where the nRMSD values stayed below 0.5 (the case of the t i parameters is exceptional, since the exponential distribution of its range disproportionately affects the error measures). For instance, Figure 2.6a shows the results for a model with three components, in which there is a nearly perfect match (nRMSD close to or below 0.1 for the non-t parameters) between the simulated and inferred values for the model parameters. For five-and six-component scenarios the results are still nearly perfect for most of the simulated scenarios but we observe an increasing number of cases (i.e., simulated scenarios) where the parameters are poorly estimated, with of scenarios with 500 rounds increased with the number of components. We found indeed that the proportion of simulations for which that maximum was reached increased with the number of components. For instance, all five-and six-component scenarios stopped their parameter search at 500 rounds, hence suggesting that at least some had not yet reached an optimum solution. For the cases with one-and two-component scenarios, all 800 independent simulations reached convergence in less than 150 rounds (see Figure 2.7). Again, the choice of the tolerance ε plays a role in these results, and selecting larger tolerances will tend to produce earlier convergence in general, but not necessarily better results. ) that the method required before converging to a scenario with a distance value smaller than the tolerance ε = 10 -10 for the unscaled case. The bottom panel shows the same information for scaled scenarios and a tolerance of ε = 10 -7 . We represent in different colors the curves corresponding to scenarios with different simulated (and inferred) number of components, ranging from c = 2 up to c = 6; the same data-set represented in figures A.1 to A.6 (top panel) and figures A.8 to A.13 (bottom panel). Higher number of rounds is not correlated with worse fit when the maximum number of rounds is not reached; and when it is reached it is not necessarily an indication of bad fit either, although all instances of bad fit stem from inferences that reached the maximum allowed number of rounds.

As a test we randomly identified a couple of scenarios with six components that had bad estimates and re-ran the algorithm with 5000 rounds. We found that the distance value consistently decreases with more rounds (see Figure 2.8), but the inferred parameter values may not converge to the simulated ones because with more components there is a higher probability that two consecutive simulated M i values are very close, thus making the corresponding event time challenging to infer.

The second possible reason for the poorly estimated parameters in Figure 2.6 may be related Inference round Optimal distance reached for scenario # 217

Figure 2.8: Convergence of the inferred IICR with the progression of the inference rounds. We measured the best distance achieved (ω = 1) as a function of the completed rounds for up to 5000 rounds. The two curves correspond to two c = 6 component scenarios that did not previously converge (ε = 10 -10 ) in 500 rounds (marked with a vertical dashed line). We note that there is a clear trend for the distance between the IICR curves to decrease with more inference rounds, but there is also diminishing returns to performing a very large number of rounds, since the issue of component misidentification may be an insurmountable difficulty for some scenarios (see main text for a more detailed discussion).

to the fact that some simulated components may have a short duration that do not leave a significant mark on the IICR curve, thus leading them to be "skipped". We refer to this issue as component misidentification, which could lead to a particular estimated parameter to be plotted in the wrong panel. For instance, the method may miss the first change in migration rate at t 1 and identify the second change in migration at t 2 . In such a case the method will assign the inferred t 2 value to the set of inferred t 1 values and plot it in the t 1 panel. This wrongly assigned t 2 value will thus appear away from the diagonal in the t 1 panel even if it was well-estimated. Such misassignment cases for one parameter will also have consequences for the M i plots, and thus will generate several misassignments across panels. They are also expected to increase in frequency as the number of components increases and as the t i values become closer to each other. This phenomenon can be observed clearly in the right panels of Figure 2.9. We also present an attempt at quantifying it for the case of c = 5 in Figure 2.5. One way to mitigate the effect of this misassignment issue in the analysis of the results is to visualize the simulated and inferred scenarios using what we call a connectivity graph. This connectivity graph represents the times at which migration changes against the values of the migration rates. Such connectivity graphs are featured in the next section.

Validation results using T-sim IICRs

The connectivity graphs and IICR plots obtained from simulated T 2 values show that again the scenarios are generally very well reconstructed (Figure 2.9 and Figures B.1 to B.15). In Figure 2.9 the connectivity graphs obtained for all the scenarios simulated with three and four components show that the inferred times at which migration rates changed (green vertical lines) and the inferred migration rates (green horizontal lines) are generally overlapping close to the simulated values (dotted vertical and horizontal gray lines). In the right panels of this figure, we show a subset of the inferred migration histogram (in red). Namely, we show the distributions of the migration values that were inferred for components with a simulated migration value of M i = 10 for panel (a) and M i = 1 for panel (b). This allows us to better appreciate the variance of the inferred migration values in the context of the simulated ones, as well as the component misassignment effect mentioned earlier. Indeed, we note here that the incorrectly inferred migration values are clustered around other simulated values, indicating a mismatch in a particular component assignment which does not affect the rest of the inferred demographic history (we present a quantification of this effect for a limited case in Figure 2.5).

For example, consider the right sub-panel of (a). We see that most repetitions correctly inferred a value close to M = 10 for the components with that simulated migration rate. However, there were cases where a given component i was simulated with a migration rate of M i = 10, but it was missed entirely (maybe because it did not generate a very different IICR or because it had a short duration), and thus the inferred migration value for component i ultimately reflected either M i-1 or M i+1 . In panel (b) we can observe the same effect with higher intensity because with more components it is more likely for them to be misaligned or misidentified during inference. See Figure 2.5 for a quantification of this effect on scenarios of c = 5 components.

These connectivity graphs (and the one obtained for five components shown in Figure B.13) also show that there are regions of the parameter space where the green lines are more widely distributed. For instance, in the recent past of Figure 2.9b (t i < 10 -3 generations) when the simulated M i value was 0.1 or 0.2 the inferred values seem to vary between 0.05 and 0.3, suggesting that the method identifies periods with low migration rates but that the exact value is difficult to estimate properly, at least in the recent times. These graphs however summarize extremely different scenarios, including scenarios in which consecutive M i values may be similar. We thus stress that the quality of the inference is dependent on the timing of the changes in migration rates and on the size of the change in M i values.

Figure 2.10 shows the results for four different scenarios. In each of the four panels, we represented the inferred and target IICR plots, connectivity graphs, N (the size of each the islands) and n (the number of islands) for the corresponding model. Panels (a) and (b) correspond to three-and four-component scenarios, whereas panels (c) and (d) show the results for two five-component scenarios, one for which we obtained very good estimates and one for which the estimates were poorer. In panels (b) and (c) the inferred and simulated M i and t i values are on top of each other as can be seen in the connectivity graphs. Similarly, N and n are also well estimated. Here, the IICR plots also overlap, although this does not always running every inference a second time, ignoring this period of possible recent expansion. This is accomplished using an option that allows to limit the interval where the distance function is computed. In this case, we restricted both this range and the bounds for the t i to be between 50 thousand and 10 million years ago, thus ignoring any population size change that may have happened in the last 50,000 years. Note that this option is also useful to ignore very ancient sections of the PSMC plots which may be difficult to trust.

Since real human PSMCs are unlikely to have been generated by an n-island model, the default value of ω used for simulated data may not be the most appropriate, and we thus performed inferences with ω ∈ {1, 0.5, 0.2}. Decreasing values of ω give increasing weight to the most ancient part of the PSMC (see the weighted distances (2.2.3)). The resulting inferred demographic scenarios are shown in section §2.4.1.

To validate the inference process using PSMC outputs, we generated 10 Seq-sim IICRs corresponding to the inferred demographic scenarios for the French, Karitiana and Yoruba individuals. We exclude the Dai and Sardinian populations from this analysis because the their corresponding inferred histories are similar to the other three (see Figures 2.18 through 2.24). For each one of the selected scenarios we simulated nreps = 30 chromosomes of length L = 10 8 base pairs, using the effective size N inferred by SNIF, a per-base per-generation mutation rate of µ = 1.25 × 10 -8 (see [START_REF] Scally | Revising the human mutation rate: implications for understanding human evolution[END_REF] and references therein). We kept for consistency the scaled recombination rate of ρ = θ/5 as in [START_REF] Li | Inference of human population history from individual whole-genome sequences[END_REF], and we ran the ms command with θ = 4µLN using: ms 2 nreps -t θ -r ρ L -p 8 -I . . . where the rest of the command follows according to the inferred demography (see Figure 2.3 for a reference). After that we prepared a .psmcfa file as input for PSMC, always choosing a bin size of s = 100. Then we ran the PSMC with the command: psmc -N25 -t15 -r5 -p "4+25*2+4+6" . . . following [START_REF] Li | Inference of human population history from individual whole-genome sequences[END_REF] on human data. We then scaled the resulting curve using the information in the generated .psmc file and used these PSMC curves as targets to determine whether we could indeed infer the parameters used for such complex scenarios.

We also applied SNIF to genomic data simulated under the scenarios used to describe recent human evolutionary history by [START_REF] Gutenkunst | Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data[END_REF] and [START_REF] Noskova | GADMA: Genetic algorithm for inferring demographic history of multiple populations from allele frequency spectrum data[END_REF]. Here, we thus ask the following two questions: if human evolution were indeed closer to such splitting models, would our method infer again an n-island model with similar parameters to those inferred from the humans PSMCs? additionally, do these models generate IICR plots that are similar to the human PSMCs? The results of these validations are shown in §2.4.2 of the Supplementary Materials.

Results

Figures 2.12 and 2.13 show the results of using SNIF on the human data.

fit to the data. The model stipulates the existence of an ancestral population that experienced an increase in size around 275 thousand years ago (Kya), and the a splitting event at about 150 Kya. This split resulted in two populations that exchanged gene flow asymmetrically: a large one that eventually became the modern African population, and a smaller one ancestral to the modern Eurasian population, whatever this terminology may mean. This ancestral lineage split about 22 Kya into the precursors of the European and Asian populations, which at this point began an exponential increase in size that continued to present day. During this period, all three lineages continued to exchange gene flow asymmetrically. The times for these resize and splitting events are represented as dotted vertical gray lines in Figure 2.14.

It is clear that the nature of this model does not lend itself to be exactly modelled by a symmetrical n-island model, but the piecewise-stationarity of our family of models should still be able to pick up the main demographic events. For example, from an n-island perspective, a merger or joining of two populations (going backwards in time) may be represented by an increase in gene flow, although this effect may be confounded by the actual changes in both the sizes of the populations and migration rates taking place during these events. Also of note is the fact that the first merger event is not visible to our method because it marks the start of the recent population expansion and is thus excluded from the distance computation.

As can be seen in panel (a) of Figure 2.14, these IICRs do not exhibit any major features past the 300 Kya mark, so they do not agree with the human PSMCs of Figure 2.11 (of which the representative ones are again shown in Figure 2.14 in dashed trace for reference), and they also do not generate significant events in the inferred demographic histories. Particularly, varying the value of the weight-shifting parameter ω did not make a great effect in this set of inferences (which is in contrast with the results shown in Figure 2.12). This inferred demographic history can be roughly summarized from panel (b) as having a period of relative high gene flow followed by a sharp decrease near the 300 Kya mark, which can be very clearly attributed to the size increase of the ancestral population in the C3PO model.

The inferred number of demes and their relatives sizes for each population can be observed in Figure 2.15. The numbers for the African population is in sharp contrast with the other two populations. We can also observe that for the three populations there is more variance (compared to the results from Figure 2.13) in the inferred values of n and N across the different values of c and ω. This may indicate a weaker link to an underlying n-island model.

In general, there is little agreement between the demographic histories inferred by our method from the PSMC data and the simulated IICRs from the C3PO model. This is expected because of how the two models have fundamentally incompatible structures, not only regarding the island versus tree aspect, but also due to the size changes in the C3PO model that affect the IICR potentially as much as gene flow does. However, we do identify the approximate timings of the two visible demographic events when using c = 5 components and the more recent-weighted value of ω = 1. These results also serve as additional validation that our method will not return the same parameter values regardless of the source of the data. They also suggest that the C3PO model is unlikely to be a good model to understand questions about ancient human structure and evolution.

Seq-Sim validations

Figure 2.17 shows the results of the validations using seq-sim IICRs. For the three chosen representative human populations (French, Karitiana and Yoruba), we selected the SNIF inferences obtained with the parameter values c = 5 components and ω = 0.2. The selection of these values as the preferred ones was supported by the fact that the d visual distance (see Figure 2.24) was best for them. For each of these inferred models we generated two independent genomic sequences of length 3 × 10 9 base pairs and applied the PSMC method to obtain two independently estimated target IICRs. These seq-sim IICRs are shown in blue on the left panels of Figure 2.17. The connectivity graphs associated with these inferred models are represented in the middle panels by the red curves. The corresponding inferred values of n and N are marked as the black reference circles in the right panels of the figure.

After obtaining the seq-sim IICRs, we applied again our inference method. The goal is to validate whether it would be able to consistently infer the same parameter values of the demographic history regardless of the origin of the source IICR curve. To this end we performed 10 independent inferences from each of the two target IICRs. The inferred IICRs are superimposed on the left panels of Figure 2.17 (transparent red curves, 20 per population). The inferred connectivity graphs are shown in the middle panels (transparent green curves) and the inferred values of n and N are presented in the right panels. We observe an agreement with the previously inferred histories, which suggests that if the real history of human evolution were close to piecewise n-island models like those used in this work, our method would be able to infer the parameters properly, and they would be similar to those maintained constant. The method can in principle infer changes in population size but this should be part of an extension of the present work. At this stage we stress that more work is needed before changes in connectivity and population size can be estimated together with confidence.

The intended use case for our method consists of several inference cycles, each preceded by adjustments to the many available parameters which include: the bounds B for the estimated variables, the number of components c, the weight distance parameter ω, the time interval where the distance is to be computed, the number of inference rounds along with the tolerance ε, some parameters of the search algorithm, and other minor options. This cycle emerges naturally from the fact that the inference process itself runs fairly quickly (a few seconds per round), so it is feasible to prepare a script that generates inferences under several combinations of these parameters and later do a general assessment of what are the most plausible scenarios for the data depending on the visual fit, the consistency of the inferred demographic histories, the distribution of the distances, etc. The SNIF software already includes a number of auxiliary scripts that may be used in this later analysis stage, including for instance the automatic generation of figures similar to the sub-panels of figures 2.12 and 2.13 in the main text, where the target and inferred IICRs can be compared, and the nature of the best fitting scenarios can be understood using the connectivity graphs or the n and N plots.

As for the number of components c, specifying a higher value will generally (but not always) result in a better fit and a lower distance, however it also incurs in longer analysis times and diminishing returns on the new information present in the resulting inferred histories. These aspects can be balanced by increasing the value of c up to the point where the inferred demographic histories start converging (for instance, as more components are added, the connectivity graphs stop revealing new major events and any additional new degrees of freedom are only used to refine already existing details).

As part of the user inference cycles mentioned above, the bounds system can be configured to infer parameters for more strict or specialized symmetrical island models. For instance, it is possible to fix the number of inferred islands to be exactly 3 by setting n min = n max = 3 in the inference bounds B. Likewise, given an independent estimate of the reference size N for the data, this information may be used to process the corresponding scaled target IICR as an unscaled one using matching inference bounds for the parameter N.

An important point concerning scaling is that the supplied value of the mutation rate µ must be accurately specified when inferring demographic parameters from a PSMC curve. This value is not used during the inference, but it is used in order to properly scale the PSMC curve and convert it into a target IICR, therefore it is important to provide the same value of µ that was used during the PSMC analysis. Otherwise, the only parameter that can be correctly inferred is n, since the rest of the parameters would only be accurate up to a scaling factor. This follows from a similar logic to that of Lemma 1.

We also suggest to use the hand-fitting Python script developed in [START_REF] Chikhi | The IICR (inverse instantaneous coalescence rate) as a summary of genomic diversity: insights into demographic inference and model choice[END_REF] as a complement to the automated inference process proposed in this work. You may want to compare the results obtained this way with the output of SNIF. Doing this will help make sense of the data, or help set the bounds for the many SNIF parameters.

Chapter 3.

Multi-sample summary statistics for demographic inference

In this chapter we explore two different multi-sample summary statistics, with the goal of using them in a future demographic inference application, similar to that of Chapter 2 with the IICR summary statistic for a sample of size two haploids. The two summary statistics in question are the IICR k of an n-island model, and the expected site frequency spectrum (SFS) of a general structured model.

The IICR k of an n-island model is the most natural generalisation of the IICR. On the other hand, the SFS is a widely used statistic that succinctly summarizes genetic diversity of a demographic sample [START_REF] Wakeley | Estimating ancestral population parameters[END_REF][START_REF] Griffiths | The age of a mutation in a general coalescent tree[END_REF][START_REF] Nielsen | Estimation of population parameters and recombination rates from single nucleotide polymorphisms[END_REF][START_REF] Gutenkunst | Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data[END_REF][START_REF] Excoffier | Robust demographic inference from genomic and snp data[END_REF]. For both cases, the methods we present in this chapter are centered on building and analyzing a Markov process that captures all necessary information to compute the given summary statistic of k sampled haploid lineages in a structured model.

Similar Markovian approaches have been used previously in the context of the coalescent [START_REF] Hobolth | Phase-type distributions in population genetics[END_REF]. Notably, [START_REF] Herbots | Stochastic models in population genetics: genealogy and genetic differentiation in structured populations[END_REF] with the study of the structured coalescent; Rodríguez et al. (2018) with an extension to piece-wise stationarity models for studying the Inverse Instantaneous Coalescence Rate (IICR); and [START_REF] Kern | Exact calculation of the joint allele frequency spectrum for isolation with migration models[END_REF] with the exact computation of the expected Joint AFS for the two-populations isolation-with-migration (IM) demographic model.

The main challenge with the exact approach is that the Markov state space grows very quickly (almost exponentially) with the size of the sample, which poses a hard limit on applicability due to computational intractability. One of the aspects we explore in this chapter to mitigate this effect is the concept of model specialization, which is a way to compress the state space by exploiting symmetries in the studied models. The n-island model is best suitable for this kind of analysis due to its high degree of symmetry and mathematical simplicity.

The majority of this chapter sections are dedicated to the study of the SFS summary statistic. The SFS is a distribution ξ = (ξ i ) k-1 i=1 where ξ j indicates the proportion of observed segregating sites with exactly i derived alleles. The sample can be of any size (although large values of k quickly become computationally prohibitive) and the demographic structure can be specified by an arbitrary migration matrix and deme-sizes vector. In §3.2 we proceed to enumerate the state space for the continuous-time Markov process that will be studied. We represent the states as k × n matrices, where k is the haploid sample size, n is the number of demes, and an entry z with coordinates (i , j) indicates that, in the corresponding state, there are z lineages in deme j, each ancestral to i lineages from the sample (see Figure 3.1). This state representation captures the necessary information to exactly compute the expected SFS of the sample under any arbitrary initial sampling configuration π 0 . Indeed, our main result (presented in §3.3) expresses the expected SFS as a conceptually simple function of the data structures derived from our representation.

During §3.2, we present algorithms that will collectively allow us to enumerate the state space and computing the corresponding transition rates. Special attention is placed to the ordering of the states, since this allows the rate matrix to have a special block structure (Figure 3.2) that can be exploited to achieve a computationally efficient method for obtaining the expected times spend under each state. This method is developed in §3.3.

In §3.4 we present a specialization of the demographic model to the symmetrical n-island model. Figure 3.3 shows how, by taking advantage of the symmetrical nature of this model (and updating the proposed algorithms accordingly), the state space can be significantly compressed, thus making the use of larger sample sizes computationally feasible. The figure also shows a comparison of the size of the state space in [START_REF] Kern | Exact calculation of the joint allele frequency spectrum for isolation with migration models[END_REF] for computing the joint SFS of a two-island isolation-with-migration model.

The structured coalescent and the IICR k

In the classical structured coalescent, the evolution of the populations can be modeled with a continuous-time Markov chain in which the states are the vectors α = (α 1 , α 2 , ..., α n ) indicating how many lineages there are in every island.

Initially, the sum i α i is k, the total number of sampled lineages. The distribution of the α i varies across the islands as migrations and coalescences occur. Rodríguez et al. (2018) goes into the details of how to build such models for several examples of demographic structures. In all cases, once we have enumerated the states and computed the transition matrix Q, then the exponential function P (t) = exp(tQ) contains information about the state of the Markov chain at time t. Specifically, we know that exp(tQ)[i , j] = P (X(t) = α j |X(0) = α i ). These probabilities are useful for computing functions like the IICR, since the IICR can be obtained given the distribution of T 2 . A special case for the n-island model and a sample size of k = 2 was used in Chapter 2.

In the presence of a larger sample (k > 2) we may be interested in computing its associated IICR k function. We know from (1.3.4) how, in a size changing population, the relative size function λ(t) is related to the distribution of T k :

λ(t) = k 2 P (T k > t) f T k (t) . (3.1.1)
But just like in the case for the IICR, the right hand side of (3.1.1) is not constrained to the model of panmixia with size change, or to any other demographic model for that matter, and as such it can be computed exactly given the distribution of T k , and can be estimated from real demographic data using methods like MSMC [START_REF] Schiffels | Inferring human population size and separation history from multiple genome sequences[END_REF] and MAGIC [START_REF] Weissman | Minimal-assumption inference from population-genomic data[END_REF]. Given the similarity of this situation to the un-identifiability problem described in [START_REF] Mazet | On the importance of being structured: instantaneous coalescence rates and human evolution-lessons for ancestral population size inference?[END_REF], it is only fitting that we call this function the IICR k :

IICR k (t) = 1 -F t k (t) f T k (t) . (3.1.2)
It is important to stress that the distribution of T k (and therefore, the IICR k ) depends on the initial state of the Markov process, i.e., the sampling vector. let call this vector π 0 , and assume that it is a probability distribution over the possible initial states of the Markov process.

We can now pose the question of what are the possible states in a structured model with n island and k lineages: "In how many different ways can k lineages migrate and coalesce in n demes?" The answer depends on whether we can distinguish, or tell apart from each other, the individual lineages and/or demes. If we represent each state as an n-tuple α = (α 1 , α 2 , . . . , α n ), where α i denotes the set of lineages in island i , then we are distinguishing the islands by their position in the tuple. If additionally we distinguish the lineages by naming them ℓ 1 , ℓ 2 , . . . , ℓ k , then we could represent the following four distinct states for n = 4 and k = 3:

({ℓ 1 , ℓ 2 }, {}, {ℓ 3 }, {}), ({ℓ 1 , ℓ 3 }, {}, {ℓ 2 }, {}), ({}, {ℓ 3 }, {ℓ 1 , ℓ 2 }, {}), ({}, {ℓ 2 }, {ℓ 1 , ℓ 3 }, {}). (3.1.3)
Now, if we were to assume that the k lineages are indistinguishable, then there would be no need to use the names ℓ j , and we could encode in each component x i simply the number of lineages it contains. The above states reduce thus to the following:

(2, 0, 1, 0), (0, 1, 2, 0).

(3.1.4)

Alternatively, if we decide we can indeed distinguish the lineages but not the individual islands, then we would need to keep the ℓ j , but also introduce a convention on how to arrange the demes. One way to do so is to consider only non-increasing tuples, i.e., encode each state as α = (α 1 , α 2 , . . . , α n ) where |α i | |α i+1 | for all valid i . Under these assumptions, the states in (3.1.3) reduce to:

({ℓ 1 , ℓ 2 }, {ℓ 3 }, {}, {}), ({ℓ 1 , ℓ 3 }, {ℓ 2 }, {}, {}). (3.1.5)
Finally, if neither islands nor lineages are distinguishable, then all of the cases considered (and more) are included in the following single, terse, state:

(2, 1, 0, 0). So, how many states are there in each case?

Ignoring coalescence events, we proceed to count in how many ways can exactly k lineages be distributed in n demes. See Stanley (2011, §1.9) for more details:

distinguishable n islands indistinguishable n islands distinguishable k lineages n k k 1 + k 2 + • • • + k n indistinguishable k lineages n k p 1 (k) + p 2 (k) + • • • + p n (k)
Here, the symbol k i denotes a Stirling number of the second kind, which counts the number of ways to partition a set of exactly k labelled objects (the lineages) into exactly i non-empty unlabeled subsets (the demes). These numbers can be computed like this:

k i = 1 i ! i j=0 (-1) i-j i j j k .
More on Stirling numbers in L Graham (1994, §6.1).

The symbol n k denotes a multiset coefficient. It counts the number of multisets of k elements which are taken from a set of n elements (multisets are sets that allow repeated elements). It can be easily shown that:

n k = n + k -1 k .
More on multisets in Stanley (2011, §1.2).

Finally, p i (k) is a restricted partition number. It counts in how may ways can the number k be expressed as the sum of exactly i non-null numbers.

In general, an integer partition of k is a multiset of positive integers that sum to k. For example, there are 7 possible partitions of k = 5, which can be represented as:

{5}, {4, 1}, {3, 2}, {3, 1, 1}, {2, 2, 1}, {2, 1, 1, 1}, {1, 1, 1, 1, 1}. (3.1.6)
The order of the terms in the multiset is unspecified, so 5 = 2 + 3 is the same partition as 5 = 3 + 2. An alternative way of representing these partitions is by using product notation,

where we specify the number of times each part is repeated, from 1 to k which is the largest possible part. The partitions in (3.1.6) would be re-written as:

(00001), ( 10010), ( 01100), ( 20100), ( 12000), ( 31000), ( 50000).

(3.1.7)

We denote by P (k), or simply P for simplicity, the total number of integer partitions of k.

There is no known closed-form expression for P , but [START_REF] Hardy | Asymptotic formulae in combinatory analysis[END_REF] provided the asymptotic approximation

P ∼ e π √ 2k/3 4k √ 3 .
Algorithm P generates all partitions of k in standard notation (3.1.6).

Algorithm P: For evaluating routine Partitions(k), which visits, in reverse-lexicographic order, all integer partitions (a 1 , a 2 , . . . , a ω ) such that

ω i=1 a i = k. 1 ω ← 1; a ← (0) 1×k 2 while True: 3 a ω ← k; q ← ω -id k-1 4 while True: 5 visit (a 1 , a 2 , . . . , a ω ) 6 if a q = 2: 7 break 8 a q ← 1; q ← q -1; ω ← ω + 1; a ω ← 1 9 if q = 0: exit x ← a q -1; a q ← x; k ← ω -q -1; ω ← q + 1 while k > x: a ω ← x; ω ← ω + 1; k ← k -x
See Knuth (2005, §7.2.1.4) for the derivation and detailed analysis of this algorithm. Integer partitions can be restricted in a multitude of ways. In this chapter, we are interested in restricting the number of parts. For instance, 5 = 2 + 3 is a partition in two parts, and 5 = 2 + 1 + 1 + 1 is a partition in four parts. The symbols p j (k) thus count the possible number of restricted partitions of k into j parts.

We may choose from these distinguishable/indistinguishable configurations depending on our needs. For instance, if the sizes of the islands are different in the general case, then the islands are inherently distinguishable. It is important to keep in mind that the choice of configuration will affect the size of the state space and the transition rates of the system. In any case, we denote the total number of states in the Markov chain by E.

Returning to the IICR k , in order to compute this function we must be able to evaluate the state of the continuous-time Markov system until the moment of first coalescence.

The migration events after the first coalescence need not be taken into account, nor any coalescences after the first, since these events do not contribute to the distribution of T k . For that reason, any coalescence event that happens in the initial states leads to the only absorbing state of the system. This produces a rate matrix Q with the structure:

Q = Q * b 0 0 E×E , where Q * is of size (E -1) × (E -1)
, and b is of size (E -1) × 1. We are interested in computing P (t) = exp(tQ). The coefficient P (t)[i , E -1] contains the probability that, given an initial sampling in state i , the system has reached the absorbing state by time t, i.e.,

P (T k t).

There are many ways to compute matrix exponentials. See for instance the extensive review in [START_REF] Moler | Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later[END_REF]. An important problem of these methods, however, is that they are based on matrix factorizations, which have a computational complexity of O(E 3 ), with E being the size of the matrix. Additionally, factorization methods generally produce dense factors even when the matrix Q is sparse, which also poses a large challenge in memory requirements. As a result, these methods scale poorly with the size of the problem (number of sampled lineages or number of islands in the model).

An alternative approach can be had with the computation of only the action of the exponential matrix. Indeed, we note that for our application, we are only interested in the last column of exp(tQ). We can express this by saying that we want to compute the vector (π 0 ) ⊤ exp(tQ) without having to compute exp(tQ) first. There are a number of methods that have focused on this in the literature (see for instance Al-Mohy and Higham ( 2011)), and with an asymptotic time complexity of O(E 2 ) they provide a great improvement over the factorization methods mentioned before. However, they still have some drawbacks: these methods are based on the truncated computation of the Taylor expansion of exp(tQ), and a significant part of their computational cost is spend towards finding the ideal cutoff point that minimizes the numerical error. Moreover, they do not take advantage of the sparse nature of the matrix Q.

Recent developments by [START_REF] Sherlock | Direct statistical inference for finite markov jump processes via the matrix exponential[END_REF] provide an improvement over these methods. Sherlock proposes an algorithm that computes the action of the exponential while taking into account the properties of Q by virtue of it being specifically the rate matrix of a Markov process. This algorithm achieves a computational cost of O(Eρ), where ρ is the maximum diagonal element of Q in absolute value.

An implementation of Sherlock's algorithm with an application to demographic inference using the IICR k will be the subject of a future study.

The structured coalescent with ancestry tracking

We proceed to describe the state space for the continuous-time Markov process that models the migrations and ancestry of a sampled population. We do this by introducing a generalization of the classical structured coalescent process.

Given our description of state matrix for a model with n demes and k samples, there is a relationship between the number of sampled lineages and the coefficients α ij of any state matrix α, which is given by the invariant:

1 i k 1 j n i α ij = k, ∀α. (3.2.1)
We refer to the factor i in (3.2.1) as the weight of the lineages that are in the row i . It represents the fact that each lineage in row i is ancestral to i lineages from the sample. The number of lineages that are present in any given state α is denoted |α|, so:

1 i k 1 j n α ij = |α|,
and we call r α the number of coalescences that have taken place up to state α, such that

k = |α| + r α , ∀α.
We note that the state matrix associated with any state α is empty for rows below r α + 1. Indeed, if we suppose α ij 1 for some i r α + 2, then the weight of that lineage is at least r α + 2, and since the remaining |α| -1 lineages have weight at least 1, then the total weight of the state is:

k = 1 i k 1 j n i α ij (r α + 2) + (|α| -1) = k + 1,
which is absurd. We denote by E the set of all possible states α, and by E the total number of states:

E = α ∈ Z k×n + : k i=1 n j=1 i α ij = k , |E| = E.
The set E is clearly defined by the sample size k and number of demes n, but in order to avoid notation clutter we keep this dependency implicit.

Another important notion related to the state matrices are the contribution vectors. The contribution vector of state α, denoted ψ α , is defined as the sum by rows of the states matrices:

ψ α = 1 j n α 1,j , 1 j n α 2,j , . . . , 1 j n α k,j .
Returning to the example from Figure 3.1, the contribution vectors are:

ψ α 1 = (4, 0, 0, 0); ψ α 2 = ψ α 3 = (2, 1, 0, 0); ψ α 4 = ψ α 5 = (1, 0, 1, 0); ψ α 6 = (0, 0, 0, 1).
such a matrix, we first need to impart an order ≺ on the set of all states E. The focus of this section will therefore be the construction of a bijective map

L : E → {1, . . . , E} α ℓ → ℓ (3.2.5) such that 1 ℓ 0 < ℓ 1 E ⇔ α ℓ 0 ≺ α ℓ 1 .
We start with the partial order given by the natural transition of the process from states with more living lineages to states with fewer living lineages. Indeed, migration events do not change the number of living lineages, and a coalescence events decrease them by one, therefore |α| (or equivalently, r α ) decreases (increases) during any given instantiation of the process. We convene then that:

r α < r β ⇒ α ≺ β. (3.2.6)
In the example from Figure 3.1, since r α 1 = 0 and r α 3 = 1, we have α 1 ≺ α 3 . By grouping together the states that share the same number of live lineages, we get a partition of the state space into k disjoint parts (one part for each possible value of r α from 0 to k -1). We denote these parts by E i , each with E i elements, such that E = ∪ i E i and E = i E i :

E i = α ∈ E : r α = i -1 , |E i | = E i , 1 i k.
The partition {E i } can be understood as the quotient set E/R, where R is the equivalence relation

α R β ⇔ r α = r β ⇔ |α| = |β|; α, β ∈ E. (3.2.7)
This partial order imposes a k × k block structure on Q, where block (i , h) has dimensions E i × E h , and contains the rates of transitioning from the states in E i to the states in E h .

Q =      [Q 11 ] E 1 ×E 1 [Q 12 ] E 1 ×E 2 • • • [Q 1k ] E 1 ×E k [Q 21 ] E 2 ×E 1 [Q 22 ] E 2 ×E 2 • • • [Q 2k ] E 2 ×E k . . . . . . . . . . . . [Q k1 ] E k ×E 1 [Q k2 ] E k ×E 2 • • • [Q kk ] E k ×E k      E×E . (3.2.8)
We note that if the process is at a non-absorbing state α ∈ E i , then it can only transition to another state in E i via migration, or to a state in E i+1 via coalescence, therefore all the coefficients of block [Q ih ] are null when h = i and

h = i + 1. Matrix Q is thus block bi-diagonal, with blocks A i = [Q ii ]
for 1 i k in the diagonal, and blocks

B i = [Q i, i+1 ] for 1 i k -1 in the upper diagonal. Q =        A 1 B 1 A 2 B 2 . . . . . . A k-1 B k-1 A k        E×E .
We call A i the migration blocks and B i the coalescence blocks of Q. We further note that migration block A k is null, since there are no more transitions once the process arrives at any of the stationary states of E k . An alternative representation of Q is thus possible:

Q = Q * B 0 0 , (3.2.9)
where

Q * =        A 1 B 1 A 2 B 2 . . . . . . A k-2 B k-2 A k-1        E * ×E * and B =        0 0 . . . 0 B k-1        E * ×E k .
Here, Q * is known as the sub-intensity matrix and contains the rates of transition between the transient states

E * = E -E k . E * is the number of transient states E -E k .
Next we focus on defining the order ≺ for the states within the sets E i . We begin by defining another equivalence relation C between states of E. We say that two states are equivalent under C if they have the same contribution vector:

α C β ⇔ ψ α = ψ β ; α, β ∈ E. (3.2.10)
This relation partitions the state space into disjoint sets identified by their common contribution vector. We can thus think of each different contribution vector as the representative of an equivalence class of states, or a macrostate. In order to count how many different macrostates there are in a process with sample size k, we need to briefly return to integer partitions. It is clear from their product representation (3.1.7) that the set of all integer partitions of k is the set P of all possible state contribution vectors defined in (3.2.2). There are thus P possible macrostates in the process. The number of parts in macrostate ψ α ∈ P is given by k i=1 ψ α (i ) = |α|. We group the macrostates by their number of parts following a similar procedure to how we group the states by their number of live lineages in (3.2.7). Overloading the name R to define a relation on P, we say that: .2.11) This partitions P into k disjoint sets P i , each with P i macrostates, such that P = ∪ i P i and P = i P i :

ψ α R ψ β ⇔ |α| = |β| ⇔ r α = r β ; ψ α , ψ β ∈ P. ( 3
P i = ψ α ∈ P : r α = i -1 , |P i | = P i , 1 i k.
We sort the macrostates of P i in reverse lexicographic order of standard notation (3.1.6) as does Algorithm P, so for k = 5 we have (4, 1) ≺ (3, 2) in P 2 , and (3, 1, 1) ≺ (2, 2, 1) in P 3 . This is a total order in P i , so for each i ∈ {1, . . . , k} we will assume the existence of a bijective map L i defined on P i such that L i (ψ) ∈ {1, . . . , P i }, and for ψ, ϕ ∈ P i , we have ψ ≺ ϕ ⇔ L i (ψ) < L i (ϕ).

Since ψ α ∈ P i ⇔ r α = i -1 ⇔ α ∈ E i , this order in P i motivates a natural way to introduce a partial order in the elements of E i . We convene that for α, β ∈ E i ,

ψ α ≺ ψ β ⇒ α ≺ β.
(3.2.12)

We can group together the uncomparable states of E i by returning to the equivalence relation C defined in (3.2.10) such that α C β ⇔ ψ α = ψ β for α, β ∈ E i . There are P i equivalence classes in the quotient set E i /C, which we will denote E ip , each with E ip states. We now have E i = ∪ P i p=1 E ip , and

E i = P i p=1 E ip E ip = α ∈ E i : L i (ψ α ) = p , |E ip | = E ip , 1 p P i .
As was the case with Q in (3.2.8), this partial order on E i imposes a P i × P i block structure on migration block A i , and a P i × P i+1 block structure on coalescence block B i . For p, q P i , sub-block (p, q) of A i has dimensions E ip × E iq , and for q ′ P i+1 , sub-block (p, q ′ ) of B i has dimensions

E ip × E i+1, q ′ .
Upon closer inspection, we see that the contribution vectors of the states in E i are invariant under migration. Indeed, if β = αe ij + e ig , then the sum of row i of both states remains unchanged, and thus ψ α = ψ β . This implies that sub-block (p, q) of migration block A i has all null entries except when p = q, and thus A i is P i × P i block-diagonal. We denote the diagonal sub-blocks of A i as A ip ∈ R E ip ×E ip for 1 p P i . There is no such simplification for coalescence block B i , but we group all its sub-blocks in each block-column as B iq ∈ R E i ×E i+1, q with 1 q P i+1 for later convenience:

A i =      A i1 A i2 . . . A i P i      E i ×E i ; B i = B i1 B i2 • • • B i P i+1 E i ×E i+1 . (3.2.13)
To make ≺ a total order, we need to define it for the states within each of the P subsub-spaces E ip . There is no further benefit in doing so in terms of the structure of the rate matrix Q in the general case, but it is still a necessary step for building it. We thus convene that the order of the states in E ip is the one obtained by applying algorithm M to state α = L -1 i (p) ⊤ , 0 k×1 , . . . , 0 k×1 and all its subsequent migration destinations repeatedly (this procedure is detailed in algorithm A). Therefore, for every i = 1, . . . , k and p = 1, . . . , P i we assume the existence of a bijective map L ip defined on E ip such that L ip (α) ∈ {1, . . . , E ip } and α ≺ β ⇔ L ip (α) < L ip (β). With this construction, we can finally map every state α ∈ E to an index ℓ ∈ {1, . . . , E}, and the definition of the global map L (3.2.5) would be

L(α) = σ ip + L ip (α), ∀α ∈ E (3.2.

14)

Algorithm A: For evaluating routine MigrationSubBlock(M, α 0 ), which given the pairwise migration rates M of the island model and a protostate α 0 (3.2.15) in E ip , builds and returns (A ip , L ip ).

1 A ip ← (0) 1×1 ; L ip ← (α 0 , 1) ; U ← (α 0 ) 1×1 2 while size(U) > 0:

3 (α, U) ← deq(U) 4 ℓ 0 ← L ip (α) 5
for each (β, m) in MigrationDestinations(M, α):

6 if β ∈ D L ip : 7 ℓ 1 ← L ip (β) 8 A ip (ℓ 0 , ℓ 1 ) ← A ip (ℓ 0 , ℓ 1 ) + m 9 else: ℓ 1 ← |L ip | + 1; L ip ← L ip ∪ (β, ℓ 1 ) A ip ← A ip 0 0 0 ℓ 1 ×ℓ 1 ; A ip (ℓ 0 , ℓ 1 ) ← m U ← enq(U, β) A ip (ℓ 0 , ℓ 0 ) ← A ip (ℓ 0 , ℓ 0 ) -m 14 return (A ip , L ip )
Since at this point all the states E ip are known, as well as the index maps L i and L ip , the procedure simply consists in visiting, for every state α, all possible states β that can be reached by coalescence. This is accomplished by using Algorithm C and accumulating the coalescence rates in coefficient B iq (σ ipσ i + L ip (α), L i+1, q (β)) for the corresponding value of q. The detailed procedures are presented in lines [9 ••• 18] of Algorithm Q. In general, Algorithm Q presents the complete construction of the Markov process, showing the details of how to combine all algorithms specified so far. Here, we omit some of the more menial details of the process, which includes how to compute the sizes E i and E ip given the existing information, and how to assemble the final rate matrix Q and index map L by arranging their respective sub-parts A ip , B ip and L ip .

Since we are interested in the normalized AFS ξ = ξ/ ξ 1 , the parameter θ is cancelled away, so we focus exclusively on the quantities E(ζ i ). These can be computed by adding, for every transient state α ∈ E\E k = E * , the number of lineages of weight i (given by ψ α (i )) multiplied by the expected total time spent in the state before absorption, which we denote by T α . We thus have:

E(ζ i ) = α∈E * ψ α (i ) E(T α ).
(3.3.1)

In order to compute the times E(T α ), we return to the sub-intensity matrix Q * defined in (3.2.9). This matrix is non-singular, and G = (-Q * ) -1 , known as the process' Green matrix, has the property that coefficient G(ℓ 0 , ℓ 1 ) is the expected time that the process spends in state α ℓ 1 before absorption, given an initialization at state α ℓ 0 (see Bladt and Nielsen (2017), Theorem 3.1.14 ). With an initial distribution for the transient states π 0 = π 0 (1), π 0 (2), . . . , π 0 (E * ) ⊤ , we can then write E(T α ℓ ) as

E(T α ℓ ) = E * ι=1 π 0 (ι) G(ι, ℓ), (3.3.2)
where

E * = |E * | = E -E k .
We can combine and rearrange equations (3.3.1) and (3.3.2) to get

E(ζ i ) = E * ι=1 π 0 (ι) E * ℓ=1 G(ι, ℓ) ψ α ℓ (i ),
which can be further simplified when considering all values of 1 i < k as (3.3.3) where Ψ, called the contribution matrix, is an E * × (k -1) matrix where coefficient Ψ(ℓ, i) = ψ α ℓ (i ) is the number of lineages of weight i in state α ℓ .

E(ζ) = (π 0 ) ⊤ • G • Ψ,

Numerical solution

We now turn to the question of how to numerically solve problem (3.3.3). The computation of matrix G involves inverting Q * , which grows exponentially in size with the number of the samples k (see Figure 3.3), so it is inevitable that the computational cost of the solution will eventually exceed the available resources regardless of the underlying method, given a large enough value of k = k critical . Due to this inherent and unavoidable complexity, our strategy simply consists in exploiting the sparsity and structure of Q ( §3.2.3) to delay k critical as much as possible.

We avoid the direct computation of the inverse G. Instead, we compute the product (π 0 ) ⊤ • G by solving the linear system (3.3.4) where A = (-Q * ) ⊤ . The unknowns of this system are the coefficients of a vector t ∈ R E * ×1 , which contains in component ℓ the time E(T α ℓ ) given initial distribution π 0 . Following the notation conventions of §3.2.3, we denote by t i and π 0 i the segments of vectors t and π 0 corresponding to states in E i , for i = 1, . . . , k -1; and by t ip and π 0 ip the sub-segments corresponding to the states in E ip for p = 1, . . . , P i . We can thus write system (3.3.4) in the following block form: .3.5) This structure immediately suggests the sequence of sub-problems:

A • t = π 0 ,
     -A ⊤ 1 -B ⊤ 1 -A ⊤ 2 . . . . . . -B ⊤ k-2 -A ⊤ k-1      •      t 1 t 2 . . . t k-1      =      π 0 1 π 0 2 . . . π 0 k-1      . ( 3
-A ⊤ 1 • t 1 = π 0 1 , -A ⊤ i • t i = b i • • = π 0 i + B ⊤ i-1 • t i-1 , i = 2, . . . , k -1. (3.3.6)
We note that in (3.3.6), the problems must be solved sequentially, since the solution segment t i depends on the previous one t i-1 . Assuming thus that i 2 and t i-1 is known, we now focus on computing t i . Exploiting the sub-structure of the system matrix A i (3.2.13), the problem -A ⊤ i • t i = b i can be decomposed into the following P i sub-problems (Figure 3.2):

-A ⊤ ip • t ip = b ip • • = π 0 ip + B ⊤ i-1, p • t i-1 , p = 1, . . . , P i . (3.3.7)
We note that the computation of the solution sub-segment t ip does not depend on any other sub-segments of t i ; only on the already-available segment t i-1 . This implies that all P i problems in (3.3.7) can be solved in parallel. For each one of them we use the Successive Over-Relaxation iterative method: given an initial guess t 0 ip , the scheme (known as SOR(ω)) generates a sequence of solutions t 1 ip , t 2 ip ,. . . ,t κ ip , where coefficient t κ ip (ℓ) is computed using information from t κ ip (ι) for ι < ℓ and t κ-1 ip (ι) for ι > ℓ, so there is no need to keep two versions of the solution sub-segment t ip in memory. We use the norm of the residual vector res κ ip = -A ⊤ ip • t κ ipb ip as the convergence metric, so for a given desired precision δ, the stopping criterion is res κ ip ∞ < δ. We denote by tip the final approximation t κ ip of the unknown vector t ip . The detailed procedure is given in lines

[7 ••• 10] of Algorithm X.
The parameter ω in the SOR method affects the rate of convergence. A necessary condition for convergence is 0 < ω < 2, and for system matrices that are M-matrices (as is the case for all the A ip ), a sufficient condition for convergence is 0 < ω 1 (see §7.2.3 in Hackbusch (1994)). In practice, values of ω larger than 1 can be used to accelerate convergence despite not satisfying the sufficient convergence condition.

Math Minute

We referenced the concept of M-matrices in the previous argumentation. An M-matrix is

Error analysis

Algorithm X guarantees that after each final iteration of the SOR scheme, the residues

res ip = -A ⊤ ip • tip -b ip satisfy res ip ∞ < δ.
It is then easy to verify that the global residue of system (3.3.4) given by res = A • tπ 0 satisfies res ∞ < δ. This is a consequence of the fact that the vectors res ip are simply sub-segments of res, and that for X ⊂ R we have max x∈X |x| < δ ⇔ ∀x ∈ X, |x| < δ. We can thus say that given δ > 0, Algorithm X computes ζ = t⊤ • Ψ, an approximation of E(ζ) (3.3.3) where

A • t -π 0 ∞ < δ.
However, the "end user" of Algorithm X is not necessarily interested in t or the residue, but rather in

E( ξ) = E(ζ) E(ζ) 1
, so it would be more meaningful to control the norm of the error err = E( ξ) -ξ, where ξ is the approximation of E( ξ) given by Algorithm X:

ξ = ζ ζ 1 .
Our goal for this section is therefore to find a way of controlling err ∞ by controlling res ∞ .

Math Minute

Some preliminary results about p-norms and their induced matrix norms:

• If 0 < p < q ∞, then for any x ∈ R n :

x p n 1 p -1 q x q and x q x p , (3.3.8)

where we convene that 1/q = 0 for the case q = ∞. This result follows form Hölder's inequality for R n vectors. No tighter bounds are possible in the general case, since equality is achieved with x = (1, 1, . . . , 1) and x = (1, 0, . . . , 0) respectively.

• For any non-null x, y ∈ R n we have

x x - y y 2 x -y x , ( 3
.3.9) regardless of the norm • . Proof:

x x - y y = x x - y x + y x - y y = 1 x (x -y ) + 1 x - 1 y y x -y x + x -y x y y 2 x -y x .
• Any norm in R n can be used to define a norm in R m×n as follows:

A = sup A • x : x ∈ R n , x = 1 = sup A • x / x : x ∈ R n , x = 0 .
Matrix norms defined this way are called induced norms. Their most useful property is a direct consequence of the definition:

A • x A x ,
In addition to the standard norm properties, induced norms also satisfy the submultiplicative property:

A • B A B .
Two notable induced norms:

A 1 = max 1 j n m i=1 |a ij |, A ∞ = max 1 i m n j=1 |a ij |. Notice that A ∞ = A ⊤ 1 .
A lower bound for err ∞ can be derived as follows:

err ∞ err 1 = E(ζ) E(ζ) 1 - ζ ζ 1 1 2 E(ζ) -ζ 1 E(ζ) 1 = 2 Ψ ⊤ • (t -t) 1 E(ζ) 1 = 2 Ψ ⊤ • (A -1 π 0 -t) 1 E(ζ) 1 = 2 Ψ ⊤ • A -1 • (π 0 -A • t) 1 E(ζ) 1 = 2 Ψ ⊤ • A -1 • res 1 E(ζ) 1 2 Ψ ⊤ 1 A -1 1 res 1 E(ζ) 1 2E * Ψ ∞ G ∞ res ∞ E(ζ) 1 .
We note that since Ψ ∞ is the largest possible row sum of Ψ, and the rows of Ψ are the contribution vectors of the states of E * , then

Ψ ∞ = max α∈E * |α| = k. We also note that G ∞ = (Q * ) -1 ∞ . With this we have err ∞ 2kE * (Q * ) -1 ∞ E(ζ) 1 res ∞ . (3.3.10)
In order to continue we need an upper bound for E(ζ) 1 . In what follows we note that since Ψ 0 and G ⊤ = A -1 0, then all of the coefficients of G • Ψ are also non-negative.

E(ζ) 1 = (π 0 ) ⊤ • G • Ψ 1 = k-1 i=1 E * ℓ=1 π 0 (ℓ) (G • Ψ)(ℓ, i) = E * ℓ=1 π 0 (ℓ) k-1 i=1 (G • Ψ)(ℓ, i ) = E * ℓ=1 π 0 (ℓ) k-1 i=1 E * ι=1 G(ℓ, ι)Ψ(ι, i ) = E * ℓ=1 π 0 (ℓ) E * ι=1 G(ℓ, ι) k-1 i=1 Ψ(ι, i ) = E * ℓ=1 π 0 (ℓ) E * ι=1 G(ℓ, ι)|α ι |.
Since ι E * , α ι is a transient state and therefore has at least two live lineages, thus:

E(ζ) 1 E * ℓ=1 2π 0 (ℓ) E * ι=1 G(ℓ, ι) E * ℓ=1 2π 0 (ℓ) G(ℓ, ℓ).
The coefficient G(ℓ, ℓ) is E(T α ℓ |X 0 = α ℓ ), the expected total time the process spends in state α ℓ prior to absorption, given that the process started in that same state. This expected time can be computed as the product of the total expected number of visits to state α ℓ E(N α ℓ |X 0 = α ℓ ), multiplied by the expected duration of each visit, which is an exponentially distributed time with parameter -Q(ℓ, ℓ). An upper bound for G(ℓ, ℓ) thus can be obtained by observing that E(N α ℓ |X 0 = α ℓ ) 1, since at least one visit is expected at state α ℓ in lieu if it being the initial state. Additionally, we note that:

-Q(ℓ, ℓ)

E * ι=1 |Q(ℓ, ι)| max 1 ℓ E * E * ι=1 |Q(ℓ, ι)| = Q * ∞ .
This implies that the expected duration of one visit to state α ℓ is greater than 1/ Q * ∞ . With this we have:

E(ζ) 1 2 Q * ∞ E * ℓ=1 π 0 (ℓ) = 2 π 0 1 Q * ∞ .
Returning then to (3.3.10): (3.3.11) where κ(Q * ) is the condition number of the matrix Q * . These results can be summarized in the following Lemma 2. For any ε > 0, it is sufficient to choose δ as

err ∞ 2kE * (Q * ) -1 ∞ E(ζ) 1 res ∞ kE * Q * ∞ (Q * ) -1 ∞ π 0 1 res ∞ = kE * κ(Q * ) π 0 1 res ∞ ,
δ = π 0 1 kE * κ(Q * ) ε (3.3.12) to have res ∞ < δ ⇒ err ∞ < ε.
Having a clear criterium for how to control the global error of the normalized expected SFS E( ξ), we can proceed to present the main and final algorithm of this chapter: Algorithm R: For evaluating routine ExpectedAFS(k, n, M, S, α 0 , ω, ε), which given the model parameters and an initial sampling state α 0 ∈ E 1 , returns an approximation of E( ξ) with global error err ∞ < ε using a SOR(ω) iterative scheme.

Early termination

The following discussion is concerned with the eventual use of Algorithm R in the context of demographic inference, where early termination is desirable.

Given a value of ε, Algorithm R computes a numerical approximation of E( ξ) with an error err ∞ < ε. However, the internal solution method we use for each of the P subsub-problems (3.3.7) is iterative, meaning that progressively-accurate approximations for which t ip are being computed, which in turn could be used to compute progressively-accurate approximations of ζ. Obtaining early approximations of ζ (i.e., approximations where err ∞ is not smaller than ε yet) could be useful for implementing more flexible stopping criteria, where we allow the method to converge to the desired high accuracy ε only if the SFS is approaching a previously designated target SFS; and otherwise terminate the process prematurely, thus saving computational resources.

The key insight for our proposal comes from the observation that Algorithm X operates by simply updating or replacing a starting solution t 0 with a more accurate one t, for which res ∞ < δ. This means that it is possible to continue refining the solution by repeatedly invoking the routine, each time passing as t 0 the final result t from the previous invocation, and a smaller value of δ.

One possible realization of this idea is as follows: suppose we are interested in observing d intermediate solutions ζ prior to the convergence with err < ε, then we may define d intermediate tolerances:

ε 1 = ε 1/d . . . ε κ = ε κ/d . . . ε d = ε, (3.3.13)
along with the corresponding δ κ = δ κ/d , and invoke Algorithm X d times: once for each of the δ κ . After each invocation we evaluate the corresponding early approximation of E( ξ), and decide whether to continue or terminate the process. The sequence (3.3.13), while certainly not unique, satisfies the desirable properties that the ε κ are in descending order (assuming ε < 1), and are uniformly log-spaced. This implies that the subsequent invocations of Algorithm X, which improve the accuracy of the solution from δ κ to δ κ+1 , should take similar times to execute.

There is a tradeoff to be balanced between execution time and level of control when choosing the value of d, since the products B i-1, p • ti-1 and Ψ ⊤ • ζ (lines [6] and [12 ••• 13] of Algorithm X, respectively) have to be computed d times each. See Figure 3.10.

The following revision of Algorithm R takes an additional parameter (d), and visits d where the columns no longer indicate the island number, but are simply an arbitrary indexing of the non-empty islands. We can define a unique class representative by sorting the columns of the equivalent state matrices. We use the reverse-lexicographic ordering on the columns (read from top to bottom). With this ordering, the class representative of the equivalence set (3.4.1) would be α 1 .

We will assume the existence of a function ClassRepresentative that maps E onto itself, and assigns to every state matrix its class representative. In this section we present modified versions for some of the previously derived algorithms in order for them to work for the n-island model specialization.

The results shown in this section were obtained with the software SISiFS (short for Symmetrical Islands Site Frequency Spectrum). The implementation was done using the D programming language, and can be found in the public repository github.com/arredondos/sisifs.

The rate matrix, revisited

We saw in §3.2.2 that the model state space is discovered by generating all possible migration and coalescence events (algorithms M and C, respectively). In order to adapt these for a symmetrical n-island specialization, all we have to do is to only visit the states classrepresentatives, which as discussed previously, are obtained by sorting the columns of the state matrices. Therefore, the adaptation will consist of, given a destination state β, generating its class representative β * and visiting it instead.

Algorithm sM: For evaluating routine SymmetricalMigrationDestinations(m, α), which given the migration rate m and a state matrix α of size k × n, visits all the tuples (β, m) where β is a class-representative state satisfying (3.2.3) for symmetrical island models, and m is the corresponding transition rate.

1 for each i in {1, 2, . . . , k -1}: 2 for each j in {1, 2, . . . , n}:

3 if α ij 1: 4 for each g in {1, 2, . . . , j -1, j + 1, . . . , n}: 5 β ← α -e ij + e ig ; β * ← ClassRepresentative(β) 6 visit (β * , α ij m)
This simple modification allows for a large compression of the state space. See Figure 3.3 for a comparison of the number of states of the Markov systems associated to the symmetrical n-island model, and the most general demographic model possible, where all the states are vector is [14,6,6], there are spikes in the expected SFS at frequencies 14, 6, 14 + 6 and 6 + 6. Another observation is that this effect is amplified when the migration value is small, since it increases the probability that coalescence happens before migration in a deme. Likewise, a smaller number of islands also amplifies this effect of structure, although to a lesser degree than M.

We can quantify the effect of structure on the SFS by measuring the deviation of the expected SFS in these n-island models from the expected SFS under panmixia (1.4.3). We present these results in Figure 3.6. We notice the same spikes near the sampling frequencies, but also that for many of the structure parameters, there is a smaller-than-expected (as compared to a panmictic model) number of unique mutations (singletons and doubletons) and a higher-than-expected number of shared mutations. These results are not entirely new (see for instance Figure 4.11 in [START_REF] Hein | Gene genealogies, variation and evolution: a primer in coalescent theory[END_REF]). However, we can now compute these figures exactly (up to machine precision), and very quickly. Figure 3.7 shows the required time to compute run Algorithm X for some parameter values of the n-island. The sizes of the associated Markov chains can be referenced in Figure 3.3. For instance, we observe that for a sample of size k = 26, the rate matrix has over 8 million states, and the implementation of Algorithm X in SISiFS iteratively solves the associated linear system to a precision of ε = 10 -6 in under 2 seconds, if the migration rate is not too large. The migration rate M affects the performance of the method because the condition number of the system matrices κ(A ip ) increases rapidly with large values of M.

The effect of the high condition number for large M can also be observed as numerical instability in Figure 3.8. This figure showcases the convergence of the expected SFS under various n-island models to the expected SFS under panmixia, for which the closed-form expression is known (1.4.3) when M → ∞. This convergence is expected since island models with a very high migration rate are known to behave like panmictic models with a population size equal to the sum of the island sizes. We verify this fact by making M increasingly larger and comparing the two expected SFSs by computing the • ∞ -based distance between them. We observe that the algorithm converges exponentially to the correct solution in all but the largest cases (k 24). This result serves as a validation of Algorithm X, but also indicates that in practice, the degradation of speed and accuracy due to very large values of M does not impose a great limitation on the method, given that the expected SFS can be well-approximated by (1.4.3) in such cases.

As an additional form of validation, we conducted simulations using the ms software and computed the empirical SFS for samples of various sizes under the n-island model. We simulated independent chromosomes continuously with a given value of θ until we achieved a certain number of observed segregating sites. After that, all the information at these sites was aggregated and an empirical SFS was computed. The goal is to compare the obtained results with the ones given by our methods. This approach validates both algorithms Q and X. Figure 3.9 shows the results of these comparisons. We find that estimating the expected SFS with this method is a poor substitute for exactly computing it using our approach for the tested sampling sizes (k 26), both from the numerical accuracy perspective as well as from the required computation time (we note that even though Figure 3.9 does not show the times required to complete these ms simulations, our method was so many orders of magnitudes faster for the same accuracies, that we felt that a direct comparison was unnecessary).

Figure3.10 showcases the early-termination feature of Algorithm R' (see §3.3.3). Here, we can see the trade-off between the speed of convergence and the flexibility afforded by the error introspection. Increasing the value of the parameter d allows for additional intermediate results of the expected SFS to be available. This could be useful for terminating the computation of the expected SFS if it is not converging towards the desired curve. For instance, in the context of a hypothetical demographic inference application, the computation of the SFS could be aborted before reaching the accuracy specified by ε if the expected SFS if not converging towards the target (observed) SFS, and the search algorithm could then focus on testing other parameters.

As expected, the fastest results are achieved with d = 1, but since the convergence is often achieved in one or two iterations, it also allows for minimum flexibility. This is the recommended value for one-time computations. In general, the optimum value of the d parameter will be heavily dependant on the application that Algorithm R' is being used in. also observed that the estimated value of n scaled well with increasing model complexity. A similar consistency can be observed with the deme size parameter N (see Figure 2.4). We give up some flexibility in the model by keeping the number of demes constant throughout the history of the population, so the timed demographic events cannot represent splits or joining of populations even though such events are likely to have taken place in the history of species. Additionally, in the n-island model we do not account for possible asymmetrical gene flow or different deme sizes, even when the theoretical framework does allow for such representations. However, it is a more challenging problem to validate due to the fact that during any given component, changing both the migration rate and the deme size have confounding effects on the IICR curve which can be hard to separate. This requires a dedicated study with a different methodology which will be explored in a future work.

A note on human evolution

An application of our IICR-based inference method to five publicly available human PSMCs suggests that the backwards long term history of the sampled individuals, when accounting for possible recent expansions and the noise introduced by the PSMC method, can be accurately modelled in the framework of a symmetrical island model of approximately 10 to 12 demes with varying levels of connectivity through time. Only one of the five samples (Yoruba) displayed less consistent evidence of this finding, which may indicate that more complex models (possibly including asymmetric gene flow, spatial modelling of the environment, or changes in deme sizes) could be needed to explain the full complexity of the data.

These findings regarding changes in connectivity and number of islands are in agreement with the results of Rodríguez et al. (2018), in which a hand-fitting approach of the IICRs was used to arrive at an estimate of 10 islands with a similar value of N and a comparable period featuring a significant increased of gene flow between 600 Kya and 2 Mya. Note that the timing in years and the deme sizes in Rodríguez et al. (2018) differ due to the change in mutation rate.

We also compared our results with the tree model for human evolution published by [START_REF] Noskova | GADMA: Genetic algorithm for inferring demographic history of multiple populations from allele frequency spectrum data[END_REF] (the C3PO model), which is a revision of the model from [START_REF] Gutenkunst | Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data[END_REF] and represents a simplified model of human evolution [START_REF] Jouganous | Inferring the joint demographic history of multiple populations: beyond the diffusion approximation[END_REF][START_REF] Kamm | Efficiently inferring the demographic history of many populations with allele count data[END_REF]. The C3PO model proposes an ancestral human population that experiences two splits: an old one that resulted in the current African "population" and another more recent one that resulted in the current European and Asian "populations". The parameters of this model include the times of these events, the population size history of these populations and their ancestral branches, and the migration rates between them. The summary statistic targeted by these methods is the AFS, and we see that a fitting AFS does not guarantee a fitting IICR and vice versa [START_REF] Chikhi | The IICR (inverse instantaneous coalescence rate) as a summary of genomic diversity: insights into demographic inference and model choice[END_REF][START_REF] Beichman | Comparison of single genome and allele frequency data reveals discordant demographic histories[END_REF]. Indeed, the IICRs of the populations from the C3PO model do not resemble those of the real humans. Likewise, when we use the C3PO model to generate IICR curves, and infer the corresponding demographic history using SNIF, these inferred n-island models do not resemble those inferred directly from the human IICRs, and the parameter values are less consistent across different runs.

On one hand, the n-island is the most symmetrical of the structured models, therefore it allows for the greatest compression on the state space, so from this perspective, other models would be expected to require larger state spaces and thus more computational resources. On the other hand, the n-island model has the densest migration matrix of all structured models (note that in the stepping stone model for instance, the migration matrix itself is sparse). This means that other models will generate even sparser rate matrices. Given that the computational cost in sparse-algebra methods are associated with the number of non-null coefficients, this could counter the effect of the increased state space. It is hard to estimate how these two opposite effects will ultimately affect the performance of the methods. This is why more specializations have to be developed and tested.

Structured demographic inference using multi-sample statistics

One immediate use case for the methods presented here is to use the SFS as a summary statistic for parameter inference in structured models.

We have shown that, at least for the case of the n-island model, we can achieve good performance and accuracy when computing the SFS. This is in stark contrast with simulationbased estimates of the SFS, which are significantly slower for all tested sample sizes and less accurate (see Figure 3.9). This opens the opportunity for an approach similar to that of Chapter 2, where the summary statistic in use would be the SFS. There are two important questions that would need to be researched here: how to manage the sampling vector and how to choose the distance function.

The sampling vector is an important aspect of any structured demographic model that incorporates multiple samples. We have seen that it can have a profound impact on the SFS of otherwise identical n-island models (Figure 3.5), and such is the case as well for the IICR k (see the differences between IICR s and IICR d in [START_REF] Chikhi | The IICR (inverse instantaneous coalescence rate) as a summary of genomic diversity: insights into demographic inference and model choice[END_REF]).

One initial approach could be to assume that the sampling vector is known, although this defeats the purpose of an inferential approach. Alternatively, we could attempt to completely infer the sampling vector. This would add a great burden of dimensionality to the parameter space which could render the method impractical. Perhaps a compromise can be achieved between these two extreme approaches, where any existing information about the samples is incorporated in the model (for instance, lineages 1 and 2 are always from the same deme, as are lineages 3 and 4, etc.). In models that are spatially aware, such as the many variants of the stepping stones model, an alternative approach to managing the sampling vector is possible. We can assume that the actual (physical) sampling locations are known, given a map of the sampling region. We can then associate each sample with the nearest model deme corresponding to that region of the map. This approach allows the sampling vector to be constructed organically from the rest of the model parameters. See Regarding the question of what distance to use for inference, this will depend on the exact nature of the method. For purely SFS based methods, log-likelihood metrics have been successfully used in the past [START_REF] Gutenkunst | Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data[END_REF][START_REF] Excoffier | Robust demographic inference from genomic and snp data[END_REF]. However, a multi-sample demographic-inference method need not be SFS-only based. In cases where full-genome information is available for the sampled individuals, it is conceivable to combine the predictive power of both the IICR and the SFS summary statistics. This could be modeled either by using several IICR curves, belonging to each of the sampled individuals. The optimization goal in these scenarios is less clearly defined, as is often the case with multiobjective optimization problems. An initial approach could be based on a linear combination of the distances used for each of the summary statistics. Another potential direction is to use multiple IICR curves simultaneously during the inference process. These multiple IICRs may come in the form of more than one IICR sampled from an asymmetrical demographic model (for which the initial sampling deme does result in different curves [START_REF] Chikhi | The IICR (inverse instantaneous coalescence rate) as a summary of genomic diversity: insights into demographic inference and model choice[END_REF], as opposed to the n-island model where demes are by definition indistinguishable). They may also be in the form of multiple IICR k curves where k is the number of sampled haploid genomes. Indeed, the IICR of [START_REF] Mazet | On the importance of being structured: instantaneous coalescence rates and human evolution-lessons for ancestral population size inference?[END_REF] was defined for k = 2, and this is the IICR that we have been studying in our previous works. However, the concept can be extended to more haploid genomes in the same way that the MSMC method [START_REF] Schiffels | Inferring human population size and separation history from multiple genome sequences[END_REF] is an extension of the PSMC to multiple genomes, which takes into consideration the distribution of the coalescent time T k . The precise concept of the IICR k is currently being developed in a separate study.

These approaches may prove beneficial in choosing between structured and non-structured models. Indeed, [START_REF] Grusea | Coalescence times for three genes provide sufficient information to distinguish population structure from population size changes[END_REF] shows that using more than one IICR curve can help discriminate between structured and non structured scenarios in the n-island model. Finally, the incorporation of larger samples not only enables exploring more complex scenarios, but it also allows using other summary statistics to complement the IICR, most notably among them the AFS, which is widely used for the purposes of demographic inference.

Closing thoughts

We have presented in this work an inference method for automatically estimating demographic parameters under a piecewise stationary symmetrical island model that uses the IICR as its summary statistic. The underlying methodology consists in quantifying the discrepancy between a target IICR and many simulated IICR curves for a large number of candidate scenarios, and using this metric to drive a global optimization process. With a large number of validations we have shown that the method works accurately and consistently for a diverse range of parameter values, and we additionally showed an application to human data that agrees with and improves upon previously published results using similar approaches. We believe that despite its current scope, our IICR-based method can be of great value during the initial exploration of the parameter space for simple models, and thus can also provide a starting point for manually fitting the IICR with models that could express spatial structure and varying N (Rodríguez et al., 2018).

With the goal of extending this method to multi-sample statistics, we have taken the first steps towards incorporating the SFS by providing a method for exactly computing it in the context of structured demographic models. The proposed method allows for an arbitrary number of demes and an arbitrary connectivity matrix between the demes. We also provide an implementation for the special case of the n-island model which was shown to be superior to simulation-based estimates of the SFS in both runtime (for the tested sample sizes) and accuracy. These developments should allow us to incorporate the SFS as a summary statistic in inference applications for demographic models of structure in which the number of islands is unknown and potentially large.

Figure 1 . 2 :

 12 Figure 1.2: Examples of models of structure.

Figure 2 . 3 :

 23 Figure 2.3:The three ways of obtaining a target IICR. The panels represent different types of discretised target IICRs for the same two-components demographic history with N = 10 3 , n = 4, t 1 = 7, (M i ) = (0.5, 1.0) and (s i ) = (1.0, 1.0). (a) Exact IICR, computed as per section 2.1.2 and discretized to 32 intervals. In gray we show a smoother discretization and keep it in the other panels for reference. This is obtained using the approach ofRodríguez et al. (2018). (b) T-sim IICR, obtained by simulating 2 × 10 4 T 2 samples using the ms command ms 2 20000 -T -L -I 4 2 0 0 0 0.5 -eM 3.5 1 -eN 3.5 1 and later scaling the IICR using the value of N. This is the approach used in[START_REF] Chikhi | The IICR (inverse instantaneous coalescence rate) as a summary of genomic diversity: insights into demographic inference and model choice[END_REF] (c) Seq-sim IICR, obtained by running PSMC on a genome simulated using the ms command ms 2 100 -t 100 -r 20 2000000 -p 8 -I 4 2 0 0 0 0.5 -eM 3.5 1 -eN 3.5 1 and later scaling with mutation rate µ = 1.25 × 10 -8 . This is obtained by running the PSMC method of[START_REF] Li | Inference of human population history from individual whole-genome sequences[END_REF]

Figure 2 . 4 :

 24 Figure 2.4: Normalised root-mean-square deviation of the inferred parameters during validation. The top part of each sub-panel corresponds to the parameters of unscaled scenarios, and the bottom part to scaled scenarios. The number of components c of the scenarios is in the horizontal axis in every case.

Figure 2 . 5 :

 25 Figure2.5: Parameter misidentification. We measure the r 2 correlation between M 2 and Mj in a sample of 100 scenarios with c = 5 components. In 'random global' we draw M j from components 0, 1, 3, and 4. In 'random far' we draw M j from components 0 and 4. In 'best far' we draw M j from components 0 and 4 again, but always keep the best fitting value of the two. Finally, in 'best near' we draw M j from components 1 and 3 and always keep the best fitting value. 'Filtered at 10%' indicates that the scenarios where the inferred value was within 10% of the simulated value were excluded.

Figure 2 . 7 :

 27 Figure 2.7: Number of performed rounds during validations with exact IICRs. The top panel shows the number of rounds (up to a maximum of 500) that the method required before converging to a scenario with a distance value smaller than the tolerance ε = 10 -10 for the unscaled case. The bottom panel shows the same information for scaled scenarios and a tolerance of ε = 10 -7 . We represent in different colors the curves corresponding to scenarios with different simulated (and inferred) number of components, ranging from c = 2 up to c = 6; the same data-set represented in figures A.1 to A.6 (top panel) and figures A.8 to A.13 (bottom panel). Higher number of rounds is not correlated with worse fit when the maximum number of rounds is not reached; and when it is reached it is not necessarily an indication of bad fit either, although all instances of bad fit stem from inferences that reached the maximum allowed number of rounds.

  58

Figure 2 . 15 :

 215 Figure 2.15: Application of our inference method to a generally-accepted tree-like human demographic scenario with three modern populations. (a) Inferred number of islands for each modern population. (b) Inferred local size of each island. Shown here are the scenarios with 4 and 5 components c, and all three values of the weight-shifting parameter ω. The bars with the darker color, marked 'restricted' in the legend, correspond to inferences realized with the option of ignoring recent population expansion.

Figure 2 . 22 :Figure 2 . 23 :

 222223 Figure 2.22: Inferred number of islands for the human populations.

  Figure 4.1 for an example.

Figure 4 . 1 :

 41 Figure 4.1: A method for parameterizing spatial structured models. A stepping-stone model is defined by overlaying a grid onto the map of the region of interest. The physical locations of the samples (the red dots in the figure) are assumed to be known, and each one is associated with the nearest deme (the association is represented by a red line connecting the sample with a deme). In this example model, it is possible to specify the number of islands, their connections, and the sampling vector by setting a single parameter: the grid spacing.

•

  Armando Arredondo, Beatriz Mourato, Khoa Nguyen, Simon Boitard, Willy Rodríguez, Camille Noûs, Olivier Mazet and Lounès Chikhi. Inferring number of populations and changes in connectivity under the n-island model. 2021. Heredity, 126(6):896-912. • Armando Arredondo, Simon Boitard, Simona Grusea, Josué Corujo, Camille Noûs, Lounès Chikhi and Olivier Mazet. Computing the exact expected Site Frequency Spectrum for structured populations. 2022. In preparation.

During my doctoral years, I also had the opportunity to pursue other collaborations. The following publications and manuscripts are the results of these exchanges. These results are not included in the present thesis. • Willy Rodríguez, Olivier Mazet, Simona Grusea, Armando Arredondo, Josué M. Corujo, Simon Boitard, Lounès Chikhi. The IICR and the non-stationary structured coalescent: towards demographic inference with arbitrary changes in population structure. 2018. Heredity, 121(6):663-678. • Simon Boitard, Armando Arredondo, Camille Noûs, Lounès Chikhi, Olivier Mazet. Heterogeneity in effective size across the genome: effects on the Inverse Instantaneous Coalescence Rate (IICR) and implications for demographic inference under linked selection. 2021. Submitted. BioRxiv preprint: doi.org/10.1101/2021.06.11.448122. • Helena Teixeira, Jordi Salmona, Armando Arredondo, Beatriz Mourato, Sophie Manzi, Romule Rakotondravony, Olivier Mazet, Lounès Chikhi, Julia Metzger, Ute Radespiel. Impact of model assumptions on demographic inferences -the case study of two sympatric mouse lemurs in northwestern Madagascar. 2021. Submitted. • Beatriz Mourato, Armando Arredondo, Pierre Lacoste, Willy Rodríguez, Khoa Nguyen, Simon Boitard, Olivier Mazet, Lounès Chikhi. Genomic data and demographic inferences under structured models applied to lemurs. 2022. In preparation.

Table 2 .1.: Results

 2

of varying the strategy parameter of the differential evolution algorithm on the speed of convergence of 10 difficult demographic scenarios with c = 4 components.

  Table 2, suggest that these parameters should be left at their default values.

	Parameter	no. 1 no. 2 no. 3 no. 4 no. 5 no. 6 no. 7 no. 8 no. 9 no. 10	Total
	popsize	100	100	100	100	100	100	100	100	100	100	1,000
	recombination	100	100	100	100	100	100	100	100	100	100	1,000
	mutation	100	100	100	100	100	100	100	100	100	100	1,000
	max-iter	100	100	100	100	100	100	100	100	100	100	1,000

Table 2 .2.:

 2 Results of varying the popsize, recombination, mutation and max-iter parameters of the differential evolution algorithm on the speed of convergence of 10 difficult demographic scenarios with c = 4 components.

Table 2 .3.:

 2 Parameter bounds used during the generation of demographic scenarios for validation.

		simulation inference
		min max min max
	n	2	40	2	50
	t i 0.1	50	0.01 100
	M i 0.1	50	0.05 60
	s i	1	1	1	1
	N 10 3 10 3	10 2 10 4

xkcd.com

(Q, L) ← RateMatrix(k, n, M, S)

π 0 ← (0) 1×E * ; π 0 L -1 (α 0 ) ← 1

δ ← ε π 0 1 / kE * κ(Q * )

t 0 ← (0) E * ×1

ζ ← ExpectedScaledAFS(Q, L, π 0 , t 0 , ω, δ)

return ζ/ ζ 1

Acknowledgements

I would like to thank the Université Fédérale Toulouse Midi Pyrénées and the Région Occitanie for providing the funding for this project, and extending it in times of duress.

To the esteemed members of the jury, thank you for such a rich and insightful discussion session, and particularly to professors Hobolth and Achaz for their very nice reports and their truly inspired questions.

We also observe that for real data it may be helpful to run the analyses for a varying number of rounds, since too few rounds may negatively affect the quality of the fit. Also, once a scenario has been inferred, it is advisable, as an additional validation step, to simulate data under the inferred scenario and use SNIF to re-infer the parameters. This is what we do with the real human data in the next section.

In general, for i < k, the coefficient ψ α (i ) represents the contribution to E(ξ i ) for each unit of time the process spends in state α on average (see §3.3). We denote by P the set of all possible contribution vectors, and by P its total number of elements:

In this case, P depends only on the number of samples k (in fact P = E when n = 1), but once again we keep this dependency implicit.

Transitions between states

From this point onwards we use the notation e ij to indicate a matrix of size k × n with null coefficients except for the one in row i and column j, which has a value of 1.

There are two possible types of state transitions: migration events and coalescence events. During a migration event, one lineage in deme j migrates to deme g = j with rate m jg . We collect these rates in the matrix M = (m jg ) with 1 j, g n; j = g.

We convene, for the sake of completeness, that the diagonal elements of M are null.

Given two states α and β, and a weight of i ∈ {1, . . . , k -1}, the process can transition from α to β via migration event if for some j, g ∈ {1, . . . , n} β = αe ij + e ig , (3.2.3) and this transition happens with rate m jg for every lineage that could potentially migrate, of which there are α ij . Algorithm M generates all possible destination states β given any departure state α.

Algorithm M: For evaluating routine MigrationDestinations(M, α), which given the migration rates M and a state matrix α of size k × n, visits all the tuples (β, m) where β is a state satisfying (3.2.3) and m (which may be null) is the corresponding transition rate.

1 for each i in {1, 2, . . . , k -1}:

for each j in {1, 2, . . . , n}:

for each g in {1, 2, . . . , j -1, j + 1, . . . , n}:

We note that Algorithm M visits all potential migration destinations, regardless of whether a transition would be possible given the migration rates M. This behaviour is desired because, in addition to generating the migration rates, we use the algorithm to build the state space E by repeatedly simulating all possible migrations (see §3.2.4 and Algorithm A).

A coalescence happens when two lineages in the same deme j are replaced with their immediate common ancestor, also in deme j. This coalescence happens with a rate that is the inverse of the deme size. We collect the deme sizes in the vector S = (s j ) n j=1 .

Given two states α and β, and two values i , h ∈ {1, . . . , k -1}, the process can transition from α to β via coalescence event if for some j ∈ {1, . . . , n} (3.2.4) and this transition happens with rate 1/s j for every pair of lineages that can coalesce, of which there are α ij 2 if h = i , and α ij α hj otherwise. Algorithm C generates all possible states to which the process can transition by a coalescence from α.

Algorithm C: For evaluating routine CoalescenceDestinations(S, α), which given the deme sizes and a state matrix α, visits all the tuples (β, c) where β is a state satisfying (3.2.4) and c = 0 is the corresponding coalescence rate.

1 for each j in {1, 2, . . . , n}: 2 for each i in {1, 2, . . . , k -1}:

Grouping and sorting the states

We now focus on defining the rate matrix of the Markov process. This is a matrix Q of size E × E, where coefficient (ℓ 0 , ℓ 1 ) contains the rate of transition from state α ℓ 0 to state α ℓ 1 when ℓ 0 = ℓ 1 ; and the diagonal elements are defined such that E ℓ=1 Q(ℓ 0 , ℓ) = 0. To define always be defined regardless of the value of n by placing their associated contribution vector in the first column:

The contribution vectors ψ are obtained by converting the partitions (a 1 , . . . , a ω ) generated by algorithm P to product notation (3.1.7). See lines

The next step is to generate (and add to E ip ) all the states that can be reached via migration from the initial state (3.2.15) using Algorithm M, then all the states that can be reached from those new states, and so on until no new states are being visited or added. This operation is accomplished using a queue data structure, or any generic data structure that implements the FIFO (first in, first out) interface, where elements are added at one end of the queue, and removed from the other end.

For our use, we may represent a queue of states as a tuple of length 0 ω < ∞ with elements in E, and two associated functions: enq and deq, which add (enqueue) and remove (dequeue) elements from the queue, respectively. More precisely, we say that U ∈ E ω is a queue of size ω; and for U = ∪ ω 0 E ω we define the maps

For a formal treatment of the queue abstract data structure (and other common data structures in computer science), see [START_REF] Lehmann | Algebraic specification of data types: A synthetic approach[END_REF].

Each time a state β is found to be reachable via migration from state α, we accumulate the corresponding migration rate m in coefficient A ip (L ip (α), L ip (β)) of the migration sub-block. The index mapping functions L ip are also progressively constructed. We represent them as subsets of the cartesian product E ip × {1, . . . , E ip }, which contains tuple (α, ℓ) if and only if L ip (α) = ℓ. Building these data structures is sufficient to fully define the state space, since the state spaces E ip are simply the domains of the maps, which we reference using the notation E ip = D L ip . Algorithm A summarizes the proposed procedure for this stage.

Implementation Notes

An efficient implementation of Algorithm A should not repeatedly reallocate and rebuild matrix A ip as indicated in step [11]. Instead, the sparsity of these matrices should be exploited by using a space-efficient representation like the DOK (dictionary of keys) format, for instance. In such cases, step [11] simply adds the key-value pair ((ℓ 0 , ℓ 1 ), m) to the dictionary. Coalescence stage. The goal for this stage is to compute the coalescence sub-blocks B ip , and consequently update the diagonal entries of sub-blocks A ip .

Algorithm Q: For evaluating routine RateMatrix(k, n, M, S), which given the model parameters, builds the rate matrix of the Markov process and the associated data structures.

for each j in {1, . . . , ω}:

for each q in {1, . . . , P i+1 }:

for each p in {1, . . . , P i }:

Computing the expected SFS

In this section we show, given the Markov process defined previously, how to compute the expected SFS for an arbitrary model structure and sample configuration. This computation is reduced to a linear algebra problem for which we propose an efficient solution in Algorithm X.

In the infinite-sites model, all mutations occur on unique sites, and as such the number of mutations is equal to the number of segregating sites. When a mutation occurs in a lineage with i descendants in the sample (a lineage of weight i ), the corresponding segregating site has i copies of the mutated allele. Additionally, since we model mutations as a Poisson process of intensity θ along each branch of the genealogy, the total number of mutations that occur in lineages of weight i has an expected value of θ times the expected total length of such lineages, which we denote as ζ i :

In the example of Figure 3.1, ζ 1 , ζ 2 and ζ 3 would be the total length of red, blue and green branches, respectively. a matrix A ∈ R n×n satisfying the conditions:

• A(i , i ) > 0 for all i = 1, . . . , n;

• A(i , j) 0 for all i = j;

• A is non-singular, with A -1 (i , j) 0 for all i and j.

It is quite straightforward to verify that all our system matrices A ip are M-matrices. Indeed, all diagonal elements of -Q * are positive, all the off-diagonal elements are non-negative, and all coefficients of (-Q * ) -1 = G are also non-negative since they can be interpreted as expected permanence times.

Algorithm X summarizes our proposed method for computing the expected SFS ξ given the rate matrix Q of the Markov process and an initial state distribution π 0 . Algorithm X: For evaluating routine ExpectedScaledAFS(Q, L, π 0 , t, ω, δ), which given the data structures associated with the Markov process and an anterior distribution of states π 0 , updates an approximate solution t to problem (3.3.4) using a SOR(ω) iterative scheme until res ∞ < δ, and returns que corresponding value of ζ.

for each ℓ in {1, 2, . . . , E ip }:

for each h in {1, 2, . . . , k -1}:

Implementation Notes

The read and write operation in line [13] must be implemented as an atomic operation due to the fact that ζ is shared across parallel threads. Otherwise, the values of ζ(h) could become corrupted as one thread reads a value, and another thread immediately writes another incompatible value to it. If atomic operations are not available, then lines [11 ••• 13] must be moved to a separate, non-parallel, p loop. increasingly-accurate approximations of E( ξ).

Algorithm R':

For evaluating an overload of routine ExpectedAFS(k, n, M, S, v , ω, d, ε) allowing early termination.

We can see in Figure 3.10 the tradeoff between computation time and error introspection granularity in an n-island model. In general, there are diminishing returns with values of d larger than 10, although the exact optimal value depends heavily on the application. In all cases, the fastest results are achieved with the fewer introspection steps.

Model specialization: Symmetrical n-island

Despite the efficient algorithms presented in the previous sections, the state space grows too quickly and becomes computationally intractable even for small values of k. For this reason, we believe that the most sensible application for these methods is through model specializations. A model specialization is a way of compressing the state space by taking advantage of the symmetries present in a specific demographic model. In this section we use the n-island model as a case study for introducing model specializations.

The n-island model represents the worst case scenario regarding the size of the state space. Indeed, since all islands are reachable from all other islands, the migration matrix M is completely dense, which in turn implies that the transition rate matrix Q is the largest and least sparse possible. On the other hand, the n-island model is the most symmetrical demographic model possible: all islands have the same sizes and connectivity, thus they are only distinguishable through the number and configuration of lineages present. Consider for instance the following states in a symmetrical 3-island model with k = 4 sampled lineages:

The states in (3.4.1) are different only because we have numbered the islands, but since the islands are indistinguishable, they can all be grouped in the same equivalence class. This equivalence class is characterized by one of the islands having two lineages of weight 1 and another island having one lineage of weight 2. This class can be represented by a state matrix present.

Algorithm sC: For evaluating routine SymmetricalCoalescenceDestinations(s, α), which given the deme size s of a symmetrical island model and a state matrix α, visits all the tuples (β, c) where β is a class-representative state satisfying (3.2.4) for this model, and c is the corresponding coalescence rate.

1 for each j in {1, 2, . . . , n}: 2 for each i in {1, 2, . . . , k -1}:

The simplification of the Markov system provided by the n-island model specialization goes beyond simply reducing the number of states. Indeed, since now all migration events happen with the same rate m, and all coalescence events with the same rate 1/s, the coefficients of the rate matrix also become greatly simplified. We can observe in Algorithm sM (line 6) that the coefficients of the migration blocks A i are of the form ℓ,ι α ℓ,ι m = m C ℓ,ι . Similarly for the coalescence blocks B i (see algorithm sC, lines 7 and 12), the transition rates are of the form 1 s C ℓ,ι . Consequently, we can decompose the rate matrix Q of a symmetrical n-island model in the following form:

The matrices Q m and Q c hold the coefficients C ℓ,ι for the migration and coalescence blocks, respectively, and they are invariant given k and n fixed. We can use this representation in order to greatly reduce the computational cost of building the rate matrix (algorithm Q) for the same values of k and n but many different values of m and s. See figure 3.4 for a comparison of the computational resources required to execute algorithm Q and those required to build the rate matrix given Q m and Q c . All time and memory benchmarks presented in this and the next sections were performed on a Ryzen 3700X desktop system with 32GB of RAM. As can be seen in Figure 3.4, the RAM capacity was the eventual limiting factor in how big a sample could be processed (k = 26 in this case). We note that larger samples could be processed by having systems with larger main memory compute the symbolic rate matrix representation remotely and then distribute this file to the local system where the analyses are being conducted. The size of this file, and

Conclusions and future work

In this final section of the thesis we summarize the main results from the previous chapters in a forward-looking context, where the natural steps for the continuation of the research are outlined.

The IICR and the inference of structure

In Chapter 2 we presented an inference framework that is able to accurately infer structure parameters (number of islands and their sizes) within a symmetrical island model given an IICR estimate like the PSMC. It is also able to date up to five events of changes in migration rate (i.e., six components) with good precision and consistency, as long as the underlying model is compatible with a symmetrical island model with constant population size.

A comprehensive validation effort using simulated data support these results. Indeed, the normalized RMSD (nRMSD) (Figure 2.4) of the simulated vs. inferred scenario parameters is zero for stationary scenarios (c = 1), and increases linearly with the number of components. The M parameters reach a value of about 0.5 at six components, and we see that the first (most recent) and last (most ancient) components are better estimated than the middle ones. It is likely that these middle components are being affected by what we called the component misidentification phenomenon, where the migration rate for a given component is being inferred for a neighboring one (previous or following), thus contributing to a negative correlation between the simulated and inferred migration values for a given component. This can happen when, for instance, a given demographic event (with an associated change of migration rate) does not have a large-enough effect on the shape of the IICR curve. Regardless of this behavior, the number of islands n and the reference effective size N are consistently well estimated, reaching an nRMSD of about 0.1 in the worst cases. The event time parameters (t i ) exhibit the worst nRMSD values, varying between 1 and 2 in the worst cases. Although in these cases, the fact that time is log-spaced and spans several orders of magnitude causes outliers to have a disproportionate contribution to this statistic.

One novel aspect of our approach is that the number of demes gets inferred as one of the model parameters, and it is in fact the best estimated parameter, which is in agreement with [START_REF] Mazet | Demographic inference using genetic data from a single individual: Separating population size variation from population structure[END_REF] that used information from the distribution of T 2 values and a likelihood approach. These authors however, only analysed stationary models. Here we found that other parameters were also well-estimated when the number of components was low, but we These findings suggest that tree models fitted with the SFS, like those considered here, do not offer a definitive answer to the question of human evolution and other families of models should be explored [START_REF] Goldstein | Human migrations and population structure: what we know and why it matters[END_REF][START_REF] Scerri | Did our species evolve in subdivided populations across Africa, and why does it matter?[END_REF][START_REF] Scerri | Beyond multiregional and simple out-of-Africa models of human evolution[END_REF]. It remains to be seen however how well models inferred with our method fit the real SFS of their respective human populations. A more general treatment of this question is a good candidate for future research. However, in §2.4.3 we compared the SFS of a sample of 216 humans from the Yoruba population [START_REF] Lapierre | Accuracy of demographic inferences from the site frequency spectrum: the case of the yoruba population[END_REF] to the one inferred by the GADMA method from [START_REF] Noskova | GADMA: Genetic algorithm for inferring demographic history of multiple populations from allele frequency spectrum data[END_REF] and the one corresponding to three variations of the inferred demographic model by our method (see Figure 2.25). These simulations suggest that existing SFSs could be easily fitted with a structured model similar to those inferred by SNIF, but in which we would allow for a recent population size change.

Computing the exact SFS for structured populations

The SFS has long been an important statistic due to its capacity to succinctly summarize a large amount of information about the genetic diversity of a sample and the ease with which it can be computed using real data.

Comparable approaches to utilize the SFS with structured populations [START_REF] Gutenkunst | Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data[END_REF][START_REF] Excoffier | Robust demographic inference from genomic and snp data[END_REF][START_REF] Kern | Exact calculation of the joint allele frequency spectrum for isolation with migration models[END_REF] focused on what is known as the joint SFS, where information regarding the origin of the lineages is preserved along with the frequency of the mutations. This approach limits the number of populations that can be modeled, since a high-dimensional joint SFS is both computationally intractable and uninformative in practice, since one would expect that most of the coefficients would be null due to lack of information.

In this work we show that by focusing on the aggregate SFS of structured populations, it is possible to have a summary statistic that allows for complex models of population structure and is sensitive to the model parameters (see figures 3.6 and 3.5). We also show that by using phase-type theory, it is possible to exactly compute this SFS for moderate sample sizes very quickly. We did this by introducing a generalization of Herbots's structured coalescent that tracks information regarding the ancestry of the lineages as well as their location in the demes. The challenge of the approach is tied to its computational requirements, since the state space of the Markov process grows almost exponentially with the sample size. To combat this effect, we recommend the usage of model specializations. These are Markov processes that are derived from the general process described in §3.2 that include state-reduction optimizations specific to the demographic model in use. These state reductions are achieved by exploiting the symmetries often exhibited in structured models. We explore the n-island model specialization in §3.4. One obvious way to continue this line of research is to explore other model specializations. A great candidate would be the stepping-stone family of models, which would enable us to draw more spatial-aware conclusions from the potential demographic inference applications.

The performance results shown in figures 3.3 to 3.7 are focused on the n-island model specialization. We do not know yet how these techniques will perform under different models. For each parameter of the demographic scenarios, we counted the instances where the inferred value was within a certain percent of the simulated one. We say that an inferred parameter p is within x% of the simulated value p if -x p < 100( pp) < xp. Note that not all parameters are present in all scenarios (for instance, the parameter M 4 only appears in scenarios with c 5 components).

Appendix B.

Additional results of validation using T-sim IICRs

In this section we show the validation results of simulating and then inferring from 100 randomly generated demographic scenarios with scaled IICRs and varying number of components c. In all cases, the scenario parameters were drawn from the following finite sets:

n ∈ {2, 5, 10, 15, 20}, t i ∈ {0.1, 0.5, 1, 2, 5, 10, 20, 50} ∀i , M i ∈ {0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50} ∀i , (B.0.1)

Unlike in the previous sections, the simulated IICRs are T-sim IICRs, meaning that the values are not exact due to the stochastic nature of the underlying ms simulation. For each value of c from c = 1 to c = 5 we show the aggregate connectivity graph for all the simulations as well as the IICR and parameters of two individual scenarios from the set.

B.1. Scenarios with 1 component