# Likelihood of social-ecological genetic model

Stéphane Dupas, Camilo Patarroyo

► To cite this version:

Stéphane Dupas, Camilo Patarroyo. Likelihood of social-ecological genetic model. 2025. hal-05167034

HAL Id: hal-05167034

https://hal.science/hal-05167034v1

Preprint submitted on 17 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1

2

3

4

5

6

7

8

# Likelihood of social-ecological genetic model

Stéphane Dupas¹*, Camilo Patarroyo²,

1 Laboratoire Evolution, Génome, Comportement et Écologie, Univ Paris-Sud, CNRS, IRD 247, 1

Avenue de la Terrasse, 91198 Gif-sur-Yvette Cedex, France

2 Department of Biological Sciences, Universidad de los Andes, Bogotá D.C., Colombia

* E-mail: Stéphane.Dupas@ird.fr

9

# Abstract

10

11

12

13

14

15

16

17

18

19

20

21

22

The genetic structure of populations depends on two parallel processes - genetic and social-ecological - providing mutual information. Models that describe species' responses within social-ecological systems are increasingly important in the context of modern environmental crises. Advances in genetic data collection now provide access to vast numbers of markers, enabling a more comprehensive understanding of species dynamics. However, current statistical inference methods do not fully integrate these processes into a single cohesive model, and rely on summary statistics derived from separate inferences and simulations.

In this work, I propose a probabilistic framework based on a coalescent model to compute the likelihood of a demogenetic model represented as connected graphs. The population graph, linking population nodes, is characterized by a backward gene adjacency matrix, which describes the probability of gene origins and is influenced by niche and dispersal functions. The genetic graph, linking allele nodes, captures mutation probabilities. A third graph represents the coalescence process of genes within the demographic and genetic graphs.

Coalescence is simulated within the demographic model, and its probability is estimated based on
the genetic model, or vice versa. Likelihood estimation is performed using a Monte Carlo algorithm.

This framework allows likelihood-based sampling and Bayesian inference, offering a robust approach for
environmental niche and meta-population modeling.

I also discuss the practical applications of this model. By combining environmental niche functions
with a coalescent framework, this approach enables probabilistic reconstructions of past species distribu-
tions based on present and historical occurrence data. It can incorporate various data sources, including
historical records, absence data, and genetic information.

By integrating niche and dispersal processes into a unified model, this framework provides a powerful
tool for improving species distribution forecasting and deepening our understanding of species' responses
to environmental change as interconnected systems.

## Author summary

Inferring biodiversity responses to social-ecosystem history using population genetic and environmental
history data is crucial in the context of modern environmental crises. However, existing models do
not fully integrate key processes such as niche modeling, coalescent modeling, and genetic modeling,
often relying on approximate Bayesian computation for inference. In this work, I propose a method to
estimate the likelihood of a fully integrated demo-genetic model that links population genetic data with
socio-ecosystem history.

The proposed method is based on a hidden coalescent process simulated within a social-ecosystem
graph, along with a sampling process represented in a genetic mutation graph. In the social-ecosystem
graph, nodes represent populations, and their effective sizes are estimated using niche models that account
for environmental variables. The edges of the graph represent effective migration, which is estimated
through connectivity models. The resulting likelihood function can be used for Bayesian or maximum
likelihood inference of genetic changes within socio-ecosystems.
This approach offers a flexible framework for inferring demo-genetic models across various fields, including phylogeography, agrobiodiversity evolution, ecology, and epidemiology. It can be applied for forecasting or hypothesis testing, providing a comprehensive tool for understanding genetic dynamics in relation to socio-ecosystem processes.

## Introduction

In the context of the modern ecological crisis, understanding species' responses to socio-ecosystem changes has become one of the most fundamental scientific challenges in ecology. Ecological niche theory is useful for predicting species' responses; however, realized niches also depend on species interactions and historical dispersal events. Genes remain the fundamental unit of information transmission, making them particularly relevant for studying species' responses in social-ecological systems [1]. The increasing availability of molecular data provides valuable complementary insights into realized niches.

One approach to integrating genetic information focuses on 'omics, examining emergent properties from individuals to populations and ecosystems [2]. Other approaches emphasize population-level processes, both adaptive and neutral. Population genetic structure is shaped by adaptive (selection) and neutral forces (drift, mutation, migration). Many genetic and genomic studies now focus on documenting adaptive processes, particularly those related to global change and niche differentiation between populations [3]. A promising research direction is to consider neutral forces, including dispersal processes and historical events [4]. The presence of a gene at a given site is determined both by its environmental niche and its ability to reach that location through spatial and socio-ecological networks.

A significant advance in recent decades has been bridging the gap between population genetics and ecological niches. Niche models provide insights into past abundance and migration patterns. Drift is influenced by reduced population size and can be linked to environmental variables through niche functions. Existing models generally separate inference into two steps: first, estimating species distribution param-
eters from occurrence data, and second, inferring genetic evolution parameters from genetic data [5].

The term Integrated Distribution Demographic and Coalescence Model (IDDC) [6] was introduced to describe these efforts. Several simulation tools now allow genetic parameter inference within discrete or continuous species distribution landscapes [6–10]. Most of these models use niche modeling to simulate species distribution, construct evolutionary force landscapes, and perform coalescence simulations using Approximate Bayesian Computation (ABC). Many spatially explicit tools are available [9], with SPLATCHE [7, 11, 12] being one of the most advanced. This software enables the simulation of genetic data based on historical demographic maps and introduction scenarios. However, these IDDC approaches still present some limitations.

The first major limitation is the lack of full integration of niche, demographic, and genetic processes into a single model. In existing IDDC models, demographic parameter maps are not inferred but rather approximated using species distribution maps. There is no inference algorithm that simultaneously integrates demographic parameters, environmental relationships, and coalescent processes. Separating demographic and genetic inference is problematic because genetic drift does not depend solely on absolute population size but also on effective population size and migration rates [4]. Additionally, niche models used as proxies for historical reconstruction are typically based on occurrence or presence-absence data, which themselves depend on historical contingencies such as migration and extinction. Thus, there is no justification for treating niche and historical processes separately in demogenetic inference. Integrating niche and genealogical modeling would address these issues. QuetzalCoalT [13] is currently the only available tool capable of performing integrated inference of genetic and niche-based processes from environmental history and genetic structure. Recognizing that demographic and genetic processes should be fully integrated, we propose the term environmental demogenetics, as the ultimate inference model should incorporate environmental, demographic, and genetic data. Furthermore, most IDDC and environmental demogenetic models fail to integrate social-ecological networks, which are increasingly crucial in the Anthropocene era (e.g., agricultural practices, ethnic groups, trade routes).
The second major challenge of IDDC approaches is performing inference in the context of complex historical and high-dimensional processes. The relationship between environmental and genetic processes is typically a high-dimensional, unknown coalescent process. Consequently, models rely on simulations, with inference performed using ABC on summary statistics. ABC inference compares observed data with data generated by random simulators, making it particularly useful for Bayesian inference when likelihoods are intractable or highly complex. This is a common situation in population biology, where ecological and genetic data are the result of multiple unobservable processes—including genealogy, ecological interactions, genetic evolution, and sampling—that are highly interdependent and contingent. In ABC, priors are sampled based on summary statistics rather than exact likelihood calculations. However, this approach has drawbacks: information may be lost when transforming data into summary statistics, and there is no universal method for selecting appropriate summary statistics. Each inference problem requires a custom selection, and it is impossible to definitively prove that the chosen summary statistics capture all relevant aspects of the modeled process [14]. These limitations do not exist in likelihood-based inference.

To address these challenges, we propose a novel approach that:

1. Integrates social-ecological history, demographic processes, and genetic processes within a single inference model.

2. Directly incorporates multiple data sources, including genetic data, ecological occurrence data, and environmental history.

3. Accounts for both spatial and non-spatial networks, including social-ecological networks.

4. Performs exact likelihood-based probabilistic modeling, rather than relying on approximations.

Our metamodel consists of three interconnected graphs:

• A genetic graph, where nodes represent alleles and edges represent mutations.
117 • A population graph, where nodes represent populations and edges represent migration events.

