# Embedding Birth-Death Processes within a Dynamic Stochastic Block Model

Gabriela Bayolo Soler, Miraine Dávila Felipe, Ghislaine Gayraud

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

## To cite this version:

Gabriela Bayolo Soler, Miraine Dávila Felipe, Ghislaine Gayraud. Embedding Birth-Death Processes within a Dynamic Stochastic Block Model. 2026. hal-05430728

HAL Id: hal-05430728

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

Preprint submitted on 26 Jan 2026

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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Embedding Birth-Death Processes within a
Dynamic Stochastic Block Model

Gabriela Bayolo Soler¹*, Miraine Dávila Felipe¹ and
Ghislaine Gayraud¹

¹*Université de technologie de Compiègne, Alliance Sorbonne
Université, LMAC, Compiègne, France.

*Corresponding author(s). E-mail(s): gbayolo@gmail.com;

## Abstract

Statistical clustering in dynamic networks aims to identify groups of nodes with similar or distinct internal connectivity patterns as the network evolves over time. While early research primarily focused on static Stochastic Block Models (SBMs), recent advancements have extended these models to handle dynamic and weighted networks, allowing for a more accurate representation of temporal variations in structure. Additional developments have introduced methods for detecting structural changes, such as shifts in community membership. However, limited attention has been paid to dynamic networks with variable population sizes, where nodes may enter or exit the network. To address this gap, we propose an extension of dynamic SBMs (dSBMs) that incorporates a birth-death process, enabling the statistical clustering of nodes in dynamic networks with evolving population sizes. This work makes three main contributions: (1) the introduction of a novel model for dSBMs with birth-death processes, (2) a framework for parameter inference and prediction of latent communities in this model, and (3) the development of an adapted Variational Expectation-Maximization (VEM) algorithm for efficient inference within this extended framework.

**Keywords:** dynamic networks, stochastic block model, variational inference, birth-death processes, latent variable models, clustering.

## 1 Introduction

Network models play a crucial role in capturing and understanding population dynamics. Unlike models based on simple random interactions, real-world networks often
exhibit complex structures that reflect underlying patterns and regularities. In the literature, there are many examples of stochastic network models that help to describe these structures, such as the random graph model introduced by Erdős and Rényi [1], the small-world network presented by [2], and the Stochastic Block Model (SBM) by [3]. These models have been applied across a broad spectrum of disciplines, including the social sciences, where they help explaining human behavior and social interactions; statistics, where they provide frameworks for analyzing relational data; physics, where they model phenomena such as percolation and diffusion; and biology, where they help making predictions about the structure and dynamics of interbreeding populations.

In this work we focus on the statistical detection of groups in dynamic networks. The problem of identifying clusters, subsets of nodes that are more strongly or just differently connected internally than with the rest of the network, has attracted considerable attention in the statistics and network science literatures in recent years. Early work on statistical community detection largely concentrated on static SBMs, which model networks with a fixed vertex set and time-invariant connectivity ([3–5]). While these models provide a powerful tool for uncovering latent structure in static networks, they do not account for the temporal evolution that is often observed in real-world systems.

This limitation has motivated a large body of work on dynamic extensions of SBMs and related models. Recent developments have extended static SBMs to weighted networks and to dynamic SBMs in which the edges evolve over time, either in discrete time (see, for example, [6–9]) or in continuous time (see [10, 11]). These extensions allow one to study how relationships and community structures change over time, providing a more realistic description of many applications in which interactions are inherently dynamic. In parallel, there has been growing interest in detecting changes in community structure, including changepoint detection methods that identify abrupt shifts in network organization ([12, 13]) and approaches designed to capture more gradual, diffuse changes in community memberships or connectivity patterns ([6, 8, 10]). Closely related block models have also been developed for bipartite or matrix data through the Latent Block Model [14].

A common limitation of most dynamic SBMs, however, is the *fixed population* assumption: the node set $V$ is taken to be constant, whereas real systems often exhibit entry and exit through birth, death, migration or isolation. Some works allow transient presence/absence or introduce a special “inactive” class within a fixed population [6, 15–17], and there is growing interest in regimes with structural shifts and phase transitions [18]. Nevertheless, to the best of our knowledge, there remains relatively little methodology that *jointly* models network evolution and genuinely *evolving-population* dynamics. This motivates the birth-death SBM developed below, in which community sizes follow explicit demographic (birth-death) processes while edges are generated according to an SBM-type mechanism.

In this work, we address this gap by extending dynamic stochastic block models (dSBMs) to account for changes in population size through the integration of a birth-death process. Birth-death processes provide a natural probabilistic description of population dynamics, representing the arrivals (births) and departures (deaths) of
individuals over time. By combining a dSBM for the network structure with a birth-death mechanism for the vertex set, we obtain a model that captures both the evolution of edges and the stochastic variation in population size.

More precisely, we introduce a *dynamic stochastic block model with birth-death population dynamics*, which we refer to as the *birth-death stochastic block model* (BD-SBM). In the BD-SBM, each node belongs to one of *K* latent communities. When a new node is born, it inherits the community of its parent and retains this community membership throughout its lifetime, so that there is no community switching over time. At each observation time, conditional on the current community memberships, the edges follow a stochastic block model. Between observation times, the population size evolves according to a continuous-time birth-death process: new individuals enter the network and existing individuals leave, while the latent community structure that drives edge formation is preserved for all living individuals. This explicit coupling between population dynamics and network structure yields a richer and more realistic description of temporal networks in which both edges and vertices evolve.

Beyond its methodological interest, the BD-SBM is particularly suited to applications where both the population and the interaction structure evolve over time, and where community memberships are naturally inherited and stable. Examples include populations structured into families, lineages or clans in ecology and population genetics, where descendants belong to the community of their parents and remain in that group until death. Social and organizational settings in which newcomers join pre-existing groups (such as departments, teams or schools) and do not change affiliation over the observation window. For example, academic research laboratories organized into thematic teams (e.g. probability, statistics, optimization, machine learning), where new PhD students or postdoctoral fellows typically join the team of their advisor and leave the system when they graduate, move to another institution or retire. In such contexts, the BD-SBM provides a coherent framework to jointly analyse the evolution of community structure and population size, while explicitly accounting for demographic turnover. In this context, the BD-SBM can be used to identify the teams of collaborations and the ways of collaborations of the researches.

From an inferential viewpoint, the BD-SBM raises non-trivial challenges. The presence of latent community labels and unobserved community sizes makes direct maximum likelihood (ML) estimation infeasible, as the observed-data likelihood would require summing over an enormous number of possible label configurations. To overcome this difficulty, we build on the mean-field variational framework introduced for static SBMs in [19] and later extended to dynamic SBMs with fixed population size in [6]. We propose an adapted Variational Expectation-Maximization (VEM) algorithm tailored to the BD-SBM. Our variational family combines a mean-field approximation for individual labels with a structured approximation for the latent community sizes, encoded through additional variational parameters that describe the evolution of community sizes over time.

In summary, our main contributions are threefold:

1. We introduce the BD-SBM, a dynamic stochastic block model in which the vertex set evolves according to a birth-death process, allowing the population size to vary over time.
2. We address the challenging problem of parameter estimation and community detection in this setting, and develop a VEM algorithm adapted to the BD-SBM that jointly estimates the model parameters and the latent community memberships. In particular, we propose a structured variational family that augments the classical mean-field approximation with latent variables describing community sizes, thereby enforcing coherence between population dynamics and community structure.

3. We assess the performance of the BD-SBM on both simulated temporal networks and a real-world collaboration network constructed from arXiv data.

These developments extend the scope of stochastic block modelling to temporal networks with evolving population size, and provide a principled framework for studying the joint dynamics of communities and populations in a wide range of applications.

# 2 The birth-death stochastic block model (BD-SBM)

In this section we formally define the BD-SBM, a dynamic stochastic block model in which the population size evolves according to a birth-death process.

## 2.1 Model description

We consider a dynamic undirected network $G(t) = (V(t), E(t))$, indexed by time $t \ge 0$, which represents the interactions between individuals at $t$. The time $t_0$ denotes the initial observation time. For each $t$, $V(t)$ is the set of active vertices (individuals), and $E(t)$ is the set of edges (interactions) between any two vertices in $V(t)$. Since we consider only binary (undirected) edges, we use a notation shorthand: we denote indistinctly by $E(t)$ the set of edges and the adjacency matrix of the graph at $t$, that is, the matrix whose coefficients are $e_{ij}(t)$, $1 \le i, j \le V(t)$, where $e_{ij}(t)$ is defined as the indicator of the presence of an edge between $i$ and $j$ at time $t$. Lastly, since the number of individuals varies over time, we write $N(t) = |V(t)|$ for the cardinality of $V(t)$.

We assume that once a node becomes active in the network, it is assigned to exactly one of the $K$ distinct communities $\{1, 2, \dots, K\}$ to which it belongs throughout its lifetime. It means that once active, a node keeps its community membership during all its lifetime. The number of communities $K$ is fixed and known in advance. We associate to any individual $i$ its community membership random vector $Z_i$ defined by $Z_i = (Z_{ik})_{k=1,\dots,K}$, where for all $k$, $Z_{ik} = 1$ if $i$ is in community $k$ and $\sum_{k=1}^K Z_{ik} = 1$.

Randomness affects the temporal edges, the population size over time and the community memberships of individuals.

1. **Population dynamics and community membership.** The BD-SBM assumes a finite population of $N$ nodes partitioned into $K$ disjoint communities as follows:

(a) For each node $i \in V_0$, the community label $Z_i = (Z_{ik})_{k=1,\dots,K}$ is drawn independently with

$$ \mathbb{P}(Z_{ik} = 1) = \beta_k, \quad k = 1, \dots, K, $$
where $\beta = (\beta_k)_{k=1,...,K}$ denotes the vector of community probabilities at time $t_0$.

(b) For each node $i \notin V_0$, the population dynamics are governed by a linear birth–death process with common birth rate $\lambda \ge 0$ and death rate $\mu \ge 0$, which are shared across communities. When a birth occurs at time $t$ and a new node $i$ is created, it inherits the community of its parent. Conditionally on a birth at time $t$, we have

$$ \mathbb{P}(Z_{ik} = 1) = \frac{\sum_{j \in V(t^{-})} Z_{jk}}{N(t^{-})}, \quad k = 1, \dots, K, $$

where $t^{-}$ denotes an instant just before $t$ and $N(t^{-}) = |V(t^{-})|$ is the population size at that time.

Labels are assigned sequentially: the label of a newly created node is defined as one unit larger than the label of the previously created node, or, initially, larger than the number of nodes present at $t_0$. Once a node dies, its label no longer appears in subsequent vertex sets.

*Remark.* The assumption of common birth and death rates across communities allows us to derive explicit maximum likelihood estimators for $(\lambda, \mu)$, as shown in Subsection 3.2. The more general case in which the birth–death parameters are community-specific is discussed in Appendix B as an extension of our VEM approach.

2. **Temporal edges.** For $t \ge t_0$, given the partition of the vertices $V(t)$ into $K$ communities, the graph $G(t)$ is generated according to a stochastic block model with connectivity matrix $\pi = (\pi_{km})_{1 \le k,m \le K}$. More precisely, for all $i < j$ such that $i, j \in V(t)$,

$$ e_{ij}(t) | \{Z_{ik} Z_{jm} = 1\} \stackrel{\text{ind.}}{\sim} \phi(\cdot, \pi_{km}), $$

where $e_{ij}(t)$ denotes the edge between $i$ and $j$ at time $t$, and $\phi(\cdot, \pi_{km})$ stands for the probability mass function of the Bernoulli distribution with success probability $\pi_{km}$.

*Notation.* In the interest of clarity, we avoid the introduction of an extra notation to distinguish between the random variables $e_{ij}(t)$ and their realisations. Whenever this notation is used, its role is whether specified or can be easily deduced from the context.

We refer to Figure 1 for a hierarchical overview of our BD-SBM. Let us define $\theta$ as the model parameters

$$ \theta = (\lambda, \mu, \pi, \beta) \in \mathbb{R}_{+} \times \mathbb{R}_{+} \times [0, 1]^{K^2} \times [0, 1]^{K}. $$
<figure><img src="image_3.png" /><figcaption>Fig. 1: Graphical representation of the BD-SBM model: community labels $Z_i$ are defined on the alive set $V(t)$, which evolves in continuous time under a birth-death process with rates $\lambda$ and $\mu$. The birth or death events occurs at ordered times $\tau = \{\tau_1, \tau_2, \dots\}$. At observation times $t_\ell$, the snapshot graph $G_{t_\ell}$ is generated from the current labels via the edge mechanism $\phi$.</figcaption></figure>

## 2.2 Observations and latent variables

Snapshots of the dynamic undirected network $G(t) = (V(t), E(t))$ are observed at discrete times $t \in T = \{t_0, t_1, \dots, t_T\}$. We denote by $V_0$ the set of vertices at time $t_0$ and by $N_0 := |V_0|$ its cardinality. The birth-death process is assumed to be observed continuously, meaning that the ordered sequence of birth-death event times, namely

$$ \tau = \{\tau_1, \tau_2, \dots, \tau_M\}, \quad t_0 \le \tau_1 \text{ and } \tau_M \le t_T, $$

the nature of these events (birth or death), and the label of the individual undergoing each birth or death event are all observed. Consequently, the number of vertices $N(t)$ is known at any $t \in [t_0, t_T]$.

The community affiliations of the individuals $i \in V := \cup_{\ell \in \tau} V(t)$ are unobserved, and hence $Z = \{Z_i\}_{i \in V}$ is a collection of latent vectors. Let $N := |V|$ be the total number of individuals present in the system over the time interval $[t_0, t_T]$.

For clarity, we summarize below the observed data and the latent variables involved in the model:

### 1. Observed variables:

(a) the set of edges at times in $T = \{t_0, \dots, t_T\}$, the so-called snapshots, defined as,

$$ e(T) := \{(e_{ij}(t_\ell))_{i,j}\}_{\ell=0,\dots,T}; $$

(b) the ordered sequence of birth-death event times $\tau = \{\tau_1, \tau_2, \dots, \tau_M\}$, with $t_0 \le \tau_1$ and $\tau_M \le t_T$; for notation convenience, we introduce an additional artificial event time $\tau_0$ which coincides with $t_0$, i.e., $\tau_0 := t_0$;
(c) the associated individuals and nature of events of the birth–death process, described by the sequence

$$b(\tau) := \{(i_{\ell}, b_{\ell})\}_{\ell=1, \dots, M},$$

where $i_\ell$ denotes the label of the individual undergoing the birth or death event at time $\tau_\ell$, and $b_\ell$ indicates whether a birth ($b_\ell = 1$) or a death ($b_\ell = -1$) occurs at time $\tau_\ell$. Consequently, the set $V(\tau_\ell)$, and the number $N(\tau_\ell)$, of alive individuals, is known for every $\ell = 0, \dots, M$.

Moreover, for each individual $i$, we denote by $\tau_i^b$ and $\tau_i^d$ her/his birth and death times, respectively, with the convention $\tau_i^b = t_0$ for $i$ present in $V_0$ and $\tau_i^d = t_T$ for $i$ present in $V(t_T)$. We also denote by $T_B$ and $T_D$ the set of birth times, respectively death times:

$$T_B = \{\tau_\ell : b_\ell = 1\}, \quad T_D = \{\tau_\ell : b_\ell = -1\}.$$

## 2. Latent variables:

(a) the community membership random vectors $Z = \{Z_i\}_{i \in V}$, where, for all $i \in V$, $Z_i = (Z_{ik})_{k=1, \dots, K}$ with $Z_{ik} = 1$ if $i$ is in community $k$ and $\sum_{k=1}^K Z_{ik} = 1$.

Figure 2 illustrates the generative mechanism of the BD-SBM, showing the evolution of active node set and community labels, and the associated observed snapshots.

For convenience, we introduce the latent community-size indicators at the jump times of the birth-death process. For each event time $\tau_\ell$ with index $\ell = 0, \dots, M$, each community $k = 1, \dots, K$ and each integer $n \ge 0$, we define,

$$L^\ell = \{L_k^\ell\}_k = (L_{k,n}^\ell)_{k,n}, \quad \text{where} \quad L_{k,n}^\ell = 1 \left\{ \sum_{i \in V(\tau_\ell)} Z_{ik} = n \right\}.$$

Thus, $L_{k,n}^\ell = 1$ if and only if community $k$ has size $n$ at time $\tau_\ell$. These variables satisfy

$$\begin{align*}
\sum_{n=0}^{N(\tau_\ell)} L_{k,n}^\ell &= 1, \\
\sum_{n=0}^{N(\tau_\ell)} n L_{k,n}^\ell &= \sum_{i \in V(\tau_\ell)} Z_{ik}, \\
\sum_{k=1}^{K} \sum_{n=0}^{N(\tau_\ell)} n L_{k,n}^\ell &= N(\tau_\ell),
\end{align*}$$

where $N(\tau_\ell) = |V(\tau_\ell)|$ denotes the number of individuals alive at time $\tau_\ell$.
<figure><img src="image_4.png" /><figcaption>Observed network snapshots (edges among currently alive nodes)</figcaption></figure>

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

snapshot at $t_0$
$|V(t)| = 3$

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

snapshot at $t_1$
$|V(t)| = 6$

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

snapshot at $t_2$
$|V(t)| = 5$

**Fig. 2:** Schematic BD-SBM with $K = 3$ communities: individuals appear and disappear following a continuous time birth–death process. For $i = 1, \dots, 8$, individual $i$ has a a community label defined as a vector $Z_i = (Z_i)_{1 \le i \le 3}$, where $Z_{ik} = 1$ if individual $i$ is in community $k$, and 0 otherwise. Communities are represented by different colors. At observation times $t_0, t_1, t_2$, we are informed by a network snapshot among the individuals alive at that time; which is represented by the adjacency matrix of the network at that time. The size of the population at snapshot time $t_m$ is $|V(t_m)|$. Diagonal blocks, corresponding to interactions from individuals in the same community, are highlighted. Interactions are symmetric and represented by dots centered in cells.

Moreover, the size of each community evolves over time according to a birth–death Markov chain and can change by at most one unit between two consecutive jump times. This is encoded, for every $\ell = 1, \dots, M$, and every $k \in \{1, \dots, K\}$, by

$$ \text{there is a birth at time } \tau_{\ell} \Leftrightarrow \begin{cases} \sum_{k,n} L_{k,n+1}^{\ell} L_{k,n}^{\ell-1} = 1, \\ \sum_{k,n' \neq n+1} L_{k,n'}^{\ell} L_{k,n}^{\ell-1} = 0, \end{cases} $$

$$ \text{there is a death at time } \tau_{\ell} \Leftrightarrow \begin{cases} \sum_{k,n} L_{k,n-1}^{\ell} L_{k,n}^{\ell-1} = 1, \\ \sum_{k,n' \neq n-1} L_{k,n'}^{\ell} L_{k,n}^{\ell-1} = 0. \end{cases} $$
# 3 Inference for the BD-SBM

From the observed data, recall that our inferential goals are twofold: (i) to recover the classes to which each individual belongs (latent community membership random vector), and (ii) to estimate the community- and population-dynamic parameters, that is, the model parameter vector $\theta = (\lambda, \mu, \pi, \beta)$.

Models with latent variables, such as the one considered here, pose specific challenges for maximum likelihood estimation because the likelihood of the observed data involves integrating over all possible configurations of the latent variables. In our setting, the observed-data likelihood of $(\tau, e(\mathcal{T}), b(\tau))$ is obtained by summing the complete-data likelihood over all possible assignments of the latent community memberships $Z$. Except in very small systems (e.g. $N$ and $M$ small), this summation is not tractable and the observed-data likelihood does not admit a closed-form expression. As a consequence, computing or maximizing the likelihood directly is computationally prohibitive, which would severely limit the practical use of the model.

To overcome this challenge, the inference is usually carried out by algorithms that either find or approximate locally maxima of the observed-data log-likelihood. A classical approach is the expectation-maximization (EM) algorithm, which alternates between computing the conditional expectation of the complete-data log-likelihood (E-step) and maximizing this expectation with respect to the parameters (M-step). In stochastic block models, however, the conditional distribution of the latent labels given the observed edges does not factorize across individuals, so the E-step is not available in closed form and becomes computationally intractable.

This motivates the use of variational approximations to EM (VEM; see [20]). For static SBMs, VEM was introduced by [19] and has since been extended in various directions; see, for instance, [21] for an overview and comparisons with alternative estimators. In parallel, Markov chain Monte Carlo (MCMC) methods have also been proposed for dynamic network models, providing accurate posterior approximations at the price of a substantial computational cost; see for example [10, 22]. In contrast, variational methods offer a scalable alternative by optimizing a tractable lower bound on the evidence (ELBO), and have been successfully applied to dynamic SBMs [6, 11]. In the setting of our BD-SBM, we adopt a Variational Expectation-Maximization (VEM) strategy.

While the birth–death parameters $(\lambda, \mu)$ admit closed-form maximum likelihood estimators, the SBM parameters $(\beta, \pi)$ together with the class at which each individual belongs do not; hence, we use a variational approximation to estimate the latter.

In the remainder of this section, we first derive the complete-data likelihood, then describe the estimation of the birth–death parameters, introduce the variational family, and finally present the VEM algorithm, model selection criterion, and initialization procedure.

## 3.1 Complete-data likelihood

We begin by deriving the complete-data likelihood of $(Z, \tau, e(\mathcal{T}), b(\tau))$ parametrized by $\theta = (\lambda, \mu, \pi, \beta)$ and associated to the probability distribution denoted $P_{\theta}$. By a slight abuse of notation, we denote this complete-data likelihood by $P_{\theta}(Z, \tau, e(\mathcal{T}), b(\tau))$.
Using the conditional independence structure of the model, this complete-data likelihood can be factorized as follows,

$$ \mathbb{P}_{\theta}(Z, \tau, e(\mathcal{T}), b(\tau)) = \mathbb{P}_{\theta}(e(\mathcal{T}) | Z, \tau, b(\tau)) \mathbb{P}_{\theta}(Z, \tau, b(\tau) | Z^{t_0}) \mathbb{P}_{\theta}(Z^{t_0}), \quad (1) $$

where $Z^{t_0} = \{Z_i\}_{i \in V_0}$ denotes the community memberships of the individuals present in the system at time $t_0$.

**Proposition 1** The complete-data log-likelihood of the BD-SBM is given by

$$ \log \mathbb{P}_{\theta}(Z, \tau, e(\mathcal{T}), b(\tau)) = \sum_{i<j} \sum_{k_1, k_2} Z_{ik_1} Z_{jk_2} \sum_{t_{\ell} \in \mathcal{Y}_{ij}} \log \phi(e_{ij}(t_{\ell}), k_1, k_2) + |T_B| \log \lambda + |T_D| \log \mu \\ -(\lambda + \mu)I_N + \sum_{\tau_{\ell} \in T_B} \sum_k \sum_{n=1}^{N(\tau_{\ell}-1)} L_{k,n+1}^{\ell} L_{k,n}^{\ell-1} \log n + \sum_{i \in V_0} \sum_k Z_{ik} \log \beta_k, \quad (2) $$

where, for any pair of distinct individuals $(i, j)$ and any pair of communities $(k_1, k_2)$,

$$ \phi(e_{ij}(t_{\ell}), k_1, k_2) = (\pi_{k_1 k_2})^{e_{ij}(t_{\ell})} (1 - \pi_{k_1 k_2})^{1 - e_{ij}(t_{\ell})}, \\ \mathcal{Y}_{ij} = \{t_{\ell} \in \mathcal{T} : t_{\ell} \in [\max(\tau_i^b, \tau_j^b), \min(\tau_i^d, \tau_j^d)]\}, \\ I_N = \sum_{\ell=1}^{M} N(\tau_{\ell-1})(\tau_{\ell} - \tau_{\ell-1}). $$

Note that by construction, $\mathcal{Y}_{ij}$ is the set of snapshot times at which both $i$ and $j$ are alive in the system.

*Proof* Under the BD-SBM setting, the three terms in the right-hand side of Equation (1) are defined as follows,

$$ \mathbb{P}_{\theta}(e(\mathcal{T}) | Z, \tau, b(\tau)) = \prod_{i<j} \prod_{k_1, k_2} \prod_{t_{\ell} \in \mathcal{Y}_{ij}} [\phi(e_{ij}(t_{\ell}), k_1, k_2)]^{Z_{ik_1} Z_{jk_2}}, \\ \mathbb{P}_{\theta}(Z, \tau, b(\tau) | Z^{t_0}) = \prod_{\ell=1}^{M} \left[ e^{-N(\tau_{\ell-1})(\lambda+\mu)\Delta_{\ell}} \prod_{k=1}^{K} \left( \lambda \sum_{i \in V(\tau_{\ell-1})} Z_{ik} \right)^{Z_{j_{\ell}k} 1\{b_{\ell}=1\}} \mu^{Z_{j_{\ell}k} 1\{b_{\ell}=-1\}} \right] \\ = \lambda^{|T_B|} \mu^{|T_D|} e^{-(\lambda+\mu)I_N} \prod_{\tau_{\ell} \in T_B} \prod_{k=1}^{K} \left( \sum_{i \in V(\tau_{\ell-1})} Z_{ik} \right)^{Z_{j_{\ell}k}} \\ = \lambda^{|T_B|} \mu^{|T_D|} e^{-(\lambda+\mu)I_N} \prod_{\tau_{\ell} \in T_B} \prod_{k=1}^{K} \prod_{n=1}^{N(\tau_{\ell}-1)} n^{L_{k,n+1}^{\ell} L_{k,n}^{\ell-1}}, \\ \mathbb{P}_{\theta}(Z^{t_0}) = \prod_{i \in V_0} \prod_k \beta_k^{Z_{ik}}, $$

with $\Delta_{\ell} = \tau_{\ell} - \tau_{\ell-1}$.

Summing the logarithm of these three terms gives the result.

□
## 3.2 Parameter estimation of the birth-death process

From Equation (2), it is apparent that the contributions of $\lambda$ and $\mu$ to the complete-data likelihood depend only on the observed birth and death times and on the total population size at each event time, but not on the latent community labels $Z$. Straight-forward maximization of the log-likelihood with respect to $(\lambda, \mu)$ yields the following maximum likelihood estimators:

$$ \hat{\lambda} = \frac{|T_B|}{I_N} \quad \text{and} \quad \hat{\mu} = \frac{|T_D|}{I_N}, $$

with

$$ I_N = \sum_{\ell=1}^{M} N(\tau_{\ell-1})\Delta_{\ell}, \quad \Delta_{\ell} = \tau_{\ell} - \tau_{\ell-1}. $$

These coincide with the classical ML estimators for a homogeneous birth-death process with constant birth and death rates that can be found, for instance, in [23, 24].

In contrast, the maximum likelihood estimation of the community proportions $\beta_k$ and block connection probabilities $\pi_{k_1k_2}$, based on the marginal distribution of the observed data alone, is neither computationally tractable nor available in closed form. We therefore resort to a VEM-type algorithm to estimate $(\pi, \beta)$ together with the latent variables (see Section 3.3).

In what follows, we denote by

$$ \hat{\theta} = (\hat{\lambda}, \hat{\mu}, \pi, \beta), $$

the parameter vector of the model composed by the ML estimators of $(\lambda, \mu)$, complemented with the remaining parameters $(\pi, \beta)$ that are still to be estimated.

## 3.3 Parameter estimation for the BD-SBM and prediction of latent variables

Classical EM algorithms rely on the posterior distribution of the latent variables given the observations, that is,

$$ \mathbb{P}_{\theta}(Z | \tau, e(\mathcal{T}), b(\tau)) = \mathbb{P}_{\theta}((Z_i)_i | \tau, e(\mathcal{T}), b(\tau)). $$

In our model, as in standard SBMs, this posterior distribution does not factorize over individuals: for any pair of distinct nodes $i$ and $j$, and any pair of communities $(k_1, k_2)$, the indicators $Z_{ik_1}$ and $Z_{jk_2}$ are no longer independent when conditioning on the data. As a result, the exact posterior is not tractable and the E-step of EM cannot be performed in closed form.

Instead, we derive a variational expectation-maximization algorithm (VEM), a particular instance of variational inference (see [25]). The idea is to replace the intractable posterior by a tractable family of distributions $Q$ on $Z$, and to optimize within this family a lower bound on the observed-data log-likelihood.
In the sequel, $q \in \mathcal{Q}$ denotes either a distribution or its associated likelihood function. The notation $\mathbb{E}_{q}[\cdot]$ refers to expectation with respect to $q$. Next, for any $q \in \mathcal{Q}$, we decompose the observed-data log-likelihood as follows,

$$ \log \mathbb{P}_{\tilde{\theta}}(\tau, e(\mathcal{T}), b(\tau)) = \mathcal{L}(q) + \mathrm{KL}(q \| \mathbb{P}_{\tilde{\theta}}(\cdot | \tau, e(\mathcal{T}), b(\tau))), $$

where

$$ \mathcal{L}(q) = \mathbb{E}_q[\log \mathbb{P}_{\tilde{\theta}}(Z, \tau, e(\mathcal{T}), b(\tau))] - \mathbb{E}_q[\log q(Z)], $$

is the evidence lower bound (ELBO), and $\mathrm{KL}(\cdot \| \cdot)$ denotes the Kullback-Leibler divergence. Both can be written as follows,

$$ \begin{aligned} \mathrm{KL}(q \| \mathbb{P}_{\tilde{\theta}}) &= -\mathbb{E}_q \left( \log \left( \frac{\mathbb{P}_{\tilde{\theta}}(Z | \tau, e(\mathcal{T}), b(\tau))}{q(Z)} \right) \right), \\ \mathcal{L}(q) &= \mathbb{E}_q \left( \log \left( \frac{\mathbb{P}_{\tilde{\theta}}(Z, \tau, e(\mathcal{T}), b(\tau))}{q(Z)} \right) \right) \\ &= \mathbb{E}_q \left( \log \mathbb{P}_{\tilde{\theta}}(\tau, e(\mathcal{T}), b(\tau), Z) \right) - \mathbb{E}_q \left( \log q(Z) \right). \end{aligned} \quad (3) $$

Since the KL divergence is non-negative, we have

$$ \log \mathbb{P}_{\tilde{\theta}}(\tau, e(\mathcal{T}), b(\tau)) \geq \mathcal{L}(q), $$

