4.1 Dispersal model

We assume that in the absence of stimuli, mosquitoes disperse following a Brownian motion (Dumont et al. 2011), which means that they fly around randomly, independently from each other, changing directions continuously, with no preferential direction of flow. This is the classical model for the dynamics of gas molecules where changes in directions are due to random collisions between molecules.

The density \(\rho\) of Brownian particles at point \(\mathbf{x} \in \mathbb{R}^2\) and at time \(t \in \mathbb{R}\) satisfies the diffusion equation

\[\begin{equation} \frac{\partial\rho}{\partial t} = D \nabla^2 \rho, \end{equation}\] where \(D\) is the diffusivity constant and represents the rate at which particles spread and \(\nabla^2 = \nabla \cdot \nabla\) is the Laplace operator (i.e., the divergence of the gradient of a scalar function).

Assuming that individuals are released from the origin of spatial coordinates at the initial time \(t = 0\), the diffusion equation has the solution

\[\begin{equation} \rho(\mathbf{x}, t) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\Big\{-\frac12 \frac{\norm{\mathbf{x}}^2}{\sigma^2}\Big\} = \frac{1}{\sqrt{8\pi Dt}} \exp\Big\{-\frac{\norm{\mathbf{x}}^2}{8Dt}\Big\} \end{equation}\]

with \(\sigma^2 = 2nDt\), and \(n\) is the dimension of the space. This is a Gaussian density function centred at the release point and with a variance equal to \(4Dt\) in \(\mathbb{R}^2\).

Figure 4.1: Evolution of the probability density over time at 4 different distances from the release point.

Working the Brownian motion with discrete time and space leads to a Random Walk, in which the walker’s position after a large number of independent steps is distributed according to a Normal distribution of total variance \(\sigma^2 = \varepsilon^2\,t/\delta_t\), where \(\varepsilon\) is the step size in space, \(t\) is the total time elapsed and \(\delta_t\) is the time step.

Thus, the diffusivity of the process in \(\mathbb{R}^2\) is given by \[\begin{equation} D = \frac{\varepsilon^2}{4\delta_t}. \end{equation}\]

4.2 Survival model

We assume that the probability that an individual survives any given day is a constant that we call Probability of Daily Survival (PDS) \(\pi_s\). As such, the probability of surviving up to day \(t\) after release is \[\begin{equation} \tag{4.1} \pi_s^t = \exp\{t\log{\pi_s}\} = \exp\{-\beta t\}, \end{equation}\] with \(\beta = - \log{\pi_s} > 0\) and \(t > 0\).

This formulation allows to work in continuous time, not only in discrete days.