118 • A coalescent graph, where nodes represent individuals and edges represent lineages.

The coalescent process serves as a fundamental link between genetic and demographic sampling processes. We propose two approaches for likelihood sampling within this metamodel using genetic status and demographic location data. First, coalescence can be simulated based on sample locations within the population graph, and the probability of the coalescent, the genetic data and graph model can be calculated. Alternatively, genealogy can be simulated within the genetic graph, and the probability of observed demographic locations can be calculated for the coalescent. This framework enables both likelihood-based and Bayesian inference. By integrating potential niche parameters (related to scenopoietic variables) with historical processes shaping the realized niche, this approach should enhance our understanding of species distributions and improve future projections (Figure 1).

## 128 Demo-genetic Model

I present below the variables and parameters of the different graph models conforming the demogenetic model. The **population graph** represents populations, the **genetic graph** represents alleles and the **coalescence history graph** represent sample coalescing lineages (Table 2, figures 1,2,3).

The population graph, noted $G_P = (V_P, E_P)$.

I first consider demographic history in a general random walk across a population graph model, representing populations as nodes (in spatial or ecological dimensions) and migration as edges (Figure 3). The population graph can represent a geographic landscape ($G_P^L$), a socio-ecological network ($G_P^{SE}$), or a combination of both $G_P^{SE+L}$ (Figure 3). The variables characterizing populations can also vary over time. Additionally, the migration function can depend on space and environment (Figure 3).

139 • Each node $V_P[i]$ represents a distinct population with an effective population size $Ne[i]$.
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

$$ \mathcal{L} = \sum_{i=1}^{N} P(H_{sim} | N) $$

**Figure 1.** Social ecological demogenetic likelihood computation. Two approaches are considered.

Approach 1: simulate coalescent from population information and calculate its probability from genetic information. Or Approach 2: simulate from genetic and calculate probability from populations. *Sim*: simulation of coalescent knowing one type of data and model. *Pr*: probability of coalescent knowing the other type of data and model.
140

141

142

143

144

* Each edge $E_P[i, j]$ is characterized by a backward gene transition probability $M_P[i, j]$, which represents the probability that a gene in node $V_P[i]$ originated from node $V_P[j]$ in the previous generation. This probability depends on the effective population sizes of both the receiving and originating demes $N_e$, and the migration rate $m[i, j]$, which is the proportion of individuals migrating from node $j$ to node $i$ in forward time.

145

* The effective population size $N_e$ is defined by a niche model $N_e(X_{i,t}, \theta_P^N)$,

146

* which relates $N_e$ to environmental variables

147

* $X_{i,t}$ and model parameters $\theta_P^N$.

148

149

* The migration rate $m[i,j]$ is defined by a migration rate model $m(E_P[i, j], \theta_P^m)$, which relates $m$ to a distance measure between nodes $i$ and $j$, incorporating migration-friction parameters $\theta_P^m$.

150

151

152

* The number of demes (the order of the graph) is denoted as $|D|$, while the number of migration connections (the degree of the graph) is denoted as $|E_P|$. The matrix $M_P$ represents the backward transition probabilities between populations.

153

154

**The genetic connected graph, denoted as $G_G = (S_G, E_G)$:** We now consider the coalescence in the genetic graph, where nodes represent alleles and edges represent mutations.

155

* Each node $V_G[i]$ represents an allele state.

156

157

158

* Each edge $E_G[i, j]$ represents a mutation event that is characterized by the probability $M_G[i, j]$, that a gene in node $V_G[i]$ have mutated from node $V_G[j]$ in the previous generation. The matrix $M_G$ represents mutation rate between alleles.

159

160

* The number of alleles (the order of the graph) is noted $|V_G|$ while the number of mutation connections (the degree of the graph) is denoted $|E_G|$.
The coalescence history graph, denoted $H_C(S)$:

$H_C(S)$ is a directed graph that traces the lineage of a genetic sample $S$ back to their common ancestor, which edges oriented opposite to those in a traditional genealogical graph (Figure 2).

* Each node $V_C[i]$ represents either a sample $S$ or a coalescence point.

* Each edge $E_C[u, v]$ represents a lineage connecting a node to its next coalescence point.

* The coalescence history graph is defined as $HC=(VC,EC)$ $H_C = (V_C, E_C)$, where $E_C = \{(u, v) \in V_C, \text{ where } v \text{ is parent of } u\}$.

* Sampled individuals $S_i = V_C[i]$ with $i = \{1..n\}$ are defined as nodes with zero in-degree. These individuals are linked to their ancestors through coalescence events, which correspond to internal nodes *In*. Internal nodes are characterized by having more than one in-degree and exactly one out-degree, except for the most recent common ancestor (MRCA), which has more than one in-degree and zero out-degree.

The individual coalescence graph, denoted $H_I(S)$ (Figure 2). The individual coalescence graph $H_I(S)$ is an extended graph that includes the coalescence history graph and all individuals along the edges of the coalescence history graph that constitute its lineages. It is defined as $H_I = (I, E_I)$. Such that $V_C \subset I$, and $E_I = \{(u, v) | v \text{ is the parent of } u \text{ in the previous generation}\}$

* Each node represents an individual kin the genealogy of the sampled genes. It is either a node from $H_C(S)$, or an individual positioned along the vertices $V_C(S)$ between two connected nodes of $H_C$.

* $H_I(S)$, as $H_C(S)$ is an acyclic graph oriented from the sampled individuals to the common ancestor. Each individual $I[i] \in I$ has a unique parent formally expressed as $\text{card}(\text{Adj}(I[i])) = 1$. The descendants of $I[i]$ are obtained by $\text{Adj}^{-1}(I[i])$.

* Each individual of $H_C$ is fully defined within the coalescent tree by one of its descendant, including sampled descendant $S_i$ and the generation time (Figure 2).
<figure><img src="image_4.png" /><figcaption>**Figure 2.** The coalescent graph $H_C$ (green nodes) and the individuals coalescent graph $H_I$ (yellow nodes). $S$: Sampled individuals. $In$: Internal coalescence nodes. $V_I[i]$: Nodes representing individuals within to the coalescent lineages. $H_C = (V_C, E_C)$ consists of vertices $V_C = S \cup In$ and edges, $E_C$, which define the lineages connecting them. The individuals along these edges belong to $H_I$. $t$: generation time.</figcaption></figure>

The probability of a coalescence history graph $H_C$ depends on both the genetic graph and the population graph. To estimate the likelihood of the demogenetic model given environmental, demographic, and genetic data, we will simulate coalescence history graphs within the population graph and compute their probabilities under the genetic graph.

All the variables and terms used in this paper are described in the table 2.

## Simulation of genealogy knowing demographic model

The first method of likelihood calculation consists of simulating a genealogy of the sample in a population graph, and then evaluating its probability in the genetic graph. The species observation data used in the model may include genetic information on the individuals reported, but this is not always required.

To analyze the movements and coalescence of gene lineages within the population graph, we first compute the **population transition matrix**, $M_P$. This matrix is a Markov transition matrix that estimates the evolution of occurrence probabilities among demes. Each cell $[i, j]$ represents the backward in-time probability, that a gene in deme (or node) $i$ originates from deme (or node) $j$.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Examples of population graphs and a coalescence history graph simulated within a population graph. (A.1) Social-Ecological Network Demographic Graph ($G_P^{SE}$) with 4 nodes, where nodes represent ecological or social niches. (A.2) Backward population transition matrix of $G_P^{SE}$.</figcaption></figure>

(B.1) Landscape Demographic Graph ($G_P^L$) with 6 nodes, where nodes correspond to raster cells or vector map contours. (B.2) Backward population transition matrix among nodes in $G_P^L$. (C.1)

Combined Social-Ecological and Landscape Demographic Graph ($G_P^{SE+L}$), integrating $G_P^{SE}$ and $G_P^L$.

