# Biology-Informed inverse problems for insect pests detection using pheromone sensors

Thibault Malou, Nicolas Parisey, Katarzyna Adamczyk-Chauvat, Elisabeta Vergu, Béatrice Laroche, Paul-André Calatayud, Philippe Lucas, Simon Labarthe

► To cite this version:

Thibault Malou, Nicolas Parisey, Katarzyna Adamczyk-Chauvat, Elisabeta Vergu, Béatrice Laroche, et al.. Biology-Informed inverse problems for insect pests detection using pheromone sensors. Peer Community Journal, 2025, 5, pp.e19. 10.24072/pcjournal.520. hal-04982331

HAL Id: hal-04982331

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

Submitted on 7 Mar 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>
# Peer Community Journal

Section: Mathematical & Computational Biology

Research article

Published
2025-02-07

**Cite as**

Thibault Malou, Nicolas Parisey, Katarzyna Adamczyk-Chauvat, Elisabeta Vergu, Béatrice Laroche, Paul-André Calatayud, Philippe Lucas and Simon Labarthe (2025)

*Biology-Informed inverse problems for insect pests detection using pheromone sensors*, Peer Community Journal, 5: e19.

**Correspondence**

simon.labarthe@inrae.fr

**Peer-review**

Peer reviewed and
recommended by
PCI Mathematical &
Computational Biology,
https://doi.org/10.24072/pci-
mcb.100313

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

BY

This article is licensed
under the Creative Commons
Attribution 4.0 License.

## Biology-Informed inverse problems for insect pests detection using pheromone sensors

Thibault Malou<sup>1,1</sup>, Nicolas Parisey<sup>1,2</sup>, Katarzyna Adamczyk-Chauvat<sup>1,1</sup>, Elisabeta Vergu<sup>1,1</sup>, Béatrice Laroche<sup>1,3</sup>, Paul-André Calatayud<sup>1,4,5</sup>, Philippe Lucas<sup>1,6</sup>, and Simon Labarthe<sup>1,7,8</sup>

Volume 5 (2025), article e19

https://doi.org/10.24072/pcjournal.520

## Abstract

Most insects have the ability to modify the odor landscape in order to communicate with their conspecies during key phases of their life cycle such as reproduction. They release pheromones in their nearby environment, volatile compounds that are detected by insects of the same species with exceptional specificity and sensitivity. Efficient pheromone detection is then an interesting lever for insect pest management in a precision agroecological culture context. A precise and early detection of pests using pheromone sensors offers a strategy for pest management before infestation. In this paper, we develop a biology-informed inverse problem framework that leverages temporal signals from a pheromone sensor network to build insect presence maps. Prior biological knowledge is introduced in the inverse problem by the mean of a specific penalty, using population dynamics PDE residuals. We benchmark the biological-informed penalty with other regularization terms such as Tikhonov, LASSO or composite penalties in a simplified toy model. We use classical comparison criteria, such as target reconstruction error, or Jaccard distance on pest presence-absence. But we also use more task-specific criteria such as the number of informative sensors during inference. Finally, the inverse problem is solved in a realistic context of pest infestation in an agricultural landscape by the fall armyworm (Spodoptera frugiperda).

<footnote n="1">¹Université Paris-Saclay, INRAE, MalAGE, 78350, Jouy-en-Josas, France, ²INRAE, Institute of Genetics, Environment and Plant Protection (IGEPP)—Joint Research Unit 1349), Le Rheu, France, ³Université Paris-Saclay, Inria, 91120, Palaiseau, France, ⁴Institut de Diversité, Ecologie et Evolution du Vivant (IDEEV), Université Paris-Saclay, CNRS, IRD, UMR Evolution, Génome, Comportement et Ecologie (EGCE), Gif-sur-Yvette, France, ⁵International Center of Insect Physiology and Ecology (icipe), Nairobi, Kenya, ⁶Institute of Ecology and Environmental Sciences of Paris (iEES-Paris - Joint Research Unit 1392 - INRAE, CNRS, IRD, Sorbonne Univ., UPEC, Univ. Paris Cité), INRAE, 78000, Versailles, France, ⁷Univ. Bordeaux, INRAE, BIOGECO, 33610, Cestas, France, ⁸Univ. Bordeaux, Inria, INRAE, 33400, Talence, France</footnote>
# 1. Introduction

Insect pests represent a major threat for agricultural systems, causing direct or indirect damages to crops (Oerke, 2006), with an even increasing burden in a warming climate (Deutsch et al., 2018). World trade, travel and climate changes increase the arrival of alien species that challenge existing pest control management for native insects (EFSA PLH Panel et al., 2018). Intensified insecticide use is not sustainable, due notably to the negative impact of pesticides on biodiversity, environment and human health. Alternative solutions must be developed: a promising perspective is the development of precision agriculture which aims in particular to leverage sensors and information systems to reduce the environmental impact of agricultural practices (Gebbers et al., 2010). In that perspective, early detection of pests before infestation settlement can help designing low-pesticide crop management and epidemiological surveillance policies.

Insect pheromones are a good target for early detection of pests (Conchou et al., 2019). Indeed, insects pheromones are species-specific, guaranteeing a good specificity of any pheromone-based detection. Pheromone specificity has long been leveraged for insect detection with the mean of pheromone traps: pheromone diffusion modifies the odor landscape, attracting the targeted insects toward traps where they can be counted (McNeil, 1991). However, if pheromone traps allow for presence/absence detections, they do not allow to accurately locate the insects in their natural habitat. Recently, new advances in sensor developments paved the way towards field-based tools to explore *in situ* pheromone plumes. These sensors have different operating principles: there are insect antenna-based biosensors (Martinez et al., 2014; Myrick et al., 2011; Pawson et al., 2020), olfactory receptor-based biosensors (Khadka et al., 2019), OBP-based biosensors (Lu et al., 2014), and e-nose technologies based on conducting polymers, metal oxide semiconductor and other devices (Moitra et al., 2016; Negri et al., 2008; Steffens et al., 2014; Wehrenfennig et al., 2012). The use of pheromone sensors could allow to directly measure the odor landscape and track back pheromone plumes towards emitters, in order to detect pests at an early stage of invasion.

From a mathematical point of view, the problem of locating the emission source of a molecule knowing measures of the emission plume belongs to the class of inverse problems (Isakov, 2017), the direct problem of which is a plume propagation model. Atmospheric dispersion, plume propagation, and plume-related inverse problems have a large literature, in particular for environmental applications, such as pollution detection (Lushi et al., 2010; Stockie, 2011). Inverse problems, and more specifically data assimilation problems, are most often formulated as an optimization problem aiming to minimize the discrepancies between the direct model outputs and in-field observations by fitting the model parameters (Isakov, 2017). In the specific case where insects are the emitters, locating the emission source amounts to searching for continuous density maps of pests, which induces a high dimensional parameter space to be explored during the inverse problem. Variational approaches are well suited for such high-dimensional optimization problem, for fast computation of gradients of the objective function (Bocquet et al., 2014). Inverse problems are most of the time ill-posed and must be regularized by additional terms in the loss function that constrain the parameter exploration, e.g. with LASSO regularization term (Adam et al., 2016). Other regularization can be searched for, that could take into account prior biological knowledge incorporated into models, like in physic or biology-informed machine learning approaches (Yazdani et al., 2020), through weak or strong constraints (Trémolet, 2007).

In conventional agricultural practices, farmers often use pesticide prophylactically: knowing the potential occurrence of an infestation, and lacking efficient ways to monitor pest populations, they apply pesticides as a preventive measure, even before a pest is observed. Consequently,
sensors capable of measuring the spatio-temporal distribution of pest populations could poten-
tially allow the farmers to reduce the frequency, dose and area of treatment by spraying when,
where and to the extent that the pest is detected. If we interpret this prophylactic strategy in
the context of a classification problem, it is as if they were classifying all their fields as "infected"
whatever their actual state of infestation, thus minimizing their miss rate, i.e. a non detected
true presence of pest. In so doing, they also maximize their false discovery rate. In order to mini-
mize pesticide use, we would like to minimize the false discovery rate, so that pesticides are only
used where and when a detected pest is actually present in the field. But, according to a classical
trade-off in classification problems, reducing the false discovery rate comes with an increase of
the miss rate, so that the detection problem could be reformulated as follows: find the minimal
false discovery rate that keeps the miss rate under an acceptable risk, from the point of view of
crop protection.

The main objective of this paper is to introduce and test a data assimilation framework,
namely the biology-informed data assimilation, which seeks to combine a pheromone propa-
gation model and pheromone sensors data, like in classical data assimilation, but also biological
information. In particular, we specifically aim to introduce biology-informed regularization terms,
on which the biology-informed data assimilation is based, in order to assess how adding biologi-
cal priors improves the accuracy of the data assimilation. Several comparison criteria will be used.
In particular, we will introduce a specific criteria designed to compare the number of active sen-
sors involved during the inference process.

## 2. Material and methods

A sketch of the global methodological framework can be found in Fig. 1. The method contains three blocks that will be detailed hereafter:

1) the direct problem that aims to model the pheromone propagation in the atmosphere over an agricultural landscape, taking into account environmental data such as land use or wind field, see Sec. 2.1;

2) the pheromone sensors that provide data of the observed variable and an observation operator that maps the state variable of the direct model to the observed variable, see Sec. 2.3.1;

3) the corresponding inverse problem that aims to infer pheromone emissions by the insects and enables the construction of pests risk maps, knowing the data and the pheromone propagation model, see Sec. 2.2.

The proposed framework is illustrated with a toy case and a more realistic case, set up in an
agricultural landscape, see Sec. 2.3.2. In the different test cases, the datasets are synthetic and
constructed from a spatio-temporal sampling of pheromone concentration and a realistic model
of the observation process.

## 2.1. 2D reaction-convection-diffusion model for the pheromone propagation

One of the most classical model to simulate propagation of chemical species in the atmo-
sphere is a 3D Chemical-Transport Model (CTM), that has been extensively reviewed (Leelössy
et al., 2014; Sportisse, 2007). The CTM is obtained from the reactive Navier-Stokes equations
by using the dilution hypothesis. This hypothesis is valid when the chemical species of interest
is present in very small amounts. In this case, the compound is assumed to be chemically pas-
sive, which results in a decoupling between atmospheric dynamics and compound concentration
dynamics. As a consequence, dynamical field such as the wind or the diffusion coefficients can
# Biology-informed Data Assimilation

## Classical Data Assimilation

<figure><img src="image_3.png" /><figcaption>**Figure 1** – Graphical abstract: summary of the proposed Biology-Informed Data Assimilation method for the pest insects localization. A classical Data Assimilation (DA) framework is first built, combining a pheromone propagation model (blue rectangle), pheromone sensors (red rectangle) into an inverse problem (green rectangle). Additionally to this classical DA framework, biological information can be supplemented through specific regularization terms of the DA problem (purple rectangle). These regularization terms are the main insights of this study.</figcaption></figure>
be computed separately, offline, and given as parameters of the CTM. It results that this (linear and offline-coupled) CTM is much less complex than the full (non-linear) reactive Navier-Stokes equations, or than a CTM coupled online with a fine scale Navier-Stokes model (Leelössy et al., 2014), while being a reliable approximation: it is able to model the main physical phenomena at play during the transportation of pheromones.

In the present context, only the lower layer of the atmosphere, where insects and sensors are located, is of interest. Hence, for sake of computation load and data availability in real applications, a 2D depth-integrated CTM is used to model the total concentration of pheromone over the depth of the layer $c(x, y, t)$ [g m$^{-2}$] on a domain $\Omega \in \mathbb{R}^2$ and during a time window $[0; T]$. This depth-integrated CTM results in a 2D reaction-convection-diffusion model, that can be derived from the classical 3D CTM, see Sec. A.1 in the annex for details, and reads:

$$ (1) \qquad \frac{\partial c}{\partial t} - \nabla \cdot (\mathbf{K} \nabla c) + \nabla \cdot (\vec{u} c) + \tau_{loss} c = s \quad \forall (x, y) \in \Omega, \forall t \in [0; T] $$

The coefficient $\mathbf{K}(x, y, t)$ [m$^2$ s$^{-1}$] is the diffusion tensor. This tensor models the turbulent diffusion at an unresolved scale (classically, scale smaller than the scale of the numerical grid) using the K theory, see e.g. Panofsky et al. (1984) and Sportisse (2007). As mentioned above, in CTM, the wind velocity field $\vec{u}(x, y, t) = (u(x, y, t), v(x, y, t))$ [m s$^{-1}$] is given as parameter, where $(u, v)$ are the velocity field components along x and y. For example, $\vec{u}$ can be computed using a large (meso-scale) meteorological model and then, taken as input data in the present CTM, see an example in Sec. 2.3.2 and on Fig. 3a. The loss coefficient $\tau_{loss}(x, y, t)$ [s$^{-1}$] models the loss of pheromone due to vertical phenomena, including pheromone settlement on the ground and vertical diffusion spreading the pheromone out of the layer of interest. Hence, $\tau_{loss}$ depends on the land cover and on the stability of the atmosphere, see estimations of this coefficient in Sec. A.1.3. Other features can also be included in this parameter, such as humidity-dependent deposition. In the following, the solution of the CTM model is also referred to as the pheromone plume.

In Eq. (1), the source term $s(x, y, t)$ [g m$^{-2}$ s$^{-1}$] is the quantity of pheromone emitted by the insects located in $(x, y, t)$ per unit of area and time. We highlight that this source term $s$ is the quantity to infer to deduce the localization of the insects. Namely, this quantity can be decomposed in

$$ (2) \qquad s(x, y, t) = q(t)p(x, y, t) $$

with $q(t)$ [g s$^{-1}$] the pheromone emission rate, and $p(x, y, t)$ [m$^{-2}$] the insect density on the ground, so that there is a direct relationship between $s$ and the local presence of pests. The emission rate $q$ can incorporate biological information, such as insect-specific emission patterns.

The 2D model (1) is closed by the following initial and boundary conditions:

* a null initial condition $c(x, y, t = 0) = 0 \quad \forall (x, y) \in \Omega$, as the initial time is set before any pheromone emission in a context of primary invasion,

* a null diffusive flux $\mathbf{K} \nabla c \cdot \vec{n} = 0 \quad \forall (x, y) \in \partial \Omega$,

* a null convective influx $\vec{u} c \cdot \vec{n} = 0 \quad \forall (x, y) \in \partial \Omega \cap \{(x, y) | \vec{u}(x, y, t) \cdot \vec{n} < 0\} \forall t \in [0; T]$ with $\vec{n}$ the outgoing normal vector.

The pheromone propagation model (1) is solved using an implicit upwind finite volume scheme on a Cartesian grid. This scheme has been implemented in the git repository of the PheroSensor project, see Malou et al. (2024c) and Sec. 5. In this repository, one can found convergence tests of this numerical scheme. The companion code (Malou et al., 2024a) is also available to reproduce the different computations and figures of this study.
## 2.2. Variational Data Assimilation for insect localization

We aim to approximate insect risk maps by inferring the spatial distribution of emitted pheromones. This inverse problem consists in estimating the source term $s$ of the pheromone propagation model (1), using time-series of pheromone concentration. To this purpose, we use Variational Data Assimilation (VDA) techniques.

### 2.2.1. Minimization problem.
VDA consists in finding the value of the control variable, $s(x, y, t)$ for all $(x, y, t) \in \Omega \times [0; T]$, that enables the prediction made by the direct model, here the pheromone propagation model (1), to be close to some observations, here the dataset of pheromone concentration time-series. Therefore, the inference problem is defined as a minimization problem of a cost function $j$ with respect to the control variable (Bocquet et al., 2014) that can be split in an observation and a regularization term:

$$ (3) \qquad \left\{ \begin{array}{l} \text{Find } s_a(x, y, t) \text{ such that:} \\ s_a(x, y, t) = \arg \min_{s(x,y,t)} j(s) \text{ with } j(s) = j_{\text{obs}}(s) + j_{\text{reg}}(s) \end{array} \right. $$

The term $j_{\text{obs}}$ aims to minimize the gap between model predictions and data (Bocquet et al., 2014; Le Dimet et al., 1986; Sasaki, 1958). It reads

$$ (4) \qquad j_{\text{obs}}(s) = \|m(c(s)) - m^{\text{obs}}\|_{\mathbb{R}^{-1}}^2 $$

where $m^{\text{obs}}$ are the observations, $c(s)$ is the direct model output obtained for a given control $s$, $m(c)$ the observation operator that maps the direct model state variable towards the observation space and $\mathbb{R}$ is the covariance operator of the observation error. Note that the equivalent statistical notation of eq. (4) reads

$$ m^{\text{obs}} = m(c(s)) + \varepsilon $$

where $\varepsilon$ is a gaussian noise with zero mean and covariance $\mathbb{R}$.

The regularization term $j_{\text{reg}}$ is commonly added in order to ensure the well-posedness of the optimization problem, by improving the convexity of $j$ (Bocquet et al., 2014). Classical regularization terms are based on a rough background estimate of the control variable, either by a model with an insufficient accuracy or by coarse data.

The novelty and the results of this study lies on the introduction of Biology-Informed (BI) regularization terms that enable to drive the DA method toward an optimum $s_a$ consistent with the biological knowledge introduced in $j_{\text{reg}}$, leading to a Biology-Informed Data Assimilation framework. In general terms, the regularization term reads:

$$ (5) \qquad j_{\text{reg}} = \sum_{i \in \mathcal{L}} \alpha_{\text{reg},i} j_{\text{reg},i} $$

with $\mathcal{L}$ the list of considered BI regularization terms $j_{\text{reg},i}$ and $\alpha_{\text{reg},i}$ the associated weight coefficient. This allows for composite regularization terms that take account of several biological inputs. The definition of these BI regularization terms is introduced in the Results in Sec. 3.1.

### 2.2.2. Adjoint model for gradient computation.
To solve the minimization problem (3), we use a gradient-based method, suitable for our convex and differentiable optimization framework. Thus, at each descent step, the gradient $\nabla_{s,j}$ must be computed and reads:

$$ \nabla_{s,j}(s) = \nabla_{s} j_{\text{obs}}(s) + \nabla_{s} j_{\text{reg}}(s) $$
Note that we denote $\nabla_s$ the gradient with respect to the control $s$ to avoid confusion with $\nabla$ the classical gradient with respect to the space variables. Since $s$ is defined in a high-dimensional space and $c(s)$ involves the computation of a PDE model, the numerical computation of this gradient is challenging, in particular for its observation term $\nabla_{s,j_{obs}}$. Noting the adjoint operator with $*$, we compute this gradient by using the adjoint model of the CTM (Lions, 1971; Pudykiewicz, 1998) that reads:

$$ (6) \quad \partial_t c^* + \nabla \cdot (\mathbf{K}^T \nabla c^*) + \nabla(\vec{u}c^*) - (\nabla \cdot \vec{u})c^* - \tau_{loss}c^* = \left(\frac{dm}{dc}(c(s))\right)^* \cdot 2\mathbf{R}^{-1} (m(c(s)) - m^{obs}) $$
with $c^*$ the adjoint state.

The adjoint model (6) is closed by the following final and boundary conditions:

* a null final condition $c^*(x, y, t = T) = 0 \quad \forall (x, y) \in \Omega$,

* a null diffusive flux $\mathbf{K}^T \nabla c^* \cdot \vec{n} = 0 \quad \forall (x, y) \in \partial\Omega$,

* a null outgoing convective flux $\vec{u}c^* \cdot \vec{n} = 0 \quad \forall (x, y) \in \partial\Omega \cap \{(x, y) | \vec{u}(x, y, t) \cdot \vec{n} > 0\} \quad \forall t \in ]0; T[$ with $\vec{n}$ the outgoing normal vector.

Further details on the derivation of the adjoint model (6) can be found in Sec. A.3.1. This adjoint model (6) can be interpreted as the retro-propagation with diffusion of the error between the prediction of the observed variable and the observations $m(c(s)) - m^{obs}$ from the final time $t = T$ to the initial time $t = 0$. The output of the adjoint model is also referred to as the retro-plume in the following. Let us note that the adjoint model (6) requires the observation operator to be differentiable.

Once the adjoint state is computed, the gradient $\nabla_{s,j}$ is computed using the expression (see more details in Sec. A.3.1):

$$ (7) \quad \nabla_{s,j}(s) = \nabla_{s,j_{reg}}(s) - c^*(s) $$

The gradient-descent algorithm is initiated with a null control $s_0 = 0$.

The adjoint model (6) is solved using an implicit downwind finite volume scheme on the same Cartesian grid as the direct problem. The whole VDA process is implemented in the git repository of the PheroSensor project, see Malou et al. (2024c) and Sec. 5, including the solver of the adjoint model and a gradient-descent algorithm. In this repository, one can find convergence tests of the numerical scheme of the adjoint model, test of the adjoint operators as well as a validation of the estimation of the gradient by comparison with finite differences.

### 2.2.3. Sensor contribution to gradient evaluation.
We stress that in our context of a network of independent sensors dispatched over the landscape, we can assume that the observation error covariance $\mathbf{R}$ is block-diagonal. Moreover, let us note that the left-hand side of Eq. (6) is linear by definition of the adjoint model. In addition, we assume that the observations operator can be split into a sum of one-sensor observation operator $m(c) = \sum_i m_i(c)$ (avoiding the degenerate case of overlapping sensors). From this assumption, we can inject $m = \sum_i m_i$ in the adjoint model (6) and split the source term to get the following one-sensor adjoint model:

$$ (8) \quad \partial_t c_i^* + \nabla \cdot (\mathbf{K}^T \nabla c_i^*) + \nabla(\vec{u}c_i^*) - (\nabla \cdot \vec{u})c_i^* - \tau_{loss}c_i^* = (d_c m_i(c(s)))^* \cdot 2\mathbf{R}^{-1} (m_i(c(s)) - m_i^{obs}) $$

with $c_i^*$ the one-sensor adjoint state related to the $i$th sensor. Since the adjoint model (6) is linear, it results that the adjoint state $c^*$ can also be expressed as the sum of the $c_i^*$:

$$ (9) \quad c^* = \sum_i c_i^* $$
It results that $c_i^*$ is the contribution of the $i^{th}$ sensor to $c^*$, and thus to the gradient of $j_{obs}$. This observation will be used later to compare the contribution of each sensor to the data assimilation. The comparison of the contribution of each sensor then enables to determine which sensors contribute the most and which contribute the less, and to propose strategies for repositioning the sensors with least contribution. Repositioning strategies based on such comparison will be used in the realistic case of this study, see more details in Sec. 2.3.2 and the results of the repositioning strategy in Sec. 3.6

More details on the derivation can be found in the annex Sec. A.3.1. Let us note that in our case, the assumption that the observations operator can be split into a sum of one-sensor observation operator is true, see Sec. 2.3.1.

2.2.4. Stopping criteria. As per Morozov's discrepancy principle, the optimization algorithm should be continued until the gap between the observations and the system state achieves the level of the accuracy of the observations. Hence, the gradient descent algorithm is terminated and the optimum $s_a$ is assumed to be reached when the cost $j_{obs}$ attains a value lower than the noise norm $||\epsilon||_{L^2(\Omega_{obs})}^2$.

## 2.3. Pheromone sensors and test cases

2.3.1. Observation operator of the pheromone sensors. The state variable $c$ is a function defined over $[0, T] \times \Omega$. We assume that we have at disposal $N_s$ sensors that are dispatched throughout the landscape and remain in the same position $(x_i, y_i)$ all along the time domain. We furthermore assume that, given $[t_1^{obs}, \dots, t_{N_T}^{obs}]$ a time sampling of the time domain, every sensor produces a data at time $t_k^{obs}$, $k = 1, \dots, N_T$, so that the observation space $\mathcal{O}$ has dimension $N_s \times N_T$. Let us note $e_{i,k}$ the vector, for sensor $i$ and time $k$, of the canonical basis of $\mathcal{O}$.

Then, we define the observation operator $c \mapsto m(c)$ mapping the state variable $c$ towards the observation space $\mathcal{O}$ with the equation

$$ (10) \qquad m(c) = \sum_{i=1}^{N_s} m_i(c) \quad \text{with} \quad m_i(c) = \sum_{k=1}^{N_T} e_{i,k} \int_{t_k^{obs}-\delta_t}^{t_k^{obs}} \int_{\Omega} \delta_{(x_i, y_i)}(x, y) c(\tau, x, y) d\tau dxdy $$

where $\delta_{(x_i, y_i)}$ is the Dirac function at the localization of the $i^{th}$ sensor and $\delta_t$ a time window. For observation $k$ of sensor $i$, this expression averages over a time window of width $\delta_t$ before the observation time the pheromone concentration at location $(x_i, y_i)$. This models a pre-concentrator that actively filters air during a time window $\delta_t$ and releases the accumulated quantity into the sensor for the measure. Let us note that the time sampling of the observations must satisfy $\delta_t < t_1^{obs}$ so that $t_1^{obs} - \delta_t > 0$. Let us note that in this case, the operator $c \mapsto m(c)$ is linear and is split in the sum of one-sensor observation operators $m_i$, as mentioned in Sec. 2.2.3.

2.3.2. Setting up test cases. All along the paper, we will consider two test cases. The first one is a simplified yet illustrative scenario, referred to as the Toy case, that will be used to test and analyze different biology-informed regularization strategies. The second one depicts a realistic and representative scenario for invasion of an alien insect in an agricultural landscape. We chose the fall armyworm (FAW), a.k.a. *Spodoptera frugiperda*, as the model insect considered in this study. The FAW is a moth that belongs to the family Noctuidae. This insect pest induces the highest yield loss for maize (Savary et al., 2019), but can grow on hundreds of plant hosts, causing severe damage in grasses (IPPC Secretariat, 2021). It is expected that global warming reduces the populations in the southern America area due to warmer and drier climates (IPPC Secretariat, 2021), but could favour the emergence of the species in the European Union (EFSA PLH Panel
et al., 2018). It is then a suitable test case for a realistic scenario in an epidemiological surveillance framework in Europe. In this section, we will detail the set up of the environmental and landscape parameters needed for solving the pheromone propagation model (1) in the two test cases. The inference method and the additional biology-informed regularization terms are tested using synthetic dataset generated from a reference source term, hereinafter referred to as the target source term $s_t$. Further details for data generation and inference will be added in this section.

**Toy case.** The computation spatial domain of the Toy case is a rectangular grid of size 40m × 60m (i.e. 0.24ha) over a time window of $T = 20$s. It is built up upon a fictive insect that constantly releases pheromone at an emission rate of $q = 2.7\text{pg}\cdot\text{s}^{-1}$ per insect. The insects density is constant $p_0^{\text{toy}} = 0.2\text{m}^{-2}$ in a fixed area of the domain located in a circular support $\Omega_{\text{supp}} = S$ centred in $(x_s, y_s) = (10, 10)$ and of radius $R = 5\text{m}$, see red circle in Fig. 2, which corresponds to approximatively 16 individuals. The wind field $u$ is constructed from a constant wind speed $|\vec{u}| = 8\text{m}\cdot\text{s}^{-1}$ and a spatially-homogeneous time-dependent wind direction that linearly varies from $\theta = \frac{\pi}{4}$ at $t = 0$ to $\theta = \frac{\pi}{2}$ at $t = T$. The loss coefficient $\tau_{\text{loss}}$ is computed to approximate realistic vertical fluxes in a 3D simplified CTM model (see more details in Sec. A.1), resulting in a constant value of $\tau_{\text{loss}} = 0.01\text{s}^{-1}$. The diffusion tensor $\mathbf{K}$ is assumed to be isotropic and homogeneous over time and space. As the diffusion value $\mathbf{K} = 50\text{m}^2\text{s}^{-1}$ is a typical value to model an unstable atmosphere (BEA Fisher et al., 1981), we use the value $\mathbf{K} = 10\text{m}^2\text{s}^{-1}$, which corresponds to a moderately unstable atmosphere. A favorable habitat, used in the Tikhonov regularization, see Sec. 3.1.1, is defined as a square with a width of 20m, overlapping the source support $S$, see zone delimited by the dotted red line on Fig. 2.

$N_s = 35$ synthetic sensors are randomly placed each at a different position in the landscapes with a uniform distribution, see Fig. 2. Moreover, the accumulation time step is set to $\delta_t = 2$s. Noise is subsequently added. We take $m_{\text{obs}}(t_i, x_j, y_j) = \max(0, m(c(s_t); t_i, x_j, y_j) + \epsilon_{i,j})$, in order to avoid negative values, with $c(s_t)$ the results of the direct model (1) computed from the source term $s_t$ and $\epsilon_{i,j} \sim N(0, \sigma_{i,j}^2)$ a Gaussian noise. Since pre-concentration uncertainties scale with pheromone density, we scale the standard deviation with the model observations: $\sigma_{i,j} = 0.01m(c(s_t); t_i, x_j, y_j)$.

To discretize the spatio-temporal domain, fixed space and time steps are used, resp. $\Delta x = \Delta y = 0.5\text{m}$ and $\Delta t = 0.1\text{s}$. This leads to a source term $s(x, y, t)$ to infer with 1920000 degree of freedom.

**Realistic set up: the FAW case.** The FAW case is based on a real landscape in central France (Adamczyk-Chauvat et al., 2024), see Fig. 3a for the computational domain (dashed lines). The same homogeneous diffusion tensor ($\mathbf{K} = 10\text{m}^2\text{s}^{-1}$) as in the Toy case is used. The wind field $\vec{u}$ is recovered from outputs of the French weather forecast model Arome, which is the regional numerical weather forecast model from MétéoFrance (French weather forecast service, data and documentation available at https://donneespubliques.meteofrance.fr/?fond=produit&id produit=131&id_rubrique=51) corresponding to the actual wind field in area of interest at noon the 2022/11/16. As the spatial sampling of the Arome model is about $0.01^\circ$ ($\approx 1.3\text{km}$), see Fig. 3a, the wind field is projected on the computation grid by linear interpolation, and kept stationary over the computation time window.

As the loss coefficient $\tau_{\text{loss}}$ depends on the land cover, the land use heterogeneities in the landscape are taken into account in the FAW case using land cover data, see Fig. 3b. The different values of $\tau_{\text{loss}}$ depending on the land cover as recapitulated in Tab. 1, see Sec. A.1 for computations. Note that in the Toy case, the constant value 0.01 corresponds to bare soil.
<figure><img src="image_4.png" /><figcaption>**Figure 2** – Map of the target source term $s_t$ circular support (constant within the red circle), the favorable habitat (square zone delimited by the red dotted line) and the sensor location (colored markers) for the Toy test case.</figcaption></figure>

<table><caption>**Table 1** – Loss coefficient $\tau_{loss}$ values depending on the land use. Sum. crop = Summer crop, wint. cereals = winter cereals.</caption><thead><tr><th>land use</th><th>sum. crop</th><th>maize</th><th>wint. cereals</th><th>grassland</th><th>forest</th><th>bare soil</th><th>urban</th><th>road</th></tr></thead><tbody><tr><td>τ<sub>loss</sub> (s<sup>-1</sup>)</td><td>0.0097</td><td>0.011</td><td>0.0105</td><td>0.0099</td><td>0.0095</td><td>0.01</td><td>0.009</td><td>0.009</td></tr></tbody></table>

FAW has a specific pattern of pheromone emission. During the mating period, the FAW females search for a suitable location and stay there until they lay their eggs. Each night, they emit a sexual pheromone that attracts the males. Pheromone emission is stopped just at the beginning of the mating process, so that at the population level, the pheromone emission exponentially decreases while females meet a mating partner and stop emitting. This pattern is repeated every night. We start the time window for pheromone recording three hours after nightfall, when emissions are at their peak (Tumlinson et al., 1986).

Recall that $s = pq$ with $p$ the density of insects and $q$ the quantity of pheromone emitted per time unit per insect. We use a constant emission rate per insect $q = 2.7pg.s^{-1}$. The reference density $p_{ref}^{FAW}$ of emitter insects is chosen to satisfy a realistic population dynamics model, see later Sec. 3.7, and is defined as

$$ (11) \qquad p_{ref}^{FAW}(t, x, y) = p_{0}^{FAW}(x, y) \exp \left( - \int_{0}^{t} \gamma(s) ds \right) $$

with $\gamma(t)$ the mating rate function representing the proportion of female that encounter a male at a given time and stop emitting (see Fig. 3c) and $p_0^{FAW}$ the initial density of females. In that regard, $t \mapsto \exp(-\int_0^t \gamma(s) ds)$ is the function representing the proportion of females still available for mating and that keep emitting at a given time. Initially, we suppose that approximately 60 females are uniformly distributed within a circular domain of radius $10m$. This corresponds to a constant initial density of $\bar{p}_0 = 0.2m^{-2}$ within $S = \{(x, y) \in \Omega | (x - x_s)^2 + (y - y_s)^2 < R^2\}$ a circular support centred in $(x_s, y_s)$ and of radius $R = 10m$, see the red circular domain on Fig. 3b:

$$ (12) \qquad p_0^{FAW}(x, y) = \bar{p}_0 \mathbb{1}_S(x, y). $$
<figure><img src="image_5.png" /></figure>

(a) Wind field as obtained from the weather forecast model, before interpolation on the computing grid (length and color of the arrows normalized by the wind speed). Nine vectors in total are available, at the center of the figure, at the corners and in the cardinal directions.

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

(b) Fragment of agricultural landscape, land cover and sensor positions (see sec. 3.6 and fig 13).

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

(c) FAW test case, mating rate $\gamma(t)$ as in Eq. (13) (blue) and total proportion of female still emitting pheromones $\exp(-\int_0^t \gamma(s)ds)$ (red, see Eq. (11)) vs time $t$.

<caption>**Figure 3** – Map of the computation domains (dashed lines, black for the total domain and red for the domain of the inference) with the wind field computed by the Arome model at a given time (left) and the division of agricultural parcels and soil occupation (right) for the FAW case.</caption>

The FAW females usually emit when laying on plants favorable for the larva. In the FAW case landscape, we assume that maize is the most favorable plant. Therefore, the support $S$ is located within a domain region covered by maize, see land cover in dark orange on Fig. 3b. Moreover, we assume that the mating rate function $\gamma$ is a logistic function, see blue curve on Fig. 3c:

$$ (13) \qquad \gamma(t) = \frac{K_{mat}}{1 + \exp(-\lambda(t - \tau))} $$

where $K_{mat}$ is the maximal mating rate, $\tau$ the midpoint time value and $\lambda$ the logistic growth rate (see Fig. 3c). Let us note that $\tau$ is related to the average time a male necessitates to reach a female while tracking back the pheromone source. Based on biological expertise, $\tau$ is taken such that 50% of the females have mated with a male within the 30 first minutes. Hence, $p_{ref}^{FAW}(t = 30 \text{min}, x, y) = \frac{1}{2}p_0^{FAW}(x, y)$. Using the mating maximal rate $K_{mat} = 0.1$ and a logistic growth rate $\lambda \approx 8.3 \times 10^{-3}$, it results that $\tau = 127.94 \text{min}$.
To begin with, 15 sensors are randomly dispatched with uniform distribution in the maize field of interest, see the land cover in dark orange on Fig. 3b. Then, the 6 sensors that contribute the less (6 last in the contribution order for the one-sensor criteria defined in Sec. 2.4.3) are randomly relocated inside the convex hull of the 9 remaining sensors. This relocation method is detailed in the results in Sec. 3.6 and illustrated in Fig. 13.

Note that the computation domain is narrowed around the 300m long and 400m wide favorable zone of maize field to reduce computation time, see domain delimited by the red dashed line on Fig. 3a. The accumulation time step is set to $\delta_t = 30min$, in the order of reasonable pre-concentration time. Noise is subsequently added with the same noise model than in the Toy case: $m_{obs}(t_i, x_j, y_j) = max(0, m(c(s_t); t_i, x_j, y_j) + \epsilon_{i,j})$ with $c(s_t)$ the results of the direct model (1) computed from the source term $s_t$ and $\epsilon_{i,j} \sim N(0, \sigma_{i,j}^2)$ a Gaussian noise with standard deviation $\sigma_{i,j} = 0.05m(c(s_t); t_i, x_j, y_j)$.

To discretized the spatio-temporal domain of this case, fixed space and time steps are used, resp. $\Delta x = \Delta y = 2m$ and $\Delta t = 30s$. This leads to a source term $s(x, y, t)$ to infer with 10800000 degree of freedom.

## 2.4. Benchmarking criteria

Once the optimization algorithm is terminated, the optimum $s_a$ is compared against the target $s_t$, that is used to generate the data, using the criteria presented below to evaluate the accuracy of the localization of the insects.

### 2.4.1. RMSE-type criteria.
A first class of criteria is based on the Root x Mean Square Error (RxMSE) between the target $s_t$ and a value of $s$, that measures the accuracy of the data as-simulation process. The x indicates here that we average over the spatial dimension only. In this study, the errors of interest are the quadratic error on the whole domain $\Omega$ and on $\Omega_{supp}$, the support of the target $s_t$. We set the total RxMSE criteria by

$$ (14) \qquad R \times \text{MSE}_{\text{tot}}(t; s) = \sqrt{\frac{1}{|\Omega|} \int_{\Omega} (s(X, t) - s_t(X, t))^2 dX} $$

and the support RMSE criteria by

$$ (15) \qquad R \times \text{MSE}_{\text{supp}}(t; s) = \sqrt{\frac{1}{|\Omega_{\text{supp}}|} \int_{\Omega_{\text{supp}}} (s(X, t) - s_t(X, t))^2 dX} $$

While the error on the whole domain gives a global assessment of the accuracy of $s_a$, by also measuring false emissions estimated outside of the target support, the error on the target support enables to quantify the error made on the specific zone where the pheromone are emitted, and thus where the insects are located.

### 2.4.2. False discovery rate and miss rate.
The presence of insects at a given time and point of the domain is assessed when the quantity of pheromone emitted is above a given threshold $\epsilon_s$. This threshold depicts the minimal quantity of pheromone emitted by a single insect over the area of a numerical mesh cell. According to this presence criteria, the subdomain where the insects are located given a source term is the superlevel set

$$ (16) \qquad \Omega_{s/s}(t; s, \epsilon_s) = \{(x, y) \in \Omega \mid s(x, y, t) > \epsilon_s\} $$

This implies that the superlevel set $\Omega_{s/s}(t; s_a, \epsilon_s)$ (set in blue dashed lines on Fig. 4) is the sub-domain where the insects are predicted to be present. On the other hand, the superlevel set $\Omega_{s/s}(t; s_t, \epsilon_s)$ (set in red dashed lines on Fig. 4) is the sub-domain where the insects are truly present.
<figure><img src="image_8.png" /><figcaption>Figure 4 - Illustration of the predicted and true superlevel sets and the false absence and presence predictions</figcaption></figure>

These two sets will be compared using dissimilarity criteria in a same way as the Jaccard similarity coefficient or the confusion matrix in supervised statistical learning, see e.g. (Stehman, 1997). The subdomain $\Omega_{fp}(t; \epsilon_s) = \Omega_{sls}(t; s_a, \epsilon_s) \cap \overline{\Omega_{sls}(t; s_t, \epsilon_s)}$ corresponds to the subdomain where the insects are wrongly predicted to be present (false presence prediction, set in red dotted lines and blue dashed lines on Fig. 4) and $\Omega_{fa}(t; \epsilon_s) = \overline{\Omega_{sls}(t; s_a, \epsilon_s)} \cap \Omega_{sls}(t; s_t, \epsilon_s)$ the subdomain where the insects are wrongly predicted to be absent (false absence prediction, set in red dashed lines and blue dotted lines on Fig. 4). Noting $|\Omega|$ the area of the set $\Omega$, we then define $\rho_{fa}(t)$ the false absence prediction rate (or miss rate, i.e. the proportion of insect true presence that is predicted as absence):

$$ (17) \qquad \rho_{fa}(t) = \frac{|\Omega_{fa}(t)|}{|\Omega_{sls}(t; s_t, \epsilon_s)|} = \frac{|\Omega_{fa}(t)|}{|\Omega_{fa}(t) \cup \Omega_{tp}(t)|} $$

and $\rho_{fp}(t)$ the false presence prediction rate (or false discovery rate, i.e. the proportion of predicted presence that is incorrect):

$$ (18) \qquad \rho_{fp}(t) = \frac{|\Omega_{fp}(t)|}{|\Omega_{sls}(t; s_a, \epsilon_s)|} = \frac{|\Omega_{fp}(t)|}{|\Omega_{fp}(t) \cup \Omega_{tp}(t)|} $$

In a context where we aim to prevent and control the infestation, this miss rate must be kept the smallest possible, or even 0. On the other hand, in a context of precision agriculture, the false discovery rate will lead to useless application of pesticide. Thus, the objectif is to have the smallest possible false discovery rate while keeping the miss rate close to 0. In the following, we use the threshold value $\epsilon_s = 0.065$.

### 2.4.3. One-sensor adjoint state-based criteria.
The purpose of this criteria is to assess how the regularization strategies influence indirectly the contribution and importance of the sensors and their data, and if some strategies require less sensors to achieve the inference with a similar accuracy on $s_a$.

Let $(s_k)_{k \in \mathbb{N}}$ denotes the sequence of source terms given the trajectory taken by the optimization algorithm in the control space. The gradient descent algorithm iterations can then be expressed as $\forall k > 0$:

$$ (19) \qquad s_{k+1} = s_k - \eta_k \nabla_s j(s_k) $$

with $(\eta_k)_{k \in \mathbb{N}}$ the sequence of step sizes. From this expression of one iteration, the optimum $s_a$ can then be expressed as function of the sequence $(s_k)_{k \in \mathbb{N}}$:

$$ (20) \qquad s_a = s_0 - \sum_k \eta_k \nabla_s j(s_k) $$
Using the expression (7) of $\nabla_{sj}$ and the splitting of the adjoint state $c^*$ into the one-sensor adjoint state $c_i^*$ as in Eq. (9), the optimum $s_a$ can then be expressed as:

$$ 
\begin{aligned}
s_a &= s_0 - \sum_k \eta_k (\nabla_{sj} \text{reg}(s_k) - c^*(s_k)) \\
&= s_0 - \sum_k (\eta_k \nabla_{sj} \text{reg}(s_k)) + \sum_i \sum_k (\eta_k c_i^*(s_k)) \\
&= s_0 - \sum_k (\eta_k \nabla_{sj} \text{reg}(s_k)) + \sum_i c_{i,tot}^*
\end{aligned}
$$

where

$$ (21) \qquad c_{i,tot}^*(x, y, t) = \sum_k \eta_k c_i^*(s_k; x, y, t) $$

This illustrates that $c_{i,tot}^*$ enables to quantify the contribution of the $i^{th}$ sensor to the computation of the optimal $s_a$. The comparison of the regularization strategies with regard of the sensors contribution will be done by looking at $||c_{i,tot}^*||_{L^\infty(\Omega \times [0, T])}$, the maximal contribution of the $i^{th}$ sensor, and the sensors will be ordered by decreasing $||c_{i,tot}^*||_{L^\infty(\Omega \times [0, T])}$, hereinafter referred to as the contribution order. This representation allows to compare the contribution of the sensors with each other and to decipher if the maximal total contribution drops from a certain rank in the contribution order, thus indicating that the sensors located upwards in the ranking may have lower, or even negligible, impact on the estimation of $s_a$. Let us note that the quantity $c_{i,tot}^*$ can be directly computed using a trajectory-integrated one-sensor adjoint model in eq. (70) and does not require to compute $c_i^*(s_k)$ for each $k$ separately, see more details in Sec. A.3.1.

## 3. Results

The core methodological improvement of this study is the Biology-Informed Data Assimilation (BI-DA) framework. It is the reason why we present it in the results and not in the material and methods section. This framework is based on the classical data assimilation optimization problem (3) but with in addition regularization terms that are informed with biological knowledge of the insect. These BI-DA framework is tested in the two test cases defined in Sec. 2.3.2, i.e. the Toy case, characterized by small spatio-temporal scales, homogeneous landscape, time-varying but space-uniform wind field and constant emission pattern of the emitting insects, and the FAW case taking into account a real agricultural landscape, observed wind field and complex emission pattern mimicking the behavior of the insect pest targetted in this study: the FAW.

We first present different biology-informed regularization terms in Sec. 3.1. In Sec. 3.2, the pheromone propagation is illustrated in the Toy case, before computing a classical VDA without regularization in Sec. 3.3, to obtain an optimal pheromone emission $s_a$ of reference. In Sec. 3.4, the BI-DA is extensively benchmarked and compared to regularization-free VDA in the Toy case. Finally, BI-DA is used to track emission source in the FAW case. In Sec. 3.5, the pheromone propagation is illustrated in this realistic setting. In Sec. 3.6, a classical VDA used as reference is performed without regularization, raising the need for a better sensor placement. In Sec. 3.7, BI-DA is performed to infer the localization of the insects.

### 3.1. Biology-Informed DA: adding biological priors in pheromone emission inference

As mentioned in Sec. 2.2, regularization terms based on additional coarse estimation of the control are often considered to improve the convexity of $j$ (Kaltenbacher et al., 2008; Tarantola, 2005). In the classical VDA framework, the optimization problem (3) is strongly constrained by the pheromone propagation model (1): the pheromone distribution $c$ is explicitly obtained from
the direct problem, given the source s. Regularization terms can be seen as additional weak constraints, i.e. constraints introduced in the minimization problem as penalties. Inference problems that rely only on optimization problems weakly-constrained by a dynamical model are already treated in VDA (Trémolet, 2007) and machine learning contexts (Raissi et al., 2019).

We specifically introduce regularization terms based on biological models, and thus term the resulting VDA problem as BI-DA. Let us assume that s satisfies a biological model that can be written $M(s) = \sigma$, where $\sigma$ is a known biological prior. Then, one can consider a regularization term of the form

$$ (22) \qquad j_{\text{reg}}(s) = \|M(s) - \sigma\|_{L^2}^2 $$

The regularization term is then designed to minimize the residual of the model $M(s) = \sigma$.

The model's operator $M$ can be any kind of operator, including partial differential operator. For instance, recall that the emission source can always be written $s = qp$ as in Eq. (2), and thus, the model $M(s) = \sigma$ can be resulting from a population dynamics model on the insect density $p$. In the VDA minimization problem (3), we work in a differentiable optimization framework. Thus, the loss function $j(s)$ must be differentiable, which implies that $j_{\text{reg}}(s)$ and the underlying operator $s \mapsto M(s)$ have to be differentiable as well. However, this framework can be extended by working in a non-smooth and subdifferentiable optimization framework (Nesterov, 2018), e.g. by using the proximal-gradient algorithm (Parikh et al., 2014). This allows us to consider subdifferentiable $j_{\text{reg}}(s)$ and $M$, such as a LASSO regularization term, i.e. of the form $j_{\text{reg}}(s) = \|M(s) - \sigma\|_{L^1}$ or alike, see more details in Sec. 3.1.2. Recall from Eq. (5) that a weighted sums of multiple regularization terms can be considered in order to take into account multiple biological models.

We now introduce several biology-informed regularization terms that can be used to regularize the optimization problem (3). First, classical regularization terms will be introduced, such as Tikhonov (a.k.a. ridge), LASSO and group-LASSO regularization. For these regularization terms, we will pay a specific attention to the biological information we can add through these terms. Then, a population dynamic-informed $j_{\text{reg}}$ will be presented, in order to explicitly model biological mechanisms.

### 3.1.1. Biology-based Tikhonov regularization term.
The most classical $j_{\text{reg}}$ is the Tikhonov regularization term (Kaltenbacher et al., 2008):

$$ (23) \qquad j_{\text{reg}, \tau}(s) = \|s - s_b\|_{C^{-1}}^2 $$

with $s_b$ a background value of s and C the background error covariance operator. Classically, $s_b$ is given by a prior estimate of the true state, e.g. based on additional data that are not accurate enough and/or too sparse or based on a rough physical estimate, see e.g. Bocquet et al. (2014). From a biological point of view, a Biology-Informed Tikhonov $j_{\text{reg}}$ ($M(s) = s$ and $\sigma = s_b$) enables to take into account prior knowledge on the insects life habits, for instance their preferred habitat and insect life cycle. In this case, the term $s_b$ describes the expected spatial distribution and density of the insects with respect to the land use or the pheromone emission pattern.

### 3.1.2. LASSO-type regularization term and corresponding biological input.
LASSO-type regularization terms are often used to ensure sparsity and perform variable selection, see Tibshirani (1996). In the VDA experiments based only on the data, see Secs. 3.3 and 3.6, some spurious sources are detected, including near several sensors even though they are far away from the target source. This residual detection, particularly visible near final time (see Fig. 6), can be due to numerical error, to noise or to diffusive phenomena in the adjoint model. When computing the superlevel set
in eq. (16) for locating insects from the inferred quantity of pheromone emitted, these spurious pheromone sources can lead to false presence prediction.

**LASSO regularization.** In our context, we assume that the data are collected at the beginning of the infestation, so that we can expect that only a few insects are present. Hence, in order to select important local pheromones sources and to mitigate false presence predictions, a LASSO penalty is considered:

$$ (24) \qquad j_{reg,LASSO}(s) = \|s\|_1 $$

Recall that the optimization problem (3) with this $j_{reg}$ is no longer differentiable. A proximal gradient algorithm, namely the ISTA method, will be used to solve it, see more details in Annex A.3.3. Note that the ISTA is a Iterative Soft-Threshold Algorithm, meaning in particular that at each optimization iterations it thresholds any value under $\alpha_{reg,LASSO}\eta$, where $\alpha_{reg,LASSO}$ is the regularization weight of $j_{reg,LASSO}$ and $\eta$ is the descent step of the gradient-descent algorithm. Hence, assuming that we have a prior knowledge of a detection threshold under which pheromone emission can be discarded, a magnitude of $\alpha_{reg,LASSO}$ can be computed from $\eta$ in order to shrink the estimate of $s_a$ under this threshold at each optimization iterations.

**Group-LASSO regularization.** It is common that the insects do not move during the pheromone emission time period. It is e.g. the case for FAW, see Sec. 2.3.2. Therefore, the total quantity of pheromone emitted during this emission period can not be lower than the total amount of pheromone emitted by one insect on this time window. Conversely, if the total amount of emitted pheromone is lower, it implies that these pheromone sources are spurious and can be discarded. The following group-LASSO term is adapted to eliminate spatial positions with lowest accumulated pheromone emission:

$$ (25) \qquad j_{reg,gLASSO}(s) = \|s\|_{L^1(\Omega; L^2(0,T))} = \int_{\Omega} \left( \int_{0}^{T} s(t, x, y)^2 dt \right)^{1/2} dxdy $$

Similarly to the LASSO regularization, a proximal gradient algorithm will be used and more details are provided in A.3.3. For the group-LASSO regularization, this algorithm thresholds at each optimization iterations $s(t, x, y)$ for all $t$ if $\|s(\cdot, x, y)\|_{L^2(0,T)}$ is under $\alpha_{reg,gLASSO}\eta$, where $\alpha_{reg,gLASSO}$ is the regularization weight of $j_{reg,gLASSO}$ and $\eta$ is the descent step of the gradient-descent algorithm.

### 3.1.3. Population dynamic-informed regularization term.
Population dynamics models that model the spatio-temporal evolution of the density of insects $p$ can also be used to derive biological models for the quantity of pheromone emitted $s$. These population dynamics models can be written in a generic way as follows:

$$ (26) \qquad \partial_t p + \nabla \left( \sum_i F_i(p) \right) - \sum_i R_i(p) = 0 $$

with $F_i$ and $R_i$ are resp. flux and reaction terms that model the behavior of the insect specie of interest. Recall that $s(x, y, t) = p(x, y, t)q(t) \Leftrightarrow p(x, y, t) = \frac{s(x, y, t)}{q(t)}$ with $q(t) > 0$ the quantity of pheromone emitted per insect, which is assumed to be a biological knowledge. A population dynamics PDE equation of the form (26) can be expressed as the following PDE equation, in which $s$ is the state variable:

$$ \frac{1}{q} \partial_t s - \frac{\partial_t q}{q^2} s + \nabla \left( \sum_i \mathcal{F}_i(s) \right) + \sum_i \mathcal{R}_i(s) = 0 $$
with $F_i(p = \frac{s}{q}) = F_i(s)$ and $R_i(p = \frac{s}{q}) = R_i(s)$. From this PDE equation, the population dynamics-informed regularization term is the following:

$$ (27) \qquad j_{reg,PD}(s) = \alpha_{reg,PD} \|\partial_t s - \frac{\partial_t q}{q} s + q \nabla \left( \sum_i F_i(s) \right) - q \sum_i R_i(s)\|_{L^2}^2 $$

We note that this residual uses the PDE equation expressed in strong form. We take into consideration no boundary effects, but if needed, boundary conditions could be added to the minimization problem. For regular $F_i$ and $R_i$ operators, this residual lies in $L^2(0, T; \Omega)$.

Several potential terms $F_i$ and $R_i$ (and the associated $F_i$ and $R_i$) are possible and have different biological meaning. The flux terms $F_i$ can be used to model the motion of the insects in a random direction or in a specific direction (e.g. toward the preferred habitat). The reaction terms can be used to model birth-death processes or the fact that certain insects stop emitting pheromones. An extensive list of additional possible terms is added in the supplementary material Sec. A.2, including advection, diffusion, quorum sensing, Keller-Segel term to model pheromone sensing, and additional reaction terms such as linear birth-death process, logistic growth, Allee effect, or non-local terms. Furthermore, the corresponding gradients are explicitly computed in the supplementary Sec. A.3.2, so that the reader has all the necessary elements to implement the population dynamic-informed penalty fitting with the behavior of other insects of interest.

For instance, the females of the model insect that we consider, FAW, start emitting pheromones once found an emission site. During the emission time lapse, the insect does not move in order to efficiently attract a mating partner and stops emitting when it finds a mating partner and the mating starts. Furthermore, a proportion of female encounter a male at a given time following a mating rate function, see Sec. 2.3.2. Thus, the dynamics of the population of the insects emitting pheromone is modelled using only a simple linear birth-death process:

$$ (28) \qquad R_{bd}(p; x, y, t) = -\gamma(x, y, t)p(x, y, t) \Leftrightarrow R_{bd}(s; x, y, t) = -\frac{\gamma(x, y, t)}{q(t)}s(x, y, t) $$

with $\gamma$ [$s^{-1}$] the birth (if $\gamma > 0$) or death (if $\gamma < 0$) rate. In general manners, this term models a linear birth-death process where a proportion $\gamma(t)$ of the population dies or reproduces at each time $t$.

### 3.1.4. Choices of hyperparameter and hot start.
As mentioned previously, we will use the different regularization terms, which includes using them simultaneously. In this regard, recall that the overall regularization term can be expressed as the weighted sum of BI regularization terms (5) and $\mathcal{L}$ is the list of the regularization terms that we take into account. In the present study, we set $\mathcal{L} = \{ 'T', 'LASSO', 'gLASSO', 'PD' \}$ where T stands for Tikhonov, PD for population dynamics and gLASSO for group-LASSO. In the following, several combination of $j_{reg,i}$ will be tested and analysed, starting from only one $j_{reg,i}$ ($\alpha_j = 0 \forall j \neq i$) and up to three $j_{reg,i}$. The composite cases 'T + LASSO', 'T + PD' and 'LASSO + PD' designate the cases when, additionally to $\alpha_{gLASSO} = 0$, we have respectively $\alpha_{PD} = 0$, $\alpha_{LASSO} = 0$ and $\alpha_{T} = 0$. Finally, the case 'T + LASSO + PD' is when only $\alpha_{gLASSO} = 0$ and is designated as 'all reg' in the Toy case. Similarly, the case 'T + gLASSO + PD' is when only $\alpha_{gLASSO} = 0$ and is designated as 'all reg' in the FAW case. Moreover, in the following, 'no reg' refers to the case without regularization terms ($\alpha_i = 0 \forall i \in \mathcal{L}$).

The values of $\alpha_i$ are chosen such that they are the largest possible while satisfying the Morozov's discrepancy principle (see Sec. 2.2.4) in a reasonable computation time. In the following, in order to facilitate the solving of optimization problem, the regularization terms are progressively
added. This process is also known as a hot start. As a first step, we optimize without regular-
ization. Then, we add one by one the regularization terms initiating the optimization with the
optimum obtained previously. Similarly, in order to reach the highest values possible, the weight
coefficients αᵢ are increased step by step, restarting each time with hot start and repeating until
the Morozov's discrepancy principle is no longer satisfied.

3.2. Numerical illustration of pheromone propagation in the Toy case

<figure><img src="image_9.png" /><figcaption>**Figure 5** – Map of the concentration in pheromone c(s_t) (pg.m⁻²) computed from the depth-integrated CTM model (1) with the reference source term s_t for the Toy case presented in Sec. 2.3.2 at several times; in addition iso-lines of c are indicated with black dashed lines, and the reference source s_t support is displayed in white at final time.</figcaption></figure>

We first present numerical results for the pheromone propagation model (1) in the Toy case,
as set up in Sec. 2.3.2, in order to illustrate the pheromone plumes computed by the model
(see Fig. 5). As expected, the convection is dominant: the pheromones mainly propagate in the
direction of the wind, with little additional diffusive effect. The pheromone crosses the domain
within a couple of seconds during which a front propagates until reaching a quasi stationary
plume. As expected, time variations of the wind field, that is north-east at the beginning and
north at the end, induces deformation of the plume.

3.3. Classical VDA without regularization used as reference in the Toy case

We now illustrate the inference of emission sources in the Toy case without regularization
term in Fig. 6. We can observe that the optimum sₐ exhibits positive values, corresponding to
pheromone emission sources, in areas spreading out upwind of the sensor locations (see Fig. 2
for sensor locations. Some sensor locations are also visible in Fig. 6, indicated by spots). This de-
tection is referred to as the retro-plume. Indeed, the optimum sₐ is computed by descending the
gradient of the cost function (7), which is evaluated using the adjoint model (6). This equation
describes the retro-propagation with diffusion of the discrepancy between predictions and ob-
servations m(c(s)) – m{s}obs from final time t = T to initial time t = 0. Consequently, the gradient
descent algorithm mostly adjusts the initial prior s₀ in diffuse areas (attributable to the diffusion
term in the adjoint model) upwind (attributable to the advection term) of the sensor locations
(attributable to the source term).
<figure><img src="image_10.png" /><figcaption>**Figure 6** – Classical VDA without regularization for the Toy case. The optimal quantity of pheromone emitted $s_a$ is obtained by inference without regularization term for the Toy case at $T = 0$s (left panel) and $T = 20$s (right panel). The support of the target quantity of emitted pheromone $s_t$ (see Fig. 2) is delimited by the dashed lines. Some spurious detection spots can be seen, as sources (light dot) or sinks (dark blue dots) in the Toy case at $T = 20$s.</figcaption></figure>

Another consequence of the retro-propagation feature can be seen in Fig. 6: if at initial time, the estimated source $s_a$ covers a large area, at final time $t = 20$s, the retro-plume is mostly confined into a small vicinity of several active sensors. Indeed, near final time, as the retro-plume is only advected for a little time interval, the propagation of the adjoint state is mostly due to diffusion.

We also note that the reconstructed source $s_a$ exhibits negative non-physical values, that can be identified by dark blue dots near some sensors (see e.g. sensor locations $X = (17, 18.5)$ and $X = (14.5, 17.0)$ at time $t = 20$s, in Fig. 6, right panel). This phenomena is not due to numerical errors: it is a feature that directly comes from the noise in the data. The VDA intends to reproduce the data, including the noise even when the noise is not consistent with physical mechanisms. Therefore, during the retro-propagation of the noisy data, the optimal source $s_a$ includes pheromone sinks, i.e. negative values, in order to approximate the noisy target, even if it is not physical. If needed, the optimization problem could be solved under positive constraints in order to avoid these drawbacks.

### 3.4. Benchmark of the different regularization terms for BI-DA in the Toy case

We now aim to benchmark the different regularization terms introduced in Sec. 3.1 for BI-DA. We use the same supervised framework as in Sec. 3.3: a target source term $s_t$ is used to produce the synthetic data, and are compared to the estimated source $s_a$. Recall that $s_t = qp_{ref}^{toy}$ with $q$ the constant emission rate, and $p_{ref}^{toy}$ the pest insect density, which remains constant over time on a circular support in the Toy case, see Sec. 2.3.2. This lead to the following population
dynamics model:

$$ (29) \qquad \partial_t p_{ref}^{toy}(t, x, y) = 0 $$

which corresponds to the population dynamic-informed regularization term

$$ (30) \qquad J_{reg,PD}^{toy}(s) = \|\partial_t s\|_{L^2}^2. $$

We also recall that a favorable zone for the emitting insect is defined in a square $\Omega_{square}$ that includes the support $p_{ref}^{toy}$, see Fig. 2. Moreover, from the typical density of insects and the per capita pheromone emission rate given by the literature, we compute a surface emission rate $\bar{s}_{b}^{toy} = 0.54pg.s^{-1}m^{-2}$ that is used as biology-informed background estimate

$$ s_{b}^{toy}(t, x, y) = \bar{s}_{b}^{toy} \mathbb{1}_{\Omega_{square}}(x, y), $$

that will be taken into consideration through the Tikhonov regularization term (23), see Sec. 3.1.1.

We include in the benchmark the Tikhonov (T), LASSO, Population Dynamic (PD) regularization, and all combination of them, in order to understand how each regularization term operates and impacts the BI-DA. The optimum $s_a$ is then compared with the target source $s_t$ using the different criteria introduced in Sec. 2.4, i.e. total and support RMSE criteria (see resp. Eqs. (14) and (15)), miss and false discovery rates (see resp. Eqs. (17) and (18)), and the one-sensor adjoint state-based criteria (see Eq. (21)).

### 3.4.1. Visual comparison of estimated pheromone emissions $s_a$.
The estimated sources $s_a$ are displayed for each regularization term, including composite terms, in Fig. 7 at initial ($t = 0s$) and final times ($t = 20s$), for one single (top panel) and composite (bottom panel) regularization terms (see Sec. 3.1), together with the case without regularization ('no reg') taken as reference inference experiment.

At $t = 0s$, only slight discrepancies can be observed between the different regularization terms, either in single (top panel of Fig. 7) or composite ones (bottom panel), except for Tikhonov regularization ('T') that over predicts source emission on the favorable zone. Note that with LASSO regularization ('LASSO'), the retro-plume is less diffuse compared to the reference case ('no reg') since small values have been thresholded by the regularization.

At final time $t = 20s$, when fewer data are available for estimation, single regularization terms do not significantly improve the estimate of $s_a$ made in the reference case ('no reg'), except with the population dynamic-informed regularization term ('PD') which substantially increases prediction accuracy. In composite regularization terms, involving population dynamics is also key to improve predictions (see 'PD + LASSO', 'PD + T' and 'all reg' compared to 'T + LASSO'). Visually, coupling population dynamics with LASSO ('PD + LASSO') slightly improves the target estimate compared to population dynamics alone, while coupling it with a Tikhonov term seems to add false positive predictions outside the target in the favorable zone.

In conclusion, adding prior information on population behavior through the population-dynamics term seems to improve the reconstruction of target emissions $s_t$, but it is not clear if coupling it with other regularization terms is beneficial or not: quantitative assessment must be conducted.

### 3.4.2. RMSE criteria.
In order to have a more quantitative assessment of the discrepancies between $s_a$ and $s_t$ in the different cases (single and composite regularization terms), we introduced in Sec. 2.4.1 $RxMSE_{tot}$ (14), the error computed on the whole computational domain, and $RxMSE_{supp}$ (15), the error computed on the support of the target source only. The latter
<figure><img src="image_11.png" /><figcaption>**Figure 7** – Optimal sources estimated with different regularization terms. Optimal quantity of pheromone emitted $s_a$ obtained by inference with several regularization strategies for the Toy case. Upper panel, from left to right: with no regularization term (no reg), with the population dynamic-informed regularization term (30) (PD), with the LASSO regularization term (24) (LASSO) and with the Tikhonov regularization term (T) (23). Lower panel, from left to right: with combination of population dynamic-informed and LASSO regularization terms (PD+LASSO), with a combination of Tikhonov and LASSO regularization terms (T+LASSO), with combination of population dynamic-informed and Tikhonov regularization terms (PD+T) and a combination of all three population dynamic-informed, LASSO and Tikhonov regularization terms (all reg). In both panels, the optimal sources are indicated at $t = 0.0$s (upper row) and $t = 19.9$s (lower row).</figcaption></figure>
enables to assess the inference accuracy on the specific zone where the insects are located. We display in Fig. 8 the time evolution of $RxMSE_{tot}$ (top panel) and $RxMSE_{supp}$ (bottom panel) for the different (single and composite) regularization terms in comparison with the reference case without regularization ('no reg').

We first note that all the regularization terms provide sensibly equal estimation error inside the target support in a large part of the time domain, until $t \approx 17.5$s (see the bottom panel of Fig. 8), while estimation error on the whole domain or near final time is more widespread. This indicates that the main differences come from reconstruction outside the target support and at the end of the time window, when data are lacking.

Again, we note that adding population dynamics to the regularization maintains low error values at final time (see 'PD', 'PD + LASSO', 'PD + T', 'all reg', compared to the other terms, bottom and top panels), as the population dynamic model enables to propagate information from the initial time to the final time. However, coupling it with LASSO and Tikhonov reduces the error both inside the target and on the whole domain (see 'all reg', red curves, Fig. 8). Compared to 'PD' term alone (orange curves), this enhancement could come from complementary effects and efficient interactions between the different regularization.

In conclusion, coupling all the regularization terms outperforms all the other combinations of regularization terms for prediction accuracy all along the time window, both on the target and the whole domain. Population dynamics is mandatory to preserve low error levels near final time, when data are lacking.

### 3.4.3. Presence prediction, false discovery and miss rates.
In a practical view point, the important information for cultural practice is not the distribution of pheromone emission (which reflects quantitative distribution of insect population), but the presence of the pests that we aim to control. Therefore, we introduced in Sec. 2.4.1 a prediction of the insects presence derived from the pheromone emission map and based on the super level set of Eq. (16) of minimal emission value, see Fig. 4.

We represent in the upper panel of Fig. 9 the presence maps of insects computed with the different regularization terms and given by the target $s_t$, at three different times ($t = 0.0s, t = 10.0s$ and $t = 20.0s$). Pests are predicted as present inside the delineated color-coded zones. As already anticipated in the visual comparison (see Sec. 3.4.1), the Tikhonov regularization overestimates the presence of pests, while the 'LASSO', 'T + LASSO' and 'no reg' strongly underestimate pest presence at final time.

Cases involving population dynamics regularization ('PD', 'PD + LASSO', 'PD + T', 'all reg') keep covering the target zone all along the time domain by propagating with the population dynamics regularization the presence predicted at initial time. We can however notice a slight overestimation at final time with the case combining all three regularizations ('all reg'), probably due to the Tikhonov term as the boundary of the square favorable zone are slightly visible (Fig. 9, upper panel, $T = 20.0s$, red line). We can also notice a slight underestimation at the final time with the case 'PD + LASSO', probably due to the effect of the LASSO regularization. We also observe spurious presence predictions with the single population dynamics regularization, the single LASSO regularization and the combination of both (Fig. 9, upper panel, resp. 'PD', 'LASSO' and 'PD + LASSO') far from the target, likely because of the noise.

To get a better view of the amount of insects that are detected and missed during the time window, we use false discovery rate $\rho_{fp}$ (18), i.e. the proportion of area where the insects are predicted to be present but are actually absent, and the miss rate $\rho_{fa}$ (17), i.e. the proportion of
<figure><img src="image_12.png" /><figcaption>**Figure 8** – Evolution of total and support RMSE of $s_a$ against $s_t$. Time evolution of the root space-mean square errors between the optimum $s_a$ obtained by inference with several regularization strategies and the target $s_t$, on the whole domain ($R \times \text{MSE}_{\text{tot}}$, Eq. (14), upper panel) and on the target support ($R \times \text{MSE}_{\text{supp}}$, Eq. (15), lower panel) in the Toy case: with no regularization term (no reg, blue), with the population dynamic-informed regularization term (30) (PD, orange), with the LASSO regularization term (24) (LASSO, pink), with the group LASSO regularization term (25) (group LASSO, yellow), with the Tikhonov regularization term (23) (T, green), with combination of population dynamic-informed and LASSO regularization terms (PD+LASSO, brown), with a combination of Tikhonov and LASSO regularization terms (T+LASSO, purple), with combination of population dynamic-informed and Tikhonov regularization terms (PD+T, grey) and a combination of all three population dynamic-informed, LASSO and Tikhonov regularization terms (all reg, red).</figcaption></figure>

area where the insects are actually present but the prediction is missed. We recall that the objective here is to reduce the false discovery rate, equivalent in our context to a useless application of pesticide in the corresponding zone, while keeping the miss rate as small as possible, or even null, to control the infestation.

Miss rates $\rho_{fa}$ are comparable for all regularization at the beginning of the simulation, but show a strong miss rate increase at final times for penalties that do not involve '$PD$' or '$T$' (Fig. 9, lower panel, '$no reg$' and '$LASSO$', respectively blue and pink plain curves). When the Tikhonov regularization is involved, the miss rate value is low, but to the price of high values of false discovery rate ($'T'$, $'T+LASSO$ and $'PD+T'$, resp. green, purple abd grey plain and dashed curves). Adding the population dynamics term, the miss rate is controlled all along the time course of the simulation (Fig. 9, lower panel, '$all reg$', '$PD$', '$PD+LASSO$', '$PD+T$', respectively red, orange, brown and gray plain curves). The combination of all three regularizations ($'all reg$' , red curves) seems to be a good trade-off, the miss and false discovery rates are very stable, $\rho_{fa}$ remains zero on the whole time window while $\rho_{fp}$ is among the smallest one on the main part of the time window. With this combination, the three terms interact in an efficient way: the LASSO term seems to keep $\rho_{fp}$ small, the Tikhonov term reduces the $\rho_{fa}$ while the population dynamics term mitigates the rates evolution by propagating information from the initial time to the end of the time domain, when the data are scarce.
<figure><img src="image_13.png" /><figcaption>Figure 9 - Presence prediction and time evolution of false discovery and miss rates.</figcaption></figure>

Prediction of insects presence with several regularization strategies: with no regularization term (no reg, blue), with the population dynamic-informed regularization term (30) (PD, orange), with the LASSO regularization term (24) (LASSO, pink), with the Tikhonov regularization term (T, green), with combination of population dynamic-informed and LASSO regularization terms (PD+LASSO, brown), with a combination of Tikhonov and LASSO regularization terms (T+LASSO, purple), with combination of population dynamic-informed and Tikhonov regularization terms (PD+T, grey) and a combination of all three population dynamic-informed, LASSO and Tikhonov regularization terms (all reg, red). Top panel: map of presence prediction at three different times. The presence is obtained from the reconstructed source $s_a$ by computing a level set with Eq. (16). The true presence (i.e. the target support) is indicated in black. Lower panel: Time evolution of the miss rate $\rho_{fa}$ (solid lines), see Eq. (17), and false discovery rate $\rho_{fp}$ (dashed lines), Eq. (18).
In conclusion, 'PD' term is very interesting for controlling the miss rate on the whole time window and can be combined to 'T' term for enhanced reduction of $\rho_{fa}$, while coupling them with 'LASSO' reduces the false discovery rate.

### 3.4.4. Estimating active sensors with the one-sensor adjoint state criteria.
The visual comparison of optimal pheromone emission $s_a$ (see Sec. 3.4.1) indicates that sensors near the emission source seem to have greater impact on the retro-plume than sensors located further away. More generally, the number of sensors needed for an accurate prediction is an hyperparameter of the experimental design that could be evaluated. In order to obtain a proxy of this number, we developed a metric enabling to detect and number sensors that most influence the optimization process. To this end, we introduced in Sec. 2.4.3 the maximal total contribution $||c_{i,tot}^*||_\infty$, which is the total contribution of the data obtained with the $i^{th}$ sensor to the decrease of the cost function along the optimization algorithm. This metric enables to understand how a specific sensor impacts the estimation of $s_a$, and to rank the sensors accordingly for each regularizations, see Fig. 10. It allows in particular to understand if a regularization is able to supplement or amplify the information collected by the sensors. Note that the sensors are not involved in the regularization terms, so that the impact of the BI-DA regularization terms on $||c_{i,tot}^*||_\infty$ is indirect: the regularization changes the trajectory taken by the optimization algorithm in the control space, hence modifying the points where the one-sensor adjoint state $c_i^*$ is evaluated along the optimization path.

We can observe that in the reference case (Fig. 10, 'no reg', blue line), the 7 first sensors are particularly impactful. Then, a plateau value is attained until the 32nd sensor, where a drop is observed, indicating that the remaining sensors do not contribute significantly to the estimate of $s_a$. With the Tikhonov regularization, the contribution of the sensors located in the favorable zone is amplified, increasing the impact of the first sensors compared to 'no reg'. In comparison, the population dynamics regularization (Fig. 10, orange line) concentrates most of the information in the first 3 sensors: it enhances their contribution while the remaining sensors reach the plateau value. The combination of 'PD' and 'T' (grey line) shows a similar pattern than Tikhonov. Conversely, the 'LASSO' regularization, and all composite penalties containing 'LASSO' (Fig. 10, pink, brown, purple and red lines) tend to amplify the contribution of all the sensors and mitigates the drops of contribution at the beginning and the end of the contribution order.

In conclusion, by providing prior information, population dynamics and Tikhonov regularizations concentrate the information needed for inference in a limited number of sensors while the LASSO term make use of all the available sensors during the estimation process. When the three regularizations are combined, sensor contribution is distributed among the whole sensor network, like in the 'LASSO' case. Hence, 'PD' or 'T' regularizations could be used with a reduced number of sensors, while any additional sensor will be leveraged by the 'LASSO' penalty term for thresholding.

## 3.5. Numerical illustration of pheromone propagation in the FAW case

The direct model (1) is now computed on the FAW case, as described in Sec. 2.3.2. The FAW case strongly differs from the Toy case as it relies on a realistic landscape, with real application time and space scales, taking into account a pest, the FAW, harbouring a characteristic emission pattern, see Sec. 2.3.2 for details. Namely, only the females emit pheromones to attract males, but stop emitting at the start of mating. The emission model hence describes the population dynamics of female and male pairing.
<figure><img src="image_14.png" /><figcaption>**Figure 10** – Ranking sensor contribution in BI-DA. Maximal total contribution $||c_{i,tot}^*||_\infty$ vs sensor rank when sorted by contribution order for several regularization strategies in the Toy case: with no regularization (no reg, blue), population dynamic-informed (PD, orange), LASSO (LASSO, pink), Tikhonov (T, green), combination of population dynamic-informed and LASSO (PD+LASSO, brown), combination of Tikhonov and LASSO (T+LASSO, purple), combination of population dynamic-informed and Tikhonov (PD+T, grey) and a combination of all three population dynamic-informed, LASSO and Tikhonov regularization terms (all reg, red). The quantity $||c_{i,tot}^*||_\infty$ indicates the maximal (with respect to time and space) total contribution of sensor *i* all along the optimization path.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 11** – Map of the concentration in pheromone $c(s_t)$ ($\text{pg.m}^{-2}$) computed from the depth-integrated CTM model (1) with the reference source term $s_t$ for the FAW case presented in Sec. 2.3.2, in logarithmic scale and at several times. In addition iso-lines of $c$ are indicated with black dashed lines. Furthermore, the agricultural parcels are displayed with grey dashed lines together with the reference source $s_t$ support in white at $T = 110.5\text{min}$ (lower right panel).</figcaption></figure>
Direct model outputs are displayed in Fig. 11. As in the Toy case, convection is also largely predominant in the FAW case: the pheromones are also transported along the wind field, with little additional diffusive effect. Note that the wind field is taken from real weather broadcasting data. The pheromones spread across the domain within a couple of minutes until reaching a quasi-steady state plume within about ten minutes. After that, the decrease in time of the source $s$ which reflects the emission pattern of FAW that stops releasing pheromones at the start of mating (see Fig. 3c and Sec. 2.3.2) implies that the concentration $c$ decreases everywhere until $t \approx 110min$. Finally, $s$ becomes negligible and $c$ does not variate significantly any longer. The landscape composition (see the grey dashed lines on the bottom right panel of Fig. 11), which is taken into account through the loss rate $\tau_{loss}$, does not seem to influence significantly $c$.

### 3.6. Classical VDA and sensor positioning improvement in the FAW case

A first classical VDA without any regularization term is first conducted in the FAW case with sensors randomly dispatched in the field of interest, see left panel in Fig. 12. We highlight that this estimation problem is by far more difficult to solve than the Toy case for several reasons: first, there is a change in the spatio-temporal scales by one order of magnitude (from tens of meters and seconds in the Toy case to hundreds of meters and hours in the FAW case); next, the pheromone emission pattern mimics the biological behavior of the FAW and is no longer stationary like in the Toy case; finally, the sensor density is much higher in the Toy case (~ $10^{-2}$ sensor/m$^2$) than in the FAW case (~ $10^{-4}$ sensor/m$^2$). In this context, classical VDA fails to detect the position of the emitting insects (left panel in Fig. 12). The VDA algorithm does not sufficiently retro-propagate the discrepancies between data and predictions: during the optimization, the VDA could be stuck in a flat zone of the objective function due to insufficient signal of the sensors.

To partially counterbalance this effect, we repositioned several sensors by leveraging the one-sensor adjoint state-based criteria (see Sec. 2.4.3) for the classical VDA. Namely, we computed $c_{i,tot}^*$ in eq. (21) resulting from this first VDA without regularization, and selected the 2/5 of the sensors with the lowest contribution to the inference of $s$. These sensors were randomly repositioned within the convex envelop of the 3/5 sensors that were kept fixed, with uniform distribution (see top right panel of Fig. 13). We highlight the fact that no knowledge of the actual position of the target pheromone emission (i.e. of the emitting insects) has been used during repositioning, but only a criteria derived from the signal received by the sensors: in that sense, the repositioning method is unsupervised. On the top right panel of Fig. 13, we can observe that the sensors that have been repositioned were those located far or upwind from the targeted emitting insects, indicating that the criteria $c_{i,tot}^*$ is efficient in selecting sensors that would have been a priori identified as potentially uninformative in a supervised context where the target is known.

After sensor repositioning, the reconstructed $s_a$ by classical VDA now overlaps the target position at $T = 0min$, see right panel of Fig. 12. Furthermore, the VDA without regularization now enables to detect the insects in most of the target $s_t$ support at initial time (blue curves on top left panel of Fig. 13), unlike before relocation where no insects were correctly located (green curves). These observations were reflected in the miss and false discovery curves, that were equal to 1 before repositioning, and dropped to about 0.2 and 0.6 respectively after sensor displacement (13 lower panel).

We note that, at $T = 30min$, insects presence prediction vanishes (see Fig. 13 upper panel), and the miss and false discovery rates go up to 1, indicating inoperative presence prediction (see
<figure><img src="image_16.png" /><figcaption>**Figure 12** – VDA without regularization for the FAW case. The optimal quantity of pheromones emitted $s_a$ is obtained by inference without regularization term for the FAW case at $T = 0\min$ (upper panels) and $T = 30\min$ (lower panels). The support of the target quantity of emitted pheromone $s_t$ (see Fig. 3b) is delimited by the dashed lines. The left panel represents the source reconstruction with the initial random sensor distribution (uniform law over the maize field), whereas the right panel represents the same estimation after sensor re-positioning (see Sec. 3.6)</figcaption></figure>

Fig. 13). This feature is completely different than the drop of presence prediction accuracy observed at final time in the Toy case (see Sec. 3.4 and Fig. 9), that was inherent to boundary effects near final time and due to the retro-propagation feature of the adjoint model. Here, the collapse of discovery accuracy comes from the emission pattern of the FAW and from the accumulation time window of sensors rather than from a limitation of the data assimilation method.

## 3.7. Pheromone emission inference and pest localization with BI-DA in the FAW case

We now test the BI-DA framework in the FAW case. We first characterize the biology-informed regularization terms and benchmark them in this realistic framework.
<figure><img src="image_17.png" /><figcaption>**Figure 13** – Presence prediction and time evolution of false discovery and miss rates in the FAW case before and after sensor repositioning. The prediction of insects presence is obtained by classical VDA without regularization before and after sensors repositioning. Before repositioning, the sensors have been randomly dispatched in the maize field. After computation of $c_{i,tot}^*$, the two-fifth of the sensors with lower contribution to the source reconstruction were randomly re-allocated in the convex envelop of the three-fifth of the sensors with higher contribution. *Upper panel:* presence predictions at $T = 0$min and $T = 30$min before (green line) and after (blue line) sensor repositioning. The position of the sensors before and after relocation is indicated in the right figure. The presence is obtained from the reconstructed source $s_a$ by computing a level set with Eq. (16) (see Sec. 2.4.2). The true presence (i.e. the target support) is indicated in black. *Lower panel.* Time evolution of the miss rate $\rho_{fa}$ (solid lines), see Eq. (17), and false discovery rate $\rho_{fp}$ (dashed lines), in eq. (18).</figcaption></figure>
### 3.7.1. Biological priors in the FAW case.
First, we define the population dynamics regularization term, that must reflect the emission pattern of the FAW. Let us note $p_{ref}^{FAW}(t, x, y)$ the proportion of females that have not mated yet, and therefore keep emitting pheromones. The dynamics of the $p_{ref}^{FAW}$ population is the following:

$$ (31) \qquad \partial_t p_{ref}^{FAW}(t, x, y) = -\gamma(t) p_{ref}^{FAW}(t, x, y) $$

with $\gamma$ the mating rate at a given time that is represented by the logistic function (13). Let us note that $p_{ref}^{FAW}$ defined in eq. (11) is an analytic solution of the population dynamics model (31). Thus, in the FAW case, the following population dynamic-informed regularization term will be considered:

$$ (32) \qquad j_{reg,PD}^{FAW}(s) = \|\partial_t s + \gamma s\|_{L^2}^2 $$

As mentioned in Sec. 2.3.2, the FAW females usually emit pheromone on plants favorable for the larvae, which in this case is assumed to be maize. However, maize covers most of the domain considered and so, it does not provide relevant information. From the biological knowledge of the insects, one can deduce an expected density of insects and a quantity of pheromone emitted per female insect, and thus a local expected value of pheromone emission $\bar{s}_{b,loc}^{FAW}$. Without any additional information, this value will be used on the whole domain at the initial time. This leads to the following Tikhonov regularization:

$$ (33) \qquad s_b^{FAW}(x, y) = \bar{s}_{b,loc}^{FAW} \quad \text{so that} \quad j_{reg,T}^{FAW}(s) = \|s(t=0) - s_b^{FAW}\|_{L^2(\Omega)}^2 $$

Finally, we note that the emission pattern of FAW is ought to be particularly impactful on the classical LASSO term. The soft-threshold will preponderantly mitigate the source reconstruction when the emission is low, accelerating miss rate at the end of the simulation when the number of emitting females decreases. Conversely, we can expect that the LASSO penalty will not filter out spurious reconstructions at the beginning of the simulation when levels of emitted pheromones are still high since these spurious zones could represent higher quantities than the shrinking threshold. To counterbalance these effects, we decided to switch the LASSO penalty for a group-LASSO penalty by filtering out spatial areas that received a total amount of pheromones, during the whole time domain, under the shrinking threshold (see eq. (25)), so that

$$ j_{reg,gLASSO}^{FAW}(s) = \int_{\Omega} \left( \int_{0}^{T} s(t, x, y)^2 \right)^{1/2} dxdy. $$

A short comparison of LASSO and group-LASSO on the Toy case is made in the supplementary Sec. A.4.

### 3.7.2. Visual comparison of estimated pheromone emissions $s_a$.
We display the map of the reconstructed pheromone emission $s_a$ obtained by BI-DA in fig. 14 for different regularization terms, namely population dynamic-informed ('PD'), group-LASSO, a combination of Tikhonov and group-LASSO ('T + group LASSO'), and a combination of all three terms ('all reg'), that are compared with the reference case without regularization term ('no reg').

We can observe that $s_a$ obtained with the population dynamics regularization (top middle panel on Fig. 14) slightly enhances the one obtained in the reference case. The group LASSO regularization (top right panel on Fig. 14) thresholds pheromone emissions at location where the total emission over time is too small. The combination of 'T + group LASSO' regularizations (bottom left panel on Fig. 14) enables to reconstruct emissions in a large area, since the Tikhonov regularization pulls the emissions on the whole domain toward the right order of magnitude, while the group-LASSO regularization counterbalances that by filtering out emission in areas with low
<figure><img src="image_18.png" /><figcaption>**Figure 14** – Optimal quantity of pheromone emitted $s_a$ obtained by inference with several regularization strategies for the FAW case. Upper panel, from left to right: no regularization term (no reg), population dynamic-informed (PD), group-LASSO. Lower panel, from left to right: combination of Tikhonov and group-LASSO (T+group LASSO), combination of all three population dynamic-informed, group-LASSO and Tikhonov regularization terms (all reg).</figcaption></figure>

total emission. Finally, the combination of the three regularizations (bottom middle panel on Fig. 14) proposes a mix between the population dynamics regularization and the combination of group-LASSO and Tikhonov regularizations.

### 3.7.3. RMSE criteria between the optimum $s_a$ and the target $s_t$. 
The RMSE$_{tot}$ and RMSE$_{supp}$ curves show comparable error magnitude and trends for the different regularization terms.

The 'group LASSO' and '$T + group LASSO$' curves (pink and purple curves, Fig. 15), like the 'no reg' case (blue curve), present bumps around $t = 30\text{min}$ and $t = 60\text{min}$, due to the effect of the pre-concentration window in the observation operator in eq. (10). When adding the population dynamics regularization, this pattern is smoothed out, like in 'PD' (orange) and 'all reg' (red) curves. The '$T + group LASSO$' curve of RMSE$_{supp}$ present a bump at $t = 0\text{s}$ due to the Tikhonov regularization term, that is also smoothed on the curve of the 'all reg' case.
<figure><img src="image_19.png" /><figcaption>**Figure 15** – Root space-mean square errors between the optimum $s_a$ and $s_t$ for several regularization strategies. We display the $RxMSE_{tot}$ introduced in Eq. (14) (upper panel) and the $RxMSE_{supp}$ introduced in Eq. (15) (lower panel), versus $t$ for the FAW case: with no regularization term (no reg, blue), with the population dynamic-informed regularization term (32) (PD, orange), with the group-LASSO regularization term (25) (pink), with the Tikhonov regularization term and group-LASSO (T+group LASSO, brown) and with a combination of all three population dynamic-informed, group-LASSO and Tikhonov regularization terms (all reg, red).</figcaption></figure>

Hence, the different regularization terms have low impact on the overall error, but the prior knowledge included by the population dynamics term filters out discontinuities in the observation operator due to pre-concentration.

3.7.4. Prediction of insect presence. We now focus on insects presence predictions provided by the different regularization terms together with the corresponding miss rate $\rho_{fa}$ and false discovery rate $\rho_{fp}$ (respectively upper and lower panels of Fig. 16). Without any regularization ('no reg', blue curve, Fig. 16), false discovery and miss rates stay stable until $t = 30$min, when pheromone emission vanishes. As expected, adding 'group LASSO', with or without Tikhonov (pink and purple curves), reduces false discoveries. Meanwhile, 'T + group LASSO' regularization increases the miss rate while 'group LASSO' alone induces a slight reduction. Adding population dynamics regularization (orange curves) allows to reduce target misses compared to 'no reg', especially during the first 20 minutes, to the price of a higher $\rho_{fp}$. As observed previously on the RMSE criteria, the population dynamics regularization enables to smooth the curves, especially around $t = 30$min where a discontinuity is visible in the reference case (blue curve). The three regularizations interact between each other in an efficient way (Fig. 16, 'all reg', red curves) since the Tikhonov regularization enables to achieve one of the lowest $\rho_{fa}$ at $t = 0$min, the population dynamics regularization enables to propagate this good initial insect presence prediction and the group-LASSO regularization enables to threshold spurious pheromone emissions and reduces $\rho_{fp}$.

Consequently, combining the different regularizations represents a good trade-off to reduce both false discovery and miss rates, which is the final operational objective for minimizing pesticide use while keeping the amount of undetected pests low.
<figure><img src="image_20.png" /><figcaption>**Figure 16** – Prediction of insects presence with several regularization strategies: with no regularization term (no reg, blue), with the population dynamic-informed regularization term (32) (PD, orange), with the group-LASSO regularization term (25) (group LASSO, pink), with Tikhonov and group-LASSO regularization terms (T+group LASSO, purple), and a combination of all three population dynamic-informed, LASSO and Tikhonov regularization terms (all reg, red). Upper panel: map of presence prediction at three different times, true presence in black. Lower panel, miss rate $\rho_{fa}$ (17) (solid lines) and false discovery rate $\rho_{fp}$ (18) (dashed lines) vs t.</figcaption></figure>

3.7.5. Estimating active sensors with the one-sensor adjoint state criteria. We now look at the maximal total contributions $||c_{i,tot}^*||_\infty$ of the sensors and how the regularizations influence these sensors contributions in the FAW case, see Sec. 2.4.3. Similar trends than those observed in the Toy case are visible without regularization (blue curve on Fig. 17): a few sensors mainly contribute, then a plateau is attained and finally, a few sensors without significant contribution. In the FAW case, adding population dynamics has low impact on sensor contribution, except for the last ones. As in the Toy case, adding shrinkage through the 'group LASSO' penalty (pink, purple and red curves) amplifies and standardizes the contributions of the sensors.
<figure><img src="image_21.png" /><figcaption>**Figure 17** – Sensor importance in BI-DA. Maximal total contribution $||c_{i,tot}^*||_\infty$ vs sensor rank when sorted by contribution order for several regularization strategies in the FAW case. It includes no regularization (no reg, blue), population dynamic-informed (PD, orange), group-LASSO (group LASSO, pink), Tikhonov and group-LASSO (T+group LASSO, purple), and a combination of all three population dynamic-informed, LASSO and Tikhonov regularization terms (all reg, red). The quantity $||c_{i,tot}^*||_\infty$ indicates the maximal (with respect to time and space) total contribution of sensor *i* all along the optimization path</figcaption></figure>

In conclusion, the different penalties tend to act in the same way than in the Toy case. We however note that sensor contributions are less contrasted compared to the Toy case since the less informative sensors have been relocated in the convex hull of the most informative ones.

# 4. Discussion

## 4.1. Incorporating biological knowledge with BI-DA

In this study, we introduced a methodology to infer insect pest presence maps from signals of a network of pheromone sensors dispatched over a field that must be protected against a pest. The inference method is built upon variational data assimilation techniques: the aim is to reconstruct a control variable or a hidden state (here the quantity of pheromone emitted by the insects) by combining a direct model (here the pheromone propagation model derived from the atmospheric fluid dynamics) and data (here from the pheromone sensors). This combination is done by minimizing the distance between the data and an estimate of the observed variable. Once the pheromone emission reconstructed, they are taken as a proxy of the pest presence map. The numerical experiments of classical Data Assimilation (DA) problem highlight that classical DA enables to reconstruct the quantity of pheromone emitted almost only in a vicinity of the sensors. To tackle this issue, we proposed several regularization terms specifically designed to add prior Biological Information (BI) of the insects to improve the accuracy of the optimal solution. In particular, we introduced a population dynamics regularization that explicitly models the behavior of the pest of interest, particularly for pheromone emission. Thus, the BI-DA reconstruction is pulled towards a solution which is consistent with the physiology of the pest. The different regularization terms have been rigorously tested in a toy case, using different criteria, such as RMSE or false discovery and miss rates. The population dynamics regularization reduces the
RMSE, especially at the end of the time window when the data are missing, and reduces the miss rates without too much enhancing the false discovery rate, compared to the reference case without regularization. But using the population dynamics term conjointly with Tikhonov and LASSO terms provide the best results regarding presence/absence prediction, pulling the miss rates towards zero with lower rate of false discovery. However, if few sensors have a strong contribution to the reconstruction with the population dynamics regularization, which allows to expect to get comparable estimates with less sensors, coupling the different regularizations increase the number of sensors that have a strong contribution to the estimate.

## 4.2. Strong vs weak constraint in biology-informed data assimilation

One of the most classical method, referred in this study as classical VDA, consists in considering the direct model as strong constraint of the minimization problem: the estimate of the observed variable is the solution of the direct model, knowing a given estimate of the control variable. On the other hand, if the hidden state to reconstruct is the state variable of the direct model, one can use methods where the direct model acts as weak constraint of the minimization problem: the minimization problem aims to minimize the distance between the data and the estimate of the observed variable computed from the hidden state, and the residual of the direct model, see e.g. Trémolet (2007).

In the BI-DA method we propose, the pheromone propagation model and the pheromone sensor data are combined with biological information, especially model of the dynamics of the insect population. Whilst the pheromone propagation model is used as strong constraint of the minimization problem, the biological information are incorporated in the minimization problem as a weak constraint through the biology-informed regularization terms, see Sec. 3.1. Thus, the combination of data, pheromone propagation model and biological information is based on a mix of strong and weak constraints of the minimization problem.

One could propose a method where the minimization problem is strongly constrained by both the pheromone propagation model and the population dynamics model. For instance, in such situation, the control to infer would be the initial density of insects $p(t = 0, x, y)$ and the observed variable would be estimated through the coupling of the population dynamics model and the pheromone propagation model in the following way. For an estimate of $p(t = 0)$, the density of insects $p(t, x, y) \ \forall (t, x, y)$ is propagated with the population dynamics model and thereby the quantity of pheromone emitted $s(t, x, y) = q(t)p(t, x, y)$ can be computed. Then, the pheromone propagation model (given $s(t, x, y) \ \forall (t, x, y)$) is used to evaluate the concentration of pheromone $c(t, x, y) \ \forall (t, x, y)$ and the observed variable can be computed. One could go even further by coupling the population dynamics model of the female insect population emitting pheromone (depending on the density of male insects), the population dynamics model of the male population (depending on the pheromone concentration and its gradient) and the pheromone propagation model (depending on the density of female insects emitting pheromone). However, these methods based only on strong constraints require a good knowledge of the population dynamics. Actually, adding a strong constraint remains to assume that the observed data are perfectly explained by the model included in the strong constraint, with small modelling error.

On the contrary, the weak constraints resulting from the biological informations included in the BI-DA intrinsically take into account the modelling errors of the underlying biological processes. As the weak constraint tries to minimize the residual of the population dynamics model, an error model can be added to this residual to control its spatio-temporal correlation structure. Therefore, the BI-DA method does not require as much precise knowledge of the pest biology as strong constraints. It is then better suited for situations where the biological informations is
only partially known or may be subject to errors or uncertainties. Moreover, strongly constrain-
ing by both the pheromone propagation model and the population dynamics model could make
the minimization problem more difficult to solve and lead to identifiability issues. Indeed, the
minimum of the cost function may not be unique if the population dynamics is non-linear, and
adding multiple strong constraints may worsen the conditioning of the cost function.

4.3. Toward an effective use of BI-DA for pheromone sensor analysis

The proposed BI-DA method has been developped with consideration of practical constraints encountered in real applications. However, several points still require improvements before an effective use in an agricultural context. These improvements must be carried out in the framework of a global benefit-risks analysis of replacing pesticide-based prophylaxis by BI-DA-guided treatments: a life cycle analysis (LCA) of the sensors should be conducted, including manufacture, transport, deployment, use, recycling and computations and compared to the LCA of pesticides. Beside carbon emissions, benefit-risks analysis should also compare the impact of the two practices on health, environment and biodiversity.

First, computation time can be redhibiting for real applications, since solving BI-DA problem presently requires tens of hours in HPC facilities. This computational load is mostly due to solving the direct and adjoint pheromone propagation models, and to the number of iterations needed by the gradient descent algorithm. To speed-up direct and adjoint models computations, adapted metamodel technics, such as POD-based methods, can be studied. In this framework, PDE computations are replaced by fast approximations obtained by the surrogate model. Moreover, such metamodels can be learned from more complex models than the 2D pheromone propagation model used in this study. For instance, a 3D pheromone propagation model that accounts for vertical variability could be considered. While such models may be computationally infeasible for real-time applications, they can be solved and used offline to train a that could be used in-line for BI-DA. Speeding-up the gradient descent algorithm can be obtained by improving the convergence rate, thus reducing the number of iterations needed for convergence, using preconditionning technics.

Second, the selectivity of pheromone sensors, which is key for precise detection, has as the counterpart that a new sensor network should be deployed for a new target species. From a modelling point of view, increasing the number of target species will not change the CTM direct model which will be identical to all the species, so that a sub-linear increase of DA computation load can be expected. However, as the new species may have its proper behavior, BI-DA regularisation terms must be adapted to the new target, inducing model development time.

Finally, we chose for this study a data assimilation framework where the whole time series
of sensor signals are obtained at once and available for inference. A more sequential frame-
work where presence predictions are updated when new data are available would be usefull for
an adaptative and sequential sensor placement and data acquisition. This framework is particu-
larly suitable for epidemio-surveillance where a dynamic screening is key to adapt the treatment
to the pest dynamics. We note that our variational approach is mathematically equivalent to
Kalman smoothing on the whole time window whenever the biology informed-regularization
terms are differentiable and quadratic. However, we note that the LASSO and group-LASSO
terms are not differentiable and that a non-linear population dynamics model would result in a
non-quadratic PD term. Furthermore, the variational approach avoid to build and store the gain
and covariance matrices needed for Kalman filtering, thus reducing the computational load.
## 4.4. Uncertainty sources in the model

The proposed BI-DA method and the inference of pheromone sources are subject to various source of uncertainties. These uncertainties are of several types: uncertainties related to the environmental parameters, to the pheromone propagation model, to the test cases, to the sensors and the data.

**Uncertainties related to the environmental parameters.** The pheromone propagation is modelled using the 2D CTM (1). Thus, the resulting pheromone concentration depends on environmental parameters such as the wind velocity field $\vec{u}$ and the loss coefficient $\tau_{loss}$. However, the estimation of these two environmental parameters is itself uncertain.

In the realistic FAW case, the wind velocity is estimated using the data provided by the Arome model, see Sec. 2.3.2. This model have a spatial sampling of 0.01°(≈ 1.3km) and a temporal sampling of 1h. Thus, the smaller scale are not taken into account in the wind data. On the other hand, the spatial and temporal scale considered in realistic application are smaller: of the order of magnitude of a couple of minutes and of tens of meters. Therefore, in this context, uncertainties related to the small spatio-temporal scale of the wind field could be taken into account.

The loss coefficient $\tau_{loss}$ enables to parametrize all the vertical fluxes not directly modelled in the 2D model, see Sec. A.1.2. This coefficient depends in particular on the vegetation that absorbs a part of the pheromone that settles on the ground, and on the atmosphere stability. In this study and up to now, this coefficient is calibrated so the 2D CTM fits to a particular solution of the 3D CTM, see Sec. A.1.3. However, the calibration of $\tau_{loss}$ with respect to the land use remains rather empirical and source of uncertainties.

In order to understand the impact of these uncertainties on the pheromone propagation and on the pheromone source inference, an uncertainty quantification and a sensitivity analysis of the model with respect to these environmental parameters could be conducted, taking into account these different scales.

**Uncertainties related to the pheromone propagation model.** The 2D CTM model is derived from a 3D CTM model by averaging the 3D model over an atmospheric layer depth. This approximation is source of uncertainties since potential vertical variability of the wind in the layer, updrafts or vertical turbulence are roughly taken into account through the loss term. Additional modelling errors can come from inaccurate model parametrization or inaccurate modelling hypothesis. These uncertainties could be tackled by devoting more effort to modelling pheromone propagation. This includes the use of models that takes into account the random nature of the uncertainties and of metamodels that potentially take into account the vertical variability of the atmosphere.

**Uncertainties related to the test cases.** To illustrate the performance of the BI-DA method compared to a classical DA method, the data are generated from a target pheromone emission that we aim to retrieve. In these test cases, the computation domain, the location of the target emission and of the sensors and more generally the landscape design or the meteorological information are sources of uncertainties. A systematic quantification of these uncertainties can be performed using sensitivity analysis to assess the impact of these parameters on the BI-DA problem. In particular, the importance of some simplificatory assumption can be explored.

For instance, in both Toy and FAW cases, it is assumed that all the pheromone emissions are inside the computation domain. If some pheromone emission are located out of the computation domain, one might expect that small spurious pheromone emissions will be retrieved inside the
domain by the inference to compensate influx through the domain boundaries, that are not considered in the inverse model. However, as done in the FAW case with the sensors reallocation, see Sec. 3.6, this type of uncertainty can be bypassed using the following steps. First, we start over a large domain and perform the inference. Then, we perform the relocation process based on the one-sensor adjoint state-based criteria as done in the FAW case and reduce the computation domain (and potentially refine the mesh) to a narrower domain around the new sensors location. These steps can be repeated several times.

*Uncertainties related to the sensors and the data.* The data used in this study are synthetic: they are generated from a target pheromone emission and an observation model that includes a pre-concentration phase, see Sec. 2.3.2. As soon as suitable sensors are available, the classical DA and new BI-DA methods should be tested in progressively less controlled environment and conditions. These tests should start with controlled laboratory experiment with pheromone spray and controlled wind condition, progressing to tests conducted in open fields with uncontrolled environment. The technological developments of pheromone sensors will aim to obtain compact devices, making it possible to move the sensors in order to implement adaptative strategies such as sensor repositioning. This will enable us to test the proposed methods in real condition and improve them by taking into account the properties of the sensors and the statistics of their data. Efficient data acquisition strategies should be considered in order to mitigate data acquisition uncertainties, such as optimal experimental design or dynamical sensor repositioning strategies.

## 4.5. Sensor placement and reallocation

In BI-DA, we try to leverage spatio-temporal data obtained through a network of sensors that measure spatio-temporal concentration maps of pheromone. The time and space sampling of these map is key for the accuracy of the inference, as illustrated in the FAW case, see Sec. 3.6. The time sampling is given by the intrinsic functioning of the sensors. Due to the small amount of pheromone in the landscape, pre-concentration chambers accumulating the signal to detect are necessary, which determine the minimal time step of the time sampling. The space sampling is given by the spatial repartition of the sensor over the landscape.

In this study, we did not investigate the sensitivity of the inverse problem to the spatial placement of the sensors. Optimal design for building sampling design minimizing the variance of the inference can be constructed in comparable spatio-temporal inference problems (Parisey et al., 2022). This optimal design indicates, given environmental and landscape characteristics, where the sensors should be placed in order to get the best possible confidence on the estimate. It also provides, after hyperparameter selection, the number of sensors necessary to achieve a given accuracy.

Optimal placement is beyond the scope of this study and will be explored in future works. However, we introduced a specific metric measuring the maximal total contribution of a given sensor. This metric allows to order the sensors by contribution in the inference, hence distinguishing between key and accessory sensors. One could then imagine an iterative method for a dynamics placement of the sensors where the marginally contributing sensors are reallocated and moved to another place where they could enhance the accuracy of pest presence/absence estimates. We investigated how this criteria was impacted for varying BI-DA regularization terms. We found that some regularization terms supplement the information taken out of most sensors and only a few key sensors still significantly contribute to the inference. This is the case for population dynamics and Tikhonov regularization terms. As these regularization terms highlight the
sensors that do not significantly contribute in the light of both the propagation of the pheromone and the biological information, they can be used to study the placement and reallocation of sensors with regards to the pheromone propagation, including the atmosphere dynamics and environmental parameter, and to the biological knowledge of the insects, including the dynamics of its population. On the contrary, other regularization terms amplify the information taken out of the sensors and the contribution of every sensors. This is particularly the case for regularizations involving a LASSO term. This stresses the fact that optimal placement strategies will depend on the selected regularization terms.

## 5. Conclusion

We introduced in this study a numerical framework to solve the inverse problem for pest detection from the assimilation of pheromone concentration data. We specifically worked on adding prior biological knowledge on the population dynamics and the behavior of the pest to the inference problem. Diverse regularizations to biologically inform the data assimilation have been proposed and benchmarked. The population dynamics model enables to propagate emissions well estimated by the data to places and times that are not covered by the assimilation of the data set.

To reduce the pest miss rate, synonymous of a better control of the pest infestation, without increasing too much the false discovery rate, synonymous of unnecessary pesticide application, the following conclusions can be formulated:

(1) Some regularization terms, such as LASSO or group-LASSO, filter out spurious signals, hence enhancing the information provided by the sensors. This leads to lower false discovery rate but higher miss rate.

(2) Some regularization terms, such as the population dynamic-informed or Tikhonov regularization terms, provide prior information that can partially supplement the signal provided by the sensors, leading to lower miss rate, sometimes to the cost of a higher false discovery rate.

(3) The sensor contribution to the inference can be measured with the one-sensor adjoint state criteria. Regularization terms can decrease the contribution of some sensors and enable to efficiently discriminate sensors with high and low contributions. This implies that less sensors are required to achieve equivalent accuracy in the inference. Moreover, sensors with the lowest contribution can be moved to increase inference accuracy of emission sources. In this respect, population-dynamic or Tikhonov regularization terms are convenient to benchmark sensor position.

(4) Once sensor location is fixed, cumulating regularization terms can reduce miss rate while keeping false discovery rates in reasonable values. In particular, it is worth using a composite regularization mixing knowledge on favorable habitat and population dynamics and thresholding spurious pheromone emissions.

## Acknowledgments

Preprint version 2 of this article has been peer-reviewed and recommended by Peer Community In Mathematical & Computational Biology (https://doi.org/10.24072/pci.mcb.100313; Tannier (2025)).
## Funding

The authors acknowledge the support of the French National Research Agency (ANR) under the grant ANR-20-PCPA-0007 (PheroSensor).

## Conflict of interest disclosure

The authors declare that they comply with the PCI rule of having no financial conflicts of interest in relation to the content of the article.

## Data, script, code, and supplementary information availability

The toolbox containing the core functions for BI-DA is available at https://forgemia.inra.fr/pherosensor/pherosensor-toolbox and Malou et al. (2024b). The companion codes needed to reproduce the different computation and figures of this paper are available at https://forgemia.inra.fr/pherosensor/companion-code-bi-da and Malou et al. (2024a).

## References

Adam L, Branda M (2016). *Sparse Optimization for Inverse Problems in Atmospheric Modelling*. Environmental Modelling & Software **79**, 256-266. https://doi.org/10.1016/j.envsoft.2016.02.002.

Adamczyk-Chauvat K, Kassa M, Papaïx J, Kiêu K, Stoica RS (2024). *Statistical Inference for Random T-tessellations Models*. Application to Agricultural Landscape Modeling. Annals of the Institute of Statistical Mathematics **76**, 447-479. https://doi.org/10.1007/s10463-023-00893-3.

Bocquet M et al. (2014). *Introduction to the principles and methods of data assimilation in geosciences*. Notes de cours, École des Ponts ParisTech.

Boukal DS, Berec L (2009). *Modelling Mate-Finding Allee Effects and Populations Dynamics, with Applications in Pest Control*. Population Ecology **51**, 445-458. https://doi.org/10.1007/s10144-009-0154-4.

Britton N (1989). *Aggregation and the Competitive Exclusion Principle*. Journal of Theoretical Biology. https://doi.org/10.1016/S0022-5193(89)80189-4.

Conchou L, Lucas P, Meslin C, Proffit M, Staudt M, Renou M (2019). *Insect Odorscapes: From Plant Volatiles to Natural Olfactory Scenes*. Frontiers in Physiology **10**, 972. https://doi.org/10.3389/fphys.2019.00972.

Deutsch CA, Tewksbury JJ, Tigchelaar M, Battisti DS, Merrill SC, Huey RB, Naylor RL (2018). *Increase in Crop Losses to Insect Pests in a Warming Climate*. Science **361**, 916-919. https://doi.org/10.1126/science.aat3466.

EFSA PLH Panel et al. (2018). *Pest Risk Assessment of Spodoptera Frugiperda for the European Union*. EFSA Journal **16**. https://doi.org/10.2903/j.efsa.2018.5351.

Fisher BEA, MacQueen JF (1981). *A Theoretical Model for Particulate Transport from an Elevated Source in the Atmosphere*. IMA Journal of Applied Mathematics **27**, 359-371. https://doi.org/10.1093/imamat/27.3.359. eprint: https://academic.oup.com/imamat/article-pdf/27/3/359/1735525/27-3-359.pdf.

Fisher RA (1937). *The Wave of Advance of Advantageous Genes*. Annals of eugenics **7**, 355-369.

Furter J, Grinfeld M (1989). *Local vs. Non-Local Interactions in Population Dynamics*. Journal of Mathematical Biology. https://doi.org/10.1007/BF00276081.

Ganga Prasath S, Mahadevan L (2023). *Rheomergy: Collective Behaviour Mediated by Active Flow-Based Recruitment*. Proceedings of the Royal Society A **479**, 20220470. https://doi.org/10.1098/rspa.2022.0470.
Gebbers R, Adamchuk VI (2010). *Precision Agriculture and Food Security*. Science **327**, 828-831.
https://doi.org/10.1126/science.1183899.

Genieys S, Volpert V, Auger P (2006). *Pattern and Waves for a Model in Population Dynamics with Nonlocal Consumption of Resources*. Mathematical Modelling of Natural Phenomena. https://doi.org/10.1051/mmnp:2006004.

IPPC Secretariat. (2021). *Scientific Review of the Impact of Climate Change on Plant Pests*. FAO on behalf of the IPPC Secretariat. https://doi.org/10.4060/cb4769en.

Isakov V (2017). *Inverse Problems for Partial Differential Equations*. Vol. 127. Applied Mathematical Sciences. Cham: Springer International Publishing. https://doi.org/10.1007/978-3-319-51658-5.

Kaltenbacher B, Neubauer A, Scherzer O (2008). *Iterative Regularization Methods for Nonlinear Ill-Posed Problems*. Vol. 6. Walter de Gruyter. https://doi.org/10.1515/9783110208276.

Kareiva PM (1983). *Local Movement in Herbivorous Insects: Applying a Passive Diffusion Model to Mark-Recapture Field Experiments*. Oecologia **57**, 322-327. https://doi.org/10.1007/bf00377175. JSTOR: 4216967.

Keller EF, Segel LA (1970). *Initiation of Slime Mold Aggregation Viewed as an Instability*. Journal of Theoretical Biology **26**, 399-415. https://doi.org/10.1016/0022-5193(70)90092-5.

Keller EF, Segel LA (1971). *Model for Chemotaxis*. Journal of Theoretical Biology **30**, 225-234.
https://doi.org/10.1016/0022-5193(71)90050-6.

Khadka R, Aydemir N, Carraher C, Hamiaux C, Colbert D, Cheema J, Malmström J, Kralicek A, Travas-Sejdic J (2019). *An Ultrasensitive Electrochemical Impedance-Based Biosensor Using Insect Odorant Receptors to Detect Odorants*. Biosensors and Bioelectronics **126**, 207-213. https://doi.org/10.1016/j.bios.2018.10.043.

Le Dimet FX, Talagrand O (1986). *Variational Algorithms for Analysis and Assimilation of Meteorological Observations: Theoretical Aspects*. Tellus A: Dynamic Meteorology and Oceanography **38**, 97-110. https://doi.org/10.3402/tellusa.v38i2.11706.

Leelössy Á, Molnár F, Izsák F, Havasi Á, Lagzi I, Mészáros R (2014). *Dispersion Modeling of Air Pollutants in the Atmosphere: A Review*. Open Geosciences. https://doi.org/10.2478/s13533-012-0188-6.

Lions JL (1971). *Optimal Control of Systems Governed by Partial Differential Equations*. Vol. 170. Springer Verlag.

Lu Y, Li H, Zhuang S, Zhang D, Zhang Q, Zhou J, Dong S, Liu Q, Wang P (2014). *Olfactory Biosensor Using Odorant-Binding Proteins from Honeybee: Ligands of Floral Odors and Pheromones Detection by Electrochemical Impedance*. Sensors and Actuators B: Chemical **193**, 420-427. https://doi.org/10.1016/j.snb.2013.11.045.

Lushi E, Stockie JM (2010). *An Inverse Gaussian Plume Approach for Estimating Atmospheric Pollutant Emissions from Multiple Point Sources*. Atmospheric Environment **44**, 1097-1107. https://doi.org/10.1016/j.atmosenv.2009.11.039.

Malou T, Labarthe S (2024a). *BI-DA-companion-code*. https://doi.org/10.5281/ZENODO.11175422.

Malou T, Labarthe S (2024b). *Pherosensor-Toolbox*. https://doi.org/10.5281/ZENODO.11110047.

Malou T, Labarthe S (2024c). *Pherosensor-toolbox: a Python package for Biology-Informed Data Assimilation*. Journal of Open Source Software **9**, 6863. https://doi.org/10.21105/joss.06863.
Martinez D, Arhidi L, Demondion E, Masson JB, Lucas P (2014). Using Insect Electroantennogram Sensors on Autonomous Robots for Olfactory Searches. Journal of Visualized Experiments, 51704. https://doi.org/10.3791/51704.

McMahon T, Denison P (1979). Empirical Atmospheric Deposition Parameters—a Survey. Atmospheric Environment (1967) **13**, 571-585. https://doi.org/10.1016/0004-6981(79)90186-0.

McNeil JN (1991). Behavioral Ecology of Pheromone-Mediated Communication in Moths and Its Importance in the Use of Pheromone Traps. Annual Review of Entomology **36**, 407-430. https://doi.org/10.1146/annurev.en.36.010191.002203.

Moitra P, Bhagat D, Pratap R, Bhattacharya S (2016). A Novel Bio-Engineering Approach to Generate an Eminent Surface-Functionalized Template for Selective Detection of Female Sex Pheromone of Helicoverpa Armigera. Scientific Reports **6**, 37355. https://doi.org/10.1038/srep37355.

Murray JD (2007). Mathematical biology: I. An introduction. Vol. 17. Springer Science & Business Media. https://doi.org/10.1007/b98868.

Myrick AJ, Baker TC (2011). Locating a Compact Odor Source Using a Four-Channel Insect Electroantennogram Sensor. Bioinspiration & Biomimetics **6**, 016002. https://doi.org/10.1088/1748-3182/6/1/016002.

Negri RM, Bernik D (2008). Tracking the Sex Pheromone of Codling Moth against a Background of Host Volatiles with an Electronic Nose. Crop Protection **27**, 1295-1302. https://doi.org/10.1016/j.cropro.2008.02.008.

Nesterov Y (2018). Lectures on Convex Optimization. Vol. 137. Springer. https://doi.org/10.1007/978-1-4419-8853-9.

Nisbet RM, Gurney W (1982). Modelling Fluctuating Populations. John Wiley.

Oerke EC (2006). Crop Losses to Pests. The Journal of Agricultural Science **144**, 31-43. https://doi.org/10.1017/S0021859605005708.

Okubo A (1986). Dynamical Aspects of Animal Grouping: Swarms, Schools, Flocks, and Herds. Advances in Biophysics **22**, 1-94. https://doi.org/10.1016/0065-227X(86)90003-1.

Panofsky HA, Dutton JA (1984). Atmospheric Turbulence. Models and Methods for Engineering Applications. New York: Wiley.

Parikh N, Boyd S, et al. (2014). Proximal Algorithms. Foundations and trends in Optimization **1**, 127-239. https://doi.org/10.1561/2400000003.

Parisey N, Leclerc M, Adamczyk-Chauvat K (2022). Optimal Spatial Monitoring of Populations Described by Reaction-Diffusion Models. Journal of Theoretical Biology **534**, 110976. https://doi.org/10.1016/j.jtbi.2021.110976.

Pawson SM, Kerr JL, O'Connor BC, Lucas P, Martinez D, Allison JD, Strand TM (2020). Light-Weight Portable Electroantennography Device as a Future Field-Based Tool for Applied Chemical Ecology. Journal of Chemical Ecology **46**, 557-566. https://doi.org/10.1007/s10886-020-01190-6.

Perthame B (2006). Transport Equations in Biology. Birkhäuser Basel. https://doi.org/10.1007/978-3-7643-7842-4.

Pudykiewicz JA (1998). Application of Adjoint Tracer Transport Equations for Evaluating Source Parameters. Atmospheric Environment **32**, 3039-3050. https://doi.org/10.1016/S1352-2310(97)00480-9.

Raissi M, Perdikaris P, Karniadakis G (2019). Physics-Informed Neural Networks: A Deep Learning Framework for Solving Forward and Inverse Problems Involving Nonlinear Partial Differential Equations. Journal of Computational Physics **378**, 686-707. https://doi.org/10.1016/j.jcp.2018.10.045.
Sasaki Y (1958). *An Objective Analysis Based on the Variational Method*. Journal of the Meteorological Society of Japan. Ser. II **36**, 77-88. https://doi.org/10.2151/jmsj1923.36.3\_77.

Savary S, Willocquet L, Pethybridge SJ, Esker P, McRoberts N, Nelson A (2019). *The Global Burden of Pathogens and Pests on Major Food Crops*. Nature Ecology & Evolution **3**, 430-439. https://doi.org/10.1038/s41559-018-0793-y.

Shigesada N (1980). *Spatial Distribution of Dispersing Animals*. Journal of mathematical biology **9**, 85-96. https://doi.org/10.1007/BF00276037.

Shigesada N, Kawasaki K, Teramoto E (1979). *Spatial Segregation of Interacting Species*. Journal of Theoretical Biology **79**, 83-99. https://doi.org/10.1016/0022-5193(79)90258-3.

Skellam JG (1951). *Random Dispersal in Theoretical Populations*. Biometrika **38**, 196-218. https://doi.org/10.2307/2332328.

Sportisse B (2007). *A Review of Current Issues in Air Pollution Modeling and Simulation*. Computational Geosciences **11**, 159-181. https://doi.org/10.1007/s10596-006-9036-4.

Steffens C, Manzoli A, Oliveira JE, Leite FL, Correa DS, Herrmann PSP (2014). *Bio-Inspired Sensor for Insect Pheromone Analysis Based on Polyaniline Functionalized AFM Cantilever Sensor*. Sensors and Actuators B: Chemical **191**, 643-649. https://doi.org/10.1016/j.snb.2013.10.053.

Stehman SV (1997). *Selecting and Interpreting Measures of Thematic Classification Accuracy*. Remote Sensing of Environment **62**, 77-89. https://doi.org/10.1016/S0034-4257(97)00083-7.

Stockie JM (2011). *The Mathematics of Atmospheric Dispersion Modeling*. SIAM Review **53**, 349-372. https://doi.org/10.1137/10080991X.

Sun GQ (2016). *Mathematical Modeling of Population Dynamics with Allee Effect*. Nonlinear Dynamics. https://doi.org/10.1007/s11071-016-2671-y.

Takahashi LT, Maidana NA, Ferreira WC, Pulino P, Yang HM (2005). *Mathematical Models for the Aedes Aegypti Dispersal Dynamics: Travelling Waves by Wing and Wind*. Bulletin of mathematical Biology **67**, 509-528. https://doi.org/10.1016/j.bulm.2004.08.005.

Tannier E (2025). *Towards Accurate Inference of Insect Presence Landscapes from Pheromone Sensor Networks*. Peer Community in Mathematical and Computational Biology, 100313. https://doi.org/10.24072/pci.mcb.100313.

Tarantola A (2005). *Inverse Problem Theory and Methods for Model Parameter Estimation*. Vol. 89. SIAM. https://doi.org/10.1137/1.9780898717921.

Tibshirani R (1996). *Regression Shrinkage and Selection via the Lasso*. Journal of the Royal Statistical Society: Series B (Methodological) **58**, 267-288. https://doi.org/10.1111/j.2517-6161.1996.tb02080.x. eprint: https://rss.onlinelibrary.wiley.com/doi/pdf/10.1111/j.2517-6161.1996.tb02080.x.

Trémolet Y (2007). *Model-Error Estimation in 4D-Var*. Quarterly Journal of the Royal Meteorological Society **133**, 1267-1280. https://doi.org/10.1002/qj.94.

Tumlinson JH, Mitchell ER, Teal PEA, Heath RR, Mengelkoch LJ (1986). *Sex Pheromone of Fall Armyworm, Spodoptera Frugiperda (J.E. Smith): Identification of Components Critical to Attraction in the Field*. Journal of Chemical Ecology **12**, 1909-1926. https://doi.org/10.1007/BF01041855.

Wehrenfennig C, Schott M, Gasch T, Sauerwald T, Düring RA, Vilcinskas A, Kohl CD (2012). *Laboratory Characterization of Metal-oxide Sensors Intended for in Situ Analyses of Pheromones - SOMMSA Approach*. physica status solidi (a) **209**, 935-939. https://doi.org/10.1002/pssa.201100586.
Yazdani A, Lu L, Raissi M, Karniadakis GE (2020). Systems Biology Informed Deep Learning for In-ferring Parameters and Hidden Dynamics. PLOS Computational Biology 16. Ed. by Vassily Hatzi-manikatis, e1007575. https://doi.org/10.1371/journal.pcbi.1007575.

# Appendix A. Supplementary material

## A.1. Derivation of the depth-integrated chemical transport model

This section aims to present the derivation of the 2D model of pheromone propagation (1). In a first part, a classical 3D model of the pheromone propagation in the atmosphere, also called 3D CTM (Chemical Transport Model), is presented. Then, the 2D model, presented in Sec. 2.1 and used in this study, is derived from this 3D model. In the last part, the calibration of the 2D model using the 3D model is presented.

A.1.1. 3D chemical transport model. In this section, $C(\vec{x}, t) (g.m^{-3})$ denotes the density of chemical (here the pheromone) at a given point $\vec{x} = (x, y, z) \in \mathbb{R}^2 \times \mathbb{R}^+$ and at a time $t \in \mathbb{R}^+$. The equation of mass conservation of C can be written as the following PDE, see e.g. reviews in Sportisse (2007) or in Leelössy et al. (2014):

$$ (34) \qquad \frac{\partial C}{\partial t}(\vec{x}, t) + \nabla_{3D} \cdot (-\mathbf{K}_{3D}(\vec{x}, t)\nabla_{3D}C(\vec{x}, t) + \vec{u}_{3D}(\vec{x}, t)C(\vec{x}, t)) = S(\vec{x}, t) $$

with $-\mathbf{K}_{3D}\nabla_{3D}C + \vec{u}_{3D}C (g.m^{-3}.m.s^{-1})$ the mass flux of pheromone due to the diffusion and the advection of the pheromone in the atmosphere, and $S(\vec{x}, t) (g.m^{-3}.s^{-1})$ the source term, i.e. the quantity of pheromone emitted by the insects per seconds and per air volume in the present context. In this section, $\nabla_{3D}$ denotes the 3D gradient operator $\nabla_{3D} = (\partial_x \partial_y \partial_z)$, in order to distinguish with the 2D gradient operator used in the rest of the article.

This equation can be obtained from the reactant Navier-Stokes equations using the dilution hypothesis. This hypothesis consists in assuming that the chemical of interest is present in the atmosphere in a very small quantity and does not significantly impact the dynamic of the atmosphere. It results that the dynamics of the atmosphere and the dynamics of the pheromone are decoupled.

As classically done in the K-theory, see e.g. Panofsky et al. (1984) or Sportisse (2007), the diffusion flux due to fluctuations of the air motion caused by the turbulence of the atmosphere at non-resolved scales is closed using the Fick's law. This law relates the diffusion flux to the gradient of the concentration and leads to the diffusion term $-\mathbf{K}_{3D}\nabla_{3D}C$ with $\mathbf{K}_{3D}(\vec{x}, t)$ the diffusion tensor. In the present case, the diffusion tensor, a priori non-homogeneous, is assumed anisotropic in the $Oxy$ plane but isotropic along the $z$ axis. Thus, the diffusion tensor can be

$$ \text{written as: } \mathbf{K}_{3D} = \begin{pmatrix} K_{xx} & K_{xy} & 0 \\ K_{yx} & K_{yy} & 0 \\ 0 & 0 & K_{zz} \end{pmatrix}. $$

The advection flux of mass $\vec{u}_{3D}C$ is due to the wind and the gravity, and $\vec{u}_{3D}(\vec{x}, t)$ is the velocity of the pheromones in the atmosphere. Thus, the velocity of the pheromones can be written as $\vec{u}_{3D} = (u, v, -w_{set})$. The components $u$ and $v$ of the velocity are given by the wind velocity. Since the vertical wind velocity is assumed to be negligible and the pheromones are heavier than the air, the vertical component $-w_{set}$ is the settling velocity, i.e. the velocity at which the pheromones fall due to the gravity.

The dilution hypothesis and decoupling of the pheromone and atmosphere dynamics imply that the diffusion tensor $\mathbf{K}_{3D}$ and the wind field $\vec{u}_{3D}$ can be computed separately by a meteorological model and given as environmental parameters of the CTM model (34).
As the pheromones are settled under the effect of the gravity, they end up on the ground and a part of the pheromones are absorbed by the ground. To model this phenomenon, the following boundary condition on the vertical mass flux is commonly used:

$$ (35) \qquad (K_{zz} \partial_z C + w_{set} C)|_{z=0} = w_{dep} C|_{z=0} $$

with $w_{dep}$ the deposition velocity.

### A.1.2. Depth integration of the 3D CTM model.
The 3D CTM (34) is now integrated over a layer of the atmosphere $z \in [h_{-}; h_{+}]$ in order to obtain the 2D model used and presented in Sec. 2.1:

$$ (36) \qquad \int_{h_{-}}^{h_{+}} \frac{\partial C}{\partial t} - \nabla_{3D} \cdot (\mathbf{K}_{3D} \nabla_{3D} C) + \nabla_{3D} \cdot (\vec{u}_{3D} C) \, dz = \int_{h_{-}}^{h_{+}} S \, dz $$

In the following, the total concentration of pheromone in the layer is denoted

$$ c(x, y, t) = \int_{h_{-}}^{h_{+}} C(x, y, z, t) dz. $$

The integral of the different terms are now computed in order to obtain a 2D model.

The integral of the time derivative becomes simply:

$$ (37) \qquad \int_{h_{-}}^{h_{+}} \partial_t C dz = \partial_t c $$

The integral of the advection term can be written as:

$$ \int_{h_{-}}^{h_{+}} \nabla_{3D} \cdot (\vec{u}_{3D} C) dz = \underbrace{\int_{h_{-}}^{h_{+}} \partial_x (uC) + \partial_y (vC) dz}_{=l_1} + \underbrace{\int_{h_{-}}^{h_{+}} \partial_z (-w_{set} C) dz}_{=l_2} $$

Assuming that the velocity components $u$ and $v$ are constant or near the average value over the layer of interest, the integral $l_1$ becomes:

$$ (38) \qquad l_1 = \partial_x \left( \int_{h_{-}}^{h_{+}} u C dz \right) + \partial_y \left( \int_{h_{-}}^{h_{+}} v C dz \right) = \partial_x (u c) + \partial_y (v c) $$

The integral of the diffusion term can be written as:

$$ (39) \qquad - \int_{h_{-}}^{h_{+}} \nabla_{3D} \cdot (\mathbf{K}_{3D} \nabla_{3D} C) \, dz = - \underbrace{\int_{h_{-}}^{h_{+}} \partial_x (K_{xx} \partial_x C + K_{xy} \partial_y C) + \partial_y (K_{yx} \partial_x C + K_{yy} \partial_y C) \, dz}_{=l_3} - \underbrace{\int_{h_{-}}^{h_{+}} \partial_z (K_{zz} \partial_z C) \, dz}_{=l_4} $$

In the present case, the diffusion coefficients in the plane $Oxy$ $K_{xx}$, $K_{xy}$, $K_{yx}$ and $K_{yy}$ is assumed to not depend on the depth $z$ which is a common assumption. Thus, the integral $l_3$ becomes:

$$ (40) \qquad \begin{aligned} l_3 &= \partial_x \left( \int_{h_{-}}^{h_{+}} K_{xx} \partial_x C + K_{xy} \partial_y C dz \right) + \partial_y \left( \int_{h_{-}}^{h_{+}} K_{yx} \partial_x C + K_{yy} \partial_y C dz \right) \\ &= \partial_x (K_{xx} \partial_x c + K_{xy} \partial_y c) + \partial_y (K_{yx} \partial_x c + K_{yy} \partial_y c) \end{aligned} $$
From the integrals of the advection and the diffusion terms, it remains the vertical mass fluxes $l_2 - l_4$ that become:

$$ (41) \qquad l_2 - l_4 = - \int_{h_{-}}^{h_{+}} \partial_z (K_{zz} \partial_z C + w_{set} C) dz = -[K_{zz} \partial_z C + w_{set} C]_{z=h_{-}}^{z=h_{+}} $$

In the present study, these vertical fluxes, which model the quantity of pheromones leaving the layer of interest, are closed using a 1st order term, i.e. by a reaction term weighted by a loss rate $\tau_{loss}(x, y)$:

$$ (42) \qquad l_2 - l_4 = -[K_{zz} \partial_z C + w_{set} C]_{z=h_{-}}^{z=h_{+}} \approx \tau_{loss} c $$

At last, the total amount of pheromone emitted by the insects in the layer of interest is denoted:

$$ (43) \qquad \int_{h_{-}}^{h_{+}} S(x, y, z, t) dz = s(x, y, t) $$

By introducing the expression of the integrals (37), (38), (40), (42) and (43) in the 3D CTM integrated over a layer of the atmosphere (36), we obtain the 2D CTM:

$$ (44) \qquad \frac{\partial c}{\partial t} - \nabla \cdot (\mathbf{K} \nabla c) + \nabla \cdot (\vec{u} c) + \tau_{loss} c = s $$

with the $\vec{u}$ and $\mathbf{K}$ resp. the wind field and diffusion tensor in the $Oxy$ plane.

### A.1.3. Calibration of the loss coefficient $\tau_{loss}$.
In the derivation of the 2D CTM from the 3D CTM, the vertical fluxes are closed by a reaction term weighted by a loss coefficient $\tau_{loss}$. This coefficient depends on the vertical fluxes. Thus, it depends on the atmosphere (through the vertical diffusion), on the pheromone's molecules (through the settling velocity) and on the vegetation (through the deposition velocity). However, the estimation of $\tau_{loss}$ is not straightforward. In this section, a method of estimation of $\tau_{loss}$ is proposed. In the present study, $\tau_{loss}$ is estimated by comparing the solution of the 2D CTM derived in the previous section and the integral over the layer of the Gaussian plume, the solution of the 3D CTM (34) in a particular case.

The following assumptions are made. It is assumed that the diffusion tensor is isotropic and constant ($\mathbf{K}_{3D} = \mathbf{K}_{ld}$) and that the wind is unidirectional and constant ($\vec{u}_{3D} = (u, 0, -w_{set})$). Moreover, it is assumed that the source term is of the form $S(\vec{x}) = Q\delta(x)\delta(y)\delta(z - H)$ with Q a constant emission rate and H the height at which the source of emission is located. The Gaussian plume is the solution of the 3D CTM under these assumptions and is of the following form, see more details and a review on the topic in Stockie (2011):

$$ (45) \qquad C_{gp}(r, y, z) = \frac{Q}{u} \frac{1}{\sqrt{4\pi r}} \mathcal{G}_h(r, y, z) \mathcal{G}_v(r, y, z) \mathcal{F}_{bc}(r, y, z) $$

with:

$$ \mathcal{G}_h(r, y, z) = \exp \left( -\frac{y^2}{4r} \right), $$

$$ \mathcal{G}_v(r, y, z) = \exp \left( -\frac{w_{set}(z-H)}{2K} - \frac{w_{set}^2 r}{4K^2} \right) \left( \exp \left( -\frac{(z-H)^2}{4r} \right) + \exp \left( -\frac{(z+H)^2}{4r} \right) \right), $$

$$ \mathcal{F}_{bc}(r, y, z) = \frac{4w_0\pi r}{K} \exp \left( -\frac{w_{set}(z-H)}{2K} - \frac{w_{set}^2 r}{4K^2} \right) \exp \left( \frac{w_0(z+H)}{K} + \frac{w_0^2 r}{K^2} \right) \text{erfc} \left( \frac{z+H}{2\sqrt{r}} + \frac{w_0\sqrt{r}}{K} \right), $$

the change of variable $r = \frac{1}{u} \int_0^x K(\xi) d\xi$, erfc the complementary error function and $w_0 = w_{dep} - \frac{w_{set}}{2}$.

In order to get realistic estimations of the loss coefficient $\tau_{loss}$, the solution of the 2D CTM (44) is empirically compared with the integral over the layer of the Gaussian plume $\int_{h_{-}}^{h_{+}} C_{gp}(r, y, z) dz$ under the assumptions presented above in a realistic case.
Based on the average of the wind speed of the meteorological data, the value $u = 12m.s^{-1}$ is used in this section. The settling velocity is set as $w_{set} = 2.7 \times 10^{-3}m.s^{-1}$ as used in (Stockie, 2011). It is assumed that the source of pheromone is located at $H = 2m$ and the layer of interest is 2m wide centred at the height of the pheromone source ($h_{\pm} = H \pm 1m$). The values of deposition velocity $w_{dep}$ with respect to the vegetation used here are found in McMahon et al. (1979). From this realistic case, it results the value of $\tau_{loss}$ given in Tab. 1.

## A.2. List of additional population dynamics-informed regularization terms

We detail here additional regularization terms that could be used during the inference process. These regularization terms are not included in the results, since they were not included in the numerical exploration of the BI-DA problems, in order to limit the number of possibilities and to focus on the known behavior of the FallArmy Worm (FAW), our insect of interest. However, we indicate here a large list of features that could be introduced in the population dynamic-informed regularization term (27). The corresponding gradient of the terms are indicated in section A.3.2.

### A.2.1. Logarithmic barrier term for exclusion zone.
As mentioned in Sec. 3.1.1, the biology-based Tikhonov enables to take into account the favorite habitat of the insects. Nonetheless, this approach does not preclude the possibility of locating insects in areas that they will avoid.

To do so, one can consider a $j_{reg}$ based on a logarithmic barrier function, see e.g. Nesterov (2018). This approach, and barrier function in general, is in particular used in penalty method that replaces an inequality constraint (e.g. $h(s; x, y, t) < b \forall t \in [0; T] \forall (x, y) \in \Omega_c$) by a penalizing term (e.g. $\int_0^T \int_{\Omega_c} -\epsilon \log(b - h(s; x, y, t))dxdydt$ with $\Omega_c$ the constraint domain $\epsilon$ a very small constant) easier to handle in optimization problems.

In the present case, this penalizing term aims to preclude locating insects in areas that they will avoid or where they will not emit pheromones, i.e. $s = 0$ or very small on the exclusion zone $\Omega_{exc}$ derived from the biological knowledge of the insect. For instance, it is known that the FAW females emit pheromones in a vertical position on vegetation that are favorable for the growth of the larva. Hence, the FAW females will not emit pheromones in some areas that are not suitable such as roads.

Therefore, the following regularization/penalization term is considered:

$$ (46) \qquad j_{reg,exc}(s) = \int_{0}^{T} \int_{\Omega_{exc}} -\epsilon \log(\epsilon - s(x, y, t)^2) dxdydt $$

Let us note that this logarithmic barrier function is smooth (and hence differentiable) and tends to the indicator function (0 if $s = 0$ locally, $+\infty$ otherwise) when $\epsilon \to 0$.

### A.2.2. Population-dynamics models.
Recall that the quantity of pheromone emitted $s$ can be written as $s = pq$ with $q$ the quantity of pheromone emitted per insect, which is assumed to be a biological knowledge, and $p$ the density of insect. Moreover, the density $p$, and consequently $s$, follows a population dynamic model that can be written in the generic way as following:

$$ (47) \qquad \begin{aligned} & \partial_t p + \nabla \left( \sum_i F_i(p) \right) - \sum_i R_i(p) = 0 \\ & \Leftrightarrow \frac{1}{q} \partial_t s - \frac{\partial_t q}{q^2} s + \nabla \left( \sum_i \mathcal{F}_i(s) \right) + \sum_i R_i(s) = 0 \end{aligned} $$

with $F_i$ and $R_i$ are resp. flux and reaction terms that model the behavior of the insect species of interest, and $F_i(p = \frac{s}{q}) = \mathcal{F}_i(s)$ and $R_i(p = \frac{s}{q}) = \mathcal{R}_i(s)$.
As mentioned earlier, only a little number of terms have been introduced in this study. In this section, we aim to indicate a large list of features and resulting population dynamics terms that could be introduced in the population dynamic-informed regularization term (27).

**Diffusion term.** The most classical term modelling population dynamics is the diffusion term:

$$ (48) \qquad F_d(p; x, y, t) = -\sigma(x, y)\nabla p(x, y, t) \Leftrightarrow \mathcal{F}_d(s) = -\frac{1}{q(t)}\sigma(x, y)\nabla s(x, y, t) $$

with $\sigma [m^2.s^{-1}]$ the diffusion coefficient of the population, depending notably on the specie.

This term models the spread of the population at a macroscopic scale caused by a motion of each individual in a random direction with a probability independent of other individuals or the environment. This diffusion term can be understood as a continuous and large population limit of a discrete microscopic unbiased random walk process with a fixed probability of jump, see e.g. Nisbet et al. (1982) or Murray (2007).

Simple models based on the diffusion term (48) are used to model insects population dynamics in Kareiva (1983) where the author provides constant values of $\sigma$ for 12 species based on experimental data.

**Advection term.** This term models the advection motion of the population toward a determined direction, see e.g. Murray (2007):

$$ (49) \qquad F_a(p; x, y, t) = \vec{v}_i(x, y, t)p(x, y, t) \Leftrightarrow \mathcal{F}_a(s) = \frac{1}{q(t)}\vec{v}_i(x, y, t)s(x, y, t) $$

with $\vec{v}_i [m.s^{-1}]$ the velocity field of the population (assumed here independent of $p$).

For instance, one can model the attraction of the population toward the preferred habitat by choosing the velocity field $\vec{v}_i(x, y) = \nabla\phi(x, y)$, where $\phi$ denotes the iso-distance line function from this habitat. This term can also be understood as an attractive force produced by the potential function $\phi$ representing how favorable the environment is, see e.g. Shigesada et al. (1979) and in the review Okubo (1986). The two following terms are specific cases of advection terms.

**Quorum sensing term.** This term models the advection motion of the population but toward a direction depending on the population density:

$$ (50) \qquad F_{qs}(p; x, y, t) = \vec{v}_i(p; x, y, t)p(x, y, t) $$

with $\vec{v}_i(p)$ the density-dependent velocity field of the population.

This term enables to model the advection and spread of the population where the individuals interact with each other. The most classical density-dependent velocity field is

$$ (51) \qquad \vec{v}_i(p) = -\sigma_{qs} \nabla p $$

with $\sigma_{qs} [m^2.m^2.s^{-1}]$ the diffusion coefficient per density unit. This velocity field enables to model a motion due to avoidance of overcrowded places and population pressure, see e.g. Murray (2007). It can also be understood as the diffusion of the population at a macroscopic scale caused by a motion of each individual in a biased random direction with a probability increasing with the density of individuals, see e.g. Nisbet et al. (1982) and Murray (2007). This density-dependent velocity field (51) leads to the term in variable $s$

$$ (52) \qquad \mathcal{F}_{qs}(s) = -\frac{\sigma_{qs}}{q^2}s\nabla s $$

Avoidance of overcrowded places and of population pressure is a common behavior among insects, see e.g. Murray (2007), Shigesada (1980) and references therein. For instance, Shigesada (1980) proposes and applies a model with the term (50) and the velocity field (51) to an
ant lion population. This model is based on experiments highlighting population pressure in this population.

**Pheromone sensing term.** This term models the advection motion of the population but toward a direction depending on the concentration of pheromone:

$$ (53) \quad F_{ps}(p; x, y, t) = \chi p(x, y, t) \nabla c(s; x, y, t) \Leftrightarrow \mathcal{F}_{ps}(s; x, y, t) = \frac{\chi}{q(t)} s(x, y, t) \nabla c(s; x, y, t) $$

with $\chi [m^2.g^{-1}.m^2.s^{-1}]$ the chemotactic gain constant ($\chi \in \mathbb{R}^{+*}$).

The insects employ pheromones to attract individuals of the same sex, such as in the case of sex pheromones for FAW, or of both sexes, as is the case with aggregation pheromones for the red palm weevils. To reach the source of the pheromone, the insects moves in the direction of the gradient of concentration of pheromones. This motion induced by a chemical signal is commonly known as chemotaxis, especially in cellular biology. A Keller-Segel-type system (Keller et al., 1970, 1971) can be used to model at the same time both the population dynamics (here of the insects) and the dynamics of the chemical attractant (here the pheromone), see e.g. Murray (2007) and Perthame (2006).

For instant, Ganga Prasath et al. (2023) proposes a Keller-Segel-type system to model the dynamics of swarm of bees communicating using a pheromone.

**Logistic term.** Introduced by Verhulst in the 19th century, the logistic growth is used to model a self-limiting birth-death process, see e.g. Murray (2007), Perthame (2006) and reference therein: (54)

$$ R_l(p; x, y, t) = \frac{r}{\kappa} p(x, y, t)(\kappa - p(x, y, t)) \Leftrightarrow R_l(s; x, y, t) = \frac{1}{q(t)^2} \frac{r}{\kappa} s(x, y, t)(\kappa q(t) - s(x, y, t)) $$

with $r [s^{-1}]$ the linear reproduction rate ($r \in \mathbb{R}^{+*}$) and $\kappa [m^{-2}]$ the carrying capacity of the environment ($\kappa \in \mathbb{R}^{+*}$).

This term enables to model a saturation phenomena that arises when the maximal occupation density (the carrying capacity) is exceeded ($p > \kappa$, the death process overcome the birth process), and invasion front. This saturation is for instance due to intra-specific competition for a limited resource that is consumed faster than it can renewed.

The equation obtained when considering diffusion and logistic source is the well-known Fisher-KPP equation, introduced in RA Fisher (1937) to model the spread of favorable gene and first used for population dynamics problems in Skellam (1951). A Fisher-KPP-type system is for instance used in Takahashi et al. (2005) to model the dynamics of mosquitoes populations where the clear water (required for egg-laying) and the amount of nutrient (human blood) are limited, potentially leading to saturation phenomenon.

**Allee effect term.** The Allee effect is characterized by the decrease of the population due to a small or sparse population. The following term is an extension of the logistic term (54) which additionally models the Allee effect, see e.g. a review in Sun (2016) and reference therein:

$$ (55) \quad R_{Ae}(p; x, y, t) = \frac{r}{\kappa_+ \kappa_-} p(x, y, t)(\kappa_+ - p(x, y, t))(p(x, y, t) - \kappa_-) \\ \Leftrightarrow R_{Ae}(s; x, y, t) = \frac{1}{q(t)^3} \frac{r}{\kappa_+ \kappa_-} s(x, y, t)(q(t)\kappa_+ - s(x, y, t))(s(x, y, t) - q(t)\kappa_-) $$

with $r [s^{-1}]$ the linear reproduction rate ($r \in \mathbb{R}^{+*}$), $\kappa_+ [m^{-2}]$ the carrying capacity of the environment ($\kappa_+ \in \mathbb{R}^{+*}$) and $\kappa_- [m^{-2}]$ the Allee threshold ($\kappa_- \in \mathbb{R}^{+*}$ and $\kappa_- < \kappa_+$) that is the critical density under which the population decrease.
This effect is e.g. due to the lack of interactions between individuals of the population. In insect populations, the Allee effect commonly occurs since lower densities reduce the chance of an individual to locate a mate for purpose of reproduction.

Since the Allee effect can lead to the extinction of the insect population, it is often used in pest control strategy, see e.g. Boukal et al. (2009).

### A.2.3. Non-local terms.
The non-local terms are generalization of the previously presented terms. These non-local terms enable to model non-local interaction and long range diffusion, leading for instance to spatial patterns, see e.g. Furter et al. (1989), Genieys et al. (2006) or Murray (2007) and reference therein. Non-local terms arise especially when the individuals of the population communicate with each other remotely, e.g. using chemical or visual signals, or when the individuals interact with the environment not local but surrounding.

For example, one can extend the logistic term (54) into a non-local version:

$$ (56) \qquad \begin{aligned} R_{nll}(p; x, y, t) &= \frac{r}{\kappa} p(x, y, t)(\kappa - k * p(x, y, t)) \\ \Leftrightarrow R_{nll}(s; x, y, t) &= \frac{1}{q(t)^2} \frac{r}{\kappa} s(x, y, t)(\kappa q(t) - k * s(x, y, t)) \end{aligned} $$

with the convolution product $k * s(x, y, t) = \int_{\Omega} k(x - x', y - y')s(x', y', t)dx'dy'$ and k the convolution kernel.

This implies that individuals are not competing for the resource locally available, as with the term (54), but for the resource available in a vicinity. In this term, the convolution product can be seen as an average weighted by the kernel k that defines how the individuals interact with the vicinity and how large is this vicinity. Let us note that the local term (54) is retrieved if k is the Dirac kernel $\delta$. Similarly, one can consider a non-local Allee effect term to model the decrease of the population due to a small or sparse population in a vicinity. For instance, Britton (1989) proposed a non-local version of the Fisher-KPP equation based on this non-local logistic term.

As mentioned before, a Keller-Segel type system can be used to model both the dynamics of the population and the dynamics of the chemical attractant, in which the quantity of chemical attractant emitted depends on the population density. In special cases, the equation of the chemical attractant dynamics can be analytically solved. The concentration of the attractant can then be expressed as a convolutional product between the density of the population and a kernel describing the attractant dynamics, which leads to a non-local flux term in the population dynamics model, see e.g. Perthame (2006).

For example, in the present case, one can approximate the pheromone dynamics, described in the direct model (1), by a Gaussian plume $c_{gp}(s) = g_{gp} * s$ with $g_{gp}$ the Gaussian plume kernel. Injecting this approximation in the pheromone sensing term (53) leads to a non-local pheromone sensing term:

$$ (57) \qquad \mathcal{F}_{nlps}(s; x, y, t) = \frac{\chi}{q(t)} s(x, y, t) \nabla (g_{gp} * s(x, y, t)) $$

Non-local terms of the shape $R_{lr}(p) \propto k * p$ can be used to model long range diffusion processes. The range of the diffusion processes depends on the choice of the kernel k and its moments, see e.g. Murray (2007).
### A.3. On solving the optimization problem (3)

In this section, we present all the calculation details for solving the optimization problem (3) presented in Sec. 2.2. In the present context, the cost function $j$ is optimized using gradient descent or proximal gradient methods. Recall that the gradient of $j$, or at least of its differentiable part, is given by

$$ (58) \qquad \nabla_{s}j = \nabla_{s}j_{obs} + \nabla_{s}j_{reg} $$

with $j_{reg} = \sum_{i} \alpha_{reg,i}j_{reg,i} \Rightarrow \nabla_{s}j_{reg} = \sum_{i} \alpha_{reg,i} \nabla_{s}j_{reg,i}$.

In Sec. 2.2.2, it is mentioned that $\nabla_{s}j_{obs} = -c^*$ with $c^*$ the adjoint state, solution of the adjoint model (6). In a first part, the derivation of the adjoint model and of the expression (7) are detailed. Moreover, related models used in this study, such that the one-sensor adjoint model, will be detailed.

In Secs. 3.1 and A.2, several Biology-Informed differentiable $j_{reg,i}$ are proposed. In a second part, the gradients of these $j_{reg,i}$ are calculated.

In Sec. 3.1.2, LASSO and group-LASSO regularization terms are proposed. However, these $j_{reg}$ are not differentiable. Thus, a proximal gradient algorithm is used in this study. In a third part, we present the proximal-gradient algorithm and its application to the LASSO and group-LASSO regularization terms, called ISTA algorithm.

#### A.3.1. Derivation of the adjoint model and related.
In this section, $\mathcal{K}_T$ denotes the control space and $\mathcal{K}_T = \mathcal{L}^2([0; T]; \mathcal{K})$ with $\mathcal{K}$ the space of functions from $\Omega \subset \mathbb{R}^2$ to $\mathbb{R}$ which regularity in space is required in order to solve the direct model (1). We assume that the spaces $\mathcal{K}_T$ and $\mathcal{K}$ are Hilbert spaces and we have

$$ < f, g >_{\mathcal{K}_T} = \int_{0}^{T} < f(t, \cdot), g(t, \cdot) >_{\mathcal{K}} dt. $$

In a similar way, $\mathcal{V}_T = \{f; f \in \mathcal{L}^2([0; T]; \mathcal{V})$ and $\partial_t f \in \mathcal{L}^2([0; T]; \mathcal{V})\}$ denotes the space of the state variable, solution of the direct model and we have $< f, g >_{\mathcal{V}_T} = \int_{0}^{T} < f(t, \cdot), g(t, \cdot) >_{\mathcal{V}} dt$. Moreover, $\mathcal{D}$ denotes the observations space.

As mentioned previously, the gradient of the cost function is $\nabla_{s}j(s) = \nabla_{s}j_{obs}(s) + \alpha_{reg} \nabla_{s}j_{reg}(s)$ with the gradient of the observation term (4) that is: $\forall \delta s \in \mathcal{K}_T$

$$ (59) \qquad \begin{aligned} < \nabla_{s}j_{obs}(s), \delta s >_{\mathcal{K}_T} &= 2 < \mathbf{R}^{-1} (m(c(s)) - m^{\text{obs}}), \frac{dm}{dc}(c(s)) \cdot \frac{dc}{ds}(s) \cdot \delta s >_{\mathcal{D}} \\ &= 2 < \left( \frac{dm}{dc}(c(s)) \right)^{*} \cdot \mathbf{R}^{-1} (m(c(s)) - m^{\text{obs}}), \frac{dc}{ds}(s) \cdot \delta s >_{\mathcal{V}_T} \\ &= 2 < \left( \frac{dc}{ds}(s) \right)^{*} \cdot \left( \frac{dm}{dc}(c(s)) \right)^{*} \cdot \mathbf{R}^{-1} (m(c(s)) - m^{\text{obs}}), \delta s >_{\mathcal{K}_T} \end{aligned} $$

with $s \mapsto c(s)$ the operator that given $s$ returns the solution of the direct model (1) and $d_s c(s)$ the derivative of this operator with respect to $s$. However, the operator $(d_s c(s))^*$ may not be easy to calculate. The adjoint model provides a method to estimate the results of this operator.

**Linear tangent model.** In this section, the direct model is written as :

$$ (60) \qquad \frac{\partial c}{\partial t} + \mathcal{A}(c) = \mathcal{B}(s) $$

In our context, the differentiable operators are:

$$ (61) \qquad \mathcal{A}(c) = \nabla \cdot (\bar{u}c - \mathbf{K}\nabla c) + \tau_{\text{loss}}c \text{ and } \mathcal{B}(s) = s $$

with the boundary conditions:
* a null diffusive flux $K \nabla c \cdot \vec{n} = 0 \forall (x, y) \in \partial \Omega$,

* a null convective influx $\vec{u}c \cdot \vec{n} = 0 \forall (x, y) \in \partial \Omega \cap \{(x, y)|\vec{u}(x, y, t) \cdot \vec{n} < 0\} \forall t \in ]0; T]$ with $\vec{n}$ the outgoing normal vector.

Denoting $\omega = d_s c \cdot \delta s \in \mathcal{V}_T$, the linear tangent model is then:

$$ (62) \qquad \frac{\partial \omega}{\partial t} + \frac{dA}{dc}(c(s)) \cdot \omega = \frac{dB}{ds}(s) \cdot \delta s $$

with the initial condition $\omega(t=0) = 0$ (since $c(t=0) = 0 \forall x, y \forall s$).

Let us note that, in the present case, $d_c A = A$ since $A$ is a linear operator. By applying the scalar product with a test function $z \in \mathcal{V}_T$ and integration by part, we then have:

$$ (63) \qquad < \frac{\partial z}{\partial t} - \left( \frac{dA}{dc}(c(s)) \right)^* \cdot z, \omega >_{\mathcal{V}_T} - \left[ < z, \omega >_{\nu} \right]_{t=0}^{t=T} = - < \frac{dB}{ds}(s) \cdot \delta s, z >_{\mathcal{V}_T} \forall z \in \mathcal{V}_T $$

**Adjoint model.** On the other hand, based on the expression of the gradient of $j$ (58) and of $j_{obs}$ (59), the gradient of $j$ satisfies: $\forall \delta s \in \mathcal{K}_T$

$$ (64) \qquad < \nabla_{sj}, \delta s >_{\mathcal{K}_T} = < \nabla_{sj_{reg}}, \delta s >_{\mathcal{K}_T} + 2 < \left( \frac{dm}{dc}(c(s)) \right)^* \cdot \mathbf{R}^{-1} (m(c(s)) - m^{\text{obs}}), \omega >_{\mathcal{V}_T} $$

By introducing the linear tangent model (63) in the equation of gradient (64), we then have: $\forall \delta s \in \mathcal{K}_T$ and $\forall z \in \mathcal{V}_T$

$$ (65) \qquad < \nabla_{sj}, \delta s >_{\mathcal{K}_T} = < \nabla_{sj_{reg}}, \delta s >_{\mathcal{K}_T} - < \frac{dB}{ds}(s) \cdot \delta s, z >_{\mathcal{V}_T} $$

$$ (66) \qquad - < \frac{\partial z}{\partial t} - \left( \frac{dA}{dc}(c(s)) \right)^* \cdot z - 2 \left( \frac{dm}{dc}(c(s)) \right)^* \cdot \mathbf{R}^{-1} (m(c(s)) - m^{\text{obs}}), \omega >_{\mathcal{V}_T} \\ + \left[ < z, \omega >_{\nu} \right]_{t=0}^{t=T} $$

Thus, the adjoint model aims to find a particular test function $c^*$ of $\mathcal{V}_T$, so-called adjoint state, that satisfies:

$$ (67) \qquad \frac{\partial c^*}{\partial t} - \left( \frac{dA}{dc}(c(s)) \right)^* \cdot c^* = 2 \left( \frac{dm}{dc}(c(s)) \right)^* \cdot \mathbf{R}^{-1} (m(c(s)) - m^{\text{obs}}) $$

with a null final condition $c^*(t=T) = 0 \forall (x, y) \in \Omega$.

For the differential operator $A$ given by Eq. (61), the adjoint operator $(d_c A)^*$ is simply $(d_c A)^* \cdot z = -\nabla \cdot (K^T \nabla z) - \nabla(\vec{u}z) + (\nabla \cdot \vec{u})z + \tau_{loss}z$ and the boundary conditions of the adjoint model are:

* a null diffusive flux $K^T \nabla c^* \cdot \vec{n} = 0 \forall (x, y) \in \partial \Omega$,

* a null outgoing convective flux $\vec{u}c^* \cdot \vec{n} = 0 \forall (x, y) \in \partial \Omega \cap \{(x, y)|\vec{u}(x, y, t) \cdot \vec{n} > 0\} \forall t \in ]0; T]$ with $\vec{n}$ the outgoing normal vector.