From a survival perspective, let \(T_0\) the time-to-natural-death in the absence of traps. Equation (4.1) gives the survival function \(S(t) = P(T_0 > t)\) which decreases exponentially at an unknown constant rate. This corresponds to a constant hazard rate equal to: \[ h_0(t) = -S'(t) / S(t) = \beta. \]

Thus, while we don’t observe any deaths directly, we observe captures in traps, which work as censoring for this time-to-event process. Indeed, if individuals were captured at some point, they must have been alive at that point.

However, we don’t observe all individuals alive at some point (only the captured fraction). Meaning that we can’t analyse this survival process independently from the capture process. The standard approach was to consider the fraction of captured individuals as a constant proportion of the total. But it is unclear whether this is a sensible or valid hypothesis.

4.3 Probability of capture

The dispersal process is only partially observed through the capture of a fraction of the individuals in a set of traps located at fixed positions in the area.

We assume that, in the absence of prior death or capture in the same or other traps, the probability of capture of a released individual is proportional to the density \(\rho\) at the location \(x_i\) of the trap \(i\) at day \(t\) and the survival probability.

Why did I state the proportionality to the survival probability? Furthermore, shouldn’t we take into account a range of influence of the trap, rather than only the density at the specific point? Or is this a acceptable approximation?

Let \(T \in \{0, 1, \ldots, I\}\) a categorical variable describing the fate of an individual, which can be either \(T = 0\) if the individual is not captured during the experiment or \(T = i\) if the individual is captured in trap \(i\).

Let \(X_i > 0 \,|\, T = i\) be the capture time of an individual, conditional to capture on trap \(i\)

\[\begin{align} \tilde\pi_i(t) & = P(X_i = t \,|\, X_i \geq t, T = i) \propto \pi_s^t\,\rho(x_i, t) \\ & = \alpha_i t^{-1/2} \exp\big\{ - \beta t - \gamma_i / t \big\},\quad t > 0 \end{align}\] with \(\alpha_i = k_i/\sqrt{\pi 8D} > 0\), \(\beta = - \log{\pi_s} > 0\) and \(\gamma_i = \norm{\mathbf{x}}^2/8D > 0\).

The normalising constant \(k_i\) ensures that the unconditional distribution sums to 1.

Here I was thinking in terms of daily probabilities in order to avoid integrating continous times. But it may be better to think in terms of hazards.

Again, this is conditionally to the absence of prior capture in the same or other traps. The quantity and location of other traps, as well as the time elapsed, would necessarily influence this conditional capture probability.

The unconditional probability is then

\[\begin{align} \pi_i(t) & = P(X_i = t) = P(X_i = t \,|\, X_i \geq t, T = i) (1 - P(\text{prior capture})) \\ & = \tilde\pi_i(t) (1 - P(\text{prior capture})) \\ & = \end{align}\]

Essentially, I think that this problem can be formulated in terms of a survival analysis with competing risks.

Since the capture of an individual on different days are mutually exclusive events, the cumulative probability of capture \(\pi_i^c(t) = \sum_{t=0}^T \pi_i(t)\) is the probability of capture at any time since release and up to time \(t\).

\[\begin{equation} \tag{4.2} \pi_i^c(t) = \alpha_i \sum_{x=0}^t x^{-1/2} \exp\big\{ - \beta x - \gamma_i / x \big\}, \end{equation}\] with \(\alpha_i = k_i/\sqrt{\pi 8D} > 0\) and \(\gamma_i = \norm{\mathbf{x}}^2/8D > 0\).

Note that in continuous time we would get an integral, which has a analytical indefinite result as:

\[\begin{multline} \int x^{-1/2} \exp(-b x - c/x)\, \mathrm{d}x = \\ \frac{ \sqrt{\pi} e^{-2 \sqrt{b} \sqrt{c}} \big( -\text{erf}(\frac{\sqrt{c} - \sqrt{b} x}{\sqrt{x}}) + e^{4 \sqrt{b} \sqrt{c}} (\text{erf}(\frac{\sqrt{b} x + \sqrt{c}}{\sqrt{x}}) - 1) + 1 \big) }{ 2 \sqrt{b} } + \text{constant}, \end{multline}\] with \(b>0, c>0\), where \(\text{erf}(x) = \frac2{\sqrt\pi}\int_0^x e^{-x^2}\,\mathrm{d}x, \; x > 0\) is the error function.

This yields the following analytical formula for the cumulative capture rate for a trap \(i\) in the absence of other captures: \[\begin{multline} \tag{4.3} \pi_i^c(t) = \frac{\alpha_i}{2} \sqrt{\frac\pi{\beta}} e^{-2 \sqrt{\beta \gamma_i}} \Big[ -\text{erf}\big(\frac{\sqrt{\gamma_i} - \sqrt{\beta} t}{\sqrt{t}}\big) + e^{4 \sqrt{\beta \gamma_i}} \Big( \text{erf}\big(\frac{\sqrt{\gamma_i} + \sqrt{\beta} t}{\sqrt{t}}\big) - 1 \Big) + 1 \Big], \end{multline}\]

Figure 4.2: Daily probability of capture for each trap using the simulated parameters and uniform trap attractiveness. The scale parameter \(k\) is arbitrary.

Figure 4.2 shows that each trap has a different capture rate with a maximum at a different moment in time, depending on the distance from the release point. The number, location and dates of captures will hopefully be informative enough to identify the diffusivity, survival and attractiveness parameters.

4.4 Observation model

We observe the numbers \(Y\) of captured individuals in traps \(i = 1, \ldots, I\) installed at survey days \(j = 1, \ldots, J_i\). The survey periods for trap \(i\) span from day \(t_{ij}^0\) to \(t_{ij}^1\).

Assuming that \(N\) individuals were released at day \(t = 0\), and that each of them has a cumulative capture probability on trap \(i\) at survey \(j\) conditional to the absence of other captures described by (4.3), the number of captures \(Y_{ij}\) follows a Binomial distribution with parameters \(N\) and \(\tilde{p}_{ij} = \pi_i^c(t_{ij}^1) - \pi_i^c(t_{ij}^0)\).

However, each of the \(N\) individuals can either end up in only one of the surveys \(i, j\) or not being captured in the course of its life with probability \(p_0\). Thus, the outcome vector \((N - \sum{Y_{ij}}, Y_{ij})\) follows a Multinomial distribution with \(N\) trials, \(K = 1 + \sum_{i = 1}^I J_i\) mutually exclusive events and probabilities \(p_k\) where, \[ p_0 + \sum_{i = 1}^I \sum_{j = 1}^{J_i} p_{ij} = 1 \]

This allows normalising the conditional probabilities \(\tilde{p}\). I can do this across traps, but not over time, because time is not exchangeable.