$V_P[i]$: Node (deme) i. $Xe_i$: Environmental variables associated with node $V_P[i]$. $Ne(X_i)$: Effective population size of node $V_P[i]$. (B.3) Construction of the 24×24 backward gene population transition matrix $M_P$, characterizing the edges $E_P^{SE+L}$ of $G_P^{SE+L}$ for coalescence simulations. (C) Movement and coalescence of a sample $S = \{S_1, S_2, \dots\}$ within $G_P^{SE+L}$. $H_C$ represents the coalescence history graph on $G_P^{SE+L}$, where $H_C = (V_C, E_C)$ consists of vertices $V_C = S \cup I_C$ (with $I_C$ as internal nodes of $V_C$) and edges $E_C$, representing the lineages connecting them.
<table><caption>**Table 1.** Definition of variables used in the model</caption><thead><tr><th>Category</th><th>Symbol</th><th>Definition</th></tr></thead><tbody><tr><td colspan="3" style="text-align:center; font-weight:bold;">Demographic Graph (G<sub>P</sub>)</td></tr><tr><td></td><td>V<sub>P</sub></td><td>Population (node) vector</td></tr><tr><td></td><td>V<sub>P</sub>[i]</td><td>Population i</td></tr><tr><td></td><td>t</td><td>Generation time</td></tr><tr><td></td><td>Ne[i, t]</td><td>Effective population size of population i at generation time t</td></tr><tr><td></td><td>Xe[i, t]</td><td>Vector of environmental variables for population i at generation time t</td></tr><tr><td></td><td>θ<sup>N</sup><sub>P</sub></td><td>Parameters of the effective population size model</td></tr><tr><td></td><td>E<sub>P</sub></td><td>Backward population transition (edge)</td></tr><tr><td></td><td>M<sub>P</sub></td><td>Backward population transition matrix: probability that the father of an individual in population i was from population j</td></tr><tr><td></td><td>M<sub>P</sub>[i, j]</td><td>Transition probability that an individual in population i originated from population j</td></tr><tr><td></td><td>m</td><td>Migration rate vector</td></tr><tr><td></td><td>m[i, j]</td><td>Probability that an individual from population j migrates to population i in the next generation</td></tr><tr><td></td><td>m(E<sub>P</sub>[i, j], θ<sup>m</sup><sub>P</sub>)</td><td>Migration rate model</td></tr><tr><td></td><td>θ<sup>m</sup><sub>P</sub></td><td>Migration rate model parameters</td></tr><tr><td colspan="3" style="text-align:center; font-weight:bold;">Genetic Graph (G<sub>G</sub>)</td></tr><tr><td></td><td>V<sub>G</sub></td><td>Allele vector</td></tr><tr><td></td><td>V<sub>G</sub>[i]</td><td>Allele i</td></tr><tr><td></td><td>E<sub>G</sub></td><td>Backward allele transition data (edge)</td></tr><tr><td></td><td>M<sub>G</sub></td><td>Backward allele transition matrix: probability that an individual in allele i originated from allele j</td></tr><tr><td></td><td>M<sub>G</sub>[i, j]</td><td>Probability that an individual in allele j mutates into allele i in the next generation</td></tr><tr><td colspan="3" style="text-align:center; font-weight:bold;">Coalescent Graph (H<sub>C</sub>)</td></tr><tr><td></td><td>V<sub>C</sub></td><td>Vector of nodes in the coalescence graph</td></tr><tr><td></td><td>S</td><td>Sample vector</td></tr><tr><td></td><td>S<sub>P</sub></td><td>Deme of collection for S</td></tr><tr><td></td><td>S<sub>G</sub></td><td>Alleles of S</td></tr><tr><td></td><td>IN</td><td>Internal node vector</td></tr><tr><td></td><td>E<sub>C</sub>[u, v]</td><td>Lineage edge from node u to its coalescence ancestor v = Adj(u)</td></tr><tr><td></td><td>p(IN[i] | θ<sub>P</sub>)</td><td>Probability of a coalescence event given demographic models and data</td></tr><tr><td colspan="3" style="text-align:center; font-weight:bold;">Full Coalescent Graph (H<sub>I</sub>)</td></tr><tr><td></td><td>I</td><td>Vector of individuals in the full coalescence history</td></tr><tr><td></td><td>t(I)</td><td>Generation times of individuals I in H<sub>I</sub> or H<sub>C</sub></td></tr><tr><td></td><td>An(I, t)</td><td>Set of ancestors of individuals I at time t (duplicates removed)</td></tr><tr><td></td><td>P<sub>P</sub>(I, t)</td><td>Matrix of probability distribution of An(I, t) across demes in G<sub>P</sub> at time t</td></tr><tr><td></td><td>P<sub>P</sub>(I)</td><td>Matrix of demographic distribution of I at t(I)</td></tr><tr><td></td><td>P<sub>G</sub>(I, t)</td><td>Matrix of probability distribution of An(I, t) across alleles in G<sub>G</sub> at time t</td></tr><tr><td></td><td>P<sub>G</sub>(I)</td><td>Matrix of genetic distribution of I at t(I)</td></tr></tbody></table>
<table><caption>**Table 2.** Signification of the variables</caption><thead><tr><th>Demographic</th><th>G<sub>P</sub></th><th>Demographic graph</th></tr></thead><tbody><tr><td></td><td>V<sub>P</sub></td><td>Population or node vector</td></tr><tr><td></td><td>V<sub>P</sub>[i]</td><td>Population i</td></tr><tr><td></td><td>Ne</td><td>Effective population size vector</td></tr><tr><td></td><td>t</td><td>Generation time</td></tr><tr><td></td><td>Ne[i, t]</td><td>Effective population size of population i at generation time t</td></tr><tr><td></td><td>Xe[i, t]</td><td>Vector of environmental variables of population i at generation time t</td></tr><tr><td></td><td>θ<sup>N</sup><sub>P</sub></td><td>Effective population size model parameters</td></tr><tr><td></td><td>E<sub>P</sub></td><td>Backward population transition or edge data variable</td></tr><tr><td></td><td>M<sub>P</sub></td><td>Edge Backward population transition matrix that individual in population i comes from population j</td></tr><tr><td></td><td>M<sub>P</sub>[i, j]</td><td>Backward population transition probability that individual in population i comes from population j</td></tr><tr><td></td><td>m</td><td>Migration rate edge vector</td></tr><tr><td></td><td>m[i, j]</td><td>Probability that an individual belonging to population j migrates to population i in the next generation</td></tr><tr><td></td><td>m(E<sub>P</sub>[i, j], θ<sup>m</sup><sub>P</sub>)</td><td>Migration rate model</td></tr><tr><td></td><td>θ<sup>m</sup><sub>P</sub></td><td>Migration rate model parameter</td></tr><tr><th>Genetic</th><th>G<sub>G</sub></th><th>Genetic graph</th></tr><tr><td></td><td>V<sub>G</sub></td><td>Allele vector</td></tr><tr><td></td><td>V<sub>G</sub>[i]</td><td>Allele i</td></tr><tr><td></td><td>E<sub>G</sub></td><td>Backward population transition or edge data variable</td></tr><tr><td></td><td>M<sub>G</sub></td><td>Edge Backward allele transition matrix that individual in allele i comes from allele j</td></tr><tr><td></td><td>M<sub>G</sub>[i, j]</td><td>Probability that an individual belonging to allele j mutates to allele i in the next generation</td></tr><tr><th>Coalescent</th><th>H<sub>C</sub></th><th>Sample coalescence history graph</th></tr><tr><td></td><td>V<sub>C</sub></td><td>Coalescence graph node vector</td></tr><tr><td></td><td>S</td><td>Sample vector</td></tr><tr><td></td><td>S<sub>P</sub></td><td>Deme of collection of S</td></tr><tr><td></td><td>S<sub>G</sub></td><td>Alleles of S</td></tr><tr><td></td><td>IN</td><td>Internal nodes vector</td></tr><tr><td></td><td>E<sub>C</sub>[u, v]</td><td>lineage edge from a node u to its coalescence ancestor v = Adj(u)</td></tr><tr><td></td><td>p(Ind[i] | θ<sub>P</sub>)</td><td>the probability of a coalescence event knowing demographic models and data.</td></tr><tr><th>Full coalescent</th><th>H<sub>I</sub></th><th>Individual sample coalescence history graph</th></tr><tr><td></td><td>I</td><td>Vector of individual nodes in full coalescence history</td></tr><tr><td></td><td>t(I)</td><td>generation times of individuaks I in the graph H<sub>I</sub> or H<sub>C</sub></td></tr><tr><td></td><td>An(I, t)</td><td>The set of ancestors at time t of the set of individuals I (doubloons being removed)</td></tr><tr><td></td><td>P<sub>P</sub>(I, t)</td><td>Matrix of probability distribution of An(I, t) in line among demes in G<sub>P</sub> in column, at ancestral time t.</td></tr><tr><td></td><td>P<sub>P</sub>(I)</td><td>Matrix of demographic distribution of I at t(I).</td></tr><tr><td></td><td>P<sub>G</sub>(I, t)</td><td>Matrix of probability distribution of An(I, t) in line among alleles in G<sub>G</sub> in column, at ancestral time t.</td></tr><tr><td></td><td>P<sub>G</sub>(I)</td><td>Matrix of genetic distribution of I at t(I).</td></tr></tbody></table>
$$M_P[i, j](X_{:,t}, \theta_P) = \frac{m(i, j, \theta_P^m)Ne(X_{i,t}, \theta_P^N)}{\sum_{k \in \Gamma(j)} m(k, j, \theta_P^m)Ne(X_{k,t}, \theta_P^N)} \quad (1)$$