Let us note that by definition, the adjoint model is a linear model.

**Computation of the gradient.** By introducing the adjoint state as test function in the equation (65), it returns that $\nabla_{sj}$ satisfies the following equation:

$$ (68) \qquad < \nabla_{sj}(s), \delta s >_{\mathcal{K}_T} = < \nabla_{sj_{reg}}(s), \delta s >_{\mathcal{K}_T} - < \left( \frac{dB}{ds}(s) \right)^* \cdot c^*(c(s)), \delta s >_{\mathcal{V}_T} \forall \delta s \in \mathcal{K}_T $$

with $(d_s B)^* = \text{Id}$ for the differential operator $B$ given by (61). Thus, it results the expression of the gradient:

$$ (69) \qquad \nabla_{sj}(s) = \nabla_{sj_{reg}}(s) - \left( \frac{dB}{ds}(s) \right)^* \cdot c^*(c(s)) $$
To estimate the gradient at each optimization iteration, it requires to solve in order the direct model (1) and then the adjoint model (6).

**One-sensor and trajectory-integrated one-sensor adjoint models.** Recall that the one-sensor adjoint model is mostly obtained by splitting the observation operator $c \mapsto m(c)$ into the sum of one-sensor observation operator $m = \sum_i m_i$. By introducing this observation operator splitting in the adjoint model (6), the adjoint model can be written the following way:

$$\partial_t c^* + \nabla \cdot (\mathbf{K}^T \nabla c^*) + \nabla (\vec{u} c^*) - (\nabla \cdot \vec{u}) c^* - \tau_{loss} c^* = \sum_i \sum_k \left( \frac{dm_k}{dc}(c(s)) \right)^* \cdot 2\mathbf{R}^{-1} \cdot (m_i(c(s)) - m_i^{obs})$$

Assuming that $(\frac{dm_k}{dc}(c(s)))^* \cdot 2\mathbf{R}^{-1} \cdot (m_i(c(s)) - m_i^{obs}) = 0$ for $i \neq k$ and as the left hand side is linear with respect to $c^*$, the adjoint model written (A.3.1) can be split into the one-sensor adjoint models (8). Let us note that, as mentioned in Sec. 2.2.3, the hypothesis $(\frac{dm_k}{dc}(c(s)))^* \cdot 2\mathbf{R}^{-1} \cdot (m_i(c(s)) - m_i^{obs}) = 0$ for $i \neq k$ is satisfied in our case if the sensors are independent ($\mathbf{R}^{-1}$ is block-diagonal) and if two sensors are not located in the same place.

Recall from Sec. 2.2.3 that the resulting one-sensor adjoint state is the contribution of a sensor to the gradient $\nabla_{s_{jobs}}$.

In this study, the sequence $(s_j)_{j \in \mathbb{N}}$ given by the trajectory taken by the optimization algorithm using different DA strategies are compared through one-sensor adjoint state-based criteria. The main one-sensor adjoint state-based criteria used, introduced in Sec. 2.4.3, is the maximal total contribution $\| \sum_j \eta_j c_i^*(s_j) \|_{L^\infty(\Omega \times [0, T])}$, with $(\eta_j)_{j \in \mathbb{N}}$ the sequence of the step sizes and $c_i^*(s_j)$ is the solution of the one-sensor adjoint model:

$$\partial_t c_{i,j}^* + \nabla \cdot (\mathbf{K}^T \nabla c_{i,j}^*) + \nabla (\vec{u} c_{i,j}^*) - (\nabla \cdot \vec{u}) c_{i,j}^* - \tau_{loss} c_{i,j}^* = (d_c m_i(c(s_j)))^* \cdot 2\mathbf{R}^{-1} (m_i(c(s_j)) - m_i^{obs})$$

Thus, by multiplying with the step size and summing over the sequences the previous one-sensor adjoint model and due the linearity of the left hand side (that is always true by definition of the adjoint model), $c_{i,tot}^* = \sum_j \eta_j c_i^*(s_j)$ is solution of the following trajectory-integrated one-sensor adjoint model:

$$(70) \qquad \partial_t c_{i,tot}^* + \nabla \cdot (\mathbf{K}^T \nabla c_{i,tot}^*) + \nabla (\vec{u} c_{i,tot}^*) - (\nabla \cdot \vec{u}) c_{i,tot}^* - \tau_{loss} c_{i,tot}^* \\ = \sum_j (\eta_j (d_c m_i(c(s_j)))^* \cdot 2\mathbf{R}^{-1} (m_i(c(s_j)) - m_i^{obs}))$$

This trajectory-integrated one-sensor adjoint model enable to compute the criteria without having to compute $c_i^*(s_j)$ for each $s_j$.

Let us note that, if the one-sensor observation operator $c \mapsto m$ is linear (which is the case in the present study), the right hand side of the trajectory-integrated one-sensor adjoint model

$$\text{RHS} = \sum_j (\eta_j (d_c m_i(c(s_j)))^* \cdot 2\mathbf{R}^{-1} (m_i(c(s_j)) - m_i^{obs})) \text{ becomes}$$

$$\text{RHS} = (d_c m_i)^* \cdot 2\mathbf{R}^{-1} \cdot (\sum_j \eta_j (m_i(c(s_j)) - m_i^{obs})).$$

Moreover, if the direct model is linear with respect to $s$ (which is the case in the present study), the right hand side of the trajectory-integrated one-sensor adjoint model even becomes

$$\text{RHS} = (d_c m_i)^* \cdot 2\mathbf{R}^{-1} \cdot (m_i(c(\sum_j \eta_j s_j)) - (\sum_j \eta_j) m_i^{obs}).$$

### A.3.2. Computation of the gradients $\nabla_{s_{jreg,i}}$
In this section, we calculate the gradient of the differentiable regularization terms $j_{reg,i}$ introduced in Sec. 3.1 and in Sec. A.2. Recall that $\nabla_s$ is the gradient operator with respect to $s$ and $\nabla_{s_{jreg,i}}(s)$ is the function given by $\frac{d_{j_{reg,i}}(s)}{ds} \cdot \delta s = \langle \nabla_{s_{jreg,i}}(s), \delta s \rangle \forall \delta s$.
**Gradient of the Tikhonov regularization term.** The gradient of the Tikhonov regularization term (23), presented in Sec. 3.1.1 is simply:

$$ \nabla_s j_{reg,T}(s) = 2\mathbf{C}^{-1}(s - s_b) $$

**Gradient of the logarithmic barrier regularization term.** The gradient of the logarithmic barrier regularization term (46), presented in Sec. A.2.1 is simply:

$$ \nabla_s j_{reg,exc}(s) = \frac{2\epsilon}{\epsilon - s^2} \mathbb{1}_{\Omega_{exc}} $$

**Gradient of the population dynamics-informed regularization term.** Recall that the population dynamics-informed regularization term (27), presented in Sec. 3.1.3, can be written as:

$$ \begin{aligned} j_{reg,PD}(s) &= \|\partial_t s - \frac{\partial_t q}{q}s + q\nabla\left(\sum_i \mathcal{F}_i(s)\right) + q\sum_i \mathcal{R}_i(s)\|_{L^2}^2 \\ &= \langle\partial_t s - \frac{\partial_t q}{q}s + q\nabla\left(\sum_i \mathcal{F}_i(s)\right) + q\sum_i \mathcal{R}_i(s), \partial_t s - \frac{\partial_t q}{q}s \\ &\quad + q\nabla\left(\sum_i \mathcal{F}_i(s)\right) + q\sum_i \mathcal{R}_i(s)\rangle \end{aligned} $$

and hence, its differentiable is

$$ \frac{d j_{reg,PD}(s)}{ds} \cdot \delta s = 2 \langle \partial_t s - \frac{\partial_t q}{q} s + q \nabla \left( \sum_i \mathcal{F}_i(s) \right) + q \sum_i \mathcal{R}_i(s) , \partial_t \delta s - \frac{\partial_t q}{q} \delta s + q \nabla \left( \sum_i \frac{d \mathcal{F}_i(s)}{ds} \cdot \delta s \right) + q \sum_i \frac{d \mathcal{R}_i(s)}{ds} \cdot \delta s \rangle $$

Hence $\frac{d j_{reg,PD}(s)}{ds} \cdot \delta s = \langle \nabla_s j_{reg,PD}(s), \delta s \rangle$ with :

$$ \nabla_s j_{reg,PD}(s) = 2 \left( \partial_t - \frac{\partial_t q}{q} + q \sum_i \frac{d \nabla \mathcal{F}_i(s)}{ds} + q \sum_i \frac{d \mathcal{R}_i(s)}{ds} \right)^* \left( \partial_t s - \frac{\partial_t q}{q}s + q \nabla \left( \sum_i \mathcal{F}_i(s) \right) + q \sum_i \mathcal{R}_i(s) \right) $$

Recall that the adjoint operator is a linear operator: $(F_1 + F_2)^* = F_1^* + F_2^*$
We know have to calculate the adjoint operators $(\partial_t - \frac{\partial_t q}{q})^*$, $(\frac{d \nabla \mathcal{F}_i(s)}{ds})^*$ and $(\frac{d \mathcal{R}_i(s)}{ds})^*$ for the $\mathcal{F}_i$ and $\mathcal{R}_i$ presented in Sec. 3.1.3 and in Sec. A.2.2. Let us note that no boundary conditions are imposed in the estimation of the population dynamics model. Hence, all the boundary terms are taken into account in the adjoint operators.

**Temporal dynamics terms.** Here, we calculate the operator $(\partial_t - \frac{\partial_t q}{q})^*$. By applying the scalar product with a test function $\varphi$, we have:

$$ \langle \partial_t \delta s - \frac{\partial_t q}{q} \delta s, \varphi \rangle = \langle -\partial_t \varphi - \frac{\partial_t q}{q} \varphi, \delta s \rangle + \int_{\Omega} [\varphi \delta s]_{t=0}^{t=T} dX $$

Therefore, the operator $(\partial_t - \frac{\partial_t q}{q})^*$ is the operator such that

$$ \varphi \mapsto -\partial_t \varphi - \frac{\partial_t q}{q} \varphi \quad \forall (x, y, t) \in \Omega \times ]0; T[ \\ \text{for } \varphi \text{ such that } \varphi = 0 \quad \forall (x, y, t) \in \Omega \times \{0; T\} $$
**Diffusion term.** Here, we calculate the operator $(\frac{d\nabla F_i(s)}{ds})^*$ for the diffusion operator $F_i = F_d$ (48). By applying the scalar product with a test function $\varphi$, we have:

$$ \begin{aligned} \langle \frac{d\nabla F_d(s)}{ds} \cdot \delta s, \varphi \rangle &= \langle -\nabla \cdot (\frac{1}{q}\sigma\nabla\delta s), \varphi \rangle \\ &= \langle -\nabla \cdot (\frac{1}{q}\sigma^T\nabla\varphi), \delta s \rangle + \int_0^T \int_{\partial\Omega} \left(\frac{1}{q}\sigma^T\nabla\varphi\delta s - \frac{1}{q}\varphi\sigma\nabla\delta s\right) \cdot \vec{n}d\Gamma dt \end{aligned} $$

Therefore, the operator $(\frac{d\nabla F_d(s)}{ds})^*$ is the operator such that

$$ \varphi \mapsto -\nabla \cdot \left(\frac{1}{q}\sigma^T\nabla\varphi\right) \quad \forall (x, y, t) \in \Omega \times [0; T] $$

for $\varphi$ such that $\frac{1}{q}\sigma^T\nabla\varphi \cdot \vec{n} = 0$ and $\varphi = 0 \quad \forall (x, y, t) \in \partial\Omega \times [0; T]$

**Advection term.** Here, we calculate the operator $(\frac{d\nabla F_i(s)}{ds})^*$ for the advection operator $F_i = F_a$ (49). By applying the scalar product with a test function $\varphi$, we have:

$$ \begin{aligned} \langle \frac{d\nabla F_a(s)}{ds} \cdot \delta s, \varphi \rangle &= \langle \nabla \cdot (\frac{1}{q}\vec{v}_i\delta s), \varphi \rangle \\ &= \langle \nabla \cdot (\frac{1}{q}\vec{v}_i)\varphi - \nabla \cdot (\frac{1}{q}\vec{v}_i\varphi), \delta s \rangle + \int_0^T \int_{\partial\Omega} \frac{1}{q}\vec{v}_i\varphi\delta s \cdot \vec{n}d\Gamma dt \end{aligned} $$

Therefore, the operator $(\frac{d\nabla F_a(s)}{ds})^*$ is the operator such that

$$ \varphi \mapsto \nabla \cdot (\frac{1}{q}\vec{v}_i)\varphi - \nabla \cdot (\frac{1}{q}\vec{v}_i\varphi) \quad \forall (x, y, t) \in \Omega \times [0; T] $$

for $\varphi$ such that $\varphi = 0 \quad \forall (x, y, t) \in \partial\Omega \times [0; T]$

**Quorum sensing term.** Here, we calculate the operator $(\frac{d\nabla F_i(s)}{ds})^*$ for the quorum sensing operator $F_i = F_{qs}$ (52). By applying the scalar product with a test function $\varphi$, we have:

$$ \begin{aligned} \langle \frac{d\nabla F_{qs}(s)}{ds} \cdot \delta s, \varphi \rangle &= \langle -\nabla \cdot \left(\frac{\sigma_{qs}}{q^2}(s\nabla\delta s + \nabla s\delta s)\right), \varphi \rangle \\ &= \langle -\nabla \cdot \left(\frac{\sigma_{qs}}{q^2}(s\nabla\varphi - \nabla s\varphi)\right) - \nabla \cdot \left(\frac{\sigma_{qs}}{q^2}\nabla s\right)\varphi, \delta s \rangle \\ &\quad + \int_0^T \int_{\partial\Omega} \frac{\sigma_{qs}}{q^2}(s\nabla\varphi\delta s - s\varphi\nabla\delta s + \nabla s\varphi\delta s) \cdot \vec{n}d\Gamma dt \end{aligned} $$

Therefore, the operator $(\frac{d\nabla F_{qs}(s)}{ds})^*$ is the operator such that

$$ \varphi \mapsto -\nabla \cdot \left(\frac{\sigma_{qs}}{q^2}(s\nabla\varphi - \nabla s\varphi)\right) - \nabla \cdot \left(\frac{\sigma_{qs}}{q^2}\nabla s\right)\varphi \quad \forall (x, y, t) \in \Omega \times [0; T] $$

for $\varphi$ such that $\varphi = 0$ and $\frac{\sigma_{qs}}{q^2}s\nabla\varphi \cdot \vec{n} = 0 \quad \forall (x, y, t) \in \partial\Omega \times [0; T]$

**Non local pheromone sensing term.** Here, we calculate the operator $(\frac{d\nabla F_i(s)}{ds})^*$ for the non-local pheromone sensing operator $F_i = F_{nlps}$ (53). By applying the scalar product with a test function
φ, we have:

$$ 
\begin{aligned}
\left\langle \frac{d\nabla \mathcal{F}_{nlps}(s)}{ds} \cdot \delta s, \varphi \right\rangle &= \left\langle \nabla \cdot \left( \frac{\chi}{q} (s\nabla(g_{gp} * \delta s) + \nabla(g_{gp} * s)\delta s) \right), \varphi \right\rangle \\
&= \left\langle \bar{g}_{pg} * \nabla \cdot \left( \frac{\chi}{q} s\nabla\varphi \right) - \nabla \cdot \left( \frac{\chi}{q} \nabla(g_{gp} * s)\varphi \right) + \nabla \cdot \left( \frac{\chi}{q} \nabla(g_{gp} * s) \right) \varphi, \delta s \right\rangle \\
&\quad + \int_0^T \int_{\partial\Omega} \frac{\chi}{q} (s\varphi\nabla(g_{gp} * \delta s) - s\nabla\varphi(g_{gp} * \delta s) + \nabla(g_{gp} * s)\varphi\delta s) \cdot \vec{n} d\Gamma dt
\end{aligned}
$$

with $\bar{g}_{pg}$ the function such that $\bar{g}_{pg}(x, y, t) = g_{pg}(-x, -y, -t)$. Therefore, the operator $(\frac{d\nabla \mathcal{F}_{nlps}(s)}{ds})^*$ is the operator such that

$$ \varphi \mapsto \bar{g}_{pg} * \nabla \cdot \left( \frac{\chi}{q} s \nabla \varphi \right) - \nabla \cdot \left( \frac{\chi}{q} \nabla (g_{gp} * s) \varphi \right) + \nabla \cdot \left( \frac{\chi}{q} \nabla (g_{gp} * s) \right) \varphi \quad \forall (x, y, t) \in \Omega \times [0; T] $$
for $\varphi$ such that $\varphi = 0$ and $\frac{\chi}{q} s \nabla \varphi \cdot \vec{n} = 0 \quad \forall (x, y, t) \in \partial\Omega \times [0; T]$

**Logistic term.** Here, we calculate the operator $(\frac{d\mathcal{R}_i(s)}{ds})^*$ for the logistic operator $\mathcal{R}_i = \mathcal{R}_I$ (54). By applying the scalar product with a test function $\varphi$, we have:

$$ 
\begin{aligned}
\left\langle \frac{d\mathcal{R}_I(s)}{ds} \cdot \delta s, \varphi \right\rangle &= \left\langle \frac{1}{q^2} \frac{r}{\kappa} (\kappa q - 2s) \delta s, \varphi \right\rangle \\
&= \left\langle \frac{1}{q^2} \frac{r}{\kappa} (\kappa q - 2s) \varphi, \delta s \right\rangle
\end{aligned}
$$

Therefore, the operator $(\frac{d\mathcal{R}_I(s)}{ds})^*$ is the operator such that

$$ \varphi \mapsto \frac{1}{q^2} \frac{r}{\kappa} (\kappa q - 2s) \varphi \quad \forall (x, y, t) \in \Omega \times [0; T] $$

**Allee effect term.** Here, we calculate the operator $(\frac{d\mathcal{R}_i(s)}{ds})^*$ for the Allee effect operator $\mathcal{R}_i = \mathcal{R}_{Ae}$ (55). By applying the scalar product with a test function $\varphi$, we have:

$$ 
\begin{aligned}
\left\langle \frac{d\mathcal{R}_{Ae}(s)}{ds} \cdot \delta s, \varphi \right\rangle &= \left\langle \frac{1}{q^3} \frac{r}{\kappa_+ \kappa_-} (-3s^2 + 2q(\kappa_+ + \kappa_-)s - q^2\kappa_+\kappa_-) \delta s, \varphi \right\rangle \\
&= \left\langle \frac{1}{q^3} \frac{r}{\kappa_+ \kappa_-} (-3s^2 + 2q(\kappa_+ + \kappa_-)s - q^2\kappa_+\kappa_-) \varphi, \delta s \right\rangle
\end{aligned}
$$

Therefore, the operator $(\frac{d\mathcal{R}_{Ae}(s)}{ds})^*$ is the operator such that

$$ \varphi \mapsto \frac{1}{q^3} \frac{r}{\kappa_+ \kappa_-} (-3s^2 + 2q(\kappa_+ + \kappa_-)s - q^2\kappa_+\kappa_-) \varphi \quad \forall (x, y, t) \in \Omega \times [0; T] $$

### A.3.3. On the proximal gradient algorithm for the LASSO regularization term.
As mentioned in Sec. 3.1, if the regularization terms $j_{reg}$ are convex and differentiable, then the cost $j$ is convex and differentiable and the optimization problem (3) is solved using a gradient descent algorithm. However, as mentioned in Secs 3.1 and 3.1.2, the differentiable framework can be extended to a non-smooth and subdifferentiable framework to consider LASSO-type regularization terms. In this context, the gradient descent algorithm can be extended to the proximal-gradient algorithm, also called the Forward-Backward algorithm.

In this section, the proximal-gradient algorithm is briefly presented in the context of the optimization problem (3) with several $j_{reg}$ including the LASSO $j_{reg}$ (24) and group-LASSO $j_{reg}$ (25).

The proximal-gradient algorithm, see e.g. Parikh et al. (2014) is an extension of the gradient descent algorithm used when the function to minimize can be split in $j(s) = g(s) + h(s)$ with:
* g(s) a convex and differentiable function, which is in the present context the differentiable part of j, i.e. the observation term $j_{obs}$ and the differentiable $j_{reg}$,

* h(s) a convex and subdifferentiable function, which is in the present case the LASSO or group-LASSO $j_{reg}$.

For such cost function, the optimum $s_a$ is a fixed point $s_a = prox_{\lambda h}(s_a - \lambda \nabla g(s_a))$ with $prox_f$ the proximal operator applied to the function $f$ and $\lambda > 0$ the gradient descent step size.

The proximal-gradient algorithm is then a fixed point algorithm based on the iteration loop:

$$s_{k+1} = prox_{\lambda h}(s_k - \lambda \nabla g(s_k))$$

In the present case, the gradient $\nabla g$ is computed using the expression (7). For the LASSO regularization term $h = \alpha_{reg,LASSO}j_{reg,LASSO}$ given by Eq. (24), the proximal operator $prox_{\lambda h}$ is

$$prox_{\lambda h}(s) = sign(s)max(0, |s| - \lambda\alpha_{reg,LASSO}).$$

For the group-LASSO regularization term $h = \alpha_{reg,gLASSO}j_{reg,gLASSO}$ given by Eq. (25), the proximal operator $prox_{\lambda h}$ is

$$prox_{\lambda h}(s) = sign(s)max\left(0, 1 - \frac{\lambda\alpha_{reg,gLASSO}}{\|s\|_{L^2(0;T)}}\right) s$$

with $\|s\|_{L^2(0;T)} = (\int_0^T s(t)^2 dt)^{1/2}$. As highlighted by the proximal operator $prox_{\lambda h}$, this algorithm acts as an iterative soft-threshold, shrinking to 0 the value of $s$ lower than the threshold defined by $\lambda\alpha_{reg,LASSO}$ for LASSO and the value of $s$ at places where the total amount $\|s\|_{L^2(0;T)}$ is lower than the threshold $\lambda\alpha_{reg,gLASSO}$ for group-LASSO. Consequently, the adaptation of the proximal-gradient algorithm to the LASSO regularization term is commonly called *Iterative Shrinkage-Threshold Algorithm* (ISTA).

## A.4. Comparison of the LASSO and group-LASSO regularization terms in the Toy case

In this section, we briefly compare the LASSO (24) and group-LASSO (25) regularization terms in the Toy case presented in Sec. 2.3.2. Recall that the effect of the LASSO regularization term in the Toy case has already been highlighted in Sec. 3.4.

### A.4.1. Visual comparison of estimated pheromone emissions $s_a$. Recall that the LASSO regularization ('LASSO', centre of Fig. 18) implies that the retro-plume is less diffuse at both initial ($t = 0$s, top row) and final ($t = 20$s, bottom row) times. On the other hand, the group-LASSO ('group LASSO', right of Fig. 18) implies that the retro-plume is even less diffuse at initial time but more diffuse at final time. This is due to the fact that the group-LASSO regularization aims to threshold all the pheromone emissions at location where the total emissions are too low, unlike the LASSO regularization that shrinks all pheromone emissions that are too low. To compensate these retro-plumes less diffuse at initial time but more diffuse at final time, the pheromone emissions are higher at initial time but lower at final time with the group-LASSO regularization than with the LASSO regularization. Let us note that more spurious emissions far from the target support are recovered with the group-LASSO regularization than with the LASSO regularization.

### A.4.2. RMSE criteria. Recall that with the LASSO regularization (pink curves on Fig. 19), the total error $RxMSE_{tot}$ is sensibly lower and the error on the target support $RxMSE_{supp}$ similar to the ones obtained without regularization (blue curve) at the beginning of the time window and up to $t \approx 17.5$s. After that time, with the LASSO regularization, both $RxMSE_{tot}$ and $RxMSE_{supp}$ significantly increase and become higher than without regularization. With the group-LASSO regularization (yellow curves on Fig. 19), the errors follow similar trend but are nearly always
<figure><img src="image_22.png" /><figcaption>Figure 18 - Optimal sources estimated with different regularization terms. Optimal quantity of pheromone emitted $s_a$ obtained by inference without regularization (left), with the LASSO regularization (24) (centre) and with the group-LASSO regularization (25) (right) for the Toy case. The optimal sources are indicated at $t = 0.0s$ (upper row) and $t = 20.0s$ (lower row).</figcaption></figure>

slightly lower. Let us note that, at the very beginning of the time window (between $t = 0s$ and $t = 2.5s$), $RxMSE_{supp}$ is higher with group-LASSO regularization than with LASSO regularization. This comes from the fact that the group-LASSO regularization thresholds all the emissions, including at the initial time, located at certain positions in the target supports.

#### A.4.3. Presence prediction, false discovery and miss rates.
As already anticipated with the visual comparison, the presence prediction (upper row of Fig. 20) is less diffuse at initial time, similarly diffuse at middle time and more diffuse at final time with the group-LASSO regularization (yellow curves) than with the LASSO regularization (pink curve). Recall the with the LASSO regularization (pink curves), the miss rate $\rho_{fa}$ and false discovery rate $\rho_{fp}$ (bottom row of Fig. 20) are quite stable until $t \approx 17.5s$, time at which $\rho_{fa}$ significantly increases. With the group-LASSO regularization (yellow curves), the trend is very similar, but two little difference can be noted. First, as observed on the errors, at the initial times, $\rho_{fa}$ is slightly higher and $\rho_{fp}$ is slightly lower with the group-LASSO regularization than with the LASSO regularization. With the LASSO regularization, $\rho_{fa}$ is even null for a short amount of time, while it is always non-zero with the group-LASSO regularization. Also, one can note that $\rho_{fa}$ significantly increases a bit later with the group-LASSO regularization than with the LASSO regularization.
<figure><img src="image_23.png" /><figcaption>**Figure 19** – Evolution of total and support RMSE of $s_a$ against $s_t$. Time evolution of the root space-mean square errors between the optimum $s_a$ obtained by inference with several regularization strategies and the target $s_t$, on the whole domain ($R_x\text{MSE}_{\text{tot}}$, Eq. (14), upper panel) and on the target support ($R_x\text{MSE}_{\text{supp}}$, Eq. (15), lower panel) in the Toy case: with no regularization term (no reg, blue), with the LASSO regularization term (24) (LASSO, pink) and with the group-LASSO regularization term (25) (group LASSO, yellow)</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Figure 20** – Presence prediction and time evolution of false discovery and miss rates.Prediction of insects presence with several regularization strategies: with no regularization term (no reg, blue), with the LASSO regularization term (24) (LASSO, pink) and with the group-LASSO regularization term (25) (group LASSO, yellow). Top panel: map of presence prediction at three different times. The presence is obtained from the reconstructed source $s_a$ by computing a level set with Eq. (16). The true presence (i.e. the target support) is indicated in black.</figcaption></figure>

Lower panel. Time evolution of the miss rate $\rho_{fa}$ (solid lines), see Eq. (17), and false discovery rate $\rho_{fp}$ (dashed lines), Eq. (18).