with the equality holding if and only if $q(Z)$ coincides with the true posterior distribution $\mathbb{P}_{\tilde{\theta}}(Z | \tau, e(\mathcal{T}), b(\tau))$. The idea of the VEM approach is to maximize the ELBO $\mathcal{L}(q)$ with respect to $q \in \mathcal{Q}$, which is equivalent to minimize the Kullback-Leibler (KL) divergence between $q$ and the true posterior distribution with respect to $q \in \mathcal{Q}$. Without restriction on the distribution family $\mathcal{Q}$, the ELBO function is maximal for $q(Z) = \mathbb{P}_{\tilde{\theta}}(Z | \tau, e(\mathcal{T}), b(\tau))$, as the KL divergence becomes zero. However in the considered model, since the posterior distribution is not tractable, the goal is to approximate it on a given family of variational distributions $\mathcal{Q}$ which are tractable.

### 3.3.1 Variational distribution family $\mathcal{Q}$

We now specify the family of variational distributions $\mathcal{Q}$ used in our VEM algorithm. A fully factorized approximation, assuming independence among all individuals, would disregard the strong constraints that the birth-death process imposes on community sizes, and would also increase the computational cost by a factor of $N$. Instead, we adopt a structured variational approximation that explicitly accounts for the joint evolution of community sizes and, through them, induces dependencies among community memberships, while still remaining computationally tractable.

Let us define the family of distributions $\mathcal{Q}$ that we consider. We decompose the variational distribution $q$ on the latent labels into a product of two terms, one is related exclusively to the labels of the individuals present at $t_0$, the other one represents the conditional distribution of the labels of the remaining individuals (those that
arrive after $t_0$, that we call in the sequel the newborns) given the size at $t_0$ of the K communities,

$$q(Z) = q(\{Z_i\}_{i \notin V_0} | Z^{t_0}) q(Z^{t_0}) = q(\{Z_i\}_{i: \tau_i^b \in T_B} | Z^{t_0}) q(Z^{t_0}), \quad (4)$$

where $Z^{t_0} = \{Z_i\}_{i \in V_0}$ are the labels of individuals present at time $t_0$ and $\{Z_i\}_{i: \tau_i^b \in T_B}$ are the labels of individuals entering the system after $t_0$.

1. Individuals present at time $t_0$. We assume independence between the label vectors of individuals present at the initial time,

$$q(Z^{t_0}) = \prod_{i \in V_0} q(Z_i) = \prod_{i \in V_0} \prod_{k=1}^{K} \delta(i, k)^{Z_{ik}}, \quad (5)$$

where $\delta(i, k) = q(Z_{ik} = 1)$ are variational parameters in $Q$ and satisfy $\sum_k \delta(i, k) = 1$ for all $i \in V_0$.

The distribution of the number of individuals at time $t_0$ across the K communities can be viewed as a generalization of the multinomial distribution, known as the Poisson Multinomial Distribution (PMD), see [26]. Unlike the classical multinomial case, here the success probabilities are not identical across individuals: each individual $i \in V_0$ is associated with a probability vector $\delta(i, k)$, with

$$\sum_{k=1}^{K} \delta(i, k) = 1, \quad \forall i.$$

Let $N_1, N_2, \dots, N_K$ denote the community counts at time $t_0$. By construction, the random variables satisfy the constraint

$$\sum_{k=1}^{K} N_k = N_0,$$

where $N_0$ is the population size at time $t_0$.

Our interest lies in the marginal distributions of the components $N_1, N_2, \dots, N_K$. For each community $k$ and count $n$, we define the marginal distribution by

$$\gamma_{\text{mar}}(0, k, n) = q(N_k = n) = q(L_{k,n}^0 = 1),$$

which corresponds to a Poisson binomial distribution (see [27]) with success probabilities $\{\delta(i, k)\}_{i \in V_0}$. More precisely,

$$q(L_{k,n}^0 = 1) = \sum_{A_n \subset V_0 : |A_n| = n} \prod_{i \in A_n} \delta(i, k) \prod_{j \in A_n^c} (1 - \delta(j, k)), \quad (6)$$

where $A_n^c$ denotes the complementary of $A_n$ in $V_0$.
2. Individuals born after $t_0$ (newborns). The community memberships of individuals not present at $t_0$ are modelled conditionally on $Z^{t_0}$ as

$$q(\{Z_i\}_{i \notin V_0} | Z^{t_0}) = q(\{Z_i\}_{i: \tau_i^b \in T_B} | Z^{t_0}), \quad (7)$$

Given the community memberships of individuals present at $t_0$, the joint distribution of the community assignments of all newborns is fully characterised by the joint distribution of the community sizes at the event times, that is

$$\{Z_i\}_{i: \tau_i^b \in T_B} | Z^{t_0} \stackrel{D}{=} (L^1, \dots, L^M | Z^{t_0}),$$

where $L^\ell = (L_{k,n}^\ell)_{k,n}$ denotes the collection of community-size indicators at the event time $\tau_\ell$.

Since each event time corresponds to either a single birth or a single death, the sequence of community sizes $(L^0, \dots, L^M)$ forms a Markov chain in the generative model: at any event time $\tau_\ell$ the distribution of $L^\ell$ given the past depends on the history only through $L^{\ell-1}$. However, in our inferential setting we condition on the full event history $(\tau, b(\tau))$, which records not only whether a birth or a death occurs at each $\tau_\ell$, but also the identity $i_\ell$ of the individual involved in the event. This additional information breaks the Markov property at death times: when $b_\ell = -1$, $L^\ell$ depends on the past not only through $L^{\ell-1}$ but also through the latent memberships $Z_i$ of individuals alive before $\tau_\ell$.

At $\tau_\ell$, the knowledge of the past $(L^{1:\ell-1}, Z^{t_0})$ is equivalent to the knowledge of $(L^{\ell-1}, \{Z_i\}_{i \in V(\tau_0: \tau_{\ell-1})})$, where $L^{1:\ell-1} = (L^1, \dots, L^{\ell-1})$ and $V(\tau_0: \tau_{\ell-1}) = \cup_{j=0}^{\ell-1} V(\tau_j)$. Therefore, we have under $q$

$$L^\ell | L^{1:\ell-1}, Z^{t_0} \stackrel{D}{=} L^\ell | L^{\ell-1}, \{Z_i\}_{i \in V(\tau_0: \tau_{\ell-1})}.$$

At birth times $(\tau_\ell \in T_B)$, the Markov property gives rise to the following simplification

$$L^\ell | L^{1:\ell-1}, Z^{t_0} \stackrel{D}{=} L^\ell | L^{\ell-1}.$$

Under our variational approximation, Equation (7) is therefore factorized as follows,

$$\begin{align*}
q(\{Z_i\}_{i \notin V_0} | Z^{t_0}) &= q(L^1, \dots, L^M | Z^{t_0}) \\
&= \prod_{\ell=1}^{M} q(L^\ell | L^{\ell-1}, \{Z_i\}_{i \in V(\tau_0: \tau_{\ell-1})}) \\
&= \prod_{\ell=1}^{M} q(L^\ell | L^{\ell-1})^{1\{b_\ell=1\}} q(L^\ell | L^{\ell-1}, \{Z_i\}_{i \in V(\tau_0: \tau_{\ell-1})})^{1\{b_\ell=-1\}}.
\end{align*}$$

At death times, once the identity of the departing individual is known, the updated community sizes are deterministically obtained from $L^{\ell-1}$ and $\{Z_i\}_{i \in V(\tau_0: \tau_{\ell-1})}$, so that the conditional distribution reduces to,
$$q(L^{\ell} | L^{\ell-1}, \{Z_i\}_{i \in V(\tau_0: \tau_{\ell-1})})^{1\{b_{\ell}=-1\}} = 1.$$

We can then rewrite Equation (7) as follows,

$$q(\{Z_i\}_{i \notin V_0} | Z^{t_0}) = \prod_{\tau_{\ell} \in T_B} \prod_{k=1}^{K} \prod_{n=1}^{N(\tau_{\ell}-1)} \gamma(\ell, k, n, n)^{L_{k,n}^{\ell-1} L_{k,n}^{\ell}} \gamma(\ell, k, n, n+1)^{L_{k,n}^{\ell-1} L_{k,n+1}^{\ell}}, \quad (8)$$

where $\gamma(\ell, k, n, n')$ denotes the transition probability to move from $n$ to $n'$ individuals in community $k$ at time $\tau_\ell$, that is

$$\gamma(\ell, k, n, n') = q(L_{k,n'}^{\ell} = 1 | L_{k,n}^{\ell-1} = 1),$$

which are constrained by