Where $\Gamma(i)$ represents the set of all demes directly connected to deme $i$, $m(i, j, \theta_P^m)$ is the probability function of migration from deme $i$ to deme $j$ according to the migration model $\theta_P^m$, and $Ne(X_{i,t}, \theta_P^N)$ is the effective population size given by a niche model that depends on environmental variables and parameters. Note that $m(d(i, i), \theta_P^m)$ corresponds to the probability that genes remain in the same deme. Additionally, the matrix satisfies the condition $\sum_j M_P[i, j] = 1$, ensuring that transition probabilities sum to one.

We will consider coalescence events sequentially, from the most recent sampling date back to the MRCA, proceeding by pairs of lineages. To represent coalescence nodes, we observe that the function $An(S[i], t)$ is an injection from $S \times t(I)$ to $I$: each individual can be uniquely identified in $G_I$ by one of its descendant, in particular in the sample, and its generation time.

Note that if coalescence event occur at time $t$ between two lineages $An(I[i], .)$ and $An(I[j], .)$, then for any earlier time $\tau$ such that $\tau \le t$, it follows that: $\Rightarrow An(I[i], \tau) = An(I[j], \tau)$.

If no coalescence occur, the evolution of the, matrix $P_P(S, t)$ from time $t$ to $t-1$ is given by :

$$P_P(S, t-1) = P_P(S, t) \cdot M_P(t) \quad (2)$$

which describes the transition between generations

Similarly, if no coalescence occurs, the general term at time $t$ is :

$$P_P(S, t) = P_P(S[i], t(S[i])) \prod_{\tau=t(S)}^{t} M_P(\tau) \quad (3)$$

Consider two ancestors $An(S[i], t)$ and $An(S[j], t)$ that have not yet coalesced at generation $t$. The probability that they coexist in the same deme in the previous generation $t-1$ is given by the scalar product of their population distribution probabilities at generation $t$:

$$p(An(S[i], t) \text{ and } An(S[j], t) \text{ coexist in a deme at } t-1) = (P_P(S[i], t) \cdot M_P(t)) \cdot (P_P(S[j], t) \cdot M_P(t))^{\top} \quad (4)$$
The probability that they coalesce in the previous generation (t − 1) is obtained by dividing the vector of coexistence probabilities across demes by twice the effective population size predicted by the niche model in each corresponding deme N$_{e}$ ( X ( t − 1)):

$$
\begin{align}
& p(\text{An}(S[i], t) \text{ and } \text{An}(S[j], t) \text{ coalesce at } t-1) \tag{5} \\
& = (P_P(S[i], t) \cdot M_P(t)) \cdot (P_P(S[j], t) \cdot M_P(t) \oslash (2Ne(X(t-1))))^{\top} \nonumber
\end{align}
$$

Where ⊗ is the element-wise division operator.

If two lineages coalesce, the probability distribution of their common ancestor given the demographic model, is given by performing an element-wise multiplication of the distribution probabilities of the two lineages, followed by an element-wise division by the population size. The result is then normalized by the probability of coalescence.

$$
P_P(An(S[{i,j}], t) | \theta_G) = \frac{\left( P_P(S[i], t) \cdot M_P(t) \right) \odot \left( P_P(S[i], t) \cdot M_P(t) \oslash (2Ne(X(t-1))) \right)}{\left( P_P(S[j], t) \cdot M_P(t) \right) \cdot \left( P_P(S[j], t) \cdot M_P(t) \oslash (2Ne(X(t-1))) \right)^{\top}} \quad (6)
$$

Where:

bla ⊙ represents element-wise

bla ⊗ represents element-wise division.

Using these equations of genetic distribution and coalescence probability from one generation to the parent, we can simulate their genealogy using the following algorithm.
<caption>**Algorithm 1** Demographic simulation of a coalescent</caption>

1. Initialize generation time:
   Define $t = \max\{t | \text{card}(An(S, t)) \ge 2\}$, i.e., the most recent time when at least two lineages exist.
2. Compute the genetic distribution:
   Calculate the genetic distribution $P_P(S, t)$ using equation (3).
3. Simulate coalescence events in pairs:
   (a) Randomly select a pair of lineages in the set $An(S, t)$.
   (b) Simulate their coalescence at $t - 1$ with the probability defined in equation (5).
4. Update lineages if coalescence occurs:
   (a) If the selected pair coalesces, merge their lineages.
   (b) Update their demographic distribution using equation (6).
   (c) Record the event in a genealogy object, including the probability distribution of each remaining lineage.
5. Repeat for all remaining pairs:
   Return to step 3 until all pairs in $An(S, t)$ have been tested for coalescence.
6. Move to the previous generation:
   Decrease the time by setting $t = t - 1$.
7. Repeat until complete coalescence:
   Return to step 2 until all lineages have merged into a single ancestor.
## Probability of a genealogy knowing genetic model

In the first approach, we consider a coalescent simulated under a known demographic model, as described in Algorithm 1.

Let us consider a coalescent graph $H_C$ simulated under a given demographic model. We calculate the probability of each coalescence point under genetic model $p(\text{In}[k] | \theta_G)$, characterized by $t(\text{In}[k]) = t_C$ and its descendant nodes denoted as $\text{In}[k]^d = \text{Adj}^{-1}(\text{In}[k])$. The number of descendant lineages coalescing at $\text{In}[k]$ is given by $\text{In}[k]^{Nd} = \text{card}(\text{In}[k]^d) = \text{card}(\text{Adj}^{-1}(\text{In}[k]))$. This number can be 2 in the case of a dichotomy or greater in the case of a polytomy. These coalescence events occur at times $t(\text{In}[k]^d)$.

If the lineages have not coalesced at any more recent time, the probability that they co-occur in the same genetic state (i.e., in the same population) at $t | t \le \min(\text{In}[k]^d)$, given the genetic model $\theta_G$, is given by the scalar product of the probability distributions of each individual in $\text{In}[k]^d$ at time $t$ (obtained from equation (3) applied to individuals in $\text{In}[k]^d$).

$$ P(\text{In}[k]^d \text{ coexist in genetic space at } t_C \mid \text{ they did not coalesce more recently}) = \prod_{i=1}^{\text{In}[k]^{Nd}} \cdot \left( P_G(\text{In}[k]^d[i]) \prod_{\tau=t(\text{In}[k]^d[i])}^{t_C+1} M_G(\tau) \right) \quad (7) $$

Where "$\prod$" represents the scalar product of the vectors.

The probability of coalescence is proportional to this probability (as it depends on demographic data and model, and is independent on genetic data and model). If the individuals have not coalesced at a more recent generation, the probability that two individuals coalesce at generation $t_C = t(\text{In}[k])$ is given by the scalar product of their probability distributions at generation $t_C$, multiplied by the parameter $\theta_G^\alpha$ (since $\theta_G^\alpha$ is independent of generation time):
$$P(\text{In}[k]^d \text{ coalesce in genetic space at } t_C \mid \text{ they did not coalesce more recently}) = \theta_G^\alpha \times \prod_{i=1}^{\text{In}[k]^{Nd}} \cdot \left( P_G(\text{In}[k]^d[i]) \prod_{\tau=t(\text{In}[k]^d[i])}^{t_C+1} M_G(\tau) \right) \quad (8)$$

The parameter $\theta^\alpha$ will be estimated as a coefficient related to the average effective population size

$N_e$.

The probability that they have not coalesced more recently is given by Equation ??:

$$P(\text{In}[k]^d \text{ did not coalesce in genetic space more recently than } t_C) = \prod_{t=\min(t(\text{In}[k]^d))}^{t_C} \left( 1 - \theta_G^\alpha \times \left( \prod_{i=1}^{\text{In}[k]^{Nd}} \cdot \left( P_G(\text{In}[k]^d[i]) \prod_{\tau=t(\text{In}[k]^d[i])}^{t+1} M_G(\tau) \right) \right) \right) \quad (9)$$

The final probability of coalescence is:

$$P(\text{In}[k]) = P(\text{In}[k]^d \text{ coalesce in an allele at } t(\text{In}[k])) = \theta_G^\alpha \times \left( \prod_{i=1}^{\text{In}[k]^{Nd}} \cdot \left( P_G(\text{In}[k]^d[i]) \prod_{\tau=t(\text{In}[k]^d[i])}^{t(\text{In}[k]^d[i])+1} M_G(\tau) \right) \right) \times \left( 1 - \theta_G^\alpha \times \left( \prod_{i=1}^{\text{In}[k]^{Nd}} \cdot \left( P_G(\text{In}[k]^d[i]) \prod_{\tau=t(\text{In}[k]^d[i])}^{t+1} M_G(\tau) \right) \right) \right) \quad (10)$$

The allele probability distribution vector of the coalescent ancestor $\text{In}[k]$ is given by:

$$P_P(I_C[i]) = \prod_{i=1}^{N(\text{In}[k])} \odot \left( P_G(\text{In}[k]^d[i]) \prod_{\tau=t(\text{In}[k]^d[i])}^{t(\text{In}[k]^d[i])+1} M_G(\tau) \right) \quad (11)$$

Where $\prod \odot$ is the element-wise vector multiplication.

We can calculate the probability of a data sample ($P_G(S)$ or $S_G$) and genealogy $H_C$, given the genetic model $\theta_G$, by multiplying the probabilities from Equation (??) over coalescence events and updating the probability distribution of nodes across the genetic space using the corresponding equation:
$$p(H_C | \theta_G, P_G(S)) = p(H_C | S_G, \theta_G) = \prod_{k=0}^{\text{card}(In)} p(In[k]). \qquad (12)$$

**Algorithm 2** Genetic Probability of a Demographically Simulated Coalescent

1: Initialize $\mathcal{P} = 1$.

2: **while** there are remaining coalescence points in the history graph **do**

3: Select the most recent remaining coalescence point.

4: Update $\mathcal{P} \leftarrow \mathcal{P} \times$ pairwise coalescence probability from Equation (10).

5: Remove daughter nodes $Adj^{-1}(V_C[i])$ from the tree.

6: Set the demic and allelic distribution probabilities of $V_C[i]$ using Equation (11).

7: **end while**

8: **return** $\mathcal{P}$

Note that this algorithm 2 can be combined with the demographic coalescence simulation algorithm (Algorithm 1) so that the coalescence simulation yields directly the genetic probability of the demographic coalescent.

**Algorithm 3: Joint simulation of coalescent on population graph and calculation of the probability on genetic graph**

1. Set generation time $t = \{\max t | \text{card}(An(S, t)) \ge 2\}$. Set $\mathcal{P}_g = 1$

2. Calculate genetic distribution of $P_P(S, t)$ (eq. 3).

3. Select randomly a pair of lineage in the set $An(S, t)$ that remain. Simulate their coalescence or not at $t-1$ with probability defined in eq. 5.

4. if coalescence occurred, merge the lineages, update demographic distribution using eq. 6. Record the event in a genealogy object including the probability distribution location of each remaining lineage. Reset current $\mathcal{P}_g$ multiplying it by the pairwise coalescence probability of this coalescence
269

event defined in (equation 10). Set the demic and allelic distribution probabilities of $V_C[i]$ according to equation (11).

270

271

5. go to 3. until every pair of the set has been tested for coalescence.

272

6. Set $t = t - 1$.

273

7. go to 2. Until all lineages have coalesced

**Algorithm 3** Joint Simulation of Coalescent on Population Graph and Calculation of Probability on Genetic Graph

1: Set generation time $t = \{\max t \mid \text{card}(An(S, t)) \ge 2\}$. Initialize $\mathcal{P}_g = 1$.

2: **while** there are multiple lineages remaining **do**

3: Compute genetic distribution $P_P(S, t)$ using Equation (3).

4: **for** each remaining pair of lineages in $An(S, t)$ **do**

5: Simulate coalescence at $t - 1$ with probability from Equation (5).

6: **if** coalescence occurs **then**

7: Merge lineages and update demographic distribution using Equation (6).

8: Record event in genealogy object with probability distribution location.

9: Update $\mathcal{P}_g \leftarrow \mathcal{P}_g \times$ pairwise coalescence probability (Equation 10).

10: Set demic and allelic distribution probabilities of $V_C[i]$ using Equation (11).

11: **end if**

12: **end for**

13: Decrease generation time: $t \leftarrow t - 1$.

14: **end while**

15: **return** $\mathcal{P}_g$
274 Simulation of genealogy knowing genetic model

The second method of likelihood calculation consists of simulating a genealogy of the sample within the genetic graph and then evaluating its probability in the population graph. We again consider coalescence events from the most recent ancestor to the common ancestor, proceeding in pairs of lineages.

The same equations, namely Equation (2), Equation (3), and Equation (4), can be used to calculate the probability of transition from the genetic distribution at time t to time t − 1 (Equation 13), the general expression of the genetic distribution at any time t (Equation 14), and the probability that two genes, which have not coalesced more recently, share the same genetic state at ancestral time t (Equation 15), respectively.

$$
P_G(S, t-1) = P_G(S, t) \cdot M_G(t) \tag{13}
$$

$$
P_G(S, t) = P_G(S[i], t(S[i])) \prod_{\tau=t(S)}^{t} M_G(\tau) \quad (14)
$$

$$
p(An(S[i], t) \text{ and } An(S[j], t) \text{ coexist in a genetic state at } t-1) = \\
(P_G(S[i], t) \cdot M_G(t)) \cdot (P_G(S[j], t) \cdot M_G(t))^\top \tag{15}
$$

The probability that two lineages coalesce in the previous generation, $t-1$, is obtained by multiplying the vector of coexistence probability in genetic space by a parameter $\theta_G^\alpha$. This parameter characterizes the probability of coalescence for individuals sharing alleles. It reflects the average effective population size, which is independent of the other $\theta_G$ parameters and $P_G(S)$ and, therefore, remains constant. However, this parameter is unknown and must be estimated.

$$
p(An(S[i], t) \text{ and } An(S[j], t) \text{ coalesce at } t-1) = \theta_G^\alpha \times (P_G(S, t)[2] \cdot M_G(t)) \cdot (P_G(S, t)[1] \cdot M_G(t))^\top \quad (16)
$$
If they coalesce, the genetic distribution of the ancestor knowing genetic model is given by element-wise multiplication of the distribution probability of the two lineages divided by the probability of coalescence (eq. 16). We can simplify by $\theta_G^\alpha$ and obtain:

If coalescence occurs, the genetic distribution of the ancestor, given the genetic model, is obtained by performing an element-wise multiplication of the probability distributions of the two lineages and normalizing by the probability of coalescence (Equation 16). Simplifying by $\theta_G^\alpha$, we obtain:

$$P_G(An(S[(i, j)], t) | \theta_G) = \frac{(P_G(S[i], t) \cdot M_G(t)) \odot (P_G(S[j], t) \cdot M_G(t))^{\top}}{(P_G(S,t)[2] \cdot M_G(t)) \cdot (P_G(S,t)[1] \cdot M_G(t))^{\top}} \quad (17)$$

Using these equations for genetic distribution and coalescence probability from one generation to the parent, we can simulate their coalescence using the following algorithm.

296
**Algorithm 4** Simulation of the genealogy given genetic graph and data

1: Set generation time $t = \max\{t \mid \text{card}(An(S, t)) \ge 2\}$.

2: Compute the genetic distribution $P_G(S, t)$ using Equation (3).

3: **repeat**

4: **repeat**

5: Randomly select a pair of lineages from the set $An(S, t)$.

6: Simulate their coalescence (or lack thereof) at $t - 1$ with probability defined in Equation (16).

7: **if** coalescence occurs **then**

8: Merge the lineages and update the genetic distribution using Equation (17).

9: Record the event in a genealogy object, including the probability distribution of each remaining lineage.

10: **end if**

11: **until** all pairs in the set have been tested for coalescence

12: Set $t = t - 1$.

13: **until** all lineages have coalesced

## Probability of a genealogy knowing demographic model

Let us consider a coalescent graph $H_C$ simulated knowing the genetic model. We calculate the probability of the genealogy given the demographic model for each coalescence point $p(In[k] | \theta_P)$, characterized by $t(In[k]) = t_C$ and its descendants, denoted as $In[k]^d = \text{Adj}^{-1}(In[k])$. The number of descendants coalescing in $In[k]$ is given by $\text{card}(In[k]^d) = \text{card}(\text{Adj}^{-1}(In[k]))$.

The number of descendants can be **2 in case of dichotomy**, or **more in case of polytomy**.
These events occur at time $t(In[k]^d)$.

If they do not coalesce more recently, the probability that they co-occur in the same genetic state (in the same population) at $t$ (where $t \le \min(In[k]^D)$) given the demographic model $\theta_P$ is obtained by
computing the scalar product of the probability distributions of each individual at time t, as given by Equation 3.

$$P(An(In[k]^d, t) \text{ coexist in a deme, if they have not coalesced more recently}) = \prod_{i=1}^{\text{card}(In[k]^d)} \left( P_P(In[k]^d[i]) \cdot \prod_{\tau=t(In[k]^d[i])}^{t+1} M_P(\tau) \right) \quad (18)$$

Let us consider a coalescent graph $H_C$ simulated knowing genetic model. We calculate the probability knowing demographic model for each coalescence point $p(In[k] | \theta_P)$ characterized by $t(In[k]) = t_C$ and its descendants noted $In[k]^d = Adj^{-1}(In[k])$. The number of descendants coalescing in $In[k]$ is noted $\text{card}(In[k]^d) = \text{card}(Adj^{-1}(In[k]))$. They can be 2 in case of dichotomy, or more in case of polytomy. They occur a times $t(In[k]^d)$

If they do not coalesce more recently, the probability that they co-occur in the same genetic state (in the same population) at $t$ ($t \le \min(In[k]^D)$) knowing demographic model $\theta_P$ is given by the scalar product of the probability distributions of each individual at time $t$ (given by eq. 3).

$$P(An(In[k]^d, t) \text{ coexist in a deme, if they have not coalesced more recently}) = \prod_{i=1}^{\text{card}(In[k]^d)} \left( P_P(In[k]^d[i]) \cdot \prod_{\tau=t(In[k]^d[i])}^{t+1} M_P(\tau) \right) \quad (19)$$

To obtain the probability that they coalesce at time $t$, we divide for each deme the probability of coexistence by the effective population size of each deme before computing the same scalar product:

$$P(In[k]^d \text{ coalesce in a deme at } t_C \mid \text{ they did not coalesce more recently}) = \left( \prod_{i=1}^{\text{card}(In[k]^d)-1} \left( P_P(In[k]^d[i]) \cdot \prod_{\tau=t(In[k]^d[i])}^{t_C+1} M_P(\tau) \right) \right) \cdot \left( P_P(In[k]^d[\text{card}(In[k]^d)]) \otimes 2Ne(t_C) \cdot \prod_{\tau=t(In[k]^d[\text{card}(In[k]^d)])}^{t_C+1} M_P(\tau) \right)^\top \quad (20)$$
We now consider the probability that they **did not coalesce** more recently than time $t_C$:

$$ P(\text{In}[k]^d \text{ did not coalesce in a deme more recently than } t_C) = \prod_{t=\min(\text{In}[k]^d)}^{t_C} \left( 1 - \left( \prod_{i=1}^{\text{card}(\text{In}[k]^d)-1} \left( P_P(\text{In}[k]^d[i]) \cdot \prod_{\tau=t(\text{In}[k]^d[i])}^{t+1} M_P(\tau) \right) \right) \right. \\ \left. \cdot \left( P_P(\text{In}[k]^d[\text{card}(\text{In}[k]^d)]) \right) \otimes 2Ne(t_C) \cdot \prod_{\tau=t(\text{In}[k]^d[\text{card}(\text{In}[k]^d)])}^{t+1} M_P(\tau) \right)^\top \right) $$

Then, the probability that the descendants $\text{In}[k]^d$ coalesce at time $t_C = t(\text{In}[k])$ (or the probability of $\text{In}[k]$) is:

$$ p(\text{In}[k] | \theta_P) = \left( \prod_{i=1}^{\text{card}(\text{In}[k]^d)-1} \left( P_P(\text{In}[k]^d[i]) \prod_{\tau=t(\text{In}[k]^d[i])}^{t(\text{In}[k]^d)+1} M_P(\tau) \right) \right) \\ \cdot \left( P_P(\text{In}[k]^d[\text{card}(\text{In}[k]^d)]) \right) \otimes 2Ne(t(\text{In}[k])) \prod_{\tau=t(\text{In}[k]^d[\text{card}(\text{In}[k]^d)])}^{t(\text{In}[k]^d)+1} M_P(\tau)^\top \\ \times \prod_{t=\min(\text{In}[k]^d)}^{t(\text{In}[k]^d)} \left( 1 - \left( \prod_{i=1}^{\text{card}(\text{In}[k]^d)-1} \left( P_P(\text{In}[k]^d[i]) \prod_{\tau=t(\text{In}[k]^d[i])}^{t+1} M_P(\tau) \right) \right) \right) \\ \cdot \left( P_P(\text{In}[k]^d[\text{card}(\text{In}[k]^d)]) \right) \otimes 2Ne(t(\text{In}[k])) \prod_{\tau=t(\text{In}[k]^d[\text{card}(\text{In}[k]^d)])}^{t+1} M_P(\tau)^\top \right) \quad (22) $$

The allele probability distribution of coalescent ancestors as is:

$$ P_P(\text{In}[k] | \theta_P) = \prod_{i=1}^{\text{card}(\text{In}[k]^d)} \odot \left( P_P(\text{In}[k]^d[i]) \prod_{\tau=t(\text{In}[k]^d[i])}^{t(\text{In}[k]^d)+1} M_P(\tau) \right) \quad (23) $$

Where $\odot$ is the elementwise multiplication operator

As we presented before for the probability of a genealogy knowing genetic data and model, the probability of a genealogy knowing population data and model is obtained by multiplying the probabilities
of equation (22) over coalescence events and updating the probability distribution of coalescent ancestors using Equation (23) at each coalescence event:

$$p(H_C | P_P(S), \theta_P) = p(H_C | S_P, \theta_P) = \prod_{k=0}^{\text{card}(In)} p(In[k] | \theta_P) \quad (24)$$

where the demographic probability distributions of the nodes $In[k]$, used to calculate $p(In[k] | \theta_P)$ using Equation (22), are successively given by Equation (23).

The algorithm to calculate the demographic probability of a genetic-simulated genealogy is as follows:

**Algorithm 5 Calculation of the probability of a coalescent history given a genetic model**

1: $\mathcal{P} \leftarrow 1$

2: $\mathcal{G} \leftarrow H_C$

3: ▷ The simulated coalescent of the sample whose probability is to be evaluated

4: **while** not all coalescence events have been processed **do**

5: Select the earliest coalescence point (backward in time) from $\mathcal{G}$.

6: Update $\mathcal{P}$ by multiplying it with the pairwise coalescence probability of this event, as defined in Equation (22).

7: Remove the daughter nodes $In[k]^d$ from $\mathcal{G}$.

8: **end while**

This algorithm can be combined with the genetic simulation algorithm so that the coalescence simulation directly yields the demographic probability of the genetic coalescent simulation.
**Algorithm 6** Joint simulation of coalescence on the genetic graph and probability calculation in the population graph

1: Set generation time $t \leftarrow \max\{t | \text{card}(An(S, t)) \ge 2\}$.

2: Compute genetic distribution $P_G(S, t)$ using Equation (3).

3: **while** all lineages have not coalesced **do**

4: **while** all pairs in $An(S, t)$ have not been tested for coalescence **do**

5: Randomly select a pair of lineages from $An(S, t)$.

6: Simulate their coalescence (or lack thereof) at $t - 1$ with probability given in Equation (16).

7: **if** coalescence occurs **then**

8: Merge the lineages and update the genetic distribution using Equation (17).

9: Record the event in a genealogy object, including the probability distribution location of each remaining lineage.

10: Update $\mathcal{P}$ by multiplying it with the pairwise coalescence probability of this event, as defined in Equation (22).

11: Remove the daughter nodes $In[k]^d$ from the tree.

12: Set the demographic and allelic distribution probabilities of the new ancestral node $In[k]$ according to Equation (23).

13: **end if**

14: **end while**

15: Set $t \leftarrow t - 1$.

16: **end while**

332 **Likelihood**

334 The likelihood function is not analytically tractable. The coalescent is considered as latent data that can be simulated from one component of the metamodel, and its probability can be calculated using the other
component. If we simulate the coalescence graph from demographic data (denoted as $H_C^D$) and calculate the likelihood from genetic data, we obtain:

$$ 
\begin{aligned}
\mathcal{L}(\theta) &= p(S_G, S_P | \theta_G, \theta_P) \\
&\propto p(S_G | \theta_G, \theta_P, S_P) \\
&= \int p(S_G | H_C^D, \theta_G) \times p(H_C^D | S_P, \theta_P) d(H_C^G) \\
&\propto \int p(H_C^D | S_G, \theta_G) \times p(H_C^D | S_P, \theta_P) d(H_C^G)
\end{aligned} 
\quad (25) $$

Reciprocally if we simulate coalescence graph from genetic data (we note this graph $H_C^G$) and calculate likelihood from demographic data:

$$ 
\begin{aligned}
\mathcal{L}(\theta) &= p(S_G, S_P | \theta_P, \theta_G) \\
&\propto p(S_P | \theta_P, \theta_G, S_G) \\
&= \int p(S_P | H_C^G, \theta_P) \times p(H_C^G | S_G, \theta_G) d(H_C^G) \\
&\propto \int p(H_C^G | S_P, \theta_P) \times p(H_C^G | S_G, \theta_G) d(H_C^G)
\end{aligned} 
\quad (26) $$

The likelihood can be estimated using a Monte Carlo algorithm, by simulating from demographic data:

$$ \hat{\mathcal{L}}(\theta) \propto \frac{1}{N} \sum_{i}^{N} p(H_C^D | S_G, \theta_G), \quad (27) $$

or from genetic data.

$$ \hat{\mathcal{L}}(\theta) \propto \frac{1}{N} \sum_{i}^{N} p(H_C^G | S_P, \theta_P) \quad (28) $$

## 340 Inference

The Monte Carlo estimated likelihood can be used in Bayesian inference or maximum likelihood estimation. For the latter, and in such a process divided into an underlying process that cannot be directly observed (in our case, the coalescent) and a sampling model (in our case, the probability of the data given the coalescent), Park [15] developed a meta-model assuming a Gaussian distribution of the log-likelihood, which allows for obtaining maximum likelihood statistics.—
Grouping generations for computational tractability

346

347

348

349

350

351

At each generation, genes migrate across demographic and genetic space, eventually coalescing when they coexist in the same demographic group (population) and genetic group (allele). However, these events occur with low probabilities, making direct simulation computationally expensive. To reduce computation time, we can raise the transition matrix to a power, allowing for faster simulation and genealogy probability calculations.

352

353

354

355

356

Since information is extracted from the joint genetic and demographic distribution of the species, the effective population size $N_{e_i}$ can be averaged within a deme over a time period corresponding to the mutation timescale (typically more than $10^4$ generations for most genetic markers). If the environment is changing, the potential effective population size $N_{e_i}$ should be approximated using the harmonic mean of potential effective sizes over that period.

357

358

To determine an appropriate time-scaling factor, simulations should be conducted while accounting for population size, migration rates, and mutation rates.

Discussion

359

360

361

362

363

364

In this paper, I propose a method to efficiently calculate a metamodel likelihood function for estimating genetic parameters (such as mutation rates) and environmental demographic parameters (such as the relationships between demogenetic factors—effective migration and effective population size—and social-ecological environmental variables) from population genetic data and measurable social-ecological history variables.

365

366

367

368

This work represents the first full integration of fundamental niche processes and realized niche pro- cesses related to historical movement in demogenetics. Additionally, it introduces the first likelihood- based approach in environmental demogenetics, opening pathways for maximum likelihood and exact Bayesian inference.
Demogenetic models are inherently complex, and integrating niche and historical processes adds further challenges. Likelihood-based inference has been restricted to simple coalescence models where the likelihood function is tractable. To address the computational challenges of complex "real-life" demographic models, Approximate Bayesian Computation (ABC) was introduced [16] and has become the primary approach for spatial demogenetic modeling [17]. Yet, ABC-based inference relies on summary statistics, whose selection is ambiguous [14, 18], and where acceptance rates decrease with an increasing number of statistics, making estimator consistency difficult to maintain [18, 19].

Instead, we employ direct likelihood-based inference. To reduce computational burden, we integrate the likelihood calculation within the simulation process—either incorporating genetic probability calculations into demographic simulations or vice versa. This approach can be further optimized by grouping demographic generations using a population transition matrix at time scales comparable to mutation time scales.

Studies must be performed to determine the identifiability of demographic parameters from genetic distributions. Since effective population size depends significantly on fluctuations, the rate of recovery may provide signals on population structure [20]. An exhaustive model may incorporate growth and carrying capacity parameters as functions of environmental data. Alternatively, to reduce computational burden we may assume a constant intrinsic growth rate (R) while allowing the carrying capacity (K) to depend on environmental data. Finally, a further simplification can be made by assuming that the species always maintains its effective population size at the carrying capacity.

From a population genetics perspective, genetic structure was historically related to simple drift, migration, mutation, and selection—forces model. With the advance of genetic sequencing and computation technologies, these later forces were progressively linked to complex landscapes [21]. Demogenetics has emerged as a field to infer historical demography models from genetic data [17]. The integration of demogenetics with environmental niche models, termed iDDC (integrated Distribution, Demographic, and Coalescent modeling) [6], has been attempted. However, previous efforts primarily coupled independent
distribution, demographic, and coalescent models rather than forming a truly integrated framework. The
novel approach we propose jointly infers niche and dispersal parameters within a single comprehensive
model. In this framework, the niche model predicts effective population size and migration, the key
variables of interest in population genetics.

From an ecological perspective, modern studies increasingly combine biogeography, paleontology, ge-
netics, and ecological modeling to explain species distributions [22]. Species presence is determined by
a complex interplay of historical constraints (evolutionary history, dispersal, past climate) and current
ecological conditions (abiotic and biotic factors). Predicting the realized niche—the fundamental goal
of ecology—requires considering fundamental niche constraints, species interactions, historical dispersal,
and adaptation. Our approach improves prediction by incorporating genealogical factors and dispersal
constraints into distribution modeling using molecular genetic data.

Traditional species distribution models infer niche parameters from occurrence data, but this inference
is biased because distribution data reflects not only niche constraints but also dispersal limitations and
historical processes. By integrating niche inference with genealogical and dispersal models, we can better
distinguish deterministic and stochastic processes underlying species occurrence. This integration removes
biases in niche parameter estimation: in predicting effective population size, we account for coalescent-
based probabilities of species presence. Although environmental niche potential is a necessary but not
sufficient condition for species presence, our coalescent-based approach helps validate whether this niche
was realized as predicted by migration models. This integration of dispersal and niche processes could
also improve forecasting of future distributions under environmental change scenarios, including climate
and land-use shifts.

A key limitation of our model is that it does not yet account for species interactions, niche evolution,
or adaptation. We assume that the relationship between species density and environmental variables
remains constant over time. However, since our niche models are inferred through genealogical pro-
cesses, they may capture niche range across evolution in the niche model. Social-ecological data could
incorporate anthropological and biotic factors beyond abiotic constraints, enabling a more comprehen-
sive ecological modeling framework. A natural next step would be to infer multiple species’ demogenetic
models simultaneously, allowing for better representation of species interactions.

We present a method to efficiently compute a likelihood function for any social-ecological demogenetic
model that can be represented as a connected graph. This advancement enhances the integration of ecol-
ogy and genetics, leveraging genetic and environmental data to infer population demographic scenarios.
Our contributions include:

• The first likelihood-based inference method for demogenetic models in heterogeneous landscapes, previously limited to ABC approaches.

• A method that directly infers niche parameters rather than merely mapping demographic parame-
ters, overcoming previous limitations.

• An algorithm capable of modeling diverse social and ecological network structures (spatial grids, social groups, ecological interactions), broadening the applicability of genetic inference models.

This algorithm is currently being implemented in the R library graphPOP [23]. The availability of
large genetic datasets—particularly for agricultural species, pests, diseases, and endangered species—combined
with hypotheses on social-ecological drivers of genetic variation, presents a significant opportunity. Our
method has the potential to enhance population genetics as a tool for biodiversity and epidemiology
studies in agriculture ecology and conservation by enabling more precise demographic and ecological
inferences.

By fully integrating potential and realized niches within a genealogical framework, and considering
that genetic data, in this respect, contain more information on ecological parameters than simple occur-
rence data, this approach enables a more accurate understanding of species distributions and their future
projections.
442

## Acknowledgments

443

This work has been supported by IRD EGCE research unit project and BASC labex. I thank to Camille Coron and LMO Orsay France Laboratory, and Arnaud Becheler and Lacey Knowles laboratory members at of Ecology and evolutionary Biology department of Michigan State University for fruit full discussions.

444

445

446

## References

447

1. Barve N, Barve V, Jiménez-Valverde A, Lira-Noriega A, Maher SP, Peterson aT, et al. The crucial role of the accessible area in ecological niche modeling and species distribution modeling. Ecological Modelling. 2011;222(11):1810-1819. doi:10.1016/j.ecolmodel.2011.02.011.

448

449

2. Matos-Maraví P. An "omics" approach to bridge community ecology and island biogeography. Molecular Ecology. 2020;29(9):1592-1595. doi:10.1111/MEC.15426.

450

451

3. Gotelli NJ, Stanton-Geddes J. Climate change, genetic markers and species distribution modelling. Journal of Biogeography. 2015;42(9):1577-1585. doi:https://doi.org/10.1111/jbi.12562.

452

453

4. Fordham DA, Brook BW, Moritz C, Nogués-Bravo D. Better forecasts of range dynamics using genetic data. Trends in Ecology & Evolution. 2014;29(8):436-443. doi:10.1016/j.tree.2014.05.007.

454

455

5. Naughtin SR, Castilla AR, Smith AB, Strand AE, Dawson A, Hoban S, et al. Integrating genetic data and simulations to evaluate alternative species distribution models and improve predictions of glacial refugia and future responses to climate change. Ecography. 2024; p. e07196. doi:10.1111/ecog.07196.

456

457

458

459

6. He Q, Edwards DL, Knowles LL. INTEGRATIVE TESTING OF HOW ENVIRONMENTS FROM THE PAST TO THE PRESENT SHAPE GENETIC STRUCTURE ACROSS LANDSCAPES. Evolution. 2013;67(12):3386-3402. doi:10.1111/evo.12159.

460

461

462
463

7. Currat M, Ray N, al E. SPLATCHE: a program to simulate genetic diversity taking into account environmental heterogeneity. Molecular Ecology Notes. 2004;4:139-142.

464

465

8. Rebaudo F, Le Rouzic A, Dupas S, Silvain J, Harry M, Dangles O. Sim<span style="font-variant:small-caps;">>A</span>dapt: an individual-based genetic model for simulating landscape management impacts on populations. Methods in Ecology and Evolution. 2013;4(6):595-600. doi:10.1111/2041-210x.12041.

466

467

468

469

9. Landguth EL, Bearlin A, Day CC, Dunham J. CDMetaPOP: an individual-based, eco-evolutionary model for spatially explicit simulation of landscape demogenetics. Methods in Ecology and Evolution. 2017;8(1):4-11. doi:10.1111/2041-210x.12608.

470

471

472

10. Yannic G, Hagen O, Leugger F, Karger DN, Pellissier L. Harnessing paleo-environmental modeling and genetic data to predict intraspecific genetic structure. Evolutionary Applications. 2020;13(6):1526-1542. doi:10.1111/eva.12986.

473

474

475

11. Ray N, Currat M, Foll M, Excoffier L. SPLATCHE2: a spatially explicit simulation framework for complex demography, genetic admixture and recombination. Bioinformatics (Oxford, England). 2010;26(23):2993-4. doi:10.1093/bioinformatics/btq579.

476

477

478

12. Currat M, Arenas M, Quilodran CS, Excoffier L, Ray N. SPLATCHE3: simulation of serial genetic data under spatially explicit evolutionary scenarios including long-distance dispersal. Bioinformatics. 2019;35(21):4480-4483. doi:10.1093/bioinformatics/btz311.

479

480

481

13. Becheler A, Coron C, Dupas S. Quetzal - an open source C++ template library for coalescence-based environmental demogenetic models inference. bioRxiv. 2017; p. 1-27. doi:10.1101/214767.

482

483

14. De Visscher L, De Baets B, Baetens JM. A critical review of common pitfalls and guidelines to effectively infer parameters of agent-based models using Approximate Bayesian Computation. Environmental Modelling & Software. 2024;172:105905. doi:https://doi.org/10.1016/j.envsoft.2023.105905.

484

485
486

487

15. Park J. Scalable simulation-based inference for implicitly defined models using a metamodel for Monte Carlo log-likelihood estimator; 2024. Available from: http://arxiv.org/abs/2311.09446.

488

489

490

16. Pritchard JK, Seielstad MT, Perez-Lezaun A, Feldman MW. Population growth of human Y chromosomes: a study of Y chromosome microsatellites. Molecular Biology and Evolution. 1999;16(12):1791-1798. doi:10.1093/oxfordjournals.molbev.a026091.

491

492

17. Csilléry K, Blum MGB, Gaggiotti OE, François O. Approximate Bayesian Computation (ABC) in practice. Trends in Ecology & Evolution. 2010;25(7):410-418. doi:10.1016/j.tree.2010.04.001.

493

494

18. Karabatsos G, Leisen F. An Approximate Likelihood Perspective on ABC Methods; 2018. Available from: http://arxiv.org/abs/1708.05341.

495

496

19. Nakagome S, Fukumizu K, Mano S. Kernel Approximate Bayesian Computation for Population Genetic Inferences; 2013. Available from: https://arxiv.org/abs/1205.3246.

497

498

499

20. Crosato E, Philippson JN, Thutupalli S, Morris RG. Dynamical demographic phases explain how population growth and mutation control the evolutionary impact of bottlenecks. Physical Review Research. 2023;5(1):013093. doi:10.1103/PhysRevResearch.5.013093.

500

501

21. Storfer A, Murphy MA, Spear SF, Holderegger R, Waits LP. Landscape genetics: where are we now? Molecular Ecology. 2010;19(17):3496-3514. doi:10.1111/j.1365-294X.2010.04691.x.

502

503

504

22. Morales-Castilla I, Davies TJ, Pearse WD, Peres-Neto PR. Combining phylogeny and co-occurrence to improve single species distribution models. Global Ecology and Biogeography. 2017;26:740-752. doi:10.1111/geb.12580.

505

506

23. Patarroyo C, Dupas S. graphPOP; 2024. Available from: https://github.com/stdupas/GraphPOP.