# Spinal study of a population model for colonial species with interactions and environmental noise

Sylvain Billiard, Charles Medous, Charline Smadi

► To cite this version:

Sylvain Billiard, Charles Medous, Charline Smadi. Spinal study of a population model for colonial species with interactions and environmental noise. 2025. hal-04921079

HAL Id: hal-04921079

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

Preprint submitted on 30 Jan 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
SPINAL STUDY OF A POPULATION MODEL FOR COLONIAL SPECIES
WITH INTERACTIONS AND ENVIRONMENTAL NOISE

SYLVAIN BILLIARD, CHARLES MEDOUS, AND CHARLINE SMADI

ABSTRACT. We introduce and study a stochastic model for the dynamics of colonial species, which reproduce through fission or fragmentation. The fission rate depends on the relative sizes of colonies in the population, and the growth rate of colonies is influenced by intrinsic and environmental stochasticities. Our setting thus captures the effect of an external noise, correlating the trait dynamics of all colonies. In particular, we study the effect of the strength of this correlation on the distribution of resources between colonies. We then extend this model to a large class of structured branching processes with interactions in which the particle type evolves according to a diffusion. The branching and death rates are general functions of the whole population. In this framework, we derive a ψ-spine construction and a Many-to-One formula, extending previous works on interacting branching processes. Using this spinal construction, we also propose an alternative simulation method and illustrate its efficiency on the colonial population model. The extended framework we propose can model various ecological systems with interactions, and individual and environmental noises.

KEY WORDS AND PHRASES: fusion-fission population dynamics, interacting populations, dif-
fusion process, spinal constructions, Many-to-One formula, simulation method, environmental
noise, ecology

MSC 2000 subject classifications: 92D25, 60J80, 60J85, 60h10, 60J60, 60K37

# 1. INTRODUCTION

The life cycle of colonial animal species often includes a mode of reproduction through fission or fragmentation after a colony has grown to a sufficiently large size (e.g. social spiders [41, 4, 25], ants [28, 48], corals [24, 40] or colonial tunicates [12]). Many plants (e.g. bryophytes [29]) or microbial species [50] show similar life cycles patterns. Describing the dynamics of such species needs to consider the dynamics of two nested levels: the size of each colony which varies because of growth (ontogeny) and fission, and the number of colonies which is governed by colony birth (fission) or death (local extinction).

Few models have considered the demographic particularities of colonial species, and most often for addressing specific questions. Some studies aimed at evaluating the effect of colonies death and fission on the extinction of populations, for instance with environmental or human perturbations, or explaining the observed structuration in size (or age) of colonies found in natural populations [40, 38, 47, 39]. Such models also generally lack widely shared processes in natural populations such as the heterogeneity between colonies that might impact their death or fission rates, the stochasticity of colonial growth, fission and death, and the effect of colony size on competition. Our goal in this paper is to provide in a single framework a general stochastic model describing the dynamics of a colonial species including all these mechanisms.

In our model, we consider a colony as an ‘individual’, itself composed of sub-units such as ‘individual’ cells, branches, ants, spiders or even birds. The trait we consider is associated to the colony and not to its sub-units. This trait can represent the size or volume of the colony, that is related to the number of sub-units, e.g. the number of workers or reproductive queens in an ant colony. The dynamics of the population of colonies is governed by three mechanisms. First, we suppose that the size of a colony can vary during its lifetime because of intrinsic factors (stochastic birth and death of sub-units) or extrinsic factors (environmental variability) [12, 16, 19, 41]. This is well represented by the evolution of the traits following a stochastic differential equation (SDE) with intrinsic and environmental components of the
fluctuations. Second, the establishment of new colonies follows fission events. Such fissions can be associated to reproductive events such as in some ant species where queens can found new nests accompanied by a variable number of workers coming from the ancestral colonies [28]. Fissions can also be due to fortuitous random events such as when corals are fragmented or spider webs are split in daughter webs [41, 25, 4]. In any case, the newly established colonies can randomly vary in their initial size, a phenomenon captured by our model as when a birth occurs, the size (or trait) of the colony is split into two random size daughter colonies. Third, we consider competitive interactions between colonies as affecting the fission rate. The fission rate is supposed to be proportional to the relative size of a colony compared to the sum of all other colonies' size. This can represent colonies competing for the resource of their environment as foraging is directly related to the size of the colony, for example the number of foraging ants, the volume of water filtered by branches of coral, the surface of spider webs [53, 51, 28, 55, 19, 5]. To sum up, our model allows to study the population dynamics of colonial species accounting for within and between colonies heterogeneity, in a time varying environment, with competition. We will in particular study the relative roles of the different sources of variability on the maintenance and composition of the population.

More specifically, we consider a structured birth and death process of a population of colonies each characterized by a trait $x$, e.g. their size. We suppose that competition between colonies affects their fission rate such as it is proportional to their relative size $x/R_t$ where $R_t$ is the sum of the sizes of all colonies. At a fission event of a colony of size $x$, resources are shared by the two daughter colonies. One colony gets a fraction $\Theta$ of the size and the other a fraction $(1-\Theta)$ with $\Theta$ a random variable on $(0,1)$. Between two branching events, the colony sizes evolve according to a diffusion; if at time $t$, $N_t$ denotes the number of colonies in the whole population, $u$ the label of a single colony and $X_t^u$ its size,

$$dX_t^u := aX_t^u dt + \frac{\sigma}{\sqrt{1+\delta^2}} X_t^u (dB_t^u + \delta dW_t) \quad (1.1)$$

where $a > 0$ is the growth rate of the colony, $\sigma > 0$ is the intensity of the total noise, $(B_t^u, u \in \mathbb{G}(t))$ denotes a $N_t$-dimensional intrinsic Wiener process and $W_t$ an independent extrinsic Wiener process. The parameter $\delta > 0$ quantifies the relative intensity of the extrinsic and intrinsic noises. The processes $B^v$ for $v$ describing the set of colonies alive are independent and all independent of $W$. The latter process corresponds to the environment and is common to all colonies alive in the same environment. Note that the value of the correlation factor $\delta$ does not change the law of a colony trait, but the correlation between the dynamics of traits of the different colonies of the population. Depending on $\delta > 0$ the noise due to the common environment may be more or less important compared to the independent intrinsic noises in the colonies.

We establish results about the long term behavior of those interacting colonies, using spine techniques. Spine techniques have been developed to study Branching diffusions since the work of Chauvin and Rouault [18] who established limit theorems for the law of the position of the right-utmost particle. Such methods are based on a Girsanov change of probability measure associated with a distinguished individual called the spine. Branching Brownian motion with constant division rate is related to the solution of Kolmogorov-Petrovskii-Piskounov (KPP) diffusion equations [44] and spine techniques have been used to establish the existence, asymptotics, and uniqueness of traveling wave solutions of the KPP equation [34, 36]. We also refer to Hardy and Harris [31] for an intuitive, spine-based proof of a large deviations result in non-structured branching Brownian motion. This framework has been extended to structured branching diffusion, where the trait follows an Ornstein-Uhlenbeck process, with results that include the long-term behavior of the system, such as the speed of colonization in space and type dimensions and population growth rates within this asymptotic shape [30]. The general case of structured branching diffusion has been studied by Hardy and Harris [32], who used this spinal decomposition to prove the $L^p$-convergence of some key martingales, which was later used to establish strong laws of large numbers [27].

We refer to Bansaye and coauthors [9, 7] for a spine-based asymptotic study of the behavior of parasite infection in cells, which follows a Feller diffusion with trait-dependent multiplicative
jumps. Later, Cloez [20] studied a structured branching process where the individual traits follow a general Markov process and determined the asymptotic behavior of the particle system with spine methods. In a similar framework, we also refer to Marguet [42], who studied a sampled individual in the asymptotic population.

The martingale change of measure associated with the spinal construction involves either an eigenfunction $\psi$ of some operator $\mathcal{G}$ related to the branching mechanisms [13, 26, 20], or a Feynman-Kac path-integral formula [33, 42, 6, 46]. In the case of Galton-Watson trees and non-structured branching processes, the constant function $\psi \equiv 1$ is an eigenfunction of the operator $\mathcal{G}$. For structured branching diffusion, Cloez [20] exhibited eigenfunctions of $\mathcal{G}$ for various branching mechanisms. In the case of branching processes with interactions, the interplay between the individuals is more complex, and it is arduous to find an explicit eigenfunction. Bansaye [6] established that, for a non-structured branching process with size-dependent interactions, the function $\psi = 1/N_t$ where $N_t$ denotes the number of individuals in the population is always an eigenfunction of $\mathcal{G}$. However, incorporating a trait dependency into the branching mechanisms compromises this property, resulting in the absence of general eigenfunctions [10]. Our motivating model of colonial species does not fit within the framework of the aforementioned studies, which must be extended to encompass a broader class of structured branching processes with interactions.

The article is organized as in the following. We first introduce our stochastic model of a population of a colonial species. We then establish a spinal construction of the latter, that enables us to study long-term properties of the population, thanks to the use of Many-to-One formulas and the study of the (simpler) spinal process. In particular we show that if the intrinsic variability is large while the environmental variability is of a smaller order, most of the resources are held by a unique colony. Otherwise, the resources are more evenly shared among colonies. We then use the spine process to propose an alternative simulation method of this colonial dynamics, that we compare to the classical method of birth and death process simulation. Finally we extend the generality of this method by derivating a spine construction with a general $\psi$ function- that will be referred to as a $\psi$-spine construction- in the context of branching diffusion with interactions. In this extended setting, we consider a wide class of continuous-time structured branching diffusions with general interactions where each individual (instead of 'colonies' in the extended model to be as general as possible) is characterized by its trait, taking values in a subset of $\mathbb{R}^d$. The lifespan of each individual is exponentially distributed with a time-inhomogeneous rate that depends on the traits of all individuals. Upon an individual's death, a random number of children are generated, each inheriting random traits at birth. Between these branching events, the evolution of the traits of all individuals in the population follows a $N_t$-dimensional SDE where $N_t$ is the population size. This SDE depends on the traits of all the individuals in the population, extending previous works [6, 46] where trait dynamics were deterministic, to include many physical models based on Brownian motion, see e.g. [22] and the references therein. In this general framework, Many-to-One formulas linking the behavior of the original and the spinal processes thus appear as a promising tool to model and study various ecological systems that include interactions, as well as individual and environmental noises.

## 2. A COLONIAL POPULATION CORRELATED THROUGH A COMMON ENVIRONMENT.

### 2.1. Description of the model.
We consider a population $\nu_t$ at time $t$, composed of $N_t \ge 0$ colonies. Each colony is associated to a unique label $u \in \mathcal{U}$, where the set $\mathcal{U}$ is defined according to the Ulam-Harris-Neveu notation:

$$ \mathcal{U} := \{\emptyset\} \cup \bigcup_{k \ge 0} (\mathbb{N}^*)^{k+1}. $$

The set of labels of the colonies present in the population at time $t$ is denoted $\mathcal{G}(t) \subset \mathcal{U}$. We denote $X_t^u \in \mathbb{R}_+^*$ the trait of a colony labeled by $u \in \mathcal{G}(t)$. $X_t^u$ can represent the mass or size
of colony $u$, or its total quantity of resources. Using these notations, we can write

$$ \nu_t := \sum_{u \in \mathbb{G}(t)} \delta_{X_t^u} \in \mathbb{U} $$

where $\mathbb{U}$ denotes the set composed of all finite point measures on $\mathbb{R}_+^*$ defined as

$$ \mathbb{U} := \left\{ \sum_{i=1}^{N} \delta_{x^i}, N \in \mathbb{N}, (x^i, 1 \le i \le N) \in (\mathbb{R}_+^*)^N \right\}. $$

We consider a dyadic birth and death process such that a colony of resource $x \in \mathbb{R}_+^*$ dies at rate $\mu$ and splits at rate $\lambda x/R_t$ where $R_t := \langle \nu_t, I_d \rangle$ denotes the total amount of resources in the population and $\lambda$ is a positive parameter. When a colony of size $x$ splits, the newly founded colony inherits $\Theta x$ resources from the parent colony, while the latter keeps $(1 - \Theta)x$ resources, with $\Theta$ a random variable on $[0, 1]$ of probability measure $q$. Between two branching events at times $T_i$ and $T_{i+1}$, the resources evolve according to a SDE (1.1).

One of the interesting features of the model we have just introduced is that we take into account the competition for space or resources in the reproduction rate. This is technically more difficult to handle than a model taking competition into account in the death term, (which corresponds to the most usual choice, the prototypical example being the competitive Lotka-Volterra model [2, Chapter 7.8]) but is much more realistic for many species (for instance the ones in which larvae compete for resources and have little chance of reaching maturity). A second interesting feature of the model is that it allows taking into account any law of resource sharing between two offspring colonies at a fission event. Recent works [43], considering structured branching processes without interaction, have shown that it has a large impact on the long time behavior of the population. We will present some results in the same direction. Finally, we consider noise arising from two distinct sources: independent internal sources within individual colonies and a shared environmental source affecting all colonies. The question of interaction and relative effects of intrinsic and environmental stochasticities is a major question in ecology [37, 15, 54]. By letting vary the parameter $\delta$ we can modify the level of correlation of the noises impacting the dynamics of the colonies, and study the effect of this correlation.

## 2.2. Spinal Construction.
In this section, we introduce the spine process $(\bar{\nu}_t, t \ge 0)$ associated with the population process and give the relationship between the law of this spine process (that does not describe any real population) and the law of the population process. In the spine process, we will also refer to the elements in the population as colonies described by a trait in $\mathbb{R}_+^*$. The associated spinal process $((E_t, \hat{\nu}_t), t \ge 0)$ distinguishes a particular colony -the spine- of label $E_t$ from the set of labels $\hat{\mathbb{G}}(t)$ in the spinal population $\hat{\nu}_t$. In order to have a simple distinctive feature between both processes, we will denote by $\hat{\nu}_t$ the marginal spinal population, and put a hat on each value that refers to the spinal process. At all time, we will denote by $\hat{Y}_t$ the trait of the spine $E_t$ and by $\hat{X}_t^u$ the trait of the colony of label $u$ belonging to $\hat{\mathbb{G}}(t) \setminus \{E_t\}$ in the spinal population. The SDE describing the evolution of the traits between two branching events involves an additional drift term compared to the resources dynamics in the biological process:

$$ d\hat{X}_t^u := \left(a + \sigma^2 \frac{\delta^2}{1 + \delta^2}\right) \hat{X}_t^u dt + \frac{\sigma}{\sqrt{1 + \delta^2}} \hat{X}_t^u \left(d\hat{B}_t^u + \delta d\hat{W}_t\right), \quad \text{for all } u \in \hat{\mathbb{G}}(t)\setminus\{E_t\}, \quad (2.1) $$

and

$$ d\hat{Y}_t := (a + \sigma^2) \hat{Y}_t dt + \frac{\sigma}{\sqrt{1 + \delta^2}} \hat{Y}_t \left( d\hat{B}_t^* + \delta d\hat{W}_t \right) \quad (2.2) $$

where $\hat{B}^*, \hat{B}^u, \hat{W}$ are independent Wiener processes. The initial spine is chosen randomly in the initial population $\hat{\nu}_0$ among the initial labels $\hat{\mathbb{G}}(0) := \{1, \dots, \langle \hat{\nu}_0, 1 \rangle\}$, with a probability proportional to its trait. That is, for all $j \in \hat{\mathbb{G}}(0)$,

$$ \mathbb{P}(E_0 = j) = X_0^j / \langle \hat{\nu}_0, I_d \rangle. \quad (2.3) $$
Each colony of trait $x$ in the spine process divides into two children at a rate $\lambda x/\hat{R}_t$, where $\hat{R}_t := \langle \hat{v}_t, I_d \rangle$ is the sum of the traits of all colonies in the spine process. Any non-spinal colony dies at rate $\mu$ while the spine does not go extinct. If the spine divides at time $t$, the new spine is chosen between the two children with probability proportional to their traits, such that

$$ \mathbb{P}(E_t = E_{t-i}) = \frac{\hat{X}^{E_{t-i}}}{\hat{X}^{E_{t-1}} + \hat{X}^{E_{t-2}}}, \quad \text{for } i \in \{1, 2\}. \qquad (2.4) $$

We can now introduce the theorem that links the law of the non-extinct biological process of interest, that properly describes the colonies dynamics, and the law of the spine process. We denote by $\hat{N}_t$ the number of colonies at time $t \ge 0$ in the spinal process.

**Theorem 2.1.** For $U_t$ a colony picked uniformly at random in $\mathbb{G}(t)$ and for every measurable function $F$ on the set of càdlàg trajectories on $\mathbb{R}_+^* \times \mathbb{U}$, we have

$$ \mathbb{E}_{\nu_0} \left[ 1_{\{\mathbb{G}(t) \neq \emptyset\}} F \left( X_t^{U_t}, (\nu_s, s \le t) \right) \right] = R_0 e^{(a-\mu)t} \mathbb{E}_{\nu_0} \left[ \frac{1}{\hat{Y}_t \hat{N}_t} F \left( \hat{Y}_t, (\hat{v}_s, s \le t) \right) \right]. \quad (2.5) $$

Notice that this equality in law holds under the event of survival of the population. The Many-to-One formula associated to this $\psi$-spine construction can be deduced from Theorem 2.1 by taking the measurable function $F$ defined on the set of càdlàg trajectories on $\mathbb{R}_+^* \times \mathbb{U}$ as

$$ F(X_t^{U_t}, (\nu_s, s \le t)) := f((X_s^{U_t}, \nu_s), s \le t) \hat{N}_t. $$

**Corollary 2.2.** For every measurable function $f$ on the set of càdlàg trajectories on $\mathbb{R}_+^* \times \mathbb{U}$:

$$ \mathbb{E}_{\nu_0} \left[ 1_{\{\mathbb{G}(t) \neq \emptyset\}} \sum_{u \in \mathbb{G}(t)} f((X_s^u, \nu_s), s \le t) \right] = R_0 e^{(a-\mu)t} \mathbb{E}_{\nu_0} \left[ \frac{f((\hat{Y}_s, \hat{v}_s), s \le t)}{\hat{Y}_t} \right]. \quad (2.6) $$

Corollary 2.2 is a key result, at the basis of the next two sections. It allows us to derive some results on the long-time behavior of the population stated in Section 2.3. Besides, from Corollary 2.2, we can deduce an algorithmic method to efficiently estimate some quantities over the population (Section 2.4).

**2.3. Long-time behavior of the population.** Formula (2.6) allows us to link the expectation of key quantities of the dynamics of the process to expectations related to the spinal process. Let us first focus on quantities describing the global dynamics of the population process, that is the population size and the total population amount of resources.

**Proposition 2.3.** The population size $N_t$ follows the law of an $M/M/\infty$ queue with arrival rate $\lambda$ and service rate $\mu$. The mean extinction time is

$$ \mathbb{E}_{\nu_0}[T_0] = \sum_{k=0}^{N_0-1} \frac{k!}{\lambda(\lambda/\mu)^k} \sum_{j \ge k+1} \frac{(\lambda/\mu)^j}{j!} \qquad (2.7) $$

and the expected size satisfies

$$ \lambda t e^{-\mu t} + N_0 e^{-\mu t} \le \mathbb{E}_{\nu_0}[N_t] \le f(t) + N_0 e^{-\mu t} \qquad (2.8) $$

where the function $f$ satisfies $f(t) = O(1/t)$ as $t$ goes to infinity.

**Remark 2.4.** An $M/M/\infty$ queue with arrival rate $\lambda$ and service rate $\mu$ is a time changed logistic process of birth rate $\lambda$, death rate 0 and competition parameter $\mu$. Indeed, if the population is in state $N$, the next state of the population is $N+1$ with probability $\lambda/\mu N = \lambda N/\mu N^2$ for both $M/M/\infty$ queue and logistic process. Since the jump rates are higher in the case of the logistic process (multiplied by a factor $N$ when the population is in state $N$), we can deduce that the $M/M/\infty$ queue has an extinction time larger than the one of the mentioned logistic process.
A first consequence of Proposition 2.3 is that the law of $N_t$ only depends on the parameters $\lambda$ and $\mu$, and not on the noise parameters of the dynamics, $\sigma$ and $\delta$. This contrasts with classical results on the role of noise on the population persistence, which conclude that a higher noise increases the population death rate (see [21] for instance). The main difference of our model with cited works is that there is a strong interaction between colonies to have access to the resource, which translates into a competition in the birth rate. This fact illustrates that taking into account interactions between colonies has a major effect on the population dynamics, and that making the approximation that colonies behave independently may lead to erroneous conclusions.

<figure><img src="image_2.png" /><figcaption>FIGURE 1. Monte-Carlo estimation of $\mathbb{E}[N_t]$ using $M = 10^4$ trajectories, for different values of the division rate $\lambda$. Other parameters are $\mu = 1$ and $N_0 = 1$.</figcaption></figure>

From (2.8), we see that the mean population size is bounded by a quantity that decreases as $1/t$ for large $t$ whereas the lower bound decreases to 0 exponentially fast. In fact, we cannot find more precise boundaries without disjointing cases, because all the asymptotics between the two bounds are possible, as illustrated by Figure 1. This depends on the value of $\lambda$ compared to $\mu$.

**Proposition 2.5.** The expectation of the resources $R_t$ in the original population verifies

$$ \mathbb{E}_{\nu_0}[R_t] = R_0 e^{(a-\mu)t}. \qquad (2.9) $$

Besides, the variance of the total resource in the population satisfies

$$ \frac{\mu + \frac{\sigma^2}{1+\delta^2}}{R_0 \left( \mu + \frac{\sigma^2}{1+\delta^2} - 2\lambda \mathbb{E}[\Theta(1-\Theta)] \right)} \left( e^{(\mu+\sigma^2-2\lambda \mathbb{E}[\Theta(1-\Theta)])t} - e^{\frac{\sigma^2\delta^2}{1+\delta^2}t} \right) + e^{\frac{\sigma^2\delta^2}{1+\delta^2}t} \\ \le \frac{\mathrm{Var}_{\nu_0}[R_t]}{R_0^2 e^{2(a-\mu)t}} + 1 \le \frac{1}{R_0} e^{(\mu+\sigma^2)t} + \frac{R_0-1}{R_0} e^{\frac{\sigma^2\delta^2}{1+\delta^2}t}. $$

Notice that the lower bound is a continuous function of the parameters and the case $\mu + \sigma^2/(1+\delta^2) = 2\lambda \mathbb{E}[\Theta(1-\Theta)]$ is still well defined.

Equation (2.9) shows that the average resource in the population only depends on its Malthusian parameter $a$ and the individual death rate $\mu$. It is not affected by the strength of the noise. Its variance, however, strongly depends on the strength of the noise (parameter $\sigma$), as well as the correlation between individual stochasticities (parameter $\delta$). As we can see in Figure 2, the bounds found in Proposition 2.5 are a good approximation of the variance. From the simulations (we observe the same trend for both lower and upper bounds), we see that the variance of $R$ increases with $\sigma$ and decreases with $\mu$ as we would expect. Notice that the lower bound
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Monte-Carlo estimation of $\mathbb{Var}[R_t]$ using $M = 10^6$ trajectories, for different values of $\lambda, \mu$ and $\sigma$. Other parameters are $R_0 = 15, a = 5$ and $\delta = 1$.</figcaption></figure>

in Proposition 2.5 decreases with $\lambda$ while the upper bound does not depend on the fission rate. The simulations in Figure 2 show that for small division rates (the blue dots with $\lambda = 1$), both boundaries are close, while when $\lambda$ is higher (the green dots with $\lambda = 10$), these bounds are less optimal and the variance decreases.

Recall that the resources do not impact the branching mechanisms at the population level, while it does at an individual level. Hence the quantity of interest for studying the behavior of the colonies in the population is not their traits but the fraction of resources $x/R$ that they own. This is the content of the next result, which focuses on the distribution of these resources among the various colonies in the population.

**Proposition 2.6.** The average resource owned by a colony picked uniformly at random at a time $t$ satisfies, as $t$ goes to infinity

$$ \mathbb{E}_{\nu_0} \left[ 1_{\{\mathbb{G}(t) \neq \emptyset\}} X_t^{U_t} \right] \sim R_0 e^{(a-\mu)t} \frac{\mu}{\lambda} \left(1 - e^{-\lambda/\mu}\right). $$

Depending on the noise intensity and correlation, the proportion of resources owned by the most successful colonies may vary. Let $\mu, \lambda$ and $q$ (the distribution of $\Theta$) be fixed. Then:

* If $\sigma^2/(1 + \delta^2) \to \infty$, for all $\varepsilon > 0$

$$ \mathbb{E}_{\nu_0} \left[ \sum_{u \in \mathbb{G}(t)} X_t^u \mathbb{1}_{\left\{ \frac{X_t^u}{R_t} \ge 1 - \varepsilon \right\}} \right] \sim \mathbb{E}_{\nu_0}[R_t], \quad (t \to \infty). \quad (2.10) $$

* If $\sigma^2/(1 + \delta^2) \to 0$, for all $1/2 < \gamma < 1$,

$$ \limsup_{t \to \infty} \frac{\mathbb{E}_{\nu_0} \left[ \sum_{u \in \mathbb{G}(t)} X_t^u \mathbb{1}_{\left\{ \frac{X_t^u}{R_t} \ge \gamma \right\}} \right]}{\mathbb{E}_{\nu_0}[R_t]} \le \frac{\sqrt{1 + 8\lambda \mathbb{E}[\Theta(1-\Theta)]/\mu - 1}}{4\gamma\lambda \mathbb{E}[\Theta(1-\Theta)]/\mu} \xrightarrow{\lambda \mathbb{E}[\Theta(1-\Theta)]/\mu \to \infty} 0. \quad (2.11) $$

We can also obtain a lower bound on the maximal quantity of resources owned by a colony

$$ \liminf_{t \to \infty} \frac{\mathbb{E}_{\nu_0} \left[ \sup_{u \in \mathbb{G}(t)} X_t^u \right]}{\mathbb{E}_{\nu_0}[R_t]} \ge \left( \frac{\beta + \beta \frac{\mu}{\lambda} (1 - e^{-\lambda/\mu}) + 2}{2\beta - \frac{4\lambda}{\mu} \mathbb{E}[\Theta \ln \Theta] + 2} \right) \vee \frac{\mu}{\lambda} (1 - e^{-\lambda/\mu}), \quad (2.12) $$

where

$$ \alpha := \lambda \mathbb{E}[\Theta(1 - \Theta)] / \mu \quad \text{and} \quad \beta := \sigma^2 / (\mu (1 + \delta^2)). \quad (2.13) $$
The distribution of resources among colonies is thus strongly affected by the intensity of the fluctuations and their correlation (i.e. the $\sigma$ and $\delta$ parameters respectively). Equation (2.10) proves that if the total noise intensity $\sigma$ is large and the fluctuations are not too correlated ($\delta$ of a smaller order than $\sigma$), most of the resources are held by one colony. Equation (2.11) shows that whatever the noise intensity, if fluctuations are correlated, resources are more evenly shared among several colonies. Indeed, the upper bound is strictly smaller than 1 if $\gamma$ is close enough to 1. Equation (2.11) also shows that if the variance $\sigma$ of the noise in the individual traits is too low, the impact of $\delta$ is not significant. This is natural as $\delta$ represents the relative intensity of environmental noise compared to intrinsic noise. Consequently, if the total variance of the noise in trait dynamics is too small, two populations with different environmental sensitivities would still exhibit similar behaviors. In addition, when the fission rate is large, the impact of the environmental noise decreases as the population size is large. However, at a fixed fission rate, an increase in the death rate also decreases the effect of environmental noise. This is due to individuals dying too fast such as noise does not have a significant impact on their trajectories. This result illustrates the importance of distinguishing between intrinsic and environmental stochasticities to predict the long time behavior of a population. Figure 3 shows the sharing of resources between individuals, estimated with Monte-Carlo methods, at different times. The stationary distribution is reached in a short time scale independently of the parameters values. The transition between a unique 'winning' colony and a more homogeneous distribution of resources among colonies is continuous with the value of $\sigma^2/(1+\delta^2)$. However, when the individual traits are noise-sensitive, the distribution of resources colonies in two populations with different environmental noise levels is distinct.

<figure><img src="image_4.png" /><figcaption>FIGURE 3. Comparison of the mean sharing of resources at different times $T$ and different values of $\sigma^2/(1+\delta^2)$. The x axis gives the label of the individuals, sorted by trait in a decreasing order. The other parameters are $\lambda = 3, \mu = 1, a = 1, N_0 = 10$ and $R_0 = 10$.</figcaption></figure>

**2.4. Spinal algorithm.** As illustrated in Section 2.3, the spine process is a powerful tool to derive analytical results on key features of the population, that would otherwise be intractable. It can also be of great interest to consider spine processes when it comes to simulations of the population. As we have emphasized, a trajectory of the spine process does not describe a potential behavior of the colonial population, but Theorem 2.1 gives the right correction to apply to estimate expectations from spinal trajectories, using Monte-Carlo methods. This alternative simulation method classifies as an importance sampling method, see e.g. [52], as we use samples generated from a different distribution (the spinal one) to evaluate properties on the distribution of interest. The spinal Algorithm 1 gives the general method to construct a trajectory of the spine process, and the implemented methods in Python can be found in [45]. It is similar to the classical acceptance-rejection method, see e.g. [17] that can be used to draw trajectories of the colonial process, with some extra steps to handle the particular behavior of the spinal colony.

According to Algorithm 1, generating a trajectory of the spine process takes a short additional amount of time compared to the direct method, using trajectories of the colonial process.
<caption>**Algorithm 1** Spinal trajectory simulation method on $[0, T]$</caption>

Initialize time $t = 0$ and choose initial spine $E_0$, according to (2.3).

While $t < T$:

1. Propose the next branching time $t_{\text{next}}$.

2. Draw the colonies trajectories on $[t, t_{\text{next}} \wedge T)$ from (2.1) and (2.2).

3. Decide the type of branching event with rejection method.

If Spinal branching event:

Choose new spine according to (2.4).

4. Update current time $t \leftarrow t_{\text{next}} \wedge T$

However, it can be faster to use this method in two cases: when the variance of the estimated quantity is smaller in the spine process, and when the computations are easier in the spine process. The first point is illustrated in Figure 5A, where the goal is to estimate the variance of the total resources $R_t$ in the colonial process by using trajectories of $R_t^2$. In this case, the variance in the Monte-Carlo method is of the order $\mathbb{E}[R_t]^4$. However, from Corollary 2.2 with $f(x) = x^2$, we see that this value can be estimated from trajectories of $\hat{R}_t$, thus the variance in the Monte-Carlo method with spinal trajectories is of the order $\mathbb{E}[R_t]^2$. As a result, the evaluated quantity converges faster using spinal than colonial trajectories, and a smaller number of iterations $M$ is required from the spine process to reach the same precision. Overall, the spinal method converges in this case approximately 10 times faster than the direct method, despite the fact that the generation of one trajectory of the spine process takes a longer time. Another interesting point in this particular case is the fact that the trajectories of $R_t^2$ generated from the colonial process can reach really high values when the intrinsic noise $\sigma$ increases. In some trajectories, the highest possible value that can be stored in simple precision (approximately $10^{38}$, according to the norm IEEE 754) is exceeded, resulting in an underestimation of the real value, as we can see in Figure 4. A simple modification consists in using double precision floating numbers, that can reach up to $10^{308}$, but this will drastically increase both the computation time and the memory usage of the algorithm. The spinal method, on the other hand, uses trajectories of $\hat{R}_t$, that have a significantly lower probability to reach such high values.

<figure><img src="image_5.png" /><figcaption>FIGURE 4. Monte-Carlo estimations of $\mathbb{Var}[R_t]$ with $10^6$ trajectories from the direct and spinal methods. Parameter values are $a = 5$, $\lambda = 6$, $\delta = 1$, $R_0 = 30$.</figcaption></figure>

The second general case for which the spine method can be significantly faster than the direct one is when the estimated quantity needs heavy computations to generate trajectories of the colonial process, while the required trajectories from the spine process are much easier to draw. In Fig.5B, we estimated the mean size of a colony picked uniformly at random in the colonial population $\mathbb{E}[1_{G(t) \neq 0} X^{U_t}]$. Direct sampling trajectories need to store the vector of traits of all
colonies and pick uniformly one colony in this array, while the spine method only necessitates to sample trajectories of $1/\tilde{N}_t$ (Th.2.1 with $F(x) = x$). Thus the spine trajectories are drawn from a simple birth and death process with an immortal colony, that is, with the parameters chosen in Fig.5B, approximately 10 times faster to compute. Furthermore, in this example, direct trajectories where the population did not survive are thrown away, and in those who survived, a colony is uniformly picked in the population. These two mechanisms largely increase the value of the variance between trajectories, resulting in an effective convergence approximately $10^3$ times slower compared to spinal method. Employing the spinal method in this context thus yields a significant time efficiency improvement, with a gain on the order of $10^4$.

<figure><img src="image_6.png" /><figcaption>FIGURE 5. Comparison of simulation times and relative error for Monte Carlo estimates from direct and spinal methods. (A) $\mathbb{V}ar[R_t]$ estimation. Parameters are $a = 1, \lambda = 1, \mu = 0.3, \sigma = 1, \delta = 0.5, T = 2$. (B) Estimation of $\mathbb{E}[1_{G(t) \neq 0} X^{U_t}]$. Parameters are $a = 1, \lambda = 2, \mu = 0.5, \sigma = 0.5, \delta = 0.5, T = 10$.</figcaption></figure>

### 3. GENERAL POPULATION MODEL AND SPINAL CONSTRUCTION

In the following we will denote by $\mathfrak{B}(A, B)$ (resp. $\mathcal{C}^i(A, B)$) the set of measurable (resp. $i$ times continuously differentiable) $B$-valued functions defined on a set $A$. We denote by $\mathcal{U}$ the Ulam-Harris-Neveu set of labels of the individuals, and by $\mathcal{X} \subset \mathbb{R}^d$ the set of traits of the individuals in the population, for $d \ge 1$. In this section, the framework is not restricted to colonies modeling and thus we will use the term individual in all generality.

Let us introduce $\bar{\mathbb{V}}$ the set composed of all finite point measures on $\mathcal{U} \times \mathcal{X}$

$$ \bar{\mathbb{V}} := \left\{ \sum_{i=1}^{N} \delta_{(u^i, x^i)}, N \in \mathbb{N}, (u^i, 1 \le i \le N) \in \mathcal{U}^N, (x^i, 1 \le i \le N) \in \mathcal{X}^N \right\}. $$

We also define the set of marginal population measures, that is

$$ \mathbb{V} := \left\{ \sum_{i=1}^{N} \delta_{x^i}, N \in \mathbb{N}, (x^i, 1 \le i \le N) \in \mathcal{X}^N \right\}. $$

For any measure $\bar{\nu} = \sum_{i=1}^{N} \delta_{(u^i, x^i)}$ in $\bar{\mathbb{V}}$, we will write $\nu := \sum_{i=1}^{N} \delta_{x^i}$ its projection on $\mathbb{V}$. By convention, if the number of points in the measure is $N = 0$, $\bar{\nu}$ and $\nu$ are the trivial zero measures on $\mathcal{U} \times \mathcal{X}$ and $\mathcal{X}$. We introduce for every $\bar{\nu} \in \bar{\mathbb{V}}$, every $g \in \mathfrak{B}(\mathcal{U} \times \mathcal{X}, \mathbb{R})$ and every $f \in \mathfrak{B}(\mathcal{X}, \mathbb{R})$

$$ \langle \bar{\nu}, g \rangle := \int_{\mathcal{U} \times \mathcal{X}} g(u, x) \bar{\nu}(du, dx), \quad \text{and} \quad \langle \nu, f \rangle := \int_{\mathcal{X}} f(x) \nu(dx). $$

Finally, we denote by $\mathbb{D}(A, B)$ the Skorohod space of $B$-valued càdlàg functions on a subset $A$ of $\mathbb{R}_+$. For every process $(\nu_t, t \in A) \in \mathbb{D}(A, B)$ and $z \in B$, we will denote $\mathbb{E}_z[f(\nu_t)] := \mathbb{E}[f(\nu_t) | \nu_0 = z]$.
3.1. **Definition of the population.** The population is described at any time $t \in \mathbb{R}_+$ by the finite point measure $\bar{\nu}_t \in \bar{\mathbb{V}}$ giving the label and the trait of every individual living in the population at this time. We introduce

$$ \mathbb{G}(t) := \left\{ u \in \mathcal{U} : \int_{\mathcal{U} \times \mathcal{X}} \mathbb{1}_{\{v=u\}} \bar{\nu}_t (dv, dx) \geq 1 \right\} $$

the set of labels of living individuals at time $t$. For every individual labeled by $u \in \mathbb{G}(t)$ we denote by $X_t^u$ its trait, and with a slight abuse of notation, $X_s^u$ will denote the trait of its unique ancestor living at time $s \in [0, t]$. Using these notations, we can write

$$ \bar{\nu}_t = \sum_{u \in \mathbb{G}(t)} \delta_{(u, X_t^u)} \quad \text{and} \quad \nu_t := \sum_{u \in \mathbb{G}(t)} \delta_{X_t^u}. $$

The initial population is given by a measure $\bar{z} = \sum_{i=1}^N \delta_{(u^i, x^i)}$, and the branching rate by a function $B \in \mathfrak{B}(\mathcal{X} \times \mathbb{V}, \mathbb{R}_+)$. An individual with trait $x$ in a population $\nu$ dies at an instantaneous rate $B(x, \nu)$. It gives birth to $n$ offspring, where $n$ is randomly chosen with distribution $(p_k(x, \nu), k \in \mathbb{N})$ of first moment $m(x, \nu)$. Hence branching events leading to $n$ children happen at rate $B_n(\cdot) := B(\cdot)p_n(\cdot)$. If $n \ge 0$, we denote by $\mathbf{y} \in \mathcal{X}^n$ the offspring traits at birth, randomly chosen according to the law $K_n(x, \nu, \cdot)$. For all $1 \le i \le n$, the $i$-th child of a parent of label $u$ is labeled $ui$ and its trait is $y^i$, the $i$-th coordinate of $\mathbf{y}$. We denote by $T_{\text{Exp}} \in \bar{\mathbb{R}}_+$ the explosion time of the process $(\bar{\nu}_t, t \ge 0)$, defined as the limit of its jumps times $(T_k, k \ge 0)$.

Recall that $N_t$ denotes the number of individuals in the population at time $t \ge 0$. It gives the dimension of the vectorial objects. For all $n \ge 1$, we introduce a number $q(n) \ge n$ that denotes the number of sources of noise in the trait dynamics. For all $q \ge n \ge 1$, we introduce $\mathbf{a}^{(n)} \in \mathfrak{B}(\bar{\mathbb{V}}, \mathcal{X}^n)$ and $\boldsymbol{\sigma}^{(n,q)} \in \mathfrak{B}(\bar{\mathbb{V}}, M_{n,q}(\mathcal{X}))$ the drift and diffusion coefficients, where $M_{n,q}(\mathcal{X})$ denotes the set of all $n \times q$ matrices with coefficients in $\mathcal{X}$. We also denote by $\mathbf{X}_t := (X_t^u, u \in \mathbb{G}(t))$ the vector of $\mathcal{X}^{N_t}$ containing the traits of every living individuals in the population. Between two successive branching events, the traits of the individuals alive in the population evolve according to population-dependent dynamics

$$ d\mathbf{X}_t = \mathbf{a}^{(N_t)}(\bar{\nu}_t) dt + \boldsymbol{\sigma}^{(N_t, q(N_t))}(\bar{\nu}_t) d\mathbf{B}_t \quad (3.1) $$

where $\mathbf{B}_t$ is a $q(N_t)$-dimensional Wiener process. The case $q(n) > n$ corresponds to an external noise added to the trait dynamics of each individual, that can be for example, the impact of a common environment that correlates the individuals. Such a process can be seen as a generalization of a multi-type continuous-state branching process in random environment [11, 8]. The population model introduced in Section 2 and studied in this paper is an example of a process satisfying (3.1).

The formalism used in (3.1) also allows to consider models where the individual dynamics differ between individuals, and change with the population size. Notice that this formalism will be useful to derive further expressions of operators and generators of the branching process.

In Section 3.2, we establish a $\psi$-spinal construction that distinguishes an individual in the population, and derive a Many-to-One formula in the context of branching processes with complex interactions. To this end we will need additional notations. The spinal construction generates a $\bar{\mathbb{V}}$-valued process along with the label of a distinguished individual that can change with time. For convenience, we will denote by $\bar{\mathbb{W}}$ and $\mathbb{W}$ the sets such that

$$ \bar{\mathbb{W}} := \{ (\mathfrak{e}, \bar{\nu}) \in \mathcal{U} \times \bar{\mathbb{V}} : \langle \bar{\nu}, \mathbb{1}_{\{\mathfrak{e}\}} \times \mathcal{X} \rangle \ge 1 \}, \quad \text{and} \quad \mathbb{W} := \{ (x, \nu) \in \mathcal{X} \times \mathbb{V} : \langle \nu, \mathbb{1}_{\{x\}} \rangle \ge 1 \}. $$

For every $(\mathfrak{e}, \bar{\nu}) \in \bar{\mathbb{W}}$, we uniquely define $x_{\mathfrak{e}}$ the trait of the individual of label $\mathfrak{e}$ in the population $\bar{\nu}$. Thus, the spine process is a $\bar{\mathbb{W}}$-valued branching process and its marginal is a $\mathbb{W}$-valued branching process. We will propose a generalized spinal construction, where branching rates are biased with weight functions chosen in a set $\mathcal{D}$, defined by

$$ \mathcal{D} := \{ F_f \in \mathfrak{B}(\bar{\mathbb{W}}, \mathbb{R}_+^*) \text{ s.t. } (f, F) \in \mathcal{C}^2(\mathcal{X}, \mathbb{R}) \times \mathcal{C}^2(\mathcal{X} \times \mathbb{R}, \mathbb{R}_+^*) \}, \quad (3.2) $$
where for every $(\mathfrak{e}, \bar{\nu}) \in \overline{\mathbb{W}}$, $F_f(\mathfrak{e}, \bar{\nu}) := F(x_\mathfrak{e}, \langle \nu, f \rangle)$. In order to alleviate the notations when there is no ambiguity, we will omit the subscript $f$. In the following we will also omit the superscript indicating the dimension of the drift and diffusion functions, and refer to the dynamics of the individual $u$ with a subscript $u$ in those functions.

**Proposition 3.1.** The infinitesimal generator $G$ of the trait process solution of (3.1), defined on the set $\mathcal{D}$ is given, for all functions $F_f \in \mathcal{D}$ and all $(\mathfrak{e}, \bar{\nu}) \in \overline{\mathbb{W}}$ by:

$$ 
\begin{aligned}
GF_f(\mathfrak{e}, \bar{\nu}) = & \partial_1 F(x_\mathfrak{e}, \langle \nu, f \rangle) a_\mathfrak{e}(\bar{\nu}) + \frac{1}{2} \partial_{1,1}^2 F(x_\mathfrak{e}, \langle \nu, f \rangle) (\boldsymbol{\sigma}\boldsymbol{\sigma}^T)_{\mathfrak{e},\mathfrak{e}}(\bar{\nu}) \\
& + \partial_{1,2}^2 F(x_\mathfrak{e}, \langle \nu, f \rangle) \int_{\mathcal{X} \times \mathcal{U}} (\boldsymbol{\sigma}\boldsymbol{\sigma}^T)_{u,\mathfrak{e}}(\bar{\nu}) f'(x) \bar{\nu}(dx, du) \\
& + \partial_2 F(x_\mathfrak{e}, \langle \nu, f \rangle) \int_{\mathcal{X} \times \mathcal{U}} \left( a_u(\bar{\nu}) f'(x) + \frac{(\boldsymbol{\sigma}\boldsymbol{\sigma}^T)_{u,u}}{2} (\bar{\nu}) f''(x) \right) \bar{\nu}(dx, du) \\
& + \frac{1}{2} \partial_{2,2}^2 F(x_\mathfrak{e}, \langle \nu, f \rangle) \int_{\mathcal{X} \times \mathcal{U}} \int_{\mathcal{X} \times \mathcal{U}} (\boldsymbol{\sigma}\boldsymbol{\sigma}^T)_{u,v}(\bar{\nu}) f'(x) f'(y) \bar{\nu}(dx, du) \bar{\nu}(dy, dv).
\end{aligned} 
\quad (3.3) $$

This proposition is proven in Section 4.1. The first two lines are related to the dynamic of a distinguished individual in the population. The first order derivative describes the drift of its trait and the second order derivatives are the diffusion parts. Notice that, thanks to the interactions, the diffusion dynamics are separated in an autocorrelated part (in the first line) and a sum of the correlations from the other individuals in the second line. The last two lines describe the dynamics of all the individuals- also with the distinguished one- in a similar way.

**Corollary 3.2.** When the individuals follow the same dynamics and are not correlated through the Brownian motion, the drift vector becomes $\mathbf{a}(\bar{\nu}_t) = (a(X_t^u, \nu_t), u \in \mathbb{G}(t))$ and the diffusion matrix is diagonal; $\boldsymbol{\sigma}(\bar{\nu}_t) = \text{diag}(\sigma(X_t^u, \nu_t), u \in \mathbb{G}(t))$. In this case the infinitesimal generator $G$ of the trait process is given, for all functions $F_f \in \mathcal{D}$ and all $(\mathfrak{e}, \bar{\nu}) \in \overline{\mathbb{W}}$ by:

$$ 
\begin{aligned}
GF_f(\mathfrak{e}, \bar{\nu}) = & \partial_1 F(x_\mathfrak{e}, \langle \nu, f \rangle) a(x_\mathfrak{e}, \nu) + \frac{1}{2} \partial_{1,1}^2 F(x_\mathfrak{e}, \langle \nu, f \rangle) \sigma^2(x_\mathfrak{e}, \nu) \\
& + \partial_2 F(x_\mathfrak{e}, \langle \nu, f \rangle) \int_{\mathcal{X}} a(x, \nu) f'(x) + \frac{1}{2} f''(x) \sigma^2(x, \nu) \nu(dx).
\end{aligned} 
$$

This form is useful in many models without interactions or when the trait dynamics are indistinguishable between individuals [9, 11].

**3.2. ψ-spine construction.** For $(u, \bar{\nu}) \in \overline{\mathbb{W}}$ we recall that $x_u$ is the trait of the individual of label $u$ in the population $\bar{\nu}$, and for $n \ge 0$ and $\mathbf{y} = (y^i, 1 \le i \le n) \in \mathcal{X}^n$ we will denote by

$$ \bar{\nu}_+(u, \mathbf{y}) := \bar{\nu} - \delta_{(u,x_u)} + \sum_{i=1}^{n} \delta_{(u_i, y^i)}, \quad \text{and} \quad \nu^+(x, \mathbf{y}) := \nu - \delta_x + \sum_{i=1}^{n} \delta_{y^i} \quad (3.4) $$

the state of the population, just after a branching event replacing the individual $u$ by the offspring of traits $\mathbf{y}$. We now introduce the key operator $\mathcal{G}$ involved in the spinal construction. It is defined for all $F \in \mathcal{D}$ and $(\mathfrak{e}, \bar{\nu}) \in \overline{\mathbb{W}}$ by

$$ 
\begin{aligned}
\mathcal{G}F(\mathfrak{e}, \bar{\nu}) := & GF(\mathfrak{e}, \bar{\nu}) + \sum_{n \ge 0} B_n(x_\mathfrak{e}, \nu) \int_{\mathcal{X}^n} \left[ \sum_{i=1}^n F(\mathfrak{e}i, \bar{\nu}^+(x_\mathfrak{e}, \mathbf{y})) - F(\mathfrak{e}, \bar{\nu}) \right] K_n(x_\mathfrak{e}, \nu, \mathrm{d}\mathbf{y}) \\
& + \sum_{n \ge 0} \int_{\mathcal{X} \times \mathcal{U} \setminus \{\mathfrak{e}\}} B_n(x, \nu) \int_{\mathcal{X}^n} \left[ F(\mathfrak{e}, \bar{\nu}^+(u, \mathbf{y})) - F(\mathfrak{e}, \bar{\nu}) \right] K_n(x, \nu, \mathrm{d}\mathbf{y}) \bar{\nu}(dx, du),
\end{aligned} 
\quad (3.5) $$

where $G$ is the generator of the traits evolution between branching events defined in (3.3). The first line corresponds to the branching event of the distinguished individuals and the choice of the new one among its children. The second line describes the branching events for every other individuals. By convention we consider that $K_0$ is the null measure on the set $\mathcal{X}^0 = \emptyset$.
Note that the operator $\mathcal{G}$ is generally not the generator of a conservative Markov process on $\mathbb{W}$.
We recall that every function $F_f$ in the set $\mathcal{D}$ defined in (3.2) is such that $F$ is $\mathcal{C}^2(\mathcal{X} \times \mathbb{R}, \mathbb{R}_+^*)$,
thus the functions $\mathcal{G}F$ are locally bounded on $\overline{\mathbb{W}}$.

Now we can introduce the $\overline{\mathbb{W}}$-valued spine process associated with a non negative bias function $\psi_\phi \in \mathcal{D}$, that will be called the $\psi_\phi$-spine process in order to keep track of the function used in its construction. For such functions $\psi_\phi$, the generator $\hat{\mathcal{L}}_{\psi_\phi}$ of the $\psi_\phi$- spine process is defined for all $F_f \in \mathcal{D}$ and all $(\mathfrak{e}, \bar{\nu}) \in \overline{\mathbb{W}} \times \mathbb{R}_+$ by

$$ \hat{\mathcal{L}}_{\psi_\phi} F_f(\mathfrak{e}, \bar{\nu}) := \frac{\mathcal{G}[\psi_\phi F_f]}{\psi_\phi}(\mathfrak{e}, \bar{\nu}) - \frac{\mathcal{G}\psi_\phi}{\psi_\phi}(\mathfrak{e}, \bar{\nu}) F_f(\mathfrak{e}, \bar{\nu}). \quad (3.6) $$

This generator can also be expressed in the following form, which is more descriptive and allows for a better understanding of the dynamics of individuals in the spinal process.

**Proposition 3.3.** The generator of the spinal process satisfies for $F_f \in \mathcal{D}$ and $(\mathfrak{e}, \bar{\nu}) \in \overline{\mathbb{W}} \times \mathbb{R}_+$

$$ \begin{aligned} \hat{\mathcal{L}}_{\psi_\phi} F_f(\mathfrak{e}, \bar{\nu}) := & \hat{G} F_f(\mathfrak{e}, \bar{\nu}) \\ & + \sum_{n \ge 0} B_n(x_{\mathfrak{e}}, \nu) \int_{\mathcal{X}^n} \sum_{i=1}^n \left[ F_f(\mathfrak{e}i, \bar{\nu}^+(x_{\mathfrak{e}}, \mathbf{y})) - F_f(\mathfrak{e}, \bar{\nu}) \right] \frac{\psi_\phi(y^i, \nu^+(x_{\mathfrak{e}}, \mathbf{y}))}{\psi_\phi(x_{\mathfrak{e}}, \nu)} K_n(x_{\mathfrak{e}}, \nu, d\mathbf{y}) \\ & + \sum_{n \ge 0} \int_{\mathcal{X} \times \mathcal{U} \setminus \{\mathfrak{e}\}} B_n(x, \nu) \int_{\mathcal{X}^n} \left[ F_f(u, \bar{\nu}^+(\mathfrak{e}, \mathbf{y})) \frac{\psi_\phi(x_{\mathfrak{e}}, \nu^+(x_u, \mathbf{y}))}{\psi_\phi(x_{\mathfrak{e}}, \nu)} - F_f(\mathfrak{e}, \bar{\nu}) \right] \\ & \times K_n(x, \nu, d\mathbf{y}) \bar{\nu}(dx, du) \end{aligned} $$

where $\nu^+$ is defined in (3.4). The continuous part is given by (with $\mathcal{X} = \mathcal{X} \times \mathcal{U}$)

$$ \begin{aligned} \hat{G} F_f(\mathfrak{e}, \bar{\nu}) ={}& \partial_1 F \left( \mathbf{a}_{\mathfrak{e}} + \frac{\partial_1 \psi}{\psi_{\phi}} (\boldsymbol{\sigma} \boldsymbol{\sigma}^T)_{\mathfrak{e}, \mathfrak{e}} + \frac{\partial_2 \psi}{\psi_{\phi}} \int_{\mathcal{X}} (\boldsymbol{\sigma} \boldsymbol{\sigma}^T)_{\mathfrak{e}, \nu} \phi'(x) \bar{\nu}(dx, du) \right) \\ & + \frac{1}{2} \partial_{1,1}^2 F (\boldsymbol{\sigma} \boldsymbol{\sigma}^T)_{\mathfrak{e}, \mathfrak{e}} + \partial_{1,2}^2 F \int_{\mathcal{X}} (\boldsymbol{\sigma} \boldsymbol{\sigma}^T)_{u, \mathfrak{e}} f'(x) \bar{\nu}(dx, du) \\ & + \partial_2 F \int_{\mathcal{X}} f'(x) \left[ \mathbf{a}_u + \frac{\partial_1 \psi}{\psi_{\phi}} (\boldsymbol{\sigma} \boldsymbol{\sigma}^T)_{u, \mathfrak{e}} + \frac{\partial_2 \psi}{\psi_{\phi}} \int_{\mathcal{X}} (\boldsymbol{\sigma} \boldsymbol{\sigma}^T)_{u, \nu} \phi'(y) \bar{\nu}(dy, dv) \right] \bar{\nu}(dx, du) \\ & + \partial_2 F \int_{\mathcal{X}} \frac{(\boldsymbol{\sigma} \boldsymbol{\sigma}^T)_{u, u}}{2} f''(x) \bar{\nu}(dx, du) + \frac{1}{2} \partial_{2,2}^2 F \int_{\mathcal{X}} \int_{\mathcal{X}} (\boldsymbol{\sigma} \boldsymbol{\sigma}^T)_{u, \nu} f'(x) f'(y) \bar{\nu}(dx, du) \bar{\nu}(dy, dv). \end{aligned} $$

We omitted the dependencies $(x_{\mathfrak{e}}, \langle \nu, f \rangle)$ and $(\bar{\nu})$ to shorten the notations, and used the notations

$$ \partial_i \psi := \partial_i \psi(x_{\mathfrak{e}}, \langle \nu, \phi \rangle), \quad \text{for } i \in \{1, 2\}. $$

From the expression of the continuous part of the generator, we get that the traits dynamics between jumps in the spinal process are the same than those of the original process with an added drift in the direction $\boldsymbol{\sigma}\boldsymbol{\sigma}^T\nabla\psi_\phi$:

$$ \hat{G} F_f(\mathfrak{e}, \bar{\nu}) = G F_f(\mathfrak{e}, \bar{\nu}) + \left( \frac{\boldsymbol{\sigma}\boldsymbol{\sigma}^T \nabla \psi_\phi}{\psi_\phi} \right)^T \nabla F_f(\mathfrak{e}, \bar{\nu}) \quad (3.7) $$

where the gradient is defined as

$$ \nabla F_f(\mathfrak{e}, \bar{\nu}) := \partial_1 F(x_{\mathfrak{e}}, \langle \nu, f \rangle) + \partial_2 F(x_{\mathfrak{e}}, \langle \nu, f \rangle) \int_{\mathcal{X} \times \mathcal{U}} f'(x) \bar{\nu}(dx, du). $$

The proof of Proposition 3.3 and the vectorial expression of the added drift are displayed in Section 4.1. Notice that the branching rates of both spinal and non-spinal individuals are biased by the function $\psi_\phi$. Dynamics 2.3 and 2.4 in [46] provide a more detailed explanation of the branching mechanisms in the spine process associated with this generator. The spinal construction is based on a martingale change of measure, that necessitates a moment assumption on the original process and the function $\psi_\phi$ involved in this change of measure. A simple condition to ensure this moment assumption is a boundary condition on the added drift.
**Assumption 1.** There exists $M \in \mathbb{R}$ such that for all $(e, \bar{\nu}) \in \bar{\mathbb{W}} \times \mathbb{R}_+$

$$ (\partial_2 \psi)^2 (e, \bar{\nu}) \int_{\mathcal{X} \times \mathcal{U}} \int_{\mathcal{X} \times \mathcal{U}} (\boldsymbol{\sigma} \boldsymbol{\sigma}^T)_{u,v} \phi'(y) \phi'(x) \bar{\nu}(dy, dv) \bar{\nu}(dx, du) \\ + (\partial_1 \psi)^2 (e, \bar{\nu}) (\boldsymbol{\sigma} \boldsymbol{\sigma}^T)_{\varepsilon, \varepsilon} + 2 \partial_1 \psi \partial_2 \psi (e, \bar{\nu}) \int_{\mathcal{X} \times \mathcal{U}} \phi'(x) (\boldsymbol{\sigma} \boldsymbol{\sigma}^T)_{u, \varepsilon} \bar{\nu}(dx, du) \leq M \psi(e, \bar{\nu}). $$

This general assumption does not take into account the branching dynamics of the spinal process. It can be lightened by considering an assumption on a path-integral exponential moment of the spinal process, under the form of the Novikov condition. For more details, we refer the interested reader to the proof of Theorem 3.4 in Section 4.1.3. Finally, we introduce for all $t \ge 0$, the $\mathcal{U}$-valued random variable $U_t$ that picks an individual alive at time $t$. Its law is characterized by the function $p_u(\bar{\nu}_t)$ which yields the probability to choose the individual of label $u$ in the set $\mathbb{G}(t)$. We can now state the main result of this section, that is a Girsanov-type formula for the spinal change of measure. It characterizes the joint probability distribution of $(U_t, (\bar{\nu}_s, s \le t))$- that is the randomly sampled individual in the population $\bar{\nu}_t$ at time $t$ and the whole trajectory of the population until this time- and links it to the law of the spine process through a path-integral formula.

**Theorem 3.4.** Let $\psi_\phi \in \mathcal{D}$ satisfying Assumption 1, $t \ge 0$, $\bar{z} \in \bar{\mathbb{V}}$. Let $((E_t, \bar{\nu}_t), t \ge 0)$ be the $\bar{\mathbb{W}}$-valued branching process with interactions defined by the infinitesimal generator $\hat{\mathcal{L}}_{\psi_\phi}$ introduced in (3.6), and $\hat{T}_{\text{Exp}}$ its explosion time. For all measurable function $H$ on $\mathcal{U} \times \mathbb{D}([0, t], \mathbb{V})$:

$$ \mathbb{E}_{\bar{z}} \left[ 1_{\{T_{\text{Exp}} > t, \mathbb{G}(t) \neq \emptyset\}} H \left( U_t, (\bar{\nu}_s, s \le t) \right) \right] = \langle z, \psi_\phi(\cdot, z) \rangle \mathbb{E}_{\bar{z}} \left[ 1_{\{\hat{T}_{\text{Exp}} > t\}} \xi \left( E_t, (\bar{\nu}_s, s \le t) \right) H \left( E_t, (\bar{\nu}_s, s \le t) \right) \right], $$

where:

$$ \xi \left( E_t, \left( \bar{\nu}_s, s \le t \right) \right) := \frac{p_{E_t}(\bar{\nu}_t)}{\psi_\phi(E_t, \bar{\nu}_t)} \exp \left( \int_0^t \frac{\mathcal{G}\psi_\phi(E_s, \bar{\nu}_s)}{\psi_\phi(E_s, \bar{\nu}_s)} ds \right). $$

The proof of this result is an extension of the proof for the case $\sigma = 0$ in [46]. We show by induction on the successive jumps times the equality in law for the trajectories between these times, using a Girsanov theorem to take into account the Brownian trajectories.

The process $((E_t, \bar{\nu}_t), t \ge 0)$ gives at any time $t \ge 0$ the label of the spinal individual- that encodes the whole spine lineage- and the spinal population. Our result thus links, for every $\psi_\phi$, the sampling of an individual and the trajectory of the population to the trajectory of the spine process. This formula gives a change of probability that involves the function $\mathcal{G}\psi_\phi/\psi_\phi$ with a path-integral formula, related to Feynman-Kac path measures and semi-groups [23]. The path integral term that links these two terms is difficult to handle in general and finding eigenfunctions of $\mathcal{G}$ allows a dual representation of the first moment semi-group of the original population with the generator of the spinal individual [3, 27, 20, 6].

**Corollary 3.5.** Let $\psi_\phi \in \mathcal{D}$ satisfying Assumption 1, such that there exists $\lambda_0 \in \mathbb{R}$ satisfying for all $(e, \bar{z}) \in \bar{\mathbb{W}}$.

$$ \frac{\mathcal{G}\psi_\phi(e, \bar{z})}{\psi_\phi(e, \bar{z})} = \lambda_0. \qquad (3.8) $$

Let $((E_t, \bar{\nu}_t), t \ge 0)$ be the $\bar{\mathbb{W}}$-valued branching process with interactions defined by the infinitesimal generator $\hat{\mathcal{L}}_{\psi_\phi}$. Let $t \ge 0$, and $\bar{z} \in \bar{\mathbb{V}}$. If the original process and the $\psi_\phi$-spine process do not explode a.s. before time $t$, then for every measurable function $f$ on $\mathcal{U} \times \mathbb{D}([0, t], \mathcal{X})$:

$$ \mathbb{E}_{\bar{z}} \left[ \sum_{u \in \mathbb{G}(t)} f(X_s^u, s \le t) \right] = \langle z, \psi_\phi(\cdot, \bar{z}) \rangle e^{\lambda_0 t} \mathbb{E}_{\bar{z}} \left[ \frac{1}{\psi_\phi(E_t, \bar{\nu}_t)} f(\hat{Y}_t, s \le t) \right]. $$
Proof. Let $f$ be a measurable function on $\mathcal{U} \times \mathbb{D}([0, t], \mathcal{X})$. We introduce $H$ the measurable function defined for all $(u, \bar{z}_s, s \le t) \in \mathcal{U} \times \mathbb{D}([0, t], \bar{\mathbb{V}})$ by

$$H(u, \bar{z}_s, s \le t) := f(X_s^u, s \le t) \langle z_t, 1 \rangle.$$

The corollary is thus a direct application of Theorem 3.4 with the function $H$ and a uniformly sampled individual under non-explosion conditions, see i.e. [42] for sufficient conditions. $\square$

For structured branching processes without interactions in some particular cases, it is possible to exhibit functions satisfying (3.8). For instance, for Ornstein-Uhlenbeck processes, Englander [27] derived eigenfunctions in the case of constant or quadratic breeding rate, and Cloez [20] in the affine breeding case. In more general cases, it is possible to partially solve the eigenproblem and derive only the eigenvalue $\lambda_0$ that gives the exponential growth of the population [20][Section 6]. In the case of non-structured branching processes with size-dependent reproduction rates, i.e. for all $k \ge 0$, $B_k(x, \nu) = B_k(\langle \nu, 1 \rangle)$, Bansaye [6][Section 3] established that the function $\psi_\phi(x, \nu) = 1/\langle \nu, 1 \rangle$ always satisfies (3.8). However, finding such functions in the structured case with interactions is particularly complex as the operator $\mathcal{G}$ is non-linear.

## 4. PROOFS

The last section is dedicated to the proofs of the results collected in the different sections of the paper. We begin by presenting the proofs from Section 3, which are essential for later establishing the results discussed in Section 2.

### 4.1. Proofs of Section 3.

#### 4.1.1. *Proof of Proposition 3.1.* We denote by $\bar{\nu}_t$ the population at time $t$ and we choose $\mathcal{E} \in \mathcal{G}(t)$ for the label of the distinguished individual. Using Itô formula we get that

$$ 
\begin{aligned}
dF_f(\mathcal{E}, \bar{\nu}_t) &= dF(X_t^\mathcal{E}, \sum_{v \in \mathcal{G}(t)} f(X_t^v)) \\
&= \sum_{u \in \mathcal{G}(t)} \partial_u F(X_t^\mathcal{E}, \langle \nu_t, f \rangle) a_u(\bar{\nu}_t) dt + \sum_{j=1}^m \sum_{u \in \mathcal{G}(t)} \partial_u F(X_t^\mathcal{E}, \langle \nu_t, f \rangle) \sigma_{u,j}(\bar{\nu}_t) dB_t^j \\
&\quad + \frac{1}{2} \sum_{u,w \in \mathcal{G}(t)} \partial_{u,w}^2 F(X_t^\mathcal{E}, \langle \nu_t, f \rangle) (\boldsymbol{\sigma}\boldsymbol{\sigma}^T)_{u,w}(\bar{\nu}_t) dt.
\end{aligned} 
\quad (4.1) $$

To derive the generator we will need to compute the following terms. For all $u \ne \mathcal{E}$

$$ 
\begin{aligned}
\partial_u F(X_t^\mathcal{E}, \langle \nu_t, f \rangle) &= \partial_2 F(X_t^\mathcal{E}, \langle \nu_t, f \rangle) f'(X_t^u) \\
\partial_\mathcal{E} F(X_t^\mathcal{E}, \langle \nu_t, f \rangle) &= \partial_1 F(X_t^\mathcal{E}, \langle \nu_t, f \rangle) + \partial_2 F(X_t^\mathcal{E}, \langle \nu_t, f \rangle) f'(X_t^\mathcal{E}).
\end{aligned} 
\quad (4.2) $$

We deduce that

$$ \sum_{u \in \mathcal{G}(t)} \partial_u F(X_t^\mathcal{E}, \langle \nu_t, f \rangle) a_u(\bar{\nu}_t) = \partial_1 F(X_t^\mathcal{E}, \langle \nu_t, f \rangle) a_\mathcal{E}(\bar{\nu}_t) + \partial_2 F(X_t^\mathcal{E}, \langle \nu_t, f \rangle) \sum_{u \in \mathcal{G}(t)} f'(X_t^u) a_u(\bar{\nu}_t). $$

Let us denote by $\hat{\mathcal{G}}(t)$ the set of labels without the spine $\mathcal{G}(t)\setminus\{E_t\}$. Then for $u \ne v \in \hat{\mathcal{G}}(t)$, and when $E_t = \mathcal{E}$, the second order derivatives write

$$ 
\begin{aligned}
\partial_{u,u}^2 F(X_t^\mathcal{E}, \langle \nu_t, f \rangle) &= \partial_2 F(X_t^\mathcal{E}, \langle \nu_t, f \rangle) f''(X_t^u) + \partial_{2,2}^2 F(X_t^\mathcal{E}, \langle \nu_t, f \rangle) (f'(X_t^u))^2 \\
\partial_{u,v}^2 F(X_t^\mathcal{E}, \langle \nu_t, f \rangle) &= \partial_{2,2}^2 F(X_t^\mathcal{E}, \langle \nu_t, f \rangle) f'(X_t^u) f'(X_t^v) \\
\partial_{u,e}^2 F(X_t^\mathcal{E}, \langle \nu_t, f \rangle) &= \partial_{2,2}^2 F(X_t^\mathcal{E}, \langle \nu_t, f \rangle) f'(X_t^u) f'(X_t^\mathcal{E}) + \partial_{1,2}^2 F(X_t^\mathcal{E}, \langle \nu_t, f \rangle) f'(X_t^u) \\
\partial_{e,e}^2 F(X_t^\mathcal{E}, \langle \nu_t, f \rangle) &= \partial_2 F(X_t^\mathcal{E}, \langle \nu_t, f \rangle) f''(X_t^\mathcal{E}) + \partial_{2,2}^2 F(X_t^\mathcal{E}, \langle \nu_t, f \rangle) (f'(X_t^\mathcal{E}))^2 \\
&\quad + \partial_{1,1}^2 F(X_t^\mathcal{E}, \langle \nu_t, f \rangle) + 2\partial_{1,2}^2 F(X_t^\mathcal{E}, \langle \nu_t, f \rangle) f'(X_t^\mathcal{E}).
\end{aligned} 
$$
We will denote by $M$ the symmetrical matrix $\sigma\sigma^T$ and we omit the dependency in the functions when there is no ambiguity. We thus have

$$
\begin{align*}
\frac{1}{2} \sum_{u,w \in \mathbb{G}(t)} \partial_{u,w}^2 F M_{u,w} &= \frac{1}{2} \sum_{u,w \in \mathbb{G}(t)} \partial_{u,w}^2 F M_{u,w} + \sum_{u \in \mathbb{G}(t)} \partial_{u,\mathcal{E}}^2 F M_{u,\mathcal{E}} + \frac{1}{2} \partial_{\mathcal{E},\mathcal{E}}^2 F M_{\mathcal{E},\mathcal{E}} \\
&= \frac{1}{2} \partial_{2,2}^2 F \sum_{u,w \in \mathbb{G}(t)} f'(X_t^u) f'(X_t^v) M_{u,w} + \frac{1}{2} \partial_{2} F \sum_{u \in \mathbb{G}(t)} f''(X_t^u) M_{u,\mathcal{E}} \\
&\quad + \partial_{1,2}^2 F \sum_{u \in \mathbb{G}(t)} f'(X_t^u) M_{u,\mathcal{E}} + \partial_{2,2}^2 F \sum_{u \in \mathbb{G}(t)} f'(X_t^u) f'(X_t^v) M_{u,\mathcal{E}} \\
&\quad + \frac{1}{2} \partial_{2} F f''(X_t^\mathcal{E}) M_{\mathcal{E},\mathcal{E}} + \frac{1}{2} \partial_{2,2}^2 f'(X_t^\mathcal{E}) f'(X_t^v) M_{\mathcal{E},\mathcal{E}}.
\end{align*}
$$

Using the definition of the infinitesimal generator and rearranging the terms conclude the proof.

4.1.2. *Proof of Proposition 3.3.* The equality for the jumps part of the generator is Proposition 2.6 in [46]. For the continuous part $G$, we use (4.1) to get that

$$
\begin{align*}
G[\psi F](\mathfrak{e}, \bar{\nu}) &= \psi(X_t^\mathfrak{e}, \langle \nu_t, f \rangle) \sum_{u \in \mathbb{G}(t)} \partial_u F(X_t^\mathfrak{e}, \langle \nu_t, f \rangle) a_u(\bar{\nu}_t) \\
&\quad + F(X_t^\mathfrak{e}, \langle \nu_t, f \rangle) \sum_{u \in \mathbb{G}(t)} \partial_u \psi(X_t^\mathfrak{e}, \langle \nu_t, f \rangle) a_u(\bar{\nu}_t) \\
&\quad + \frac{\psi(X_t^\mathfrak{e}, \langle \nu_t, f \rangle)}{2} \sum_{u,w \in \mathbb{G}(t)} \partial_{u,w}^2 F(X_t^\mathfrak{e}, \langle \nu_t, f \rangle) (\boldsymbol{\sigma}\boldsymbol{\sigma}^T)_{u,w}(\bar{\nu}_t) \\
&\quad + \frac{F(X_t^\mathfrak{e}, \langle \nu_t, f \rangle)}{2} \sum_{u,w \in \mathbb{G}(t)} \partial_{u,w}^2 \psi(X_t^\mathfrak{e}, \langle \nu_t, f \rangle) (\boldsymbol{\sigma}\boldsymbol{\sigma}^T)_{u,w}(\bar{\nu}_t) \\
&\quad + \frac{1}{2} \sum_{u,w \in \mathbb{G}(t)} (\partial_w F \partial_u \psi + \partial_u F \partial_w \psi) (X_t^\mathfrak{e}, \langle \nu_t, f \rangle) (\boldsymbol{\sigma}\boldsymbol{\sigma}^T)_{u,w}(\bar{\nu}_t).
\end{align*}
$$

Using the symmetry of $\sigma\sigma^T$ we get that

$$
\frac{G[\psi F]}{\psi}(\mathfrak{e}, \bar{\nu}) - \frac{G\psi}{\psi}F(\mathfrak{e}, \bar{\nu}) = GF(\mathfrak{e}, \bar{\nu}) \\
+ \frac{1}{\psi(X_t^\mathfrak{e}, \langle \nu_t, f \rangle)} \sum_{u,w \in \mathbb{G}(t)} \partial_u \psi(X_t^\mathfrak{e}, \langle \nu_t, f \rangle) (\boldsymbol{\sigma}\boldsymbol{\sigma}^T)_{u,w}(\bar{\nu}_t) \partial_w F(X_t^\mathfrak{e}, \langle \nu_t, f \rangle).
$$

Notice that the last line can also be expressed in matricial form as

$$
\left(\frac{\boldsymbol{\sigma}\boldsymbol{\sigma}^T \nabla \psi}{\psi}\right)^T \nabla F.
$$

In the following we omit the dependencies in the functions to shorten the expressions. We slice
the sum by distinguishing the label $\mathfrak{e}$, and use the first order derivative expressions (4.2) to get

$$
\begin{align*}
(\boldsymbol{\sigma}\boldsymbol{\sigma}^T \nabla \psi)^T \nabla F &= \partial_2 F \partial_2 \psi \int_{\mathcal{X} \times \mathcal{U}} \int_{\mathcal{X} \times \mathcal{U}} f'(x) (\boldsymbol{\sigma}\boldsymbol{\sigma}^T)_{u,v} \phi'(y) \bar{\nu}(dy, dv) \bar{\nu}(dx, du) \\
&\quad + \partial_1 F \partial_1 \psi (\boldsymbol{\sigma}\boldsymbol{\sigma}^T)_{\mathcal{E},\mathcal{E}} + \partial_2 F \partial_1 \psi \int_{\mathcal{X} \times \mathcal{U}} f'(x) (\boldsymbol{\sigma}\boldsymbol{\sigma}^T)_{u,\mathcal{E}} \bar{\nu}(dx, du) \\
&\quad + \partial_1 F \partial_2 \psi \int_{\mathcal{X} \times \mathcal{U}} (\boldsymbol{\sigma}\boldsymbol{\sigma}^T)_{\mathcal{E},u} \phi'(x) \bar{\nu}(dx, du).
\end{align*}
$$

We end the proof by summing this last equation and expression (3.3) of the generator $G$.
4.1.3. *Proof of Theorem 3.4.* We obtain Theorem 3.4 with a modification from the proof of Theorem 2.2 in [46], to take into account the diffusive trait dynamics between jumps.

We denote $(U_k, k \ge 0)$ the sequence of $\mathcal{U}$-valued random variables giving the label of the branching individuals at the jumps times $(T_k, k \ge 0)$. Let $(N_k, k \ge 0)$ be the sequence of $\mathbb{N}$-valued random variables giving the number of children at each branching event and we denote for brevity $A_k := (U_k, N_k)$ for all $k \ge 0$. At the $k$-th branching time $T_k$, we denote $\mathcal{Y}_k$ the $\mathcal{X}^{N_k}$-valued random variable giving the vector of offspring traits. Finally we introduce for all $k \ge 0$, $\mathcal{V}_k := (T_k, A_k, \mathcal{Y}_k)$. We similarly define $((\hat{T}_k, \hat{U}_k, \hat{N}_k, \hat{\mathcal{Y}}_k), k \ge 0)$, the sequence of jumps times, labels of the branching individual, number of children and trait of these children at birth in the spinal construction. Notice that the distribution of number of children and traits depend on whether the branching individual is the spinal one or not. We will also use for all $k \ge 0$, $\hat{\mathcal{V}}_k := (\hat{T}_k, \hat{A}_k, \hat{\mathcal{Y}}_k)$ where $\hat{A}_k := (\hat{U}_k, \hat{N}_k)$. At time $s \in [\hat{T}_{k-1}, \hat{T}_k)$, the label of the spinal individual is denoted $E_k$ and its trait $\hat{\mathcal{Y}}_k$. For a given initial population $\bar{z} = \sum_{i=1}^n \delta_{(i,x^i)} \in \bar{\mathbb{V}}$, we use by convention $U_0 = \emptyset$, $N_0 = n$, $\mathcal{Y}_0 = (x^i, 1 \le i \le n)$ almost surely. The same convention holds for the spine process. For all $0 \le k$, we introduce the associated filtrations

$$ \mathcal{F}_k = \sigma(\mathcal{V}_i, 0 \le i \le k), \quad \text{and} \quad \hat{\mathcal{F}}_k = \sigma(E_k, (\hat{\mathcal{V}}_i, 0 \le i \le k)). $$

The core of the proof of Theorem 2.2 in [46], lies in the following Girsanov-type result

**Lemma 4.1** (Lemma 1.5.1 in [46]).

We denote by $\hat{\nu}$ the spine process. For any $k > 0$, $\bar{z} \in \bar{\mathbb{V}}$ and all $a = ((u_i, n_i), 0 \le i \le k)$, let $F$ be a measurable function on $\Pi_{i=1}^k (\mathbb{R}_+ \times \mathcal{U} \times \mathbb{N} \times \mathcal{X}^{n_i})$. For any $\mathfrak{c} \in \hat{\mathbb{G}}(\hat{T}_k)$,

$$ \mathbb{E}_{\bar{z}} \left[ F(\mathcal{V}_i, 0 \le i \le k) \prod_{i=0}^{k} \mathbb{1}_{\{A_i=a_i\}} \right] = \langle z, \psi(\cdot, z, 0) \rangle \\ \times \mathbb{E}_{\bar{z}} \left[ \xi_k \left( E_k, (\hat{\mathcal{V}}_i, 0 \le i \le k) \right) F(\hat{\mathcal{V}}_i, 0 \le i \le k) \mathbb{1}_{\{E_k=e\}} \prod_{i=0}^{k} \mathbb{1}_{\{\hat{A}_i=a_i\}} \right], $$

$$ \text{where } \xi_k (E_k, (\hat{\mathcal{V}}_i, 0 \le i \le k)) := \frac{1}{\psi(\hat{Y}_{\hat{T}_k}, \hat{\nu}_{\hat{T}_k}, \hat{T}_k)} \exp \left( \int_0^{\hat{T}_k} \frac{G\psi(\hat{Y}_s, \hat{\nu}_s, s)}{\psi(\hat{Y}_s, \hat{\nu}_s, s)} ds \right). $$

We will prove this Lemma in our extended framework- with diffusive trait dynamics- and the proof of Theorem 3.4 follows the same proof as in [46]. Lemma 4.1 will be proved by induction on the number of branching events $k$ as in [46], and the core Equation (1.5.10) is shown using Lemma 1.5.2 in [46] related to the fundamental theorem of calculus. However Lemma 1.5.2 in [46] does not hold anymore in this framework because of the stochastic dynamics of the traits between branching events. We derive an equivalent result that allows us to conclude the proof of Lemma 4.1 in the same way. We recall that the set $\mathcal{D}$, defined in (3.2), is composed of twice continuously differentiable and continuous functions $F_f$ such that for all $(\mathfrak{c}, \bar{\nu}) \in \bar{\mathbb{W}}$, $F_f(\mathfrak{c}, \bar{\nu}) := F(x_\mathfrak{c}, \langle \nu, f \rangle)$. Without loss of generality on the set $\mathcal{D}$, to simplify the notations and for more coherence between our notations and those in [46], we assume that $\psi_\phi$ is a function of the trait of the spine and the marginal population and will denote without distinction

$$ \psi_\phi(\mathfrak{c}, \bar{\nu}_t) = \psi_\phi(\hat{Y}_t, \hat{\nu}_t). $$

**Lemma 4.2.** If Assumption 1 holds, between two successive jumps, for all $t \in [\hat{T}_{k-1}, \hat{T}_k)$, for every measurable function $H$ on the set of trajectories on $[\hat{T}_{k-1}, \hat{T}_k)$, we get that

$$ \mathbb{E} \left[ H(\nu_s, \hat{T}_{k-1} \le s \le t) \mid \hat{\mathcal{F}}_k, \{\nu_{\hat{T}_{k-1}} = \hat{\nu}_{\hat{T}_{k-1}}\} \right] = \\ \mathbb{E} \left[ \frac{\psi_\phi(\hat{Y}_{\hat{T}_{k-1}}, \hat{\nu}_{\hat{T}_{k-1}})}{\psi_\phi(\hat{Y}_t, \hat{\nu}_t)} \exp \left\{ \int_{\hat{T}_{k-1}}^t \frac{G\psi_\phi(\hat{Y}_s, \hat{\nu}_s) ds}{\psi_\phi(\hat{Y}_s, \hat{\nu}_s) ds} \right\} H(\hat{\nu}_s, \hat{T}_{k-1} \le s \le t) \mid \hat{\mathcal{F}}_k \right]. $$
Notice that we did not consider time inhomogeneities in this framework as it would have complicated the notations, and our motivating example does not depend on time. Extending the results to such a non-homogeneous context only necessitates adding a time derivative in the expression of G and in the generator of the spine process, and the rest of the proof still holds.

*Proof.* Let $k \ge 1$ and $t \in [\hat{T}_{k-1}, \hat{T}_k)$. Using the expression of the generator $\hat{G}$ of the trait dynamics between branching events in the spine process, introduced in (3.7), we got that the vector $\hat{\mathbf{X}}_t = (\hat{\mathbf{X}}_t^u, u \in \hat{\mathbb{G}}(t))$ containing the vector of the individuals in the spine process satisfies, for a spine of label $\mathfrak{e}$

$$d\hat{\mathbf{X}}_t = (\mathbf{a}(\bar{\nu}_t) + \frac{\boldsymbol{\sigma}\boldsymbol{\sigma}^T \nabla\psi_\phi}{\psi_\phi}(\mathfrak{e}, \bar{\nu}_t))dt + \boldsymbol{\sigma}(\bar{\nu}_t) d\hat{\mathbf{B}}_t, \quad (4.3)$$

where the vector $\hat{\mathbf{B}}_t := (\hat{B}_t^u, u \in \hat{\mathbb{G}}(\hat{\nu}_t))$ is a Wiener process with i.i.d components under the law $\hat{\mathbb{P}}$ of the auxiliary process. We introduce $(\hat{\gamma}_s, s \ge 0)$ the $\hat{\mathcal{F}}_t$-adapted process via:

$$\hat{\gamma}_s := - \frac{\boldsymbol{\sigma}^T \nabla \psi_\phi}{\psi_\phi} (\mathfrak{e}, \bar{\nu}_t).$$

The SDE (4.3) can be written:

$$d\hat{\mathbf{X}}_t = \mathbf{a}(\bar{\nu}_t) dt + \boldsymbol{\sigma}(\bar{\nu}_t) d(\hat{\mathbf{B}}_t - \int_{\hat{T}_{k-1}}^t \hat{\gamma}_s ds). \quad (4.4)$$

Let $\mathbb{P}_{k-1}$ and $\hat{\mathbb{P}}_{k-1}$ denote respectively the laws of the original and the spinal branching processes conditioned on the filtration $\hat{\mathcal{F}}_{k-1}$. The Girsanov-Cameron-Martin theorem (see [35, p.72]), applied to (4.4) ensures that if the local $\hat{\mathcal{F}}_t$-martingale $\hat{L}$ defined by

$$\hat{L}_t := \exp \left\{ \int_{\hat{T}_{k-1}}^t \hat{\gamma}_s \cdot d\hat{\mathbf{B}}_s - \frac{1}{2} \int_{\hat{T}_{k-1}}^t \|\hat{\gamma}_s\|_2^2 ds \right\} \quad (4.5)$$

is a martingale, thus $(\hat{\mathbf{B}}_t - \int_{\hat{T}_{k-1}}^t \hat{\gamma}_s ds)$ is a Wiener process under the law $\mathbb{P}_{k-1}$ of the original branching process, and $\hat{L}$ is the Radon-Nykodym derivative such that

$$\left. \frac{d\hat{\mathbb{P}}_{k-1}}{d\mathbb{P}_{k-1}} \right|_{\hat{\mathcal{F}}_t} = \hat{L}_t. \quad (4.6)$$

Note that a classical sufficient condition for $\hat{L}$ to be a martingale is Novikov's condition:

$$\mathbb{E} \left[ \exp \left\{ \frac{1}{2} \int_{\hat{T}_{k-1}}^{t} \| \hat{\gamma}_s \|_2^2 ds \right\} \right] < \infty. \quad (4.7)$$

Assumption 1 ensures that the process $\hat{\gamma}$ is almost surely bounded, and that (4.7) holds. However this sufficient condition is not necessary to ensure this finite moment hypothesis, and a more detailed analysis on the process dynamics might give a less restrictive set of hypotheses. Now we will establish that the martingale $\hat{L}$ satisfies:

$$\hat{L}_t = \frac{\psi_\phi(\hat{Y}_{\hat{T}_{k-1}}, \hat{\nu}_{\hat{T}_{k-1}})}{\psi_\phi(\hat{Y}_t, \hat{\nu}_t)} \exp \left\{ \int_{\hat{T}_{k-1}}^{t} \frac{G\psi_\phi(\hat{Y}_s, \hat{\nu}_s) ds}{\psi_\phi(\hat{Y}_s, \hat{\nu}_s)} \right\}. \quad (4.8)$$

Using the vectorial formalism introduced earlier and applying Ito's formula we get that

$$\ln \psi_\phi (\mathfrak{e}, \bar{\nu}_t) = \ln \psi_\phi (\mathfrak{e}, \bar{\nu}_{\hat{T}_{k-1}}) + \int_{\hat{T}_{k-1}}^{t} \hat{G} \ln \psi_\phi (\mathfrak{e}, \bar{\nu}_s) ds + \int_{\hat{T}_{k-1}}^{t} \nabla \ln \psi_\phi^T (\mathfrak{e}, \bar{\nu}_s) \boldsymbol{\sigma} (\bar{\nu}_s) d\hat{\mathbf{B}}_s$$
Using the expression of $\hat{G}$ introduced in Proposition 3.3 we get that

$$
\begin{aligned}
\ln \psi_{\phi}(\mathbf{e}, \bar{\nu}_t) &= \ln \psi_{\phi}(\mathbf{e}, \bar{\nu}_{\hat{T}_{k-1}}) + \int_{\hat{T}_{k-1}}^{t} G \ln \psi_{\phi}(\mathbf{e}, \bar{\nu}_s) \, ds \\
&\quad + \int_{\hat{T}_{k-1}}^{t} \frac{\nabla \psi_{\phi}^{T}}{\psi_{\phi}} \boldsymbol{\sigma} \boldsymbol{\sigma}^{T} \frac{\nabla \psi_{\phi}}{\psi_{\phi}}(\mathbf{e}, \bar{\nu}_{s}) \, ds + \int_{\hat{T}_{k-1}}^{t} \frac{\nabla \psi_{\phi}^{T}}{\psi_{\phi}}(\mathbf{e}, \bar{\nu}_{s}) \boldsymbol{\sigma}(\bar{\nu}_{s}) \, d\hat{\mathbf{B}}_{s} \\
&= \ln \psi_{\phi}(\mathbf{e}, \bar{\nu}_{\hat{T}_{k-1}}) + \int_{\hat{T}_{k-1}}^{t} G \ln \psi_{\phi}(\mathbf{e}, \bar{\nu}_{s}) \, ds + \int_{\hat{T}_{k-1}}^{t} \|\hat{\gamma}_{s}\|_{2}^{2} \, ds - \int_{\hat{T}_{k-1}}^{t} \hat{\gamma}_{s} \, d\hat{\mathbf{B}}_{s}.
\end{aligned}
$$

Thus taking the exponential in the previous equation, and using (4.5), we get that

$$
\hat{L}_t = \frac{\psi_{\phi}(\mathbf{e}, \bar{\nu}_{\hat{T}_{k-1}})}{\psi_{\phi}(\mathbf{e}, \bar{\nu}_t)} \exp \left\{ \int_{\hat{T}_{k-1}}^{t} G \ln \psi_{\phi}(\mathbf{e}, \bar{\nu}_s) + \frac{1}{2} \|\hat{\gamma}_s\|_2^2 ds \right\}.
$$

We now detail the expression of the generator $G$ applied to the function $\ln \psi_{\phi}$. Using (3.3) to the function $F_f = \ln \circ \psi_{\phi}$, we have

$$
\begin{aligned}
G \ln \psi_{\phi}(\mathbf{e}, \bar{\nu}) ={}& \frac{\partial_1 \psi}{\psi_{\phi}}(x_{\mathbf{e}}, \langle \nu, \phi \rangle) \mathbf{a}_{\mathbf{e}}(\bar{\nu}) + \frac{1}{2} \left[ \frac{\partial_{1,1}^2 \psi}{\psi_{\phi}^2} - \left(\frac{\partial_1 \psi}{\psi_{\phi}}\right)^2 \right] (x_{\mathbf{e}}, \langle \nu, \phi \rangle) (\boldsymbol{\sigma} \boldsymbol{\sigma}^T)_{\mathbf{e}, \mathbf{e}}(\bar{\nu}) \\
& + \left[ \frac{\partial_{1,2}^2 \psi}{\psi_{\phi}^2} - \frac{\partial_1 \psi \partial_2 \psi}{\psi_{\phi}^2} \right] (x_{\mathbf{e}}, \langle \nu, \phi \rangle) \int_{\mathcal{X} \times \mathcal{U}} (\boldsymbol{\sigma} \boldsymbol{\sigma}^T)_{u, \mathbf{e}}(\bar{\nu}) \phi'(x) \bar{\nu}(dx, du) \\
& + \frac{\partial_2 \psi}{\psi_{\phi}}(x_{\mathbf{e}}, \langle \nu, \phi \rangle) \int_{\mathcal{X} \times \mathcal{U}} \mathbf{a}_{u}(\bar{\nu}) \phi'(x) + \frac{(\boldsymbol{\sigma} \boldsymbol{\sigma}^T)_{u, u}}{2} (\bar{\nu}) \phi''(x) \bar{\nu}(dx, du) \\
& + \frac{1}{2} \left[ \frac{\partial_{2,2}^2 \psi}{\psi_{\phi}^2} - \left(\frac{\partial_2 \psi}{\psi_{\phi}}\right)^2 \right] (x_{\mathbf{e}}, \langle \nu, \phi \rangle) \times \int_{\mathcal{X} \times \mathcal{U}} \int_{\mathcal{X} \times \mathcal{U}} (\boldsymbol{\sigma} \boldsymbol{\sigma}^T)_{u, v}(\bar{\nu}) \phi'(x) \phi'(y) \bar{\nu}(dx, du) \bar{\nu}(dy, dv).
\end{aligned}
$$

First we recognize the expression of $G\psi_{\phi}$ in this expression

$$
\begin{aligned}
G \ln \psi_{\phi}(\mathbf{e}, \bar{\nu}) ={}& \frac{G \psi_{\phi}}{\psi_{\phi}}(\mathbf{e}, \bar{\nu}) - \frac{1}{2} \left(\frac{\partial_1 \psi}{\psi_{\phi}}\right)^2 (x_{\mathbf{e}}, \langle \nu, \phi \rangle) (\boldsymbol{\sigma} \boldsymbol{\sigma}^T)_{\mathbf{e}, \mathbf{e}}(\bar{\nu}) \\
& - \frac{\partial_1 \psi \partial_2 \psi}{\psi_{\phi}^2} (x_{\mathbf{e}}, \langle \nu, \phi \rangle) \int_{\mathcal{X} \times \mathcal{U}} (\boldsymbol{\sigma} \boldsymbol{\sigma}^T)_{u, \mathbf{e}}(\bar{\nu}) \phi'(x) \bar{\nu}(dx, du) \\
& - \frac{1}{2} \left(\frac{\partial_2 \psi}{\psi_{\phi}}\right)^2 (x_{\mathbf{e}}, \langle \nu, \phi \rangle) \times \int_{\mathcal{X} \times \mathcal{U}} \int_{\mathcal{X} \times \mathcal{U}} (\boldsymbol{\sigma} \boldsymbol{\sigma}^T)_{u, v}(\bar{\nu}) \phi'(x) \phi'(y) \bar{\nu}(dx, du) \bar{\nu}(dy, dv),
\end{aligned}
$$

then we recognize the expression of $\nabla\psi_{\phi}^T \boldsymbol{\sigma} \boldsymbol{\sigma}^T \nabla\psi_{\phi}$:

$$
G \ln \psi_{\phi}(\mathbf{e}, \bar{\nu}) = \frac{G \psi_{\phi}}{\psi_{\phi}}(\mathbf{e}, \bar{\nu}) - \frac{1}{2} \frac{1}{\psi_{\phi}^2} \nabla \psi_{\phi}^T \boldsymbol{\sigma} \boldsymbol{\sigma}^T \nabla \psi_{\phi}(\mathbf{e}, \bar{\nu}).
$$

Then, using the definition of $\gamma$ and applying again Ito's formula on the previous equation yields

$$
G \psi_{\phi}(\mathbf{e}, \bar{\nu}_t) = G \ln \psi_{\phi}(\mathbf{e}, \bar{\nu}_t) + \frac{1}{2} \|\hat{\gamma}_t\|_2^2.
$$

This last equation gives (4.8), and the change of measure (4.6) concludes the proof. $\square$

## 4.2. Proofs of Section 2.
This section is devoted to the proofs concerning the colonial population model at the heart of the article.

4.2.1. *Proof of Theorem 2.1.* We introduce the sets $E_0 := \mathcal{U} \times \mathbb{R}_+$ and $E_2 := \mathcal{U} \times \mathbb{R}_+ \times [0, 1]$. We introduce the Poisson Point measures $Q_2$ and $Q_0$ respectively on the sets $[0, t] \times E_2$ and $[0, t] \times E_0$ and of respective intensity measures $d\mathit{sn}(du)d\mathit{rq}(d\theta)$ and $d\mathit{sn}(du)dr$, where $n(du)$ denotes the counting measure on $\mathcal{U}$ and $q$ is the probability measure giving the sharing of resources at a splitting event. The measure-valued process $(\nu_t, t \ge 0)$ following the branching mechanisms
previously introduced and where the colonies traits follow the SDE (1.1), is a weak-solution of
the following SDE, for all function $F \in C^2(\mathbb{R}_+, \mathbb{R})$:

$$
\begin{align*}
\langle \nu_t, F \rangle = \langle \nu_0, F \rangle + M_F(t) &+ \int_0^T \int_{\mathbb{R}_+} \left( axF'(x) + \frac{\sigma^2 x^2}{2} F''(x)\nu_s(dx) \right) ds \\
&+ \int_{[0,t] \times E_2} \mathbb{1}_{\left\{ u \in \mathbb{G}_s, r \le \lambda \frac{X_{s-}^u}{\langle \nu_{s-}, I_d \rangle} \right\}} \left( F(\theta X_{s-}^u) + F((1-\theta)X_{s-}^u) - F(X_{s-}^u) \right) Q_2(ds, du, dr, d\theta) \\
&- \int_{[0,t] \times E_0} \mathbb{1}_{\{u \in \mathbb{G}_s\}} \mathbb{1}_{\{r \le \mu\}} F(X_{s-}^u) Q_0(ds, du, dr).
\end{align*}
$$

The martingale $M_F(t)$ is given by

$$
M_F(t) := \int_0^T \dots \int_0^T \sum_{u \in \mathbb{G}(t)} \partial_x F(X_s^u) \frac{\sigma}{\sqrt{1+\delta^2}} X_s^u (dB_s^u + \delta dW_s).
$$

The first line is related to the traits dynamics of the colonies between the branching events. The second line gives the fission events and the extinction events are encompassed in the last line. Notice that this process does not explode in finite time as the branching rates are bounded and the mass is conserved at birth, see e.g. [42].

Without loss of generality, we rearrange the labels in the lexicographic order to fit the vectorial formalism introduced in Section 3. Thus we have

$$
\mathbf{X}_t := (X_t^u, u \in \mathbb{G}(t)) = (X_t^1, \dots, X_t^{N_t}) \quad \text{and} \quad \mathbf{a}(\bar{\nu}_t) := (aX_t^1, \dots, aX_t^{N_t}).
$$

The vector of Brownian processes is then $\mathbf{B}_t := (B_t^1, \dots, B_t^{N_t}, W_t)$ and the diffusion matrix

$$
\boldsymbol{\sigma}(\bar{\nu}_t) := \begin{bmatrix} \frac{\sigma}{\sqrt{1+\delta^2}} X_t^1 & 0 & \cdots & 0 & \frac{\sigma}{\sqrt{1+\delta^2}} \delta X_t^1 \\ 0 & \frac{\sigma}{\sqrt{1+\delta^2}} X_t^2 & \cdots & \vdots & \frac{\sigma}{\sqrt{1+\delta^2}} \delta X_t^2 \\ \vdots & \ddots & \ddots & 0 & \vdots \\ 0 & \cdots & 0 & \frac{\sigma}{\sqrt{1+\delta^2}} X_t^{N_t} & \frac{\sigma}{\sqrt{1+\delta^2}} \delta X_t^{N_t} \end{bmatrix}. \quad (4.9)
$$

With this formalism, the dynamics of the traits between jumps satisfy Equation (3.1). From the dynamics describing the colonial process in Section 3.2, the operator $\mathcal{G}$ introduced in (3.5) is defined for all measurable function $\psi_\phi(\mathbf{e}, \bar{\nu}) = \psi(x_\mathbf{e})$ by

$$
\mathcal{G}\psi_{\phi}(\mathbf{e}, \bar{\nu}) := \psi'(x_{\mathbf{e}})ax_{\mathbf{e}} + \frac{\lambda x_{\mathbf{e}}}{\langle \nu, Id \rangle} \int_{X^n} [\psi(\theta x_{\mathbf{e}}) + \psi((1-\theta)x_{\mathbf{e}}) - \psi(x_{\mathbf{e}})] q(d\theta) - \mu\psi(x_{\mathbf{e}}).
$$

Then, as the fission event is mass conservative, the bias function $\psi(x_\epsilon) = x_\epsilon$ is an eigenfunction
of the operator $\mathcal{G}$ for the eigenvalue $(a - \mu)$. To construct the spine process $((\hat{Y}_t, \hat{\nu}_t), t \ge 0)$, we
use the non-negative function $\psi$ defined for all $x \in \mathbb{R}_+^*$ as $\psi(x) = x$. Using (4.9), we derive the
expression of the additional drift term involved in the SDE that describes the evolution of the
traits in the $\psi$-spine process between two branching events

$$
\frac{\boldsymbol{\sigma} \boldsymbol{\sigma}^T (\bar{\nu}_t) \nabla \psi (E_t, \bar{\nu}_t)}{\hat{Y}_t} = \left( \sigma^2 \hat{Y}_t, \left( \sigma^2 \frac{\delta^2}{1 + \delta^2} \hat{X}_t^u, u \in \hat{\mathbb{G}}(t) \setminus \{E_t\} \right) \right).
$$

From this expression, we obtain the SDEs (2.1) and (2.2). The rates of jumps in the spine
process are given by the expression of the generator of the process given by Proposition 3.3,
applied with the chosen function $\psi$. We obtain that the branching rates in the spine process are
the same as the ones in the colonial process. The choice of the new spinal individual follows the
law (2.3) and (2.4). The death rates in the spine process are the same as the ones in the colonial
process for the non spinal individuals ($\mu$), and the spine does not extinct. Finally, applying
Theorem 3.4 for $U_t$ an individual picked uniformly at random in $\mathbb{G}(t)$ and for all measurable
function $F$ on the set of càdlàg trajectories on $\mathbb{R}_+^* \times \mathbb{U}$, we obtain Theorem 2.1.
4.2.2. Study of the spinal process. In order to derive results on the original population process, we first need to establish some technical lemmas on the spinal process. The first step is to use a set of autonomous processes whose law will well describe the law of the set of processes (($\hat{N}_t$, $\hat{R}_t$, $\hat{Y}_t$, ($X_t^u$, $u \in \hat{G}_t \setminus \{E_t\}$)), $t \ge 0$). To state the following result, we introduce the fraction of resources contained in the spine $\hat{Z}_t$ as well as in an individual $u \in \hat{G}_t \setminus \{E_t\}$, $\hat{Z}_t^u$ via

$$ \hat{Z}_t := \hat{Y}_t / \hat{R}_t \quad \text{and} \quad \hat{Z}_t^u := \hat{X}_t^u / \hat{R}_t. $$

**Lemma 4.3.** The process $(\mathcal{X}_t, t \ge 0) := (($\hat{N}_t$, $\hat{Z}_t$, ($\hat{Z}_t^u$, $u \in \hat{G}_t \setminus \{E_t\}$)), $t \ge 0$)

(1) is an autonomous and regenerative process

(2) does not depend on the parameter a

Recall that a regenerative process is a process $(\mathcal{X}_t, t \ge 0)$ such that there exists a renewal process with epochs $S_1, S_2, ...$ such that for all $n \ge 1$,

$$ ((\mathcal{X}_{t+S_n}, t \ge 0), (S_{m+n} - S_n, m \in \mathbb{N})) \stackrel{d}{=} ((\mathcal{X}_{t+S_1}, t \ge 0), (S_{m+1} - S_1, m \in \mathbb{N})) $$

and (($\mathcal{X}_{t+S_n}, t \ge 0$), ($S_{m+n} - S_n, m \in \mathbb{N}$)) is independent of (($\mathcal{X}_t, t \in [0, S_n]$), $S_1, S_2, ..., S_n$). Such a process admits long term properties which will be helpful in the study of the spinal process. Lemma 4.3 will be a key tool to study the long time behavior of the fraction of resources $\hat{Z}$ contained in the spinal individual.

*Proof.* From (3.3) with the bias function chosen in Section 2.2, we get that the trait $\hat{Y}_t$ of the spinal individual is the solution of the following SDE

$$ \begin{aligned} \hat{Y}_t = \hat{Y}_0 &+ \int_0^t (a + \sigma^2) \hat{Y}_s ds + \int_0^t \frac{\sigma}{\sqrt{1 + \delta^2}} \hat{Y}_s (d\hat{B}_s^* + \delta d\hat{W}_s) \quad (4.10) \\ &- \int_0^t \int_{\mathbb{R}_+} \int_0^1 \mathbb{1}_{\{r \le \frac{\lambda\theta\hat{Y}_{s-}}{\langle\hat{Y}_{s-}, I_d\rangle}\}} (1-\theta) \hat{Y}_s Q_1^*(ds, dr, d\theta) - \int_0^t \int_{\mathbb{R}_+} \int_0^1 \mathbb{1}_{\{r \le \frac{\lambda(1-\theta)\hat{Y}_{s-}}{\langle\hat{Y}_{s-}, I_d\rangle}\}} \theta \hat{Y}_s Q_2^*(ds, dr, d\theta) \end{aligned} $$

where $\hat{B}_s^* := \hat{B}_s^{E_s}$ is a Wiener process and $Q_1^*, Q_2^*$ are independent Poisson Point Measures on $\mathbb{R}_+ \times \mathbb{R}_+ \times [0, 1]$ of intensity $dsdrq(d\theta)$. For all individuals $u$ outside the spine $E_t$ of resources $\hat{Y}_t$, i.e. for all $u \in \hat{G}(t) \setminus \{E_t\}$, between the branching events, we denote $\hat{X}_t^u$ the trait of the individual of label $u$ living at time $t$ and with a slight abuse of notation it also denotes the trait of its ancestor at times $s < T_{\text{birth}}^u$. Thus, for all $u \in \hat{G}(t) \setminus \{E_t\}$ we can write the evolution of the trait of the 'lineage' of $u$ that survived until time $t$ as follows

$$ \begin{aligned} \hat{X}_t^u = \hat{X}_0^u &+ \int_0^t \left( a \hat{X}_s^u + \sigma^2 \frac{\delta^2}{1 + \delta^2} \hat{X}_t^u \right) ds + \int_0^t \frac{\sigma}{\sqrt{1 + \delta^2}} \hat{X}_s^u (d\hat{B}_s^u + \delta d\hat{W}_s) \\ &- \int_0^t \int_{\mathbb{R}_+} \int_0^1 \mathbb{1}_{\{r \le \lambda \frac{\hat{X}_{s-}^u}{\hat{R}_{s-}}\}} (1 - \theta) Q^u(ds, dr, d\theta) \quad (4.11) \end{aligned} $$

with $Q^u$ a Poisson Point Measure on $\mathbb{R}_+ \times \mathbb{R}_+ \times [0, 1]$ of intensity $dsdrq(d\theta)$. From (4.10) and (4.11) we obtain that the total resources $\hat{R}_t$ in the spinal population follows the SDE

$$ \begin{aligned} \hat{R}_t = \hat{R}_0 &+ \int_0^t \left( a \hat{R}_s + \frac{\sigma^2}{1 + \delta^2} (\hat{Y}_s + \delta^2 \hat{R}_s) \right) ds + \int_0^t \frac{\sigma}{\sqrt{1 + \delta^2}} \left( \sum_{u \in \hat{G}_s} \hat{X}_s^u d\hat{B}_s^u + \delta \hat{R}_s d\hat{W}_s \right) \\ &- \int_0^t \int_{\mathcal{U}} \int_{\mathbb{R}_+} \mathbb{1}_{\{u \in \hat{G}_{s-} \setminus \{E_t\}\}} \mathbb{1}_{\{r \le \mu\}} \hat{X}_{s-}^u \hat{Q}_0(ds, du, dr). \quad (4.12) \end{aligned} $$
We want to establish SDEs for $\hat{X}_t^u/\hat{R}_t$ for all individual $u \in \hat{\mathbb{G}}(t)$. To this end we need to differentiate an individual in (4.12). For all $\mathfrak{e} \in \mathbb{G}(t)$

$$ \hat{R}_t = \hat{R}_0 + \int_0^t \left( a\hat{R}_s + \frac{\sigma^2}{1+\delta^2} (\hat{Y}_s + \delta^2\hat{R}_s) \right) ds + \int_0^t \frac{\sigma}{\sqrt{1+\delta^2}} \left( \sum_{u \in \hat{\mathbb{G}}_s \setminus \{\mathfrak{e}\}} (\hat{X}_s^u)^2 \right)^{1/2} d\hat{B}_s \quad (4.13) \\ + \int_0^t \frac{\sigma}{\sqrt{1+\delta^2}} (\hat{X}_s^\mathfrak{e} d\hat{B}_s^\mathfrak{e} + \delta\hat{R}_s d\hat{W}_s) - \int_0^t \int_{\mathcal{U}} \int_{\mathbb{R}_+} \mathbb{1}_{\{u \in \hat{\mathbb{G}}_{s-} \setminus \{E_{s-}\}\}} \mathbb{1}_{\{r \le u\}} \hat{X}_{s-}^u \hat{Q}_0(ds, du, dr) $$

with $\hat{Q}_0$ a Poisson Point Measure on $\mathbb{R}_+ \times \mathcal{U} \times \mathbb{R}_+$ of intensity $dsn(du)dr$, and

$$ \left( \sum_{u \in \hat{\mathbb{G}}_s \setminus \{\mathfrak{e}\}} (\hat{X}_s^u)^2 \right)^{1/2} \stackrel{d}{=} \sum_{u \in \hat{\mathbb{G}}_s \setminus \{\mathfrak{e}\}} \hat{X}_s^u d\hat{B}_s^u. $$

Finally the size $\hat{N}_t$ of the spinal population is given by:

$$ \hat{N}_t = \hat{N}_0 + \int_0^t \int_{\mathbb{R}_+} \mathbb{1}_{\{r \le \lambda\}} \hat{Q}_2'(ds, dr) - \int_0^t \int_{\mathbb{R}_+} \mathbb{1}_{\{r \le \mu(\hat{N}_{s-1})\}} \hat{Q}_0'(ds, dr). \quad (4.14) $$

Note that $\hat{Q}_0'$ is equal to the measure $\hat{Q}_0$, introduced in (4.13), integrated over the individuals. To prove the first point of Proposition 4.3, we will apply Itô formula to the function $\phi(y, r) := y/r$, using (4.10) and (4.13) with $\mathfrak{e} = E_t$ to obtain the SDE satisfied by the processes $\hat{Z}$. Recall that $\hat{B}^* = \hat{B}^{E_t}$ is the Wiener process associated with the spinal individual. We get

$$ d\hat{Z}_t = \left( (a + \sigma^2) \hat{Y}_t \partial_1 \phi + \left( a\hat{R}_s + \frac{\sigma^2}{1+\delta^2} \hat{Y}_s + \sigma^2 \frac{\delta^2}{1+\delta^2} \hat{R}_s \right) \partial_2 \phi \right) dt \\ + \frac{1}{2} \left( \Sigma_{1,1}^s \partial_{1,1}^2 \phi + 2\Sigma_{1,2}^s \partial_{1,2}^2 \phi + \Sigma_{2,2}^s \partial_{2,2}^2 \phi \right) dt + \frac{\sigma\delta}{\sqrt{1+\delta^2}} \left( \hat{Y}_t \partial_1 \phi + \hat{R}_t \partial_2 \phi \right) d\hat{W}_t \\ + \frac{\sigma}{\sqrt{1+\delta^2}} \left( \hat{Y}_t (\partial_1 \phi + \partial_2 \phi) d\hat{B}_t^* + \left( \sum_{u \in \hat{\mathbb{G}}_s \setminus \{\mathfrak{e}\}} (\hat{X}_s^u)^2 \right)^{1/2} \partial_2 \phi d\hat{B}_t \right) \\ - \int_{\mathbb{R}_+} \int_0^1 \mathbb{1}_{\{r \le \lambda \theta \frac{\hat{Y}_{t-}}{\hat{R}_{t-}}\}} (1-\theta) \frac{\hat{Y}_{t-}}{\hat{R}_{t-}} Q_1^*(dt, dr, d\theta) - \int_{\mathbb{R}_+} \int_0^1 \mathbb{1}_{\{r \le \lambda(1-\theta) \frac{\hat{Y}_{t-}}{\hat{R}_{t-}}\}} \theta \frac{\hat{Y}_{t-}}{\hat{R}_{t-}} Q_2^*(dt, dr, d\theta) \\ + \int_{\mathcal{U}} \int_{\mathbb{R}_+} \mathbb{1}_{\{u \in \hat{\mathbb{G}}_{t-} \setminus \{E_{t-}\}\}} \mathbb{1}_{\{r \le \mu\}} \frac{\hat{Y}_{t-}}{\hat{R}_{t-}} \frac{\hat{X}_{t-}^u}{\hat{R}_{t-} - \hat{X}_{t-}^u} \hat{Q}_0(dt, du, dr) $$

with

$$ \Sigma^s := \begin{pmatrix} \sigma^2 \hat{Y}_s^2 & \frac{\sigma^2}{1+\delta^2} \hat{Y}_s (\hat{Y}_s + \delta^2 \hat{R}_s) \\ \frac{\sigma^2}{1+\delta^2} \hat{Y}_s (\hat{Y}_s + \delta^2 \hat{R}_s) & \frac{\sigma^2}{1+\delta^2} \left( \sum_{u \in \hat{\mathbb{G}}_s} (\hat{X}_s^u)^2 + \delta^2 \hat{R}_s^2 \right) \end{pmatrix}, $$

and

$$ \partial_1 \phi = 1/\hat{R}_s, \quad \partial_2 \phi = -\hat{Y}_s/\hat{R}_s^2 \quad \partial_{1,1}^2 \phi = 0, \quad \partial_{1,2}^2 \phi = -1/\hat{R}_s^2, \quad \partial_{2,2}^2 \phi = 2\hat{Y}_s/\hat{R}_s^3. $$

From the previous computation, we obtain

$$ \hat{Z}_t = \hat{Z}_0 + \int_0^t \frac{\sigma^2 \hat{Z}_s}{1+\delta^2} \left( 1 - 2\hat{Z}_s + \sum_{u \in \hat{\mathbb{G}}_s} (\hat{Z}_s^u)^2 \right) ds + \int_0^t \frac{\sigma \hat{Z}_s}{\sqrt{1+\delta^2}} (1-\hat{Z}_s) d\hat{B}_s^* \quad (4.15) \\ - \int_0^t \frac{\sigma \hat{Z}_s}{\sqrt{1+\delta^2}} \left( \sum_{u \in \hat{\mathbb{G}}_s \setminus \{\mathfrak{e}\}} (\hat{X}_s^u)^2 \right)^{1/2} d\hat{B}_s - \int_0^t \int_{\mathbb{R}_+} \int_0^1 \mathbb{1}_{\{r \le \lambda \theta \hat{Z}_{s-}\}} (1-\theta) \hat{Z}_{s-} Q_1^*(ds, dr, d\theta) \\ - \int_0^t \int_{\mathbb{R}_+} \int_0^1 \mathbb{1}_{\{r \le \lambda(1-\theta) \hat{Z}_{s-}\}} \theta \hat{Z}_{s-} Q_2^*(ds, dr, d\theta) + \int_0^t \int_{\mathcal{U}} \int_{\mathbb{R}_+} \mathbb{1}_{\{u \in \hat{\mathbb{G}}_{s-} \setminus \{E_{s-}\}, r \le \mu\}} \frac{\hat{Z}_{s-} \hat{Z}_s^u}{1-\hat{Z}_{s-}^u} \hat{Q}_0(ds, du, dr). $$
Similarly we obtain for all $u \in \hat{\mathbb{G}}_{t-} \setminus \{E_t\}$

$$
\begin{align}
\hat{Z}_t^u = \hat{Z}_0 &+ \int_0^t \frac{\sigma^2 \hat{Z}_s^u}{1+\delta^2} \left( \sum_{u \in \hat{\mathbb{G}}_s} (\hat{Z}_s^u)^2 - \hat{Z}_s^u - \hat{Z}_s \right) ds + \int_0^t \frac{\sigma \hat{Z}_s^u}{\sqrt{1+\delta^2}} (1-\hat{Z}_s^u) d\hat{B}_s^u \tag{4.16} \\
& - \int_0^t \frac{\sigma \hat{Z}_s^u}{\sqrt{1+\delta^2}} \left( \sum_{v \in \hat{\mathbb{G}}_s \setminus \{u\}} (Z_s^v)^2 \right)^{1/2} dB_s - \int_0^t \int_{\mathbb{R}_+} \int_0^1 1_{\{r \le \lambda \theta \hat{Z}_{s-}^u\}} (1-\theta) \hat{Z}_{s-}^u Q_1^*(ds, dr, d\theta) - \nonumber \\
& \qquad \int_0^t \int_{\mathbb{R}_+} \int_0^1 1_{\{r \le \lambda(1-\theta) \hat{Z}_{s-}^u\}} \theta \hat{Z}_{s-}^u Q_2^*(ds, dr, d\theta) + \int_0^t \int_{\mathcal{U}} \int_{\mathbb{R}_+} 1_{\{v \in \hat{\mathbb{G}}_{s-} \setminus \{E_{s-}, u\}, r \le \mu\}} \frac{\hat{Z}_{s-}^u Z_s^v}{1-\hat{Z}_{s-}^v} \hat{Q}_0(ds, du, dr). \nonumber
\end{align}
$$

Finally, the cardinal of $\hat{\mathbb{G}}$ is given by $\hat{N}$. The latter is autonomous and satisfies SDE (4.14). Hence, from SDEs (4.14), (4.15) and (4.16) we obtain that $\mathcal{X}$ is an autonomous process. To prove that it is also a regenerative process, it is enough to notice that when $\hat{N}_t = 1$, then $\hat{Z}_t = 1$ and $\hat{\mathbb{G}}_t \setminus \{E_t\} = \emptyset$, hence in this case $\mathcal{X}_t = (1, 1)$. Moreover, $S_1 < \infty$ almost surely, and the successive durations $(S_2 - S_1, S_3 - S_2, ...)$ are non-lattice iid random variables with finite expectation. It ends the proof of point 1. of Proposition 4.3. To prove point 2. it is enough to notice that the parameter $a$ does not appear in Equations (4.14), (4.15) and (4.16). $\square$

A natural question is to determine how the parameters $(\lambda, \mu, \sigma, \delta)$ influence the long time behavior of the system. Obtaining quantitative bounds is delicate. We may however get some indications, stated in Lemma 4.4. Recall (2.13) for the definition of $\alpha$ and $\beta$.

**Lemma 4.4.** The random variable $\hat{Z}_t$ converges in law to a random variable $\hat{Z}_\infty$ as $t \to \infty$.

(1) For every bounded continuous function on $\mathbb{N}^* \times (0, 1] \times [0, 1]^{\mathbb{N}^*}$,

$$ \mathbb{E}[f(\hat{Z}_{\infty})] = \lim_{t \to \infty} \frac{1}{t} \int_0^t f(\hat{Z}_s) ds = \lim_{t \to \infty} \frac{1}{t} \int_0^t \mathbb{E}[f(\hat{Z}_s)] ds. \quad (4.17) $$

(2) Let us denote $M := \left(\sqrt{1 + 8\alpha\left(1 + \frac{8}{27}\beta\right)} - 1\right) / 4\alpha$.

$$ \mathbb{E}[\hat{Z}_{\infty}] \le \frac{2\beta - 1}{4(\alpha + \beta)} \left(1 + \mathrm{sign}(2\beta - 1)\sqrt{1 + \frac{8(\alpha + \beta)}{(2\beta - 1)^2}}\right) \wedge M < 1. $$

$$ (3) \qquad \mathbb{E}[\hat{Z}_{\infty}] \ge \left( \frac{\beta + \beta \frac{\mu}{\lambda} (1 - e^{-\lambda/\mu}) + 2}{2\beta - \frac{4\lambda}{\mu} \mathbb{E}[\Theta \ln \Theta] + 2} \right) \vee \frac{\mu}{\lambda} (1 - e^{-\lambda/\mu}). $$

Notice that the upper bound has a limit $(1/\sqrt{2\alpha + 1})$ as $2\beta \to 1$.

*Proof.* The first point is a direct application of Theorems 7.1.4 and 7.1.5 of [14]. It is a consequence of the regenerative property of the process $\mathcal{X}$, stated in Proposition 4.3. Let us now focus on point (2), which states an upper bound for the mean fraction of resources contained in the spine in large time. Let us introduce $K_\theta := \mathbb{E}[\Theta(1-\Theta)]$. From (4.15), we deduce the SDE

$$ \hat{Z}_t = \hat{Z}_0 + M_t + \int_0^t \hat{Z}_s \left[ \frac{\sigma^2}{1+\delta^2} \left( 1 - 2\hat{Z}_s + \sum_{u \in \hat{\mathbb{G}}_s} (\hat{Z}_s^u)^2 \right) + \mu \sum_{u \in \hat{\mathbb{G}}_s \setminus \{E_s\}} \frac{\hat{Z}_s^u}{1-\hat{Z}_s^u} - 2\hat{Z}_s \lambda K_\theta \right] ds, \quad (4.18) $$

where $M_t$ is a square integrable martingale. Notice that almost surely, for all $u \in \hat{\mathbb{G}}_s \setminus \{E_s\}$, $\hat{Z}_s^u \le 1 - \hat{Z}_s$ and the second sum in (4.18) can be upper bounded as follows

$$ \sum_{u \in \hat{\mathbb{G}}_{s-} \setminus \{E_t\}} \frac{\hat{Z}_s^u}{1 - \hat{Z}_s^u} \le \frac{1 - \hat{Z}_s}{\hat{Z}_s}. $$
Hence,

$$ \begin{aligned} \hat{Z}_t &\le \hat{Z}_0 + \int_0^t \hat{Z}_s \left[ 2 \frac{\sigma^2}{1+\delta^2} (1-\hat{Z}_s) + \mu \frac{1-\hat{Z}_s}{\hat{Z}_s} - 2\hat{Z}_s \lambda K_\theta \right] ds + M_t \\ &= \hat{Z}_0 + \mu \int_0^t \left( 1 + \hat{Z}_s \left[ 2\beta - 1 - 2(\alpha + \beta)\hat{Z}_s \right] \right) ds + M_t \end{aligned} $$

with the definition of $\alpha$ and $\beta$ in (2.13). From Jensen inequality, we obtain

$$ \frac{d\mathbb{E}[\hat{Z}_t]}{dt} \le \mu \left( 1 + \mathbb{E}[\hat{Z}_t] \left[ 2\beta - 1 - 2(\alpha + \beta)\mathbb{E}[\hat{Z}_t] \right] \right). $$

Solving this second order polynomial we obtain the first upper bound of Lemma 4.4 point (2). The second upper bound, $M$, is obtained using that:

$$ \left( \sum_{u \in \hat{\mathbb{G}}_s} \hat{X}_s^u \right)^2 \ge \sum_{u \in \hat{\mathbb{G}}_s} \left( \hat{X}_s^u \right)^2 + 2\hat{Y}_s \sum_{u \in \hat{\mathbb{G}}_s \setminus \{E_s\}} \hat{X}_s^u $$

that gives an upper bound of the sum of ratios

$$ \sum_{u \in \hat{\mathbb{G}}_s} (\hat{Z}_s^u)^2 \le 1 - 2\hat{Z}_s (1 - \hat{Z}_s). $$

Using this refined upper bound we have that

$$ \begin{aligned} \frac{d\mathbb{E}[\hat{Z}_t]}{dt} &\le \mu + (2\sigma^2 - \mu) \mathbb{E}[\hat{Z}_t] - 2(\sigma^2 + \lambda K_\theta)\mathbb{E}[\hat{Z}_t^2] - 2\sigma^2 \mathbb{E}[\hat{Z}_s^2 (1 - \hat{Z}_s)] \\ &\le \mu + (2\sigma^2 - \mu) \mathbb{E}[\hat{Z}_t] - 2\lambda K_\theta \mathbb{E}[\hat{Z}_t^2] - 2\sigma^2 \mathbb{E}[2\hat{Z}_s^2 - \hat{Z}_s^3] \\ &\le \mu + \frac{8}{27}\sigma^2 - \mu \mathbb{E}[\hat{Z}_t] - 2\lambda K_\theta \mathbb{E}[\hat{Z}_t^2] \end{aligned} $$

where we used that for $x \in [0, 1]$, $2x^2 - x^3 \ge x - b$ with a constant $b$ such that $b := \max_{[0,1]} x(x - 1)^2 = 4/27$. This entails that $\limsup_{t \to \infty} \mathbb{E}[\hat{Z}_t] \le M$ and concludes the proof of point (2).

Let us now focus on point (3). We aim at finding a lower bound for the mean fraction of resources contained in the spine in large time. From (4.15) and Itô formula we get

$$ \begin{aligned} \ln \frac{\hat{Z}_t}{\hat{Z}_0} = & \int_0^t \frac{\sigma^2}{2(1+\delta^2)} \left[ (1-\hat{Z}_s)^2 + \sum_{u \in \hat{\mathbb{G}}_s \setminus \{E_s\}} (\hat{Z}_s^u)^2 \right] ds \\ & + \int_0^t \int_{\mathbb{R}_+} \int_0^1 \mathbb{1}_{\{r \le \lambda \theta \hat{Z}_{s-}\}} \ln \theta Q_1^*(\mathrm{d}s, \mathrm{d}r, \mathrm{d}\theta) + \int_0^t \int_{\mathbb{R}_+} \int_0^1 \mathbb{1}_{\{r \le \lambda(1-\theta) \hat{Z}_{s-}\}} \ln(1-\theta) Q_2^*(\mathrm{d}s, \mathrm{d}r, \mathrm{d}\theta) \\ & - \int_0^t \int_{\mathcal{U}} \int_{\mathbb{R}_+} \mathbb{1}_{\{u \in \hat{\mathbb{G}}_{s-} \setminus \{E_{s-}\}\}} \mathbb{1}_{\{r \le \mu\}} \ln (1 - \hat{Z}_{s-}^u) \hat{Q}_0(\mathrm{d}s, \mathrm{d}u, \mathrm{d}r) + \mathrm{Mart}_t, \end{aligned} $$

where $Mart_t$ is the martingale part due to integrals with respect to Brownian motions.

Using that $\mathbb{E}[\hat{Z}_s^2] \le \mathbb{E}[\hat{Z}_s]$, the following inequalities

$$ \sum_{u \in \hat{\mathbb{G}}_s} (\hat{Z}_s^u)^2 \ge \frac{1}{\hat{N}_t}, \quad \sum_{u \in \hat{\mathbb{G}}_{s-} \setminus \{E_{s-}\}} -\ln(1-\hat{Z}_{s-}^u) \ge \sum_{u \in \hat{\mathbb{G}}_{s-} \setminus \{E_{s-}\}} \hat{Z}_{s-}^u = 1-\hat{Z}_{s-}, $$

and taking the expectation yield

$$ \begin{aligned} \mathbb{E}[\ln \hat{Z}_t] \ge \mathbb{E}[\ln \hat{Z}_0] + & \int_0^t \frac{1}{2(1+\delta^2)} \left( 1 - 2\mathbb{E}[\hat{Z}_s] + \mathbb{E}\left[\frac{1}{\hat{N}_t}\right] \right) ds \\ & + \int_0^t \lambda \mathbb{E}[\hat{Z}_s] \mathbb{E}[\Theta \ln \Theta + (1-\Theta) \ln(1-\Theta)] ds + \int_0^t \mu(1-\mathbb{E}[\hat{Z}_s]) ds. \end{aligned} $$
Recalling that the law of $\Theta$ is symmetrical, we obtain for any $0 \le t_0 \le t$,

$$ \mathbb{E}[\ln \hat{Z}_t] \ge \mathbb{E}[\ln \hat{Z}_{t_0}] + \int_{t_0}^{t} \left[ \left( \frac{1}{2} \frac{\sigma^2}{1+\delta^2} \left( 1 + \mathbb{E}\left[\frac{1}{\hat{N}_t}\right] \right) + \mu \right) - \left( \frac{\sigma^2}{1+\delta^2} - 2\lambda \mathbb{E}[\Theta \ln \Theta] + \mu \right) \mathbb{E}[\hat{Z}_s] \right] ds. $$

We can check that the unique stationary distribution for the random variable $\hat{N}$ writes

$$ \mathbb{P}(\hat{N}_{\infty} = k) = e^{-\lambda/\mu} \frac{(\lambda/\mu)^{k-1}}{(k-1)!}, \quad \text{for } k \ge 1. $$

In particular,

$$ \mathbb{E}[1/\hat{N}_t] \underset{t \to \infty}{\implies} \mathbb{E}[1/\hat{N}_\infty] = \frac{\mu}{\lambda} (1 - e^{-\lambda/\mu}). \quad (4.19) $$

Let $\varepsilon > 0$ and $t_0$ be such that

$$ \mathbb{E}[1/\hat{N}_t] \ge \frac{\mu}{\lambda} (1 - e^{-\lambda/\mu}) - \varepsilon, \quad \forall t \ge t_0. $$

Then using Jensen inequality, we obtain that for any $t \ge t_0$,

$$ \mathbb{E}[\hat{Z}_t] \ge \exp\left\{\mathbb{E}[\ln \hat{Z}_{t_0}]\right\} \exp\left\{\left[\frac{\sigma^2}{2(1+\delta^2)}\left(1 + (1-e^{-\lambda/\mu})\frac{\mu}{\lambda} - \varepsilon\right) + \mu\right](t-t_0)\right\} \\ \times \exp\left\{-\left(\frac{\sigma^2}{1+\delta^2} - 2\lambda\mathbb{E}[\Theta \ln \Theta] + \mu\right)\int_{t_0}^{t} \mathbb{E}[\hat{Z}_s]ds\right\}. \quad (4.20) $$

For the sake of readability, let us introduce the two real numbers:

$$ \mathfrak{a} := 1 + (1 - e^{-\lambda/\mu}) \frac{\mu}{\lambda} - \varepsilon + \frac{2\mu(1 + \delta^2)}{\sigma^2} \le 2 - \frac{(1 + \delta^2)}{\sigma^2} (4\lambda\mathbb{E}[\Theta \ln \Theta] + 2\mu) =: \mathfrak{b}, $$

as well as

$$ \mathfrak{c} := \mathfrak{b} \exp \left( \mathbb{E}[\ln \hat{Z}_{t_0}] \right) \quad \text{and} \quad \mathfrak{g}_t := \mathfrak{b} \int_{t_0}^{t} \mathbb{E}[\hat{Z}_s] ds. $$

Equation (4.20) reads

$$ \mathfrak{g}'_t \ge \mathfrak{c} \exp \left\{ \frac{\sigma^2}{2(1+\delta^2)} \mathfrak{a}(t-t_0) \right\} \exp \left\{ -\frac{\sigma^2}{2(1+\delta^2)} \mathfrak{g}_t \right\}. $$

Separating variables, we get

$$ \exp\left\{\frac{\sigma^2}{2(1+\delta^2)}\mathfrak{g}_t\right\} - 1 \ge \frac{\mathfrak{c}}{\mathfrak{a}}\left(\exp\left\{\frac{\sigma^2}{2(1+\delta^2)}\mathfrak{a}(t-t_0)\right\} - 1\right), $$

which finally leads to

$$ \liminf_{t \to \infty} \frac{1}{t} \int_{t_0}^{t} \mathbb{E}[\hat{Z}_s] ds \ge \frac{1 + (1 - e^{-\lambda/\mu}) \frac{\mu}{\lambda} + \frac{2\mu(1+\delta^2)}{\sigma^2}}{2 - \frac{4\lambda(1+\delta^2)}{\sigma^2} \mathbb{E}[\Theta \ln \Theta] + \frac{2\mu(1+\delta^2)}{\sigma^2}}. $$

The second lower bound is a direct consequence of the Many-to-One formula (2.6):

$$ \mathbb{E}_{\nu_0} \left[ 1_{\{\mathfrak{G}(t) \neq \emptyset\}} \sum_{u \in \mathfrak{G}(t)} \frac{(X_t^u)^2}{R_t} \right] = R_0 e^{(a-\mu)t} \mathbb{E}_{\nu_0} \left[ \frac{Y_t^2 / \hat{R}_t}{\hat{Y}_t} \right] = R_0 e^{(a-\mu)t} \mathbb{E}_{\nu_0} \left[ \hat{Z}_t \right] $$

and

$$ \mathbb{E}_{\nu_0} \left[ 1_{\{\mathfrak{G}(t) \neq \emptyset\}} \frac{R_t}{\hat{N}_t} \right] = \mathbb{E}_{\nu_0} \left[ 1_{\{\mathfrak{G}(t) \neq \emptyset\}} \sum_{u \in \mathfrak{G}(t)} \frac{X_t^u}{\hat{N}_t} \right] = R_0 e^{(a-\mu)t} \mathbb{E}_{\nu_0} \left[ \frac{\hat{Y}_t / \hat{N}_t}{\hat{Y}_t} \right] = R_0 e^{(a-\mu)t} \mathbb{E}_{\nu_0} \left[ \frac{1}{\hat{N}_t} \right]. $$

We conclude the proof thanks to the inequality $R_t/N_t \le \sum_{u \in \mathfrak{G}(t)} (X_t^u)^2 / R_t$. $\square$

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

The next Lemma explores the role of noise in the behavior of the spinal individual. It will appear as a key result in the proof of Proposition 2.6.

**Lemma 4.5.** Let $\mu$ and $\lambda$ be fixed. Then

$$ \lim_{\sigma^2/(1+\delta^2) \to \infty} \mathbb{E}[\hat{Z}_\infty] = 1. $$
Proof. From (4.17) and (4.18) we get that

$$ \lim_{t \to \infty} \left\{ \mathbb{E} \left[ \frac{\sigma^2}{1 + \delta^2} \hat{Z}_t \left( 1 - 2\hat{Z}_t + \sum_{u \in \hat{G}_t} (\hat{Z}_t^u)^2 \right) + \mu \hat{Z}_t \sum_{u \in \hat{G}_t \setminus \{E_t\}} \frac{\hat{Z}_t^u}{1 - \hat{Z}_t^u} - 2\hat{Z}_t^2 \lambda K_\theta \right] \right\} = 0. $$

Adding (4.17) and recalling that $K_\theta := \mathbb{E}[\Theta(1 - \Theta)]$, we deduce that

$$ \lim_{t \to \infty} \mathbb{E} \left[ \frac{\sigma^2 \hat{Z}_t}{1 + \delta^2} \sum_{u \in \hat{G}_t \setminus \{E_t\}} (\hat{Z}_t^u)^2 + \mu \hat{Z}_t \sum_{u \in \hat{G}_t \setminus \{E_t\}} \frac{\hat{Z}_t^u}{1 - \hat{Z}_t^u} \right] = \\ \mathbb{E} \left[ 2\lambda K_\theta \hat{Z}_\infty^2 + \frac{\sigma^2}{1 + \delta^2} (2\hat{Z}_\infty^2 - \hat{Z}_\infty - \hat{Z}_\infty^3) \right]. \quad (4.21) $$

Noticing that the function $x \mapsto 2x^2 - x - x^3$ is non-positive on $[0, 1]$ and using that $Z_t^u \in (0, 1)$ for any $u \in \hat{G}_t \setminus \{E_t\}$, we obtain

$$ \lim_{t \to \infty} \frac{\sigma^2}{1 + \delta^2} \mathbb{E} \left[ \hat{Z}_t \sum_{u \in \hat{G}_t \setminus \{E_t\}} (\hat{Z}_t^u)^2 \right] \le 2\lambda K_\theta \mathbb{E} \left[ \hat{Z}_\infty^2 \right], $$

and thus

$$ \lim_{t \to \infty} \mathbb{E} \left[ \hat{Z}_t \sum_{u \in \hat{G}_t \setminus \{E_t\}} (\hat{Z}_t^u)^2 \right] \le (1 + \delta^2) \frac{2\lambda K_\theta}{\sigma^2}. \quad (4.22) $$

Now using that

$$ \sum_{u \in \hat{G}_t \setminus \{E_t\}} (\hat{Z}_t^u)^2 = (1 - \hat{Z}_t)^2 \sum_{u \in \hat{G}_t \setminus \{E_t\}} \left(\frac{\hat{Z}_t^u}{1 - \hat{Z}_t}\right)^2 \ge \frac{(1 - \hat{Z}_t)^2}{\hat{N}_t} $$

we deduce, using Cauchy-Schwarz inequality

$$ \lim_{t \to \infty} \left( \mathbb{E} \left[ \hat{N}_t^{1/2} \hat{Z}_t^{1/2} \frac{(1 - \hat{Z}_t)}{\hat{N}_t^{1/2}} \right] \right)^2 \le \mathbb{E} \left[ \hat{N}_\infty \right] \lim_{t \to \infty} \mathbb{E} \left[ \hat{Z}_t \frac{(1 - \hat{Z}_t)^2}{\hat{N}_t} \right] \le \frac{\lambda + \mu}{\lambda} (1 + \delta^2) \frac{2\lambda K_\theta}{\sigma^2}. $$

At the end, we thus get

$$ \lim_{\sigma^2/(1+\delta^2) \to \infty} \mathbb{E} \left[ \hat{Z}_\infty^{1/2} (1 - \hat{Z}_\infty) \right] = 0 $$

which implies that the random variable $\hat{Z}_\infty^{1/2}(1 - \hat{Z}_\infty)$ converges in probability to 0 when $\sigma^2/(1+\delta^2)$ goes to infinity. Adding (4.22) implies that

$$ \lim_{t \to \infty} \mathbb{E} \left[ \sum_{u \in \hat{G}_t \setminus \{E_t\}} \hat{Z}_t^u \sum_{u \in \hat{G}_t \setminus \{E_t\}} (\hat{Z}_t^u)^2 \right] \le (1 + \delta^2) \frac{2\lambda K_\theta}{\sigma^2}. $$

Hence $\sum_{u \in \hat{G}_t \setminus \{E_t\}} (\hat{Z}_t^u)^2$ converges in probability to 0 when $\sigma^2/(1+\delta^2)$ goes to $\infty$, which allows us to conclude that the random variable $1 - \hat{Z}_\infty$ converges in probability to 0 when $\sigma^2$ goes to infinity. This ends the proof. $\square$

Now we can establish the results presented in Section 2.3.

4.2.3. *Proof of Proposition 2.3.* Notice that the process $N_t$ has the same law as an $M/M/\infty$ queue with initial state 1, until it reaches 0. The exact law of $N_t$ seems out of reach (we may find its Laplace transform in the literature, but it is not easily invertible), but we know its mean extinction time, which is enough for our purposes. If we denote by $T_0 := \inf\{t \ge 0, N_t = 0\}$, we derive from [49] Equation (2.7)

$$ \mathbb{E}_{\nu_0}[T_0] = \frac{1}{\lambda} \sum_{j \ge 1} \left(\frac{\lambda}{\mu}\right)^j \sum_{k=0}^{N_0-1} \frac{k!}{(j+k)!} \le N_0 \frac{e^{\lambda/\mu} - 1}{\lambda}. $$
The inequality is an equality if $\langle \nu_0, 1 \rangle = 1$. Applying Markov Inequality, we deduce

$$ \mathbb{P}(\mathbb{G}(t) \neq \emptyset) = \mathbb{P}(T_0 \geq t) \leq 1 \wedge \frac{\mathbb{E}_{\nu_0}[T_0]}{t} \leq 1 \wedge N_0 \frac{e^{\lambda/\mu} - 1}{\lambda t} \quad (4.23) $$

which goes to 0 when $t$ goes to infinity. A simple lower bound on $\mathbb{P}(\mathbb{G}(t) \neq \emptyset)$ can be derived from the Many-to-One formula (2.2)

$$ \mathbb{P}(\mathbb{G}(t) \neq \emptyset) = \mathbb{E}\left[1_{\mathbb{G}(t) \neq \emptyset} \frac{\sum_{u \in \mathbb{G}(t)} X_t^u}{R_t}\right] = R_0 e^{(a-\mu)t} \mathbb{E}\left[\frac{1}{\hat{R}_t}\right]. \quad (4.24) $$

Then, from Jensen inequality we have

$$ \mathbb{E}\left[\frac{1}{\hat{R}_t}\right] \geq \exp\left(-\mathbb{E}\left[\ln \hat{R}_t\right]\right). $$

Using Itô formula, we obtain the SDE for $\ln \hat{R}_t$ and taking the expectation we get

$$ \frac{d\mathbb{E}[\ln \hat{R}_t]}{dt} \leq a - \mu + \frac{\sigma^2}{2(1+\delta^2)} \left(\delta^2 - E\left[\frac{1}{\hat{N}_t}\right]\right) + \left(\frac{\sigma^2}{1+\delta^2} + \mu\right) \mathbb{E}[\hat{Z}_t]. \quad (4.25) $$

We now use the almost sure bounds $\mathbb{E}[\hat{Z}_t] \leq 1$ and $E[1/\hat{N}_t] \geq 0$. Then (4.25) becomes

$$ \mathbb{E}[\ln \hat{R}_t] \leq \mathbb{E}[\ln \hat{R}_0] + \left(a + \frac{\sigma^2}{1+\delta^2} \left(\frac{\delta^2}{2} + 1\right)\right) t. $$

Notice that the law of the population size $N_t$ only depends on the parameters $\lambda$ and $\mu$. We can thus inject the previous inequality in (4.24) and take $\sigma = 0$. This reads

$$ \mathbb{P}(\mathbb{G}(t) \neq \emptyset) \geq e^{-\mu t}. \quad (4.26) $$

From the branching dynamics of the colonial process, we obtain the differential equation

$$ d\mathbb{E}_{\nu_0}[N_t]/dt = \lambda \mathbb{P}(\mathbb{G}(t) \neq \emptyset) - \mu \mathbb{E}_{\nu_0}[N_t]. $$

Using the previously established bounds (4.23) and (4.26) in the previous equality we get that

$$ \lambda e^{-\mu t} - \mu \mathbb{E}_{\nu_0}[N_t] \leq \frac{d}{dt} \mathbb{E}_{\nu_0}[N_t] \leq \lambda \left(1 \wedge N_0 \frac{e^{\lambda/\mu} - 1}{\lambda t}\right) - \mu \mathbb{E}_{\nu_0}[N_t]. \quad (4.27) $$

Integrating the lower bound gives that

$$ (\lambda t + N_0) e^{-\mu t} \leq \mathbb{E}_{\nu_0}[N_t]. $$

For the upper bound, we integrate on $[0, t_1] \cup (t_1, t]$, where $t_1 := N_0(e^{\lambda/\mu} - 1)/\lambda$. Integrating on $[0, t_1]$, we have that

$$ \mathbb{E}_{\nu_0}[N_t] \leq \frac{\lambda}{\mu} (1 - e^{-\mu t}) + N_0 e^{-\mu t}. $$

Integrating on $(t_1, t]$, we have that

$$ \mathbb{E}_{\nu_0}[N_t] \leq N_0 \left(e^{\frac{\lambda}{\mu}} - 1\right) \left(E_i(\mu t) - E_i(\mu t_1)\right) e^{-\mu t} + n_{t_1} e^{-\mu(t-t_1)} $$

where $n_{t_1} := (1 - e^{-\mu t_1})\lambda/\mu + N_0 e^{-\mu t_1}$ and the function $E_i$ is the exponential integral function defined for all $x \in \mathbb{R}_+^*$ as

$$ E_i(x) := \int_{-\infty}^{x} \frac{e^t}{t} dt. $$

We conclude the proof by noticing that $E_i(x)e^{-x} \sim 1/x$ as $x$ goes to infinity [1].
4.2.4. *Proof of Proposition 2.5.* First notice that the expression (2.9) of the mean resources can be obtained with classical means from the SDEs describing the dynamics of the colonies. However the result can be obtain as a direct application of the Many-to-One formula (2.6) with the function $f(x) = x$. Now we will establish boundaries on the variance of the total resources in the population. From the Many-to-One formula (2.6), with the function $f(x) = xR$ we have

$$ \mathbb{E}_{\nu_0}[R_t^2] = R_0 e^{(a-\mu)t} \mathbb{E}_{\nu_0}[\hat{R}_t]. \quad (4.28) $$

We can thus simply find bounds on $Var[R_t]$ by finding inequalities on the mean resources of the spinal process. From the SDE (4.13) we have

$$ \frac{d\mathbb{E}[\hat{R}_t]}{dt} = \mathbb{E}[\hat{R}_t] \left(a - \mu + \frac{\sigma^2 \delta^2}{1 + \delta^2}\right) + \mathbb{E}[\hat{Y}_t] \left(\mu + \frac{\sigma^2}{1 + \delta^2}\right). \quad (4.29) $$

We now take the expectation in (4.10) we have that

$$ \frac{d\mathbb{E}[\hat{Y}_t]}{dt} = \mathbb{E}\left[\hat{Y}_t \left(a + \sigma^2 - 2\lambda \mathbb{E}[\theta(1-\theta)]\frac{\hat{Y}_t}{\hat{R}_t}\right)\right]. \quad (4.30) $$

Then we use, in (4.30), the fact that $0 \le \hat{Z}_t \le 1$ a.s. to obtain the following bounds on $\mathbb{E}[\hat{Y}_t]$

$$ \mathbb{E}[\hat{Y}_t] (a + \sigma^2 - 2\lambda \mathbb{E}[\theta(1-\theta)]) \le \frac{d\mathbb{E}[\hat{Y}_t]}{dt} \le \mathbb{E}[\hat{Y}_t](a + \sigma^2). $$

Recall that $K_\theta := \mathbb{E}[\theta(1-\theta)]$. Integrating these inequalities and injecting into (4.29) yield

$$ \frac{d\mathbb{E}_{\nu_0}[\hat{R}_t]}{dt} \le \mathbb{E}[\hat{R}_t](a - \mu + \frac{\sigma^2}{1 + \delta^2}\delta^2) + \left(\mu + \frac{\sigma^2}{1 + \delta^2}\right) \mathbb{E}[\hat{Y}_0]e^{(a+\sigma^2)t}, $$

that integrates into

$$ \mathbb{E}[\hat{R}_t] \le e^{(a+\sigma^2)t} + (R_0 - 1)e^{(a-\mu+\frac{\sigma^2}{1+\delta^2}\delta^2)t}. $$

And

$$ \frac{d\mathbb{E}_{\nu_0}[\hat{R}_t]}{dt} \ge \mathbb{E}[\hat{R}_t] \left(a - \mu + \frac{\sigma^2 \delta^2}{1 + \delta^2}\right) + \left(\mu + \frac{\sigma^2}{1 + \delta^2}\right) \mathbb{E}[\hat{Y}_0]e^{(a+\sigma^2-2\lambda K_\theta)t}, $$

that integrates into

$$ \mathbb{E}[\hat{R}_t] \ge \frac{\mu + \frac{\sigma^2}{1+\delta^2}}{\mu + \frac{\sigma^2}{1+\delta^2} - 2\lambda K_\theta} \left( e^{(a+\sigma^2-2\lambda K_\theta)t} - e^{(a-\mu+\frac{\sigma^2}{1+\delta^2}\delta^2)t} \right) + R_0 e^{(a-\mu+\frac{\sigma^2}{1+\delta^2}\delta^2)t} $$

if $\mu+\sigma^2/(1+\delta^2) \ne 2\lambda K_\theta$, and into $\mathbb{E}[\hat{R}_t] \ge (2\lambda K_\theta t + R_0)e^{(a+\sigma^2-2\lambda K_\theta)t}$ if $\mu+\sigma^2/(1+\delta^2) = 2\lambda K_\theta$. We conclude the proof by using (4.28) and the definition of the variance.

4.2.5. *Proof of Proposition 2.6.* To obtain the quantity of resources owned by an individual picked uniformly at random in the population, we apply (2.5) to $F = x^u$. It yields

$$ \mathbb{E}_{\nu_0}[1_{\{\mathbb{G}(t) \ne \emptyset\}} X_t^{U_t}] = R_0 e^{(a-\mu)t} \mathbb{E}_{\nu_0}\left[\frac{1}{\hat{N}_t}\right] \sim R_0 e^{(a-\mu)t} \frac{\mu}{\lambda} \left(1 - e^{-\lambda/\mu}\right). $$

From the Many-to-One formula (2.6), with the function $f(x, \nu) := x1_{\{\frac{x}{\langle \nu, t \rangle} \ge \gamma\}}$, we get

$$ \mathbb{E}_{\nu_0}\left[\sum_{u \in \mathbb{G}(t)} X_t^u 1_{\left\{\frac{X_t^u}{\hat{R}_t} \ge \gamma\right\}}\right] = \mathbb{E}_{\nu_0}[R_t] \mathbb{P}_{\nu_0}\left(\frac{\hat{Y}_t}{\hat{R}_t} \ge \gamma\right). \quad (4.31) $$

Equation (2.10) is a direct consequence of point (1) of Lemma 4.5 and (2.11) is a direct consequence of point (2) of Lemma 4.4 with $\beta = 0$. Equation (2.12) is a consequence of the Many-to-One formula (2.6), with the function $f(x, \nu) := x^2/R$ and the inequalities

$$ \sum_{u \in \mathbb{G}(t)} \frac{(X_t^u)^2}{R_t} \le \sum_{u \in \mathbb{G}(t)} \sup_{u \in \mathbb{G}(t)} X_t^u \frac{X_t^u}{R_t} = \sup_{u \in \mathbb{G}(t)} X_t^u, $$

as well as point (3) of Lemma 4.4.
ACKNOWLEDGMENTS.

We thank J.-F. Delmas and S. C. Harris for their diligent proofreading and constructive com-
ments on this work. We also thank L. Coquille and A. Marguet for the interesting discussions
we had on this project.

This work is supported by the French National Research Agency in the framework of the
"France 2030" program (ANR-15-IDEX-0002) and by the LabEx PERSYVAL-Lab (ANR-11-
LABX-0025-01). We acknowledge partial support by the Chair "Modélisation Mathématique et
Biodiversité" of VEOLIA-Ecole Polytechnique-MNHN-F.X.

REFERENCES

[1] Milton Abramowitz and Irene A Stegun. *Handbook of mathematical functions with formulas, graphs, and mathematical tables*. Vol. 55. US Government printing office, 1968.

[2] Linda JS Allen. *An introduction to stochastic processes with applications to biology*. CRC press, 2010.

[3] Krishna B. Athreya. "Change of Measures for Markov Chains and the LlogL Theorem for Branching Processes". In: *Bernoulli* 6.2 (2000), pp. 323-338. ISSN: 13507265.

[4] Leticia Avilés. “Nomadic behaviour and colony fission in a cooperative spider: life history evolution at the level of the colony?” In: *Biological Journal of the Linnean Society* 70.2 (2000), pp. 325–339.

[5] Leticia Avilés and Paul Tufino. “Colony size and individual fitness in the social spider Anelosimus eximius”. In: *The American Naturalist* 152.3 (1998), pp. 403–418.

[6] Vincent Bansaye. “Spine for interacting populations and sampling”. In: *Bernoulli* 30.2 (2024), pp. 1555–1585.

[7] Vincent Bansaye, Juan Carlos Pardo Millan, and Charline Smadi. “On the extinction of continuous state branching processes with catastrophes”. In: *Electronic Journal of Probability* 18 none (2013), pp. 1–31.

[8] Vincent Bansaye, Juan Carlos Pardo, and Charline Smadi. “Extinction rate of continuous state branching processes in critical Lévy environments”. In: *ESAIM: Probability and Statistics* 25 (2021), pp. 346–375.

[9] Vincent Bansaye and Viet Chi Tran. “Branching Feller diffusion for cell division with parasite infection”. In: *ALEA : Latin American Journal of Probability and Mathematical Statistics* 8 (2011), pp. 95–127.

[10] Vincent Bansaye et al. “A non-conservative Harris ergodic theorem”. In: *Journal of the London Mathematical Society* 106.3 (2022), pp. 2459–2510.

[11] Mátyás Barczy, Zenghu Li, and Gyula Pap. “Stochastic differential equation with jumps for multi-type continuous state and continuous time branching processes with immigration”. In: *ALEA. Latin American Journal of Probability and Mathematical Statistics* 12.1 (2015), pp. 129–169.

[12] Oshrat Ben-Hamo et al. “Fission in a colonial marine invertebrate signifies unique life history strategies rather than being a demographic trait”. In: *Scientific Reports* 12.1 (2022), p. 15117.

[13] J. D. Biggins and A. E. Kyprianou. “Measure Change in Multitype Branching”. In: *Advances in Applied Probability* 36.2 (2004), pp. 544–581. ISSN: 00018678.

[14] Mogens Bladt and Bo Friis Nielsen. *Matrix-exponential distributions in applied probability*. Vol. 81. Springer, 2017.

[15] Michael B Bonsall and Alan Hastings. “Demographic and environmental stochasticity in predator–prey metapopulation dynamics”. In: *Journal of Animal Ecology* 73.6 (2004), pp. 1043–1055.

[16] Raphaël Boulay et al. “Intraspecific competition affects population size and resource alloca-
tion in an ant dispersing by colony fission”. In: *Ecology* 91.11 (2010), pp. 3312–3321.

[17] Nicolas Champagnat and Amaury Lambert. “Evolution of discrete populations and the canonical diffusion of adaptive dynamics”. In: *The Annals of Applied Probability* 17.1 (2007), pp. 102–155.
[18] Brigitte Chauvin and Alain Rouault. "K-P-P equation and supercritical branching Brownian motion in the subcritical speed area. Application to spatial trees". In: *Probability theory and related fields* 80.2 (1988), pp. 299-314.

[19] Ron Chen et al. "Resource limitation drives fission-fusion dynamics of group composition and size in a social bird". In: *Animal Behaviour* 191 (2022), pp. 15-32.

[20] Bertrand Cloez. "Limit theorems for some branching measure-valued processes". In: *Advances in Applied Probability* 49.2 (2017), pp. 549-580.

[21] Mary M Conner and Gary C White. "Effects of individual heterogeneity in estimating the persistence of small populations". In: *Natural Resource Modeling* 12.1 (1999), pp. 109-127.

[22] P Damos, Alexandros Rigas, and Matilda Savopoulou-Soultani. *Application of Markov Chains and Brownian motion models on insect ecology*. 2011.

[23] Pierre Del Moral. *Feynman-Kac Formulae: Genealogical and Interacting Particle Systems with Applications*. New York, NY: Springer New York, 2004, pp. 427-522. ISBN: 978-1-4684-9393-1.

[24] LM DeVantier, R Endean, et al. "Observations of colony fission following ledge formation in massive reef corals of the genus Porites." In: *Marine ecology progress series*. Oldendorf 58.1 (1989), pp. 191-195.

[25] Emily S. Durkin et al. "Population differences in the aggregation and collective foraging behavior of fragmented social spider colonies". In: *Ethology* 129.4-5 (2023), pp. 224-231.

[26] János Engländer. "Branching diffusions, superdiffusions and random media". In: *Probability Surveys* 4 (2007), pp. 303-364.

[27] János Engländer, Simon C. Harris, and Andreas E. Kyprianou. "Strong Law of Large Numbers for branching diffusions". In: *Annales de l'Institut Henri Poincaré, Probabilités et Statistiques* 46.1 (2010), pp. 279-298.

[28] Basile Finland et al. "Solitary foundation or colony fission in ants: an intraspecific study shows that worker presence and number increase colony foundation success". In: *Oecologia* 204.3 (2024), pp. 517-527.

[29] Wolfgang Frey and Harald Kürschner. "Asexual reproduction, habitat colonization and habitat maintenance in bryophytes". In: *Flora-Morphology, Distribution, Functional Ecol-ogy of Plants* 206.3 (2011), pp. 173-184.

[30] Y. Git, J. W. Harris, and S. C. Harris. "Exponential growth rates in a typed branching diffusion". In: *The Annals of Applied Probability* 17.2 (2007), pp. 609-653.

[31] Robert Hardy and Simon C. Harris. "A conceptual approach to a path result for branching Brownian motion". In: *Stochastic Processes and their Applications* 116.12 (2006), pp. 1992-2013. ISSN: 0304-4149.

[32] Robert Hardy and Simon C. Harris. "A Spine Approach to Branching Diffusions with Applications to Lp-Convergence of Martingales". In: *Séminaire de Probabilités XLII* 1979 (2009). Ed. by Catherine Donati-Martin et al., pp. 281-330.

[33] Simon C Harris and David Williams. "Large-deviations and martingales for a typed branching diffusion, 1". In: *Asterique* 236 (1996), pp. 133-154.

[34] Simon C. Harris. "Travelling-waves for the FKPP equation via probabilistic arguments". In: *Proceedings of the Royal Society of Edinburgh: Section A Mathematics* 129.3 (1999), pp. 503-517.

[35] Monique Jeanblanc, Marc Yor, and Marc Chesney. *Mathematical methods for financial markets*. Springer Science & Business Media, 2009.

[36] A.E. Kyprianou. "Travelling wave solutions to the K-P-P equation: alternatives to Simon Harris' probabilistic analysis". In: *Annales de l'Institut Henri Poincare (B) Probability and Statistics* 40.1 (2004), pp. 53-72. ISSN: 0246-0203.

[37] Russell Lande. "Demographic stochasticity and Allee effect on a scale with isotropic noise". In: *Oikos* (1998), pp. 353-358.

[38] Howard R Lasker. "Clonal propagation and population dynamics of a gorgonian coral". In: *Ecology* 71.4 (1990), pp. 1578-1589.
[39] Nicholas A Linacre and Michael J Keough. "Demographic effects of fragmentation history in modular organisms: illustrated using the bryozoan Mucropetraliella ellerii (MacGillivray)". In: *Ecological modelling* 170.1 (2003), pp. 61–71.

[40] Diego Lirman. "A simulation model of the population dynamics of the branching coral Acropora palmata Effects of storm intensity and frequency". In: *Ecological modelling* 161.3 (2003), pp. 169–182.

[41] Yael D Lubin. "Population studies of two colonial orb-weaving spiders". In: *Zoological Journal of the Linnean Society* 70.3 (1980), pp. 265–287.

[42] Aline Marguet. "Uniform sampling in a structured branching population". In: *Bernoulli* 25.4A (2019), pp. 2649–2695.

[43] Aline Marguet and Charline Smadi. "Parasite infection in a cell population: role of the partitioning kernel". In: *Journal of Mathematical Biology* 89.1 (2024), p. 14.

[44] Henry P McKean. "Application of brownian motion to the equation of kolmogorov, petrovskii, piskunov". In: *Communications on pure and applied mathematics* 28.3 (1975), pp. 323–331.

[45] Charles Medous. *Github repository Colonial-model*. 2024.

[46] Charles Medous. "Spinal constructions for continuous type-space branching processes with interactions". In: *Electronic Journal of Probability* 29 (2024), pp. 1–46.

[47] Mayuko Nakamaru et al. "Ecological conditions favoring budding in colonial organisms under environmental disturbance". In: *PLoS One* 9.3 (2014), e91210.

[48] Christian Peeters and Mathieu Molet. "Colonial reproduction and life histories". In: *Ant ecology* (2010), pp. 159–176.

[49] Frank Roijers, Michel Mandjes, and Hans Van Den Berg. "Analysis of congestion periods of an M/M/∞-queue". In: *Performance Evaluation* 64.7-8 (2007), pp. 737–754.

[50] Karin Sauer et al. "The biofilm life cycle: expanding the conceptual model of biofilm formation". In: *Nature Reviews Microbiology* 20.10 (2022), pp. 608–620.

[51] JZ Shik. "Ant colony size and the scaling of reproductive effort". In: *Functional Ecology* (2008), pp. 674–681.

[52] Peter J. Smith, Mansoor Shafi, and Hongsheng Gao. "Quick simulation: A review of importance sampling techniques in communications systems". In: *IEEE journal on selected areas in communications* 15.4 (1997), pp. 597–613.

[53] Walter R Tschinkel, Eldridge S Adams, and Thomas Macom. "Territory area and colony size in the fire ant Solenopsis invicta". In: *Journal of animal ecology* (1995), pp. 473–480.

[54] Christopher C Wilmers, Eric Post, and Alan Hastings. "A perfect storm: the combined effects on population fluctuations of autocorrelated environmental noise, age structure, and density dependence". In: *The American Naturalist* 169.5 (2007), pp. 673–683.

[55] Paul M Yoshioka. "Are large colonies a" key factor" in the dynamics of gorgonian populations?" In: *Revista de biología tropical* 46.S5 (1998), pp. 137–143.

SYLVAIN BILLIARD, UNIV. LILLE, CNRS, UMR 8198 – EVO-ECO-PALEO, F-59000 LILLE, FRANCE.

Email address: sylvain.billiard@univ-lille.fr

CHARLES MEDOUS, UNIV. GRENOBLE ALPES, CNRS, IF, 38000 GRENOBLE, FRANCE AND UNIV. GRENOBLE ALPES, INRIA, 38000 GRENOBLE, FRANCE

Email address: charles.medous@univ-grenoble-alpes.fr

CHARLINE SMADI, UNIV. GRENOBLE ALPES, INRAE, LESSEM, 38000 GRENOBLE, FRANCE AND UNIV. GRENOBLE ALPES, CNRS, INSTITUT FOURIER, 38000 GRENOBLE, FRANCE

Email address: charline.smadi@inrae.fr