$$\gamma(\ell, k, n, n + b_{\ell}) + \gamma(\ell, k, n, n) = 1 \quad \text{for all } \ell \ge 1, \\ \text{and } \gamma(\ell, k, n, n') = 0 \quad \text{for all } n' \notin \{n, n + b_{\ell}\}, \quad \text{for all } \ell \ge 1.$$

To summarize, the variational distribution $q$ in (4) is fully specified by Equations (5) and (8), and is parametrized by $\delta(i, k)$ for all $i \in V_0$, $k = 1, \dots, K$, and by the transition probabilities $\gamma(\ell, k, n, n')$ for all $k, n, n'$ and event times $\tau_\ell \in T_B$. Intuitively, the $\delta(i, k)$'s control the community memberships at the initial time $t_0$, while the $\gamma(\ell, k, n, n')$'s govern the evolution of community sizes at subsequent birth events.

Once the transition probabilities $\gamma(\ell, k, n, n')$ have been introduced, we need to relate them to the membership probabilities $\delta(i, k)$ for individuals born after $t_0$. To this end, we work with the marginal distributions of the community sizes at all event times $\tau_\ell \in \tau$. In what follows, we first define these marginal probabilities $\gamma_{\text{mar}}(\ell, k, n)$ and then express the probabilities $\delta(i_\ell, k)$ for individuals $i_\ell$, who arrive at $\tau_\ell$ in the system, in terms of $\gamma_{\text{mar}}$ and $\gamma$.

(a) Marginal distributions of community sizes. For each $\ell, k$ and $n$, we define the marginal distribution of $L_{k,n}^\ell$ by

$$\gamma_{\text{mar}}(\ell, k, n) = q(L_{k,n}^{\ell} = 1),$$

which satisfies the recursion

$$\gamma_{\text{mar}}(\ell, k, n) = \sum_{n' \in \{n, n-b_{\ell}\}} \gamma(\ell, k, n', n) \gamma_{\text{mar}}(\ell-1, k, n'), \quad \forall \ell \ge 1, k, n, \quad (9)$$

and which have been already introduced for $t_0$, by $\gamma_{\text{mar}}(0, k, n)$, defined in (6). Note that the ELBO does not explicitly depend on the transition probabilities $\gamma(\ell, k, n, n')$ at event times $\tau_\ell \in T_D$. Nevertheless, the marginals $\gamma_{\text{mar}}(\ell, k, n)$ are required at every event time (both births and deaths) in order to apply
the recursion (9) and propagate information forward in time, so the full set of transitions is needed to compute all $\gamma_{\text{mar}}(\ell, k, n)$.

(b) Community membership probabilities induced by $\gamma$ and $\gamma_{\text{mar}}$. Once probabilities $\gamma$ and $\gamma_{\text{mar}}$ are specified, we can express the variational community membership probabilities at birth times in terms of these transitions. The key observation is that the event “the newborn $i_\ell$ at time $\tau_\ell$ belongs to community $k$” exactly coincides with an increase by one unit in the size of community $k$ at $\tau_\ell$. Therefore, the quantity $\delta(i_\ell, k)$ can be written as

$$ \delta(i_\ell, k) = \sum_{n=1}^{N(\tau_\ell-1)} \gamma(\ell, k, n, n+1) \gamma_{\text{mar}}(\ell-1, k, n) \\ = \sum_{n=1}^{N(\tau_\ell-1)} [1 - \gamma(\ell, k, n, n)] \gamma_{\text{mar}}(\ell-1, k, n), $$

where $i_\ell$ denotes the newborn at time $\tau_\ell$ and $\gamma_{\text{mar}}(\ell-1, k, n)$ is the marginal probability that community $k$ has size $n$ just before $\tau_\ell$.

### 3.3.2 Variational expectation maximization inference

Proposition 1 and the family of variational distributions lead to the following expression for the ELBO function given by Equation (3):

$$ \begin{aligned} \mathcal{L}(q) &= \mathbb{E}_q (\log(\mathbb{P}_{\tilde{\theta}}(Z, \tau, e(\mathcal{T}), b(\tau))) - \mathbb{E}_q (\log q(Z)) \\ &= \sum_{i<j} \sum_{k_1,k_2} \delta(i, k_1)\delta(i, k_2) \sum_{\ell \in \mathcal{Y}_{ij}} \log \phi(e_{ij}(t_\ell), k_1, k_2) + |T_B| \log \hat{\lambda} + |T_D| \log \hat{\mu} - (\hat{\lambda} + \hat{\mu}) I_N \\ &\quad + \sum_{\tau_\ell \in T_B} \sum_k \sum_{n=1}^{N(\tau_\ell-1)} \gamma(\ell, k, n, n+1) \gamma_{\text{mar}}(\ell-1, k, n) \log n + \sum_{i \in V_0} \sum_k \delta(i, k) \log \beta_k \\ &\quad - \sum_{i \in V_0} \sum_k \delta(i, k) \log \delta(i, k) - \sum_{\tau_\ell \in T_B} \sum_k \sum_{n=1}^{N(\tau_\ell-1)} [\gamma(\ell, k, n, n) \gamma_{\text{mar}}(\ell-1, k, n) \log \gamma(\ell, k, n, n) \\ &\quad + \gamma(\ell, k, n, n+1) \gamma_{\text{mar}}(\ell-1, k, n) \log \gamma(\ell, k, n, n+1)]. \end{aligned} \quad (10) $$

maximizing $\mathcal{L}(q)$ over $q \in \mathcal{Q}$ is equivalent to optimizing over the variational parameters $\delta(i, k)$ for $i \in V_0$ and $\gamma(\ell, k, n, n')$ under the constraints

$$ \left\{ \begin{array}{ll} (a) & \sum_k \delta(i, k) = 1, \quad \forall i, \\ (b) & \sum_k \sum_{n=0}^{N(\tau_\ell-1)} \gamma(\ell, k, n, n) \gamma_{\text{mar}}(\ell-1, k, n) = K-1, \quad \forall \ell \ge 1, \\ (c) & \gamma(\ell, k, n, n+b_\ell) + \gamma(\ell, k, n, n) = 1, \quad \forall \ell \ge 1. \end{array} \right. \quad (11) $$
Note that relations (b) and (c) imply,

$$ \sum_k \sum_{n=0}^{N(\tau_{\ell}-1)} \gamma(\ell, k, n, n + b_{\ell}) \gamma_{\text{mar}}(\ell - 1, k, n) = 1, \quad \forall \ell \ge 1, $$

which expresses the fact that exactly one community size changes at each event time.

Adapting Variational Expectation Maximization algorithms to our setting leads to the pseudo-code in Algorithm 1.

**Algorithm 1 Variational Expectation Maximization**

1. **Input:** Observed data $(\tau, e(\mathcal{T}), b(\tau))$; initial parameters $(\pi^{(0)}, \beta^{(0)})$; initial variational distribution $q^{(0)}$; pre-estimated birth-death parameters $(\hat{\lambda}, \hat{\mu})$.

2. **Output:** Estimated parameters $(\pi, \beta)$ and optimal variational distribution $q \in \mathcal{Q}$.

3. Set iteration counter $m = 0$

4. **repeat**

5. **VE-Step:** ▷ Update the variational distribution $q^{(m+1)}$

6. Maximize the ELBO with respect to $q$

$$ q^{(m+1)} = \arg\max_{q \in \mathcal{Q}} \{\mathbb{E}_q[\log \mathbb{P}_{\tilde{\theta}^{(m)}}(Z, \tau, e(\mathcal{T}), b(\tau))] - \mathbb{E}_q[\log q(Z)]\} $$

7. **VM-Step:** ▷ Update the model parameters $\tilde{\theta}^{(m+1)} = (\hat{\lambda}, \hat{\mu}, \pi^{(m+1)}, \beta^{(m+1)})$

8. Maximize the ELBO with respect to $(\pi, \beta)$ a part of $\tilde{\theta} = (\hat{\lambda}, \hat{\mu}, \pi, \beta)$:

$$ (\pi^{(m+1)}, \beta^{(m+1)}) = \arg\max_{(\pi, \beta)} \mathbb{E}_{q^{(m+1)}}[\log \mathbb{P}_{\tilde{\theta}}(Z, \tau, e(\mathcal{T}), b(\tau))] $$

9. Increment iteration counter $m = m + 1$

10. **until** a convergence criterion is met.

In the VE-step, we must determine the $\delta(i, k)$ for all $i \in V_0$ and the $\gamma(\ell, k, n, n)$ for all birth times $\tau_{\ell} \in T_B$ that maximize $\mathcal{L}(q)$ under the constraints (11). Recall that under $q$, $\delta(i, k)$ is the variational community membership distribution for individual $i$, whereas $\gamma(\ell, k, n, n)$ encodes the transition probabilities for community-size changes.

**Proposition 2** For $\ell \in T_B$, we refer to $i_{\ell}$ as the newborn at time $\tau_{\ell}$, and we denote by $\hat{\gamma}(\ell, k, n, n)$ the quantity that maximizes $\mathcal{L}(q)$ with respect to $\gamma(\ell, k, n, n)$ when $\tilde{\theta} = (\hat{\lambda}, \hat{\mu}, \pi, \beta)$ is fixed. Then there exists a unique $\rho_{\ell} > 0$ such that

$$ \hat{\gamma}(\ell, k, n, n) = \frac{\rho_{\ell}}{\rho_{\ell} + n \prod_{j \neq i_{\ell}} \prod_{k'} \prod_{\ell' \in \Upsilon_{i_{\ell}j}} [\phi(e_{ij}(t_{\ell'}), k, k')]^{\tilde{\delta}(j, k')}}, \quad \text{for } n \in [0, N_{\ell-1}[ $$

where the parameter $\rho_{\ell}$ is determined by the constraint (b) in (11).
*Proof* The result follows by maximizing $\mathcal{L}(q)$ with respect to $\gamma(\ell, k, n, n)$, using Lagrange multipliers for the constraints in (11). See Appendix A.3 for details. □

*Remark 1* In order to compute $\hat{\gamma}(\ell, k, n, n)$ in practice, one must first determine $\rho_{\ell}$. Since the defining equation for $\rho_{\ell}$ has no closed-form solution in general, $\rho_{\ell}$ is obtained numerically.

Using Proposition 2 together with (c) in (11) immediately yields, for $\ell \in T_B$,

$$ \hat{\gamma}(\ell, k, n, n+1) = 1 - \hat{\gamma}(\ell, k, n, n). $$

Optimizing $\mathcal{L}(q)$ with respect to the initialization parameters $\delta(i, k)$ for $i \in V_0$ is more delicate. Following the approximation strategy of [6, 28], we neglect the dependence of certain terms in $\mathcal{L}(q)$ on $\delta(i, k)$ and update $\delta(i, k)$ by solving the fixed-point equation

$$ \hat{\delta}(i, k) \propto \beta_k \prod_{j \neq i} \prod_{k'} \prod_{\ell' \in \Upsilon_{ij}} [\phi(e_{ij}(t_{\ell'}), k, k')]^{\hat{\delta}(j, k')}, \quad (12) $$

where $\phi(e_{ij}(t_{\ell'}), k, k')$ and $\Upsilon_{ij}$ are defined in Proposition 1. In simulations (Section 4.1), this approximation performs well. For completeness, we provide in the Appendix A.3 the derivative used in this update.

For $i \notin V_0$, the computation of $\delta(i, k)$ does not entail $\mathcal{L}(q)$ optimization. Indeed, the event that the newborn $i_{\ell}$ at time $\tau_{\ell}$ belongs to community $k$ is equivalent to the size of community $k$ increasing by one at $\tau_{\ell}$, we can express the variational community membership probability $\hat{\delta}(i_{\ell}, k)$ in terms of the transition probabilities $\hat{\gamma}(\ell, k, n, n+1)$ as

$$ \hat{\delta}(i_{\ell}, k) = \sum_{n=1}^{N(\tau_{\ell}-1)} \hat{\gamma}(\ell, k, n, n+1) \hat{\gamma}_{\text{mar}}(\ell-1, k, n), \quad (13) $$

where

$$ \hat{\gamma}_{\text{mar}}(\ell - 1, k, n) = \begin{cases} \sum_{n' \in \{n, n-b_{\ell}\}} \hat{\gamma}(\ell - 1, k, n', n) \hat{\gamma}_{\text{mar}}(\ell - 2, k, n'), & \text{for } \ell - 1 > 0, \\ \sum_{A_n \subset V_0: |A_n| = n} \prod_{i \in A_n} \hat{\delta}(i, k) \prod_{j \in A_n^c} (1 - \hat{\delta}(j, k)), & \text{for } \ell - 1 = 0, \end{cases} \quad (14) $$

using the PMD expression (6). The term $n = 0$ is excluded from the sum in (13), since we cannot have any birth event on a community once this community is empty.

We optimize the ELBO explicitly only with respect to the transition probabilities at birth times, since these are the only ones appearing directly in (10). Transition probabilities at death times do not appear in the ELBO in closed form since they are already determined implicitly by the probabilities of previous events, because we observe the individual involved in each death event. Computing them exactly would require summing over all possible past configurations of community memberships. In
other words, it requires to do the deconvolution of the γ probabilities backwards in time up until the time of birth of the departing individual (so, as many times as the number of events that take place during the lifetime of this individual), which is intractable,

$$ \gamma(\ell, k, n, n-1) = q(L_{k,n-1}^{\ell} = 1 | L_{k,n-1}^{\ell-1} = 1) \\ = \sum_{Z_i \in \{0,1\}^K, \forall i \in V(\tau_0: \tau_{\ell-1})} q(L_{k,n-1}^{\ell} = 1, \{Z_i\}_{i \in V(\tau_0: \tau_{\ell-1})} | L_{k,n-1}^{\ell-1} = 1). $$

These probabilities are still needed to propagate the marginals via (14). Therefore, instead of computing them directly, for $b_{\ell} = -1$ and individual $i_{\ell}$ dying at time $\tau_{\ell}$, we approximate $\gamma(\ell, k, n+1, n)$ by enforcing the additional constraint

$$ \delta(i_{\ell}, k) = \sum_{n=0}^{N(\tau_{\ell-1})-1} \gamma(\ell, k, n+1, n) \gamma_{\text{mar}}(\ell-1, k, n+1), \quad (15) $$

which simply states that the probability that the size of community $k$ decreases by one at $\tau_{\ell}$ must coincide with the probability that $i_{\ell}$ belongs to community $k$.

Observing that $\gamma(\ell, k, n+1, n)$ should increase with $n$, we adopt the following parametric form

$$ \gamma(\ell, k, n+1, n) = \frac{\rho_{\ell k} n}{1 + \rho_{\ell k} n}, \quad (16) $$

where $\rho_{\ell k} \ge 0$ is a parameter specific to the $\ell$-th BD (which is a death) event and community $k$. We show in Appendix A.3.1 that there is a unique solution of this form satisfying (15).

Let $\hat{\rho}_{\ell k}$ being the positive numerical solution of the next equation in $\rho_{\ell k}$

$$ \sum_{n=0}^{N(\tau_{\ell-1})-1} \frac{\rho_{\ell k} n}{1 + \rho_{\ell k} n} \hat{\gamma}_{\text{mar}}(\ell - 1, k, n + 1) = \hat{\delta}(i_{\ell}, k). $$

which is the constraint (15) with $\gamma(\ell, k, n+1, n)$ replaced by its parametric expression.

In the VM-step of the algorithm, we update the parameters $\pi_{k_1 k_2}$ and $\beta_k$ by maximizing the ELBO function $\mathcal{L}(\cdot)$ written in equation (10) for fixed $q$.

**Proposition 3** Let $q \in \mathcal{Q}$ be fixed. Then the maximizers of $\mathcal{L}(q)$ with respect to $\pi_{k_1 k_2}$ and $\beta_k$ are given by

$$ \hat{\pi}_{k_1 k_2} = \frac{\sum_{i \neq j} \delta(i, k_1) \delta(j, k_2) \sum_{\ell \in \Upsilon_{ij}} e_{ij}(t_{\ell})}{\sum_{i \neq j} \delta(i, k_1) \delta(j, k_2) |\Upsilon_{ij}|}, $$

where $\Upsilon_{ij}$ is defined in Proposition 1, and

$$ \hat{\beta}_k = \frac{1}{|V_0|} \sum_{i \in V_0} \delta(i, k). $$
*Proof* For a fixed $q \in \mathcal{Q}$, the result is straightforward by solving the maximization of the ELBO function $\mathcal{L}(q)$ with respect to $\pi_{k_1 k_2}$. Similarly, maximizing $\mathcal{L}(q)$ with respect to $\beta_k$, under the constraint $\sum_k \beta_k = 1$, gives the stated formula for $\hat{\beta}_k$. Detailed calculations are given in Appendix A.4. $\square$

VEM was first introduced in the SBM context by [19], who derived the variational updates for the static binary SBM under a mean-field variational family. Our variational expressions can be viewed as natural temporal extensions of this classical mean-field framework for stochastic block models. See, for instance, [6] for a dynamic generalisation of the static SBM using mean-field variational inference in the case of a fixed population size, and related developments in other directions.

Our results recover the same functional form for the updates of $(\delta(i, k))_{i,k}$ for $i \in V_0$ and of $(\pi_{k_1 k_2}, \beta_k)$, but embed them in a richer temporal model where the population size is itself random and evolves according to a birth-death mechanism. In particular, the additional variational parameters $(\gamma, \gamma_{\text{mar}})$ explicitly encode the evolution of community sizes and enforce consistency between individual memberships and the latent population process.

From this perspective, our BD-SBM can be seen as a dynamic generalisation of the static variational SBM. When there are no births or deaths (so that the population size remains constant), our updates reduce to the standard VEM updates of a dynamic SBM with a fixed vertex set.

## 3.4 Model selection with the Integrated Completed Likelihood (ICL)

A key aspect of fitting a Stochastic Block Model is the choice of $K$, the number of latent communities, which is generally unknown. Using too few classes may oversimplify the network structure, whereas using too many may lead to overfitting. To select an appropriate value of $K$, we rely on the Integrated Completed Likelihood (ICL) criterion [19, 29].

For any number of communities $K \ge 1$, we make explicit the dependence of the model parameters on $K$ by writing $\theta_K = \{\lambda, \mu, \beta_K, \pi_K\}$. For fixed $K$, the corresponding estimator is denoted by $\hat{\theta}_K = (\hat{\lambda}, \hat{\mu}, \hat{\beta}_K, \hat{\pi}_K)$, as described in Sections 3.2 and 3.3.2.

Given $K$, we denote by $\hat{Z}$ and $\hat{L}$ the maximum a posteriori (MAP) estimates of the community memberships and community sizes, respectively, under the VEM. The complexity of our model, as a function of $K$, is driven by:

1. the $K-1$ degrees of freedom associated with the proportion parameters $\beta_K$ for the $|V_0|$ individuals present in the system at time $t_0$;

2. the $\frac{K(K+1)}{2}$ block connection probabilities $\pi_{k_1 k_2}$ governing the $\sum_{i<j} |\Upsilon_{ij}|$ Bernoulli trials;

3. the birth-death parameters $(\lambda, \mu)$, which do not depend on $K$ and are therefore not included in the penalisation term.
The ICL criterion is then defined as follows,

$$ \text{ICL}(K) = \log \mathbb{P}_{\hat{\theta}_K} \left( e(\mathcal{T}), \hat{Z}, \tau, b(\tau) \right) - \frac{1}{2}(K-1) \log N_0 - \frac{K(K+1)}{2} \log \left( \sum_{i<j} |\Upsilon_{ij}| \right). $$

Maximizing ICL($K$) over candidate values of $K$ provides a commonly use for selecting the number of communities. In practice, we fit the model for a range of values of $K$, compute the ICL for each, and select the value that maximizes the criterion. This approach balances model fit and complexity while explicitly accounting for the uncertainty in the community assignments.

## 3.5 Initialization of the Variational EM Algorithm

The initialization of the latent variables is a crucial step for the variational expectation-maximization (VEM) algorithm, since poor starting values may lead to convergence towards suboptimal local maxima. We describe here the procedure used to initialize the community membership probabilities $\delta(i, k)$ for all $i \in V$ and $k = 1, \dots, K$. Throughout this subsection, $\delta^{(0)}(i, k)$ denotes the initial value of $\delta(i, k)$.

To obtain a first partition of the nodes into $K$ communities, we build a similarity score between each pair of nodes $(i, j)$, defined by

$$ s_{ij} = \frac{2 \sum_{\ell \in \Upsilon_{ij}} e_{ij}(t_{\ell}) - |\Upsilon_{ij}|}{|\Upsilon_{ij}|}, $$

which measures the average tendency of nodes $i$ and $j$ to be connected during the time interval when they are simultaneously present in the network. Let $S = (s_{ij})$ denote the resulting similarity matrix. We apply a $K$-means clustering algorithm to $S$ and obtain an initial hard partition of the nodes into $K$ groups. The associated binary indicators are denoted by $\hat{Z}_{ik}^{(0)}$, where $\hat{Z}_{ik}^{(0)} = 1$ if node $i$ is assigned to cluster $k$, and 0 otherwise.

We then construct the initial variational membership probabilities $\delta^{(0)}(i, k)$ so as to be consistent with this hard clustering while allowing for uncertainty. For a fixed weight $\omega \in (0, 1)$, we set

$$ \delta^{(0)}(i, k) = \omega \hat{Z}_{ik}^{(0)} + \frac{1 - \omega}{K}, \quad i \in V, k = 1, \dots, K. $$

Thus each node receives a larger probability $\omega + (1 - \omega)/K$ for its assigned cluster and a smaller, uniform probability $(1 - \omega)/K$ for the remaining communities.

We decided to initialize the model parameters $\beta$ and $\pi$ according to the following formulas.

1. The initial community proportions are given by

$$ \beta_k^{(0)} = \frac{1}{|V_0|} \sum_{i \in V_0} \delta^{(0)}(i, k), $$
that is, the average assignment probability of individuals in $V_0$ to community $k$.

2. The initial block connection probabilities are set to

$$ \pi_{k_1 k_2}^{(0)} = \frac{\sum_{i<j} \delta^{(0)}(i, k_1) \delta^{(0)}(j, k_2) \sum_{\ell \in \mathcal{Y}_{ij}} e_{ij}(t_\ell)}{\sum_{i<j} \delta^{(0)}(i, k_1) \delta^{(0)}(j, k_2) |\mathcal{Y}_{ij}|} $$

This initialization scheme provides informative starting values that reflect the empirical clustering structure of the network, while preserving enough flexibility for the VEM algorithm to refine the variational distribution in subsequent iterations. In particular, the initialization does not enforce the birth-death structure on community sizes, and the constraints induced by the birth-death dynamics are only imposed through the subsequent VEM updates.

# 4 Experiments

We evaluate the proposed model through two complementary types of experiments: two controlled simulated data sets and a real-world temporal co-authorship network built from arXiv. An open-source Python implementation of the methods described in this paper is available at this GitHub repository.

With simulated data, we assess the performance of our model through clustering accuracy, the recovery of group memberships, and the accuracy of parameter estimation. The arXiv application demonstrates how the model behaves on a large, sparse network with heterogeneous communities and genuine temporal turnover, providing evidence of its practical relevance.

## 4.1 Simulated data

We generate two synthetic data sets, one according to our birth-death Stochastic Block Model (BD-SBM) and the other one simulated from a pure birth process embedded within a dynamic Stochastic Block Model (B-SBM).

For the BD-SBM, the birth-death process is simulated with birth rate $\lambda = 0.04$ and death rate $\mu = 0.02$ starting at $t_0 = 0$ and end up at $t_T = 150$. We consider $K = 4$ latent communities with initial sizes 10, 11, 7, and 12, yielding an initial population $N_0 = 40$.

At $t_T = 150$, the total number of distinct individuals who have been involved in the process reaches $N = 1208$.

For the pure birth process within a dynamic SBM ("B-SBM"), the birth process is generated according to the same parameters as for the birth-death process described above except that the death rate is set to $\mu = 0$ and the horizon is $t_{end} = 100$, while keeping $\lambda = 0.04$. In this setting, At $t_T = 100$, the total number of distinct individuals who have been involved in the process reaches $N = 2601$.

For both scenarios, we construct two dynamic SBMs by specifying different within- and between-community connection probabilities. The two regimes are encoded by the following symmetric block matrices:
$$ \pi_{\text{low}} = \begin{pmatrix} 0.05 & 0.09 & 0.05 & 0.04 \\ 0.09 & 0.10 & 0.055 & 0.06 \\ 0.05 & 0.055 & 0.20 & 0.07 \\ 0.04 & 0.06 & 0.07 & 0.06 \end{pmatrix} \qquad \pi_{\text{high}} = \begin{pmatrix} 0.75 & 0.36 & 0.20 & 0.16 \\ 0.36 & 0.91 & 0.22 & 0.24 \\ 0.20 & 0.22 & 0.82 & 0.28 \\ 0.16 & 0.24 & 0.28 & 0.66 \end{pmatrix} $$

The first matrix, $\pi_{\text{low}}$, corresponds to a low-signal regime in which within- and between-community probabilities are relatively close (resulting in a sparser and less clearly separated network connection probabilities), whereas $\pi_{\text{high}}$ represents a high-signal regime with much stronger within-community connectivity.

Network snapshots are recorded at equally spaced times over the interval $[t_0 = 0, t_T]$, with one snapshot per unit of time, resulting in 150 snapshots for the BD-SBM and 100 snapshots for the pure birth model.

Note that in our implementation, a single VEM iteration requires

$$ \mathcal{O}(KK + MKN) = \mathcal{O}(MKN) \quad \text{operations, since usually } MN \gg K, $$

where $M$ is the number of birth-death events, $N = |V| = |\cup_{t \in [t_0, t_M]} V(t)|$ is the total population size over the observation window, and $K$ is the number of communities. This cost is dominated by the updates of the pairwise edge contributions and community responsibilities across all event times.

### 4.1.1 Birth-death process within dSBM

From the simulated dataset, we first choose the optimal $K$ using the ICL criterion described in section 3.4. In particular, we perform 50 random initializations of the VEM algorithm and, for each of them, we retain the value of $K$ that maximizes the ICL. Figure 3 reports, for each regimes, how often each value of $K$ is selected across the random runs. For both connectivity regimes, the best trade-off between goodness-of-fit and model complexity is obtained for $K = 4$, which coincides with the true number of communities used to generate the data. Second, we apply the proposed VEM algorithm with $K = 4$.

As shown in Figure 4, our algorithm converges quickly, in fewer than 10 iterations. As expected, the algorithm converges more rapidly under the high-signal regime $\pi_{\text{right}}$ than under the low-signal regime $\pi_{\text{low}}$. Since our implementation relies on an approximate version of the VEM, the classical monotonicity property of the ELBO cannot be theoretically guaranteed. Nevertheless, in practice, the ELBO trajectories are nearly monotonic, rapidly reach a plateau, and lead to highly accurate parameter estimates.

Figure 5 shows the posterior community membership probabilities for all individuals in the two data sets. For both regimes, the posterior probabilities are typically very close to 0 or 1, especially in the high-signal regime, indicating that individuals are assigned to a single community with high confidence. Most exceptions involve the most recent individuals (those arriving later in the process), particularly in the low-signal case, where fewer interactions are available. Overall, the posterior distributions are very similar across the two regimes, with sharper assignments under $\pi_{\text{high}}$.

Figure 6 compares the true and inferred community memberships. The confusion matrices in Figures 6a and 6b show that, in both scenarios, the vast majority of
<figure><img src="image_8.png" /><figcaption>**Fig. 3:** Number of times the ICL criterion selects each value of *K* as the optimal number of communities over all random initializations in the BD-SBM experiment with $\pi_{\text{low}}$ (left) and $\pi_{\text{high}}$ (right).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 4:** Evolution of the ELBO over VEM iterations for the two simulated BD-SBM networks with *K* = 4 communities: $\pi_{\text{low}}$ (left) and $\pi_{\text{high}}$ (right).</figcaption></figure>

individuals are assigned to their correct communities, with only limited off-diagonal elements.

Figures 6c and 6d display the evolution of clustering accuracy as a function of time (i.e., as the population size evolves). For the low-signal case $\pi_{\text{low}}$, accuracy is relatively unstable at early times, then increases and remains high during the middle part of the process before decreasing slightly at the very end of the observation window. This behaviour is expected: information about individuals present at $t_0$ is limited, whereas for individuals arriving later, both their entry time and the population size at entry are known, and their subsequent interactions accumulate over time. Examining accuracy as a function of the proportion of initial individuals $N_0$ relative to the total population, we find that accuracy increases as this ratio decreases. In the high-signal case $\pi_{\text{high}}$, accuracy rapidly reaches a high level and remains largely stable, with only a slight decline for individuals arriving at the very end of the observation window,
<figure><img src="image_10.png" /><figcaption>Fig. 5: Posterior community membership probabilities for the two simulated BD-SBM networks with $K = 4$ communities. Left: data generated with $\pi_{low}$; right: data generated with $\pi_{high}$.</figcaption></figure>

those who have very few recorded interactions. As expected, overall performance is superior in the high-signal scenario, where intra-community connection probabilities clearly outweigh inter-community ones.

The birth and death rates $\lambda$ and $\mu$ do not depend on the latent communities, and their estimates are very close to the true values. Table 1 reports both the true and the estimated values of the birth and death rates.

<table><caption>Table 1: True vs. estimated birth-death parameter values.</caption><thead><tr><td>Parameter</td><td>True value</td><td>Estimated value</td></tr></thead><tbody><tr><td>λ</td><td>0.04</td><td>0.0400</td></tr><tr><td>μ</td><td>0.02</td><td>0.0207</td></tr></tbody></table>

The estimation of the remaining model parameters is summarized in Table 2 and Figure 7. For both regimes $\pi_{low}$ and $\pi_{high}$, the estimates of the community probabilities at $t_0$, i.e. $\beta$, are reported in Table 2. In the low-signal case, this estimation is slightly less accurate than in the high-signal case, which is consistent with the lower clustering performance at early times.

The estimated block probability matrices $\hat{\pi}$ in Figure 7 closely reproduce the true matrices $\pi_{low}$ and $\pi_{high}$ used in the simulations. While minor deviations appear in blocks with low connection probabilities, the relative ordering of intra- and inter-community connection strengths is preserved, indicating that the model effectively captures both structural heterogeneity and the impact of temporal population turnover.

Allowing for departures in our model makes the inference problem both statistically and algorithmically more challenging. In particular, when individuals are allowed to
<figure><img src="image_11.png" /><figcaption>(a) Confusion matrix for the data generated with $\pi_{low}$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(c) Evolution of clustering accuracy over time for $\pi_{low}$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(b) Confusion matrix for the data generated with $\pi_{high}$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(d) Evolution of clustering accuracy over time for $\pi_{high}$.</figcaption></figure>

<caption>**Fig. 6:** Comparison between true and inferred communities for the two simulated BD-SBM networks with $K = 4$ communities. Top: confusion matrices. Bottom: evolution of clustering accuracy over time as the population evolves.</caption>

leave the network, the process governing community sizes is no longer Markovian at departure times, rendering the exact VEM updates intractable. This motivates the approximation introduced in Section 3.3.2, which restores tractability at the cost of a slight weakening of the theoretical guarantees on ELBO monotonicity.

### 4.1.2 Pure birth process within dSBM

We now turn to the pure birth process within a dynamic SBM, where the population size may only increase over time ($\mu = 0$), while the within- and between-community connection probabilities are the same block matrices $\pi_{low}$ and $\pi_{high}$ defined in 4.1. This scenario serves as a baseline in which temporal dynamics arise exclusively from arrivals, allowing us to isolate the impact of departures on community recovery.

As for the BD-SBM case, we first choose the optimal $K$ using the ICL criterion described in section 3.4. In particular, we perform 50 random initializations of the
<figure><img src="image_15.png" /><figcaption>(a) Data generated with $\pi_{\text{low}}$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(b) Data generated with $\pi_{\text{high}}$.</figcaption></figure>

<caption>**Fig. 7:** Estimated block probability matrices $\hat{\pi}$ (right of each panel) compared to the true matrices $\pi$ (left) for the two BD-SBM scenarios.</caption>

VEM algorithm and, for each of them, we retain the value of *K* that maximizes the ICL. Figure 8 reports, separately for $\pi_{\text{low}}$ and $\pi_{\text{high}}$, how often each value of *K* is selected across the random runs. In both regimes, the ICL criterion most frequently selects *K* = 4, in agreement with the true number of communities used to generate the data.

We therefore fix *K* = 4 and apply the proposed VEM algorithm using the same initialization procedure as in the BD-SBM experiment (described in Subsection 3.5). In this pure birth setting, the algorithm converges in a small number of iterations. Figure 9 displays the evolution of the ELBO over the VEM iterations for both connectivity regimes. As in the BD-SBM case, the ELBO trajectories quickly stabilize, indicating that the algorithm reaches a local optimum after only a few updates.

Figure 10 displays the posterior community membership probabilities for all individuals in the two B-SBM data sets. We observe the same phenomena as for the BD-SBM: the posterior probabilities are typically very close to 0 or 1, indicating highly
<table><caption>**Table 2:** True vs. estimated community proportions $\beta$ for the BD-SBM.</caption><thead><tr><td>Parameter</td><td>True value</td><td>Estimated (π<sub>low</sub>)</td><td>Estimated (π<sub>high</sub>)</td></tr></thead><tbody><tr><td>β<sub>0</sub></td><td>0.25</td><td>0.2751</td><td>0.2525</td></tr><tr><td>β<sub>1</sub></td><td>0.275</td><td>0.2463</td><td>0.2725</td></tr><tr><td>β<sub>2</sub></td><td>0.175</td><td>0.1622</td><td>0.1750</td></tr><tr><td>β<sub>3</sub></td><td>0.30</td><td>0.3165</td><td>0.3000</td></tr></tbody></table>

<figure><img src="image_17.png" /><figcaption>**Fig. 8:** Number of times the ICL criterion selects each value of $K$ as the optimal number of communities over all random seeds in the BSBM experiment with $\pi_{low}$ (left) and $\pi_{high}$ (right).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 9:** Evolution of the ELBO over VEM iterations for the two simulated B-SBM networks with $K = 4$ communities: $\pi_{low}$ (left) and $\pi_{high}$ (right).</figcaption></figure>

confident assignments. The absence of departures makes the trajectory of each individual longer on average, which in turn yields more interaction data and thus facilitates community identification.

We next compare the inferred community labels to the true ones. Figure 11 reports the confusion matrices for $\pi_{low}$ and $\pi_{high}$. In both signal regimes, the model achieves
<figure><img src="image_19.png" /><figcaption>**Fig. 10:** Posterior community membership probabilities for the two simulated BSBM networks with *K* = 4 communities. Left: data generated with $\pi_{\text{low}}$; right: data generated with $\pi_{\text{high}}$.</figcaption></figure>

perfect clustering accuracy: off-diagonal entries in the confusion matrices are essentially zero, and the accuracy remains equal to one throughout the observation window. This sharp contrast with the BD-SBM results highlights the additional difficulty introduced by departures, which both require an approximation in the VEM algorithm and shorten some trajectories, thereby reducing the amount of information available for those individuals.

<figure><img src="image_20.png" /><figcaption>**Fig. 11:** Comparison between true and inferred communities for the two simulated BSBM networks with *K* = 4 communities. Top: confusion matrices. Bottom: evolution of clustering accuracy over time.</figcaption></figure>

Finally, we examine parameter estimation in the pure birth setting. Table 3 reports the true and estimated values of the birth rate $\lambda$. The estimated value is extremely
close to the true one, confirming that the model accurately recovers the underlying birth dynamics when no departures are present.

The estimated block probability matrices $\hat{\pi}$ in Figure 12 also closely match the true matrices $\pi_{low}$ and $\pi_{high}$, preserving the ordering of within- and between-community connection strengths in both regimes.

<figure><img src="image_21.png" /><figcaption>**Fig. 12:** Estimated block probability matrices $\hat{\pi}$ (right of each panel) compared to the true matrices $\pi$ (left) for the two BSBM scenarios.</figcaption></figure>

<table><caption>**Table 3:** True vs. estimated birth-death parameter values.</caption><thead><tr><td>Parameter</td><td>True value</td><td>Estimated value</td></tr></thead><tbody><tr><td>λ</td><td>0.04</td><td>0.0405</td></tr></tbody></table>
## 4.2 Data Extraction from arXiv

The dataset under study is constructed from the arXiv electronic repository. We focus on all articles categorized under Mathematics and published between January 2019 and September 2025. Metadata are retrieved through the official arXiv API, which provides structured access to publication records, including titles, author lists, subject categories, and submission dates.

Each article is treated as an undirected temporal interaction among its co-authors (distinct nodes) and makes the pair of co-authors connected by an edge dated at the time of the paper's first submission (the published field). In addition, to capture revision activity, we treat different versions of the same paper as distinct temporal interaction events between the same set of co-authors: version 1 is timestamped at published, while subsequent versions (if any) are timestamped at updated. The resulting temporal co-authorship network spans more than five years of activity and covers all mathematical subfields in arXiv (e.g., math.ST, math.PR, math.OO, etc., see table 5). To summarize, nodes represent individual authors while temporal edges represent collaboration events occurring at specific dates.

The raw scale of the Mathematics corpus over this period exceeds 200 000 works and 50 000 unique authors. Because daily-resolution snapshots are extremely sparse and full-graph inference would be computationally prohibitive, we first restrict attention to a denser working network of size N = 5000 authors, chosen to maximize edge density under this size constraint.

From this reduced dataset, we define the ancestor set $V_0$ as the authors who published at least once during the first six months (182 days) of 2019. The observation window for our analysis begins in July 2019. Starting from $V_0$, we iteratively expand the node set by adding all authors who co-author with any member of the current set at least once within the observation window, and we repeat this expansion until no new authors are added. The transitive co-authorship neighborhood grown from $V_0$, together with all co-authorship links among its members, forms the induced subgraph we analyze. Applying this construction yields $|V_0| = 2365$ ancestors and a total of N = 4622 authors (ancestors plus their retained collaborators) present in the sequence of interaction networks during the observation window. The observed total number of birth-death events is M = 3725.

To mitigate sparsity while retaining temporal structure, we discretize time into two-month bins, resulting in 38 time steps over the study period. Within each bin, two authors are linked if they co-author at least one paper (or version-event) during that bin. For authors who do not belong to $V_0$ (i.e., non-ancestors), we define their birth time as the date of their first observed publication, and their death time as min(last-publication-date+2 months, $t_T$), where $t_T$ denotes the end of the observation window. In the rare case of identical birth or death times for multiple authors, we break ties by discarding the author(s) with fewer publications, ensuring unique event times.

We fit our birth-death SBM model (BD-SBM) to this temporally aggregated co-authorship network.

*Model selection.* First, we perform model selection over the number of communities by running our VEM algorithm for $K \in \{2, 3, 4\}$ and computing the Integrated
<figure><img src="image_22.png" /><figcaption>**Fig. 13:** Adjacency structure of the arXiv Mathematics co-authorship network.</figcaption></figure>

Completed Likelihood (ICL) for each candidate. Given the large product $N \times M$, we deliberately restrict attention to small values of $K$ to keep computation tractable while still allowing for non-trivial community structure. The best trade-off between fit and complexity is obtained for $K = 4$, which attains the highest ICL and is therefore retained in the subsequent analysis.

*Initialization of the VEM algorithm.* As mentioned in section 3.5 the initialization of a VEM algorithm plays a crucial role. Given the observation of a very sparse dynamic networks, we therefore use a dynamic extension of the sparse SBM proposed in [30] to construct the initialization for the latent memberships. In practice, we first fit this dynamic sparse SBM to obtain coarse cluster assignments and then use the resulting soft memberships as initial values for our BD-SBM VEM algorithm.

*Post-inference results.* Figure 13 summarizes the global connectivity structure of the arXiv co-authorship network, showing the raw adjacency matrix in Figure 13a while displaying in Figure 13b the same matrix after permuting authors according to their inferred community labels. This highlights an overall low level of interaction among authors, with stronger interaction intensity observed within communities than between them.

The inference results are summarized in Figures 14 and 15. Figure 14a shows the posterior community membership probabilities for all authors, while Figure 14b reports the estimated block probability matrix $\hat{\pi}$. The ELBO trajectory in Figure 14c indicates that the algorithm converges and stabilizes after approximately 25 iterations.

To better understand the structural role of the inferred communities, we analyse mean connection rates by cluster (Figure 15). Figure 15a reports the average degree (mean connection rate) of authors in each cluster, while Figure 15b focuses on the mean rate of within-cluster connections. Cluster 2 exhibits a low within-cluster rate, suggesting a weakly connected community compared to the others.

Figure 16 provides a content-based interpretation of the inferred communities. For each author, we consider as “primary category” every arXiv primary subject in
<figure><img src="image_23.png" /><figcaption>(a) Posterior community membership.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(b) Estimated block probability matrix $\hat{\pi}$.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(c) ELBO over VEM iterations.</figcaption></figure>

**Fig. 14:** Inference results for the arXiv Mathematics co-authorship network.

<figure><img src="image_26.png" /><figcaption>(a) Individual connection rate (average degree) by cluster.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>(b) Individual within-cluster connection rate by cluster.</figcaption></figure>

**Fig. 15:** Boxplots of individual connection rates by inferred communities

which they have at least one paper during the study period (an author can therefore contribute to several categories). We retain the top 20 categories with the largest number of authors. For each category, we plot its proportion within each cluster. We observe that clusters differ not only in connectivity patterns but also in their predominant mathematical areas, indicating that the inferred communities capture both structural and thematic heterogeneity.

Cluster 2 gathers many authors who tend to collaborate less frequently (those with lower mean degree). This cluster is particularly enriched in Optimization and Control (math.OC), Numerical Analysis (math.NA), and authors whose primary labels combine these areas, such as [math.OC, math.NA]. Authors publishing in Probability (math.PR), or jointly in [math.PR, math.ST] and [math.PR, math.OC], also tend to fall into this low-collaboration cluster. Cluster 0 is dominated by Combinatorics (math.CO). Cluster 1 is strongly influenced by Analysis of PDEs (math.AP); it contains both authors who publish exclusively in math.AP and those whose primary labels combine math.AP with other areas, such as [math.AP, math.PR], [math.AP, math.OC] or [math.AP, math.DG].
<figure><img src="image_28.png" /><figcaption>**Fig. 16:** Top arXiv Mathematics subject categories by inferred community (proportion within each cluster).</figcaption></figure>

Cluster 3 includes a large fraction of authors in Number Theory (math.NT), Algebraic Geometry (math.AG) and their combination, e.g. [math.AG, math.CO]. See Table 5 for details on the subject-category labels.

Finally, Table 4 reports the estimated birth and death rates, together with the inferred community proportions $\beta_k$ for $K = 4$. The fitted model suggests a relatively small but non-zero departure rate, together with a dominant community (cluster 2) containing approximately 60% of the authors.

<table><caption>**Table 4:** Estimated BD-SBM parameters for the arXiv Mathematics co-authorship network.</caption><thead><tr><td>Parameter</td><td>Estimated value</td></tr></thead><tbody><tr><td>λ</td><td>0.271 · 10<sup>-3</sup></td></tr><tr><td>μ</td><td>0.176 · 10<sup>-3</sup></td></tr><tr><td>β<sub>0</sub></td><td>0.1420</td></tr><tr><td>β<sub>1</sub></td><td>0.1467</td></tr><tr><td>β<sub>2</sub></td><td>0.6012</td></tr><tr><td>β<sub>3</sub></td><td>0.1101</td></tr></tbody></table>

The clustering obtained from the BD-SBM reveals marked heterogeneity among authors: it highlights sparsely connected nodes, researchers who tend to collaborate less within all communities, and groups of authors whose collaboration patterns align
closely with the primary arXiv subject categories. Taken together, these results illustrate the ability of the BD-SBM to disentangle latent research communities while accounting for the demographic turnover induced by the arrival and departure of authors over time.

# 5 Conclusion

In this work we have introduced the birth-death stochastic block model (BD-SBM), a dynamic stochastic block model that explicitly accounts for demographic turnover. The BD-SBM combines a continuous-time birth-death process for the population size with a block-structured model for edges observed at discrete times. Individuals belong to one of *K* latent communities, inherit the community of their parent at birth, and retain this membership throughout their lifetime. This construction allows the vertex set to evolve over time while preserving an underlying community structure that governs edge formation, and recovers classical static and dynamic SBMs with fixed vertex set as special cases.

From an inferential viewpoint, we have derived the complete-data likelihood of the BD-SBM and shown how the presence of latent community labels and unobserved community sizes naturally leads to an intractable observed-data likelihood. To overcome this difficulty, we have developed a variational expectation-maximization (VEM) algorithm tailored to the BD-SBM. Our variational family combines a mean-field approximation for individual labels present at time $t_0$ with a structured approximation for the latent community-size process, encoded through the additional variational parameters of the marginal and transition probabilities. Within this framework, the block parameters $(\pi, \beta)$ are estimated via variational updates, while the birth-death parameters $(\lambda, \mu)$, which do not depend on the latent communities, admit closed-form maximum likelihood estimators that can be computed directly from the observed event history.

We have also proposed an adaptation of the Integrated Completed Likelihood (ICL) criterion to the BD-SBM for selecting the number of communities. The resulting penalisation reflects the complexity of the block structure through the community proportions and block connection parameters, while the birth-death parameters, which do not depend on *K*, are excluded from the penalty. Maximizing this ICL criterion over a range of candidate values of *K* provides a practical and principled tool for model selection in dynamic networks with evolving population size. In addition, we have detailed an initialization strategy for the VEM algorithm based on a similarity matrix and *K*-means clustering, which yields informative starting values for the variational distribution without enforcing the birth-death structure at initialization.

The numerical study illustrates the ability of the BD-SBM to recover latent communities and model parameters in synthetic temporal networks under two signal regimes. In a low-signal setting, where within- and between-community connection probabilities are relatively close and small (yielding a sparse network with weakly separated communities), the method still performs well and recovers the main community structure, although with lower clustering accuracy than in the high-signal case. In a high-signal regime, characterized by much stronger within-community connectivity, the BD-SBM
behaves as expected and achieves substantially higher accuracy in both clustering and parameter estimation. We also considered the special case of a pure birth process, for which the algorithm coincides with an exact VEM scheme and the community-size process is Markovian without additional approximations; in this setting, the procedure almost perfectly recovers the true parameters and labels.

We further applied the BD-SBM to a temporal co-authorship network constructed from arXiv mathematics articles over a range of 6 years. Each article has been treated as an undirected temporal interaction among its co-authors adding an edge the pair of its co-authors at the date of the paper’s first submission. The clustering obtained from the BD-SBM reveals marked heterogeneity among authors: it highlights sparsely connected nodes, researchers who tend to collaborate less within all communities, and groups of authors whose collaboration patterns align closely with the primary arXiv subject categories.

Looking ahead, several directions for future research could further enrich the framework developed in this work. A first natural extension is to allow community-specific birth and death parameters, so that demographic dynamics may differ across communities; this extension is discussed in Appendix B, but implementing it and evaluating its performance remains a subsequent step. It would also be interesting to relax the assumption of fixed community memberships and to combine the birth-death mechanism with models that permit community switching, overlapping or mixed-membership structures, as well as the possibility for new communities to emerge over time. Another important avenue is to relax the assumption that the birth-death process is fully observed, and instead consider settings where only discrete-time measurements of the population (e.g. presence/absence or counts at $t_0, \dots, t_T$) are available, so that the demographic process becomes partially latent and must be inferred jointly with the community structure. A further promising direction is to let edges evolve in continuous time according to a Markovian dynamics, for instance via a continuous-time Markov chain on edge states; in such settings, the VEM algorithm developed here would no longer yield closed-form updates and it would be natural to resort to MCMC-based or hybrid variational-sampling schemes to approximate the posterior distribution. Finally, applying the BD SBM to additional real datasets in ecology, epidemiology or more finely resolved academic collaboration networks, in which one can follow, over long time periods, the evolution of emerging or trending research topics, would further assess its practical relevance and may suggest additional refinements of the model.

**Supplementary information.**
<table><caption>Table 5: arXiv Mathematics subject categories and their descriptions.</caption><thead><tr><th>Category</th><th>Description</th></tr></thead><tbody><tr><td>math.AC</td><td>Commutative Algebra.</td></tr><tr><td>math.AG</td><td>Algebraic Geometry.</td></tr><tr><td>math.AP</td><td>Analysis of PDEs.</td></tr><tr><td>math.AT</td><td>Algebraic Topology.</td></tr><tr><td>math.CA</td><td>Classical Analysis and ODEs.</td></tr><tr><td>math.CO</td><td>Combinatorics.</td></tr><tr><td>math.CT</td><td>Category Theory.</td></tr><tr><td>math.CV</td><td>Complex Variables.</td></tr><tr><td>math.DG</td><td>Differential Geometry.</td></tr><tr><td>math.DS</td><td>Dynamical Systems.</td></tr><tr><td>math.FA</td><td>Functional Analysis.</td></tr><tr><td>math.GM</td><td>General Mathematics.</td></tr><tr><td>math.GN</td><td>General Topology.</td></tr><tr><td>math.GR</td><td>Group Theory.</td></tr><tr><td>math.GT</td><td>Geometric Topology.</td></tr><tr><td>math.HO</td><td>History and Overview.</td></tr><tr><td>math.IT</td><td>Information Theory.</td></tr><tr><td>math.KT</td><td>K-Theory and Homology.</td></tr><tr><td>math.LO</td><td>Logic.</td></tr><tr><td>math.MG</td><td>Metric Geometry.</td></tr><tr><td>math.MP</td><td>Mathematical Physics.</td></tr><tr><td>math.NA</td><td>Numerical Analysis.</td></tr><tr><td>math.NT</td><td>Number Theory.</td></tr><tr><td>math.OA</td><td>Operator Algebras.</td></tr><tr><td>math.OC</td><td>Optimization and Control.</td></tr><tr><td>math.PR</td><td>Probability.</td></tr><tr><td>math.QA</td><td>Quantum Algebra.</td></tr><tr><td>math.RA</td><td>Rings and Algebras.</td></tr><tr><td>math.RT</td><td>Representation Theory.</td></tr><tr><td>math.SG</td><td>Symplectic Geometry.</td></tr><tr><td>math.SP</td><td>Spectral Theory.</td></tr><tr><td>math.ST</td><td>Statistics Theory.</td></tr></tbody></table>

# A Proofs

In this appendix we collect the proofs of Propositions 2 and 3. We first recall the Lagrangian associated with the constraints introduced in Section 11.

## A.1 Lagrangian associated with the constraints

The constraints on the variational parameters $(\delta, \gamma, \gamma_{\text{mar}})$ can be written in the form:

$$ (a) \sum_k \delta(i, k) = 1, \text{ for all } i \in V_0; $$
(b) $\sum_{k,n} \gamma(\ell, k, n, n) \gamma_{\text{mar}}(\ell - 1, k, n) = K - 1$, for all $\tau_{\ell} \in \tau$;

(c) $\gamma(\ell, k, n, n + b_{\ell}) + \gamma(\ell, k, n, n) = 1$, for all $\ell, k, n$;

(d) $\sum_k \sum_{n=0}^{N(\tau_{\ell}-1)} \gamma(\ell, k, n, n + b_{\ell}) \gamma_{\text{mar}}(\ell - 1, k, n) = 1$, for all $\tau_{\ell} \in \tau$,

where $b_{\ell} \in \{-1, 1\}$ indicates whether the event at time $\tau_{\ell}$ is a birth or a death.

We associate to these constraints the following Lagrange terms:

$$ 
\left\{ 
\begin{aligned}
\text{Lag(a)} &= \sum_{i} \nu_i (\sum_{k} \delta(i, k) - 1), \\
\text{Lag(b)} &= \sum_{\ell} \sigma_{\ell} (\sum_{k,n} \gamma(\ell, k, n, n) \gamma_{\text{mar}}(\ell - 1, k, n) - (K - 1)), \\
\text{Lag(c)} &= \sum_{\ell, k, n} \eta_{\ell kn} (\gamma(\ell, k, n, n + b_{\ell}) + \gamma(\ell, k, n, n) - 1), \\
\text{Lag(d)} &= \sum_{\ell} \nu_{\ell} (\sum_{k} \sum_{n=0}^{N(\tau_{\ell}-1)} \gamma(\ell, k, n, n + b_{\ell}) \gamma_{\text{mar}}(\ell - 1, k, n) - 1).
\end{aligned}
\right. 
$$

We then consider the augmented objective function,

$$ \mathcal{J}(\delta, \gamma, \nu, \sigma, \eta, \nu) = \mathcal{L}(q) + \text{Lag(a)} + \text{Lag(b)} + \text{Lag(c)} + \text{Lag(d)}, $$

where $\mathcal{L}(q)$ is the ELBO function defined in 10.

## A.2 Update of $\delta(i, k)$ for $i \in V_0$

We maximize $\mathcal{J}(\delta, \gamma, \nu, \sigma, \eta, \nu)$ with respect to $\delta(i, k)$ for $i \in V_0$, under constraint (a). Following the approximation used in [6, 28], we neglect the dependence on $\delta(i, k)$ of certain terms in $\mathcal{L}(q)$. Differentiating the resulting objective with respect to $\delta(i, k)$ gives

$$ \frac{\partial \mathcal{J}}{\partial \delta(i, k)} = \sum_{j \neq i} \sum_{k'} \sum_{\ell \in \Upsilon_{ij}} \delta(j, k') \log \phi(e_{ij}(t_{\ell}), k, k') + \log \beta_k - \log \delta(i, k) - 1 + \nu_i. $$

Setting this derivative equal to zero leads to

$$ \log \delta(i, k) = \sum_{j \neq i} \sum_{k'} \sum_{\ell \in \Upsilon_{ij}} \delta(j, k') \log \phi(e_{ij}(t_{\ell}), k, k') + \log \beta_k - 1 + \nu_i. $$

Hence,

$$ \delta(i, k) = \beta_k e^{\nu_i - 1} \prod_{j \neq i} \prod_{k'} \prod_{\ell \in \Upsilon_{ij}} \phi(e_{ij}(t_{\ell}), k, k')^{\delta(j, k')}, \quad \text{for all } i \in V_0. $$

Finally, the normalisation $\sum_k \delta(i, k) = 1$ for each $i$ determines $e^{\nu_i - 1}$ as the appropriate normalizing constant and yields the expression proposed.
### A.3 Proof of Proposition 2

We now maximize $\mathcal{J}(\delta, \gamma, \nu, \sigma, \eta, \nu)$ with respect to $\gamma(\ell, k, n, n+1)$ and $\gamma(\ell, k, n, n)$ for birth times $\ell \in T_B$. Recall that, if $i_\ell$ is the individual born at time $\tau_\ell$, then

$$ \delta(i_\ell, k) = \sum_{n \ge 1} \gamma(\ell, k, n, n+1) \gamma_{\text{mar}}(\ell-1, k, n), $$

so that

$$ \frac{\partial \delta(i_\ell, k)}{\partial \gamma(\ell, k, n, n+1)} = \gamma_{\text{mar}}(\ell-1, k, n). \quad (17) $$

* Using the above relation and differentiating $\mathcal{J}(\delta, \gamma, \nu, \sigma, \eta, \nu)$ with respect to $\gamma(\ell, k, n, n+1)$, we obtain

$$ \frac{\partial \mathcal{J}}{\partial \gamma(\ell, k, n, n+1)} = \sum_{j \neq i_\ell} \sum_{k'} \sum_{\ell' \in \mathcal{Y}_{i_\ell j}} \gamma_{\text{mar}}(\ell-1, k, n) \delta(j, k') \log \phi(e_{ij}(t_{\ell'}), k, k') \\ + \gamma_{\text{mar}}(\ell-1, k, n) \log n - \gamma_{\text{mar}}(\ell-1, k, n) \log \gamma(\ell, k, n, n+1) - \gamma_{\text{mar}}(\ell-1, k, n) \\ + v_\ell \gamma_{\text{mar}}(\ell-1, k, n) + \eta_{\ell kn} \gamma_{\text{mar}}(\ell-1, k, n). $$

Setting this derivative equal to zero and dividing by $\gamma_{\text{mar}}(\ell-1, k, n)$ leads to

$$ \log \gamma(\ell, k, n, n+1) = \sum_{j \neq i_\ell} \sum_{k'} \sum_{\ell' \in \mathcal{Y}_{i_\ell j}} \delta(j, k') \log \phi(e_{ij}(t_{\ell'}), k, k') + \log n - 1 + v_\ell + \eta_{\ell kn}. $$

Equivalently,

$$ \gamma(\ell, k, n, n+1) = e^{v_\ell + \eta_{\ell kn} - 1} n \prod_{j \neq i_\ell} \prod_{k'} \prod_{\ell' \in \mathcal{Y}_{i_\ell j}} \phi(e_{ij}(t_{\ell'}), k, k')^{\delta(j, k')}, \quad \text{for } \ell \in T_B. $$

* Similarly, differentiating with respect to $\gamma(\ell, k, n, n)$ gives

$$ \frac{\partial \mathcal{J}}{\partial \gamma(\ell, k, n, n)} = -\gamma_{\text{mar}}(\ell-1, k, n) \log \gamma(\ell, k, n, n) - \gamma_{\text{mar}}(\ell-1, k, n) \\ + \sigma_\ell \gamma_{\text{mar}}(\ell-1, k, n) + \eta_{\ell kn} \gamma_{\text{mar}}(\ell-1, k, n), $$

so that

$$ \log \gamma(\ell, k, n, n) = -1 + \sigma_\ell + \eta_{\ell kn}. $$

For notational convenience, define

$$ p(\ell, k) = \prod_{j \neq i_\ell} \prod_{k'} \prod_{\ell' \in \mathcal{Y}_{i_\ell j}} \phi(e_{ij}(t_{\ell'}), k, k')^{\delta(j, k')}. $$
Then we can rewrite

$$ \gamma(\ell, k, n, n+1) = e^{v_{\ell} + \eta_{\ell k n} - 1} n p(\ell, k) \quad \text{and} \quad \gamma(\ell, k, n, n) = e^{\sigma_{\ell} + \eta_{\ell k n} - 1}. $$

Constraint (c) imposes $\gamma(\ell, k, n, n+1) + \gamma(\ell, k, n, n) = 1$, hence

$$ e^{\eta_{\ell k n}} [e^{v_{\ell}-1} n p(\ell, k) + e^{\sigma_{\ell}-1}] = 1, $$

which gives

$$ e^{\eta_{\ell k n}} = \frac{1}{e^{v_{\ell}-1} n p(\ell, k) + e^{\sigma_{\ell}-1}}. $$

Substituting back, we obtain

$$ \gamma(\ell, k, n, n+1) = \begin{cases} \frac{e^{v_{\ell}} n p(\ell, k)}{e^{\sigma_{\ell}} + n e^{v_{\ell}} p(\ell, k)} & \text{if } n \in [0, N_{\ell-1}], \\ 1 & \text{if } n = N_{\ell-1}, \end{cases} $$

and

$$ \gamma(\ell, k, n, n) = \begin{cases} \frac{e^{\sigma_{\ell}}}{e^{\sigma_{\ell}} + n e^{v_{\ell}} p(\ell, k)} & \text{if } n \in [0, N_{\ell-1}], \\ 0 & \text{if } n = N_{\ell-1}. \end{cases} $$

Finally, setting $\rho_{\ell} = e^{\sigma_{\ell} - v_{\ell}} > 0$, we can rewrite

$$ \gamma(\ell, k, n, n+1) = \frac{n p(\ell, k)}{\rho_{\ell} + n p(\ell, k)} \quad \text{and} \quad \gamma(\ell, k, n, n) = \frac{\rho_{\ell}}{\rho_{\ell} + n p(\ell, k)}, \quad 0 \le n < N_{\ell-1}. $$

The value of $\rho_{\ell}$ is determined from constraint (d) (or equivalently from (b)), which requires

$$ \sum_{k,n} \gamma(\ell, k, n, n+1) \gamma_{\text{mar}}(\ell-1, k, n) = 1, $$

that is,

$$ \sum_{k,n} \frac{n p(\ell, k)}{\rho_{\ell} + n p(\ell, k)} \gamma_{\text{mar}}(\ell - 1, k, n) = 1. $$

First, we need to prove the existence of this solution. Define

$$ f_{\ell}(\rho) = \sum_{k,n} \frac{n p(\ell, k)}{\rho + n p(\ell, k)} \gamma_{\text{mar}}(\ell - 1, k, n). $$

For each fixed $(k, n)$ with $n \ge 1$ and $p(\ell, k) > 0$, the map

$$ \rho \mapsto \frac{n p(\ell, k)}{\rho + n p(\ell, k)} \gamma_{\text{mar}}(\ell - 1, k, n) $$
is continuous and strictly decreasing on $(0, \infty)$, with

$$ 
\begin{aligned}
\lim_{\rho \to 0^+} \frac{n p(\ell, k)}{\rho + n p(\ell, k)} \gamma_{\text{mar}}(\ell - 1, k, n) &= \gamma_{\text{mar}}(\ell - 1, k, n), \\
\lim_{\rho \to +\infty} \frac{n p(\ell, k)}{\rho + n p(\ell, k)} \gamma_{\text{mar}}(\ell - 1, k, n) &= 0.
\end{aligned} 
\quad (18) $$

Hence $f_\ell$ is continuous and strictly decreasing on $(0, \infty)$, and

$$ \lim_{\rho \to 0^+} f_\ell(\rho) = \sum_k \sum_{n \ge 1} \gamma_{\text{mar}}(\ell - 1, k, n), $$

is the expected number of non-empty communities at time $\tau_{\ell-1}$. If $N(\tau_{\ell-1}) \ge 1$ there is at least one individual alive, so at least one community is non-empty with positive probability, and then the number of non-empty communities is almost surely greater or equal than 1. This implies in particular that $f_\ell(0^+) \ge 1$, while

$$ \lim_{\rho \to +\infty} f_\ell(\rho) = 0. $$

By continuity and monotonicity, there exists a unique $\rho_\ell > 0$ such that $f_\ell(\rho_\ell) = 1$. Therefore the equation

$$ \sum_{k,n} \frac{n p(\ell, k)}{\rho_\ell + n p(\ell, k)} \gamma_{\text{mar}}(\ell - 1, k, n) = 1, $$

admits a unique solution $\rho_\ell > 0$. Since this equation has no closed-form solution in general, $\rho_\ell$ must be found numerically. This concludes the proof of Proposition 2.

### A.3.1 Approximation of $\gamma$ for death times

We show here that the parametrization adopted in (16) admits a unique solution. For each fixed $(\ell, k)$, define

$$ f_{\ell k}(\rho) = \sum_{n=0}^{N(\tau_{\ell-1})-1} \frac{\rho n}{1 + \rho n} \gamma_{\text{mar}}(\ell - 1, k, n + 1), \quad \rho \ge 0. $$

For $n \ge 1$, the map

$$ \rho \mapsto \frac{\rho n}{1 + \rho n} $$

is continuous and strictly increasing on $[0, \infty)$, with limits 0 as $\rho \to 0^+$ and 1 as $\rho \to \infty$, and derivative $n/(1 + \rho n)^2 > 0$.
Then, $f_{\ell k}$ is continuous and strictly increasing on $[0, \infty)$, and

$$ \lim_{\rho \to 0^+} f_{\ell k}(\rho) = 0, \quad \lim_{\rho \to \infty} f_{\ell k}(\rho) = \sum_{n=1}^{N(\tau_{\ell-1})} \hat{\gamma}_{\text{mar}}(\ell-1, k, n) = 1 - \hat{\gamma}_{\text{mar}}(\ell-1, k, 0). $$

In our model, a death in community $k$ cannot occur when the community is empty, so it is natural to have

$$ \delta(i_{\ell}, k) \le 1 - \gamma_{\text{mar}}(\ell - 1, k, 0). $$

By continuity and strict monotonicity of $f_{\ell k}$, there then exists a unique $\rho_{\ell k} \ge 0$ such that

$$ f_{\ell k}(\rho_{\ell k}) = \delta(i_{\ell}, k), $$

which justifies the result.

### A.4 Proof of Proposition 3

For the proofs of the expressions obtained for the model parameters, we maximize the ELBO $\mathcal{L}(q)$ with respect to $\pi_{k_1 k_2}$ and $\beta_k$, under the constraint $\sum_{k=1}^K \beta_k = 1$.

Then we maximize $\mathcal{L}(q)$ in function of $\pi$

$$ \frac{\partial \mathcal{L}(q)}{\partial \pi_{k_1 k_2}} = \sum_{i \neq j} \sum_{k_1, k_2} \delta(i, k_1) \delta(j, k_2) \left[ \left( \sum_{\ell \in \Upsilon_{ij}} e_{ij}(t_{\ell}) \right) \log \pi_{k_1 k_2} + \left( |\Upsilon_{ij}| - \sum_{\ell \in \Upsilon_{ij}} e_{ij}(t_{\ell}) \right) \log(1 - \pi_{k_1 k_2}) \right], $$

and setting this derivative equal to zero we obtain the result.

To maximize $\mathcal{L}(q)$ with respect to $\beta$ under the constraint $\sum_k \beta_k = 1$, we introduce a Lagrange multiplier $\xi$ and consider

$$ \mathcal{J}(\beta, \xi) = \mathcal{L}(q) + \xi \left( \sum_{k=1}^{K} \beta_k - 1 \right) = \sum_{i \in V_0} \sum_{k=1}^{K} \delta(i, k) \log \beta_k + \xi \left( \sum_{k=1}^{K} \beta_k - 1 \right). $$

We differentiate with respect to $\beta_k$:

$$ \frac{\partial \mathcal{J}}{\partial \beta_k} = \frac{1}{\beta_k} \sum_{i \in V_0} \delta(i, k) + \xi. $$

Setting this equal to zero yields,

$$ \beta_k = -\frac{1}{\xi} \sum_{i \in V_0} \delta(i, k). $$

Using the constraint $\sum_k \beta_k = 1$, we obtain the result.
# B Community-specific birth and death rates in the BD-SBM

In the main part of this work, we assumed shared demographic rates across communities, that is $\lambda_k \equiv \lambda$, $\mu_k \equiv \mu$ for all $k \in \{1, \dots, K\}$. This appendix relaxes this assumption and allows each community to have its own birth and death rates, leading to the following vector parameters

$$ \lambda = (\lambda_1, \dots, \lambda_K), \quad \mu = (\mu_1, \dots, \mu_K), \quad \text{with } \lambda_k \ge 0, \mu_k \ge 0, $$

where $\lambda_k$ and $\mu_k$ correspond respectively to the birth and death parameters of community $k$. For simplicity, we denote this vector parameters by $\lambda$ and $\mu$ throughout this appendix. All other components of the BD-SBM (SBM edge model, variational family, and optimization scheme) remain unchanged; only the birth-death contribution to the likelihood/ELBO and the corresponding VM-step updates are modified.

## B.1 Birth-death likelihood with ($\lambda_k, \mu_k$)

Let $Z_i$ be the (latent) community labels of individual $i$. Conditionally on the labels, each alive individual of class $k$ gives birth at rate $\lambda_k$ and dies at rate $\mu_k$. Let $N_k(t)$ be the size of community $k$ at time $t$, and denote the event times by $\tau_1 < \dots < \tau_M$, with inter-event durations $\Delta_\ell = \tau_\ell - \tau_{\ell-1}$ (and $\tau_0 = t_0$). Recall that $T_B$ (resp. $T_D$) denotes the set of birth (resp. death) event indices. As in the main text, let $i_\ell$ be the newborn at $\tau_\ell$ if $\tau_\ell \in T_B$, and the dying individual at $\tau_\ell$ if $\tau_\ell \in T_D$.

The total event rate at $\tau_{\ell-1}$ is

$$ R(\tau_{\ell-1}) = \sum_{k=1}^{K} (\lambda_k + \mu_k) N_k(\tau_{\ell-1}). $$

Under this parametrization, the second term in the likelihood (1) can be written as

$$ \begin{aligned} \mathbb{P}_{\theta}(Z, \tau, b(\tau) | Z^{t_0}) &= \prod_{\ell=1}^{M} \left[ e^{-R(\tau_{\ell-1})\Delta_{\ell}} \prod_{k=1}^{K} \left( \lambda_k \sum_{i \in V(\tau_{\ell-1})} Z_{ik} \right)^{Z_{j_{\ell}k} 1\{b_{\ell}=1\}} \mu_k^{Z_{j_{\ell}k} 1\{b_{\ell}=-1\}} \right] \\ &= \prod_{\ell=1}^{M} \left[ e^{-R(\tau_{\ell-1})\Delta_{\ell}} \prod_{k=1}^{K} \prod_{n=1}^{N(\tau_{\ell-1})} \left( \lambda_k n \right)^{L_{k,n+1}^{\ell} L_{k,n}^{\ell-1} L_{k,n-1}^{\ell} L_{k,n}^{\ell-1}} \mu_k^{L_{k,n-1}^{\ell}} \right]. \end{aligned} $$

Taking logs yields the log-likelihood component

$$ \log \mathbb{P}_{\theta}(Z, \tau, b(\tau) | Z^{t_0}) = \sum_{\ell=1}^{M+1} \left[ -R(\tau_{\ell-1})\Delta_{\ell} + \sum_{k=1}^{K} \sum_{n=1}^{N(\tau_{\ell-1})} L_{k,n+1}^{\ell} L_{k,n}^{\ell-1} \log(\lambda_k n) + L_{k,n-1}^{\ell} L_{k,n}^{\ell-1} \log \mu_k \right]. $$
## B.2 Contribution to the ELBO

Recall that $\gamma_{\text{mar}}(\ell, k, n)$ denotes the variational marginal pmf of $N_k(\tau_\ell)$, and $\gamma(\ell, k, n, n')$ the variational transition probability of passing in community $k$ from $n$ to $n'$ at a time event $\tau_\ell \in \tau$ (see Section 3.3.1). Under the variational distribution $q$, the expected total rate satisfies

$$ \mathbb{E}_q [R(\tau_{\ell-1})] = \sum_{k=1}^{K} (\lambda_k + \mu_k) \sum_{n=1}^{N(\tau_{\ell-1})} n \gamma_{\text{mar}}(\ell-1, k, n). $$

Therefore, the BD contribution to the ELBO, that is, the sum of the second, third, fourth and fifth terms in (10) becomes,

$$ \begin{aligned} & \sum_{\tau_\ell \in \tau} \sum_{k=1}^{K} \sum_{n=1}^{N(\tau_{\ell-1})} \left[ - (\lambda_k + \mu_k) \Delta_\ell n \gamma_{\text{mar}}(\ell-1, k, n) \right. \\ & \qquad \left. + \gamma(\ell, k, n, n+1) \gamma_{\text{mar}}(\ell-1, k, n) \log(\lambda_k n) + \gamma(\ell, k, n, n-1) \gamma_{\text{mar}}(\ell-1, k, n) \log(\mu_k) \right] \\ & = - \sum_{\tau_\ell \in \tau} \sum_{k=1}^{K} \sum_{n=1}^{N(\tau_{\ell-1})} (\lambda_k + \mu_k) \Delta_\ell n \gamma_{\text{mar}}(\ell-1, k, n) + \sum_{\tau_\ell \in T_D} \sum_{k=1}^{K} \delta(i_\ell, k) \log \mu_k \\ & \qquad + \sum_{\tau_\ell \in T_B} \sum_{k=1}^{K} \sum_{n=1}^{N(\tau_{\ell-1})} \gamma(\ell, k, n, n+1) \gamma_{\text{mar}}(\ell-1, k, n) \log(\lambda_k n). \end{aligned} $$

Then, the ELBO becomes,

$$ \begin{aligned} \mathcal{L}(q) = & \sum_{i<j} \sum_{k_1, k_2} \delta(i, k_1)\delta(i, k_2) \sum_{\ell \in \Upsilon_{ij}} \log \phi(e_{ij}(t_\ell), k_1, k_2) - \sum_{\tau_\ell \in \tau} \sum_{k=1}^{K} \sum_{n=1}^{N(\tau_{\ell-1})} (\lambda_k + \mu_k) \Delta_\ell n \gamma_{\text{mar}}(\ell-1, k, n) \\ & + \sum_{\tau_\ell \in T_D} \sum_{k=1}^{K} \delta(i_\ell, k) \log \mu_k + \sum_{\tau_\ell \in T_B} \sum_{k=1}^{K} \sum_{n=1}^{N(\tau_{\ell-1})} \gamma(\ell, k, n, n+1) \gamma_{\text{mar}}(\ell-1, k, n) \log(\lambda_k n) \\ & + \sum_{i \in V_0} \sum_{k=1}^{K} \delta(i, k) \log \beta_k - \sum_{i \in V_0} \sum_{k=1}^{K} \delta(i, k) \log \delta(i, k) \\ & + \sum_{\tau_\ell \in T_B} \sum_{k=1}^{K} \sum_{n=1}^{N(\tau_{\ell-1})} \left[ \gamma(\ell, k, n, n) \gamma_{\text{mar}}(\ell-1, k, n) \log \gamma(\ell, k, n, n) \right. \\ & \qquad \left. + \gamma(\ell, k, n, n+1) \gamma_{\text{mar}}(\ell-1, k, n) \log \gamma(\ell, k, n, n+1) \right]. \end{aligned} $$
## B.3 VM-step updates for $\lambda_k$ and $\mu_k$

We now summarize the coordinate updates used in the VEM algorithm.

In analogy to equation (12), we choose to update $\delta(i, k)$ for $i \in V_0$ by solving the fixed-point equation

$$ \hat{\delta}(i, k) \propto \beta_k \mu_k^{1\{\tau_i^d \in T_D\}} \prod_{j \neq i} \prod_{k'} \prod_{\ell' \in \Upsilon_{ij}} [\phi(e_{ij}(t_{\ell'}), k, k')]^{\hat{\delta}(j, k')} $$

where $\phi(e_{ij}(t_{\ell'}), k, k')$ and $\Upsilon_{ij}$ are given in Proposition 1.

**Proposition 4** For $\ell \in T_B$, we refer to $i_\ell$ as the newborn at time $\tau_\ell$, and we denote by $\hat{\gamma}(\ell, k, n, n)$ the quantity that maximizes $\mathcal{L}(q)$ with respect to $\gamma(\ell, k, n, n)$ when $\theta = (\lambda, \mu, \pi, \beta)$ is fixed. Then there exists a unique $\rho_\ell > 0$ such that

$$ \hat{\gamma}(\ell, k, n, n) = \frac{\rho_{\ell}}{\rho_{\ell} + \mu_k^{1\{\tau_{i_{\ell}}^d \in T_D\}}} \lambda_k n \prod_{j \neq i_{\ell}} \prod_{k'} \prod_{\ell' \in \Upsilon_{i_{\ell}j}} [\phi(e_{ij}(t_{\ell'}), k, k')]^{\hat{\delta}(j, k')}, \quad \text{for } n \in [0, N_{\ell-1}], $$

where the parameter $\rho_\ell$ is determined by the constraint (b) in (11).

*Proof* The result follows by maximizing $\mathcal{L}(q)$ with respect to $\gamma(\ell, k, n, n)$, using Lagrange multipliers for the constraints in (11), in analogy to the derivation provided in Appendix A. $\square$

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

*Remark 2* In order to compute $\hat{\gamma}(\ell, k, n, n)$ in practice, one must first determine $\rho_\ell$. Since the defining equation for $\rho_\ell$ has no closed-form solution in general, $\rho_\ell$ is obtained numerically.

**Proposition 5** Let $q \in \mathcal{Q}$ be fixed. Then the maximizers of $\mathcal{L}(q)$ with respect to $\lambda_k$ and $\mu_k$ are given by

$$ \hat{\lambda}_k = \frac{\sum_{\ell \in T_B} \delta(i_\ell, k)}{\sum_{\tau_\ell \in \tau} \sum_n \Delta_\ell n \gamma_{\text{mar}}(\ell - 1, k, n)}, \quad \hat{\mu}_k = \frac{\sum_{\ell \in T_D} \delta(i_\ell, k)}{\sum_{\tau_\ell \in \tau} \sum_n \Delta_\ell n \gamma_{\text{mar}}(\ell - 1, k, n)}. $$

The updates of $\pi$ and $\beta$ are unchanged and are given in Proposition 3.

*Proof* The result is straightforward by solving the maximization of the ELBO function $\mathcal{L}(q)$ with respect to the parameters. $\square$

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

## References

[1] Erdős, P., Rényi, A., et al.: On the evolution of random graphs. Publ. Math. Inst. Hung. Acad. Sci 5(1), 17–60 (1960)

[2] Watts, D.J., Strogatz, S.H.: Collective dynamics of ‘small-world’networks. Nature 393(6684), 440–442 (1998)
[3] Holland, P.W., Laskey, K.B., Leinhardt, S.: Stochastic blockmodels: First steps. Social Networks 5(2), 109–137 (1983) https://doi.org/10.1016/0378-8733(83) 90021-7

[4] Snijders, T.A.B., Nowicki, K.: Estimation and prediction for stochastic blockmodels for graphs with latent block structure. Journal of Classification 14(1), 75–100 (1997) https://doi.org/10.1007/s003579900031

[5] Aicher, C., Jacobs, A.Z., Clauset, A.: Learning latent block structure in weighted networks. Journal of Complex Networks 3(2), 221–248 (2014) https://doi.org/10. 1093/comnet/cnu026

[6] Matias, C., Miele, V.: Statistical clustering of temporal networks through a dynamic stochastic block model. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 79(4), 1119–1141 (2017) https://doi.org/10.1111/rssb. 12200

[7] Tang, X., Yang, C.C.: Dynamic community detection with temporal dirichlet process. In: 2011 IEEE Third International Conference on Privacy, Security, Risk and Trust and 2011 IEEE Third International Conference on Social Computing, pp. 603–608 (2011). IEEE

[8] Xu, K.S., Hero, A.O.: Dynamic stochastic blockmodels: Statistical models for time-evolving networks. CoRR abs/1403.1358 (2014)

[9] Jiang, B., Li, J., Yao, Q.: Autoregressive networks. Journal of Machine Learning Research 24(227), 1–69 (2023)

[10] Ludkin, M., Eckley, I., Neal, P.: Dynamic stochastic block models: parameter estimation and detection of changes in community structure. Statistics and Computing 28, 1201–1213 (2018)

[11] Corneli, M.: Dynamic stochastic block models, clustering and segmentations in dynamic graphs. PhD thesis, Université Paris 1-Panthéon Sorbonne (2017)

[12] Peel, L., Clauset, A.: Detecting change points in the large-scale structure of evolving networks. Proceedings of the AAAI Conference on Artificial Intelligence 29(1), 2914–2920 (2015)

[13] Wang, T., Aggarwal, C.C., Wang, W.: A survey of algorithms for change detection in dynamic graphs. Knowledge and Information Systems 51(3), 635–669 (2017) https://doi.org/10.1007/s10115-016-0987-z

[14] Brault, V., Mariadassou, M.: Co-clustering through latent bloc model: a review. Journal de la Société Française de Statistique 156(3), 120–139 (2015)

[15] Yang, T., Chi, Y., Zhu, S., Gong, Y., Jin, R.: Detecting communities and their
evolutions in dynamic social networks—a bayesian approach. In: Proceedings of the 2011 ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pp. 990–998 (2011)

[16] Xu, K.: Stochastic Block Transition Models for Dynamic Networks. In: Lebanon, G., Vishwanathan, S.V.N. (eds.) Proceedings of the Eighteenth International Conference on Artificial Intelligence and Statistics. Proceedings of Machine Learning Research, vol. 38, pp. 1079–1087. PMLR, San Diego, California, USA (2015). https://proceedings.mlr.press/v38/xu15.html

[17] Pensky, M., Zhang, T.: Dynamic network models and graphon estimation. Electronic Journal of Statistics 13(1), 678–709 (2019)

[18] Brânzei, S., Kumar, N., Ranade, G.: Phase transitions of diversity in stochastic block model dynamics. In: 2023 59th Annual Allerton Conference on Communication, Control, and Computing (Allerton), pp. 1–8 (2023). IEEE

[19] Daudin, J.-J., Picard, F., Robin, S.: A mixture model for random graphs. Statistics and computing 18(2), 173–183 (2008)

[20] Jordan, M.I., Ghahramani, Z., Jaakkola, T.S., Saul, L.K.: An introduction to variational methods for graphical models. Machine learning 37(2), 183–233 (1999)

[21] Matias, C., Robin, S.: Modeling heterogeneity in random graphs through latent space models: a selective review. ESAIM: Proceedings and Surveys 47, 55–74 (2014)

[22] DuBois, C., Butts, C., Smyth, P.: Stochastic blockmodeling of relational event dynamics. In: Carvalho, C.M., Ravikumar, P. (eds.) Proceedings of the Sixteenth International Conference on Artificial Intelligence and Statistics. Proceedings of Machine Learning Research, vol. 31, pp. 238–246. PMLR, Scottsdale, Arizona, USA (2013). https://proceedings.mlr.press/v31/dubois13a.html

[23] Reynolds, J.F.: On estimating the parameters of a birth-death process. Australian Journal of Statistics 15(1), 35–43 (1973)

[24] Keiding, N.: Maximum likelihood estimation in the birth-and-death process. The Annals of Statistics 3(2), 363–372 (1975)

[25] Bishop, C.M.: Pattern recognition and machine learning. Springer google schola 2, 1122–1128 (2006)

[26] Lin, Z., Wang, Y., Hong, Y.: The computing of the poisson multinomial distribution and applications in ecological inference and machine learning. Computational Statistics 38(4), 1851–1877 (2023)
[27] Tang, W., Tang, F.: The poisson binomial distribution—old & new. Statistical Science **38**(1), 108–119 (2023)

[28] Agarwal, A., Lee, K.H., Xue, L.: Clustering time-evolving networks using temporal exponential-family random graph models with conditional dyadic independence and dynamic latent blocks. Journal of Computational and Graphical Statistics, 1–13 (2025)

[29] Biernacki, C., Celeux, G., Govaert, G.: Assessing a mixture model for clustering with the integrated completed likelihood. IEEE Transactions on Pattern Analysis & Machine Intelligence **22**(07), 719–725 (2000)

[30] Frisch, G., Leger, J.-B., Grandvalet, Y.: Sparsebm: A python module for handling sparse graphs with block models (2021)