Loading [MathJax]/jax/output/SVG/jax.js
Case report

Laparoscopic uterosacral nerve block: A fertility preserving option in chronic pelvic pain

  • Chronic pelvic pain (CPP) can cause extreme physical distress in women and has widespread socio-economic consequences. Nerve root blocks have become a safe and effective treatment modality in multiple specialties in both the diagnosis and treatment of pain. We describe a novel technique of a laparoscopic uterosacral nerve block (USNB) and demonstrate its effectiveness in the treatment of a complex case of CPP. USNB has potential diagnostic, prognostic and therapeutic implications. It should therefore be considered as part of the multi-disciplinary management of women with CPP of suspected uterine origin such as adenomyosis, degenerating fibroids or following myomectomy.

    Citation: Benjamin P Jones, Srdjan Saso, Timothy Bracewell-Milnes, Jen Barcroft, Jane Borley, Teodor Goroszeniuk, Kostas Lathouras, Joseph Yazbek, J Richard Smith. Laparoscopic uterosacral nerve block: A fertility preserving option in chronic pelvic pain[J]. AIMS Medical Science, 2019, 6(4): 260-267. doi: 10.3934/medsci.2019.4.260

    Related Papers:

    [1] Yunfeng Liu, Guowei Sun, Lin Wang, Zhiming Guo . Establishing Wolbachia in the wild mosquito population: The effects of wind and critical patch size. Mathematical Biosciences and Engineering, 2019, 16(5): 4399-4414. doi: 10.3934/mbe.2019219
    [2] Luis Almeida, Michel Duprez, Yannick Privat, Nicolas Vauchelet . Mosquito population control strategies for fighting against arboviruses. Mathematical Biosciences and Engineering, 2019, 16(6): 6274-6297. doi: 10.3934/mbe.2019313
    [3] Bo Zheng, Wenliang Guo, Linchao Hu, Mugen Huang, Jianshe Yu . Complex wolbachia infection dynamics in mosquitoes with imperfect maternal transmission. Mathematical Biosciences and Engineering, 2018, 15(2): 523-541. doi: 10.3934/mbe.2018024
    [4] Bo Zheng, Lihong Chen, Qiwen Sun . Analyzing the control of dengue by releasing Wolbachia-infected male mosquitoes through a delay differential equation model. Mathematical Biosciences and Engineering, 2019, 16(5): 5531-5550. doi: 10.3934/mbe.2019275
    [5] Yun Li, Hongyong Zhao, Kai Wang . Dynamics of an impulsive reaction-diffusion mosquitoes model with multiple control measures. Mathematical Biosciences and Engineering, 2023, 20(1): 775-806. doi: 10.3934/mbe.2023036
    [6] Daiver Cardona-Salgado, Doris Elena Campo-Duarte, Lilian Sofia Sepulveda-Salcedo, Olga Vasilieva, Mikhail Svinin . Optimal release programs for dengue prevention using Aedes aegypti mosquitoes transinfected with wMel or wMelPop Wolbachia strains. Mathematical Biosciences and Engineering, 2021, 18(3): 2952-2990. doi: 10.3934/mbe.2021149
    [7] Mugen Huang, Moxun Tang, Jianshe Yu, Bo Zheng . The impact of mating competitiveness and incomplete cytoplasmic incompatibility on Wolbachia-driven mosquito population suppressio. Mathematical Biosciences and Engineering, 2019, 16(5): 4741-4757. doi: 10.3934/mbe.2019238
    [8] Diego Vicencio, Olga Vasilieva, Pedro Gajardo . Monotonicity properties arising in a simple model of Wolbachia invasion for wild mosquito populations. Mathematical Biosciences and Engineering, 2023, 20(1): 1148-1175. doi: 10.3934/mbe.2023053
    [9] Rajivganthi Chinnathambi, Fathalla A. Rihan . Analysis and control of Aedes Aegypti mosquitoes using sterile-insect techniques with Wolbachia. Mathematical Biosciences and Engineering, 2022, 19(11): 11154-11171. doi: 10.3934/mbe.2022520
    [10] Chad Westphal, Shelby Stanhope, William Cooper, Cihang Wang . A mathematical model for Zika virus disease: Intervention methods and control of affected pregnancies. Mathematical Biosciences and Engineering, 2025, 22(8): 1956-1979. doi: 10.3934/mbe.2025071
  • Chronic pelvic pain (CPP) can cause extreme physical distress in women and has widespread socio-economic consequences. Nerve root blocks have become a safe and effective treatment modality in multiple specialties in both the diagnosis and treatment of pain. We describe a novel technique of a laparoscopic uterosacral nerve block (USNB) and demonstrate its effectiveness in the treatment of a complex case of CPP. USNB has potential diagnostic, prognostic and therapeutic implications. It should therefore be considered as part of the multi-disciplinary management of women with CPP of suspected uterine origin such as adenomyosis, degenerating fibroids or following myomectomy.


    1. Introduction

    In recent years, the spread of chikungunya, dengue, and zika has become a major public health issue, especially in tropical areas of the planet [1, 7]. All those diseases are caused by arboviruses whose main transmission vector is the Aedes aegypti. One of the most important and innovative ways of vector control is the artificial introduction of a maternally transmitted bacterium of genus Wolbachia in the mosquito population (see [8, 22, 36]). This process has been successfully implemented on the field (see [19]). It requires the release of Wolbachia-infected mosquitoes on the field and ultimately depends on the prevalence of one sub-population over the other. Other human interventions on mosquito populations may require such spatial release protocols (see [2, 3] for a review of past and current field trials for genetic mosquito population modification). Designing and optimizing these protocols remains a challenging problem for today (see [17, 34]), and may be enriched by the lessons learned from previous release experiments (see [18, 26, 38])

    This article studies a spatially distributed model for the spread of Wolbachia-infected mosquitoes in a population and its success as far as non-extinction probabilities are concerned. We address the question of the release protocol to guarantee a high probability of invasion. More precisely, what quantity of mosquitoes need to be released to ensure invasion, if we have only one release point? What if we have multiple release points and if there is some uncertainty in the release protocol? We obtain lower bounds so as to quantify the success probability of spatial spread of the introduced population according to a mathematical model.

    We define here an ad hoc framework for the computation of this success probability. As a totally new feature added to the previous works on this topic (see [10, 15, 16, 21, 33, 37]), it involves space variable as a key ingredient. In this paper we provide quantitative estimate and numerical results in dimension $1$.

    It is well accepted that stochasticity plays a significant role in biological modeling. Probabilities of introduction success have already been investigated for genes or other agents into a wild biological population. The recent work [4] makes use of reaction-diffusion PDEs to describe the biological phenomena underlying sucessful introduction as cytoplasmic analogues of the Allee effect. The infection of the mosquito population by Wolbachia is seen as an "alternative trait", spreading across a population having initially a homogeneous regular trait. Other recent models have been proposed either to compute the invasion speed ([9]), or get an insight into the induced time dynamics of more complex systems, including humans or pathogens (see [14, 20]). In the mosquito part, models usually feature two stable steady states: invasion (the regular trait disappears) and extinction (the alternative trait disappears). Since this phenomenon is currently being investigated as a tool to fight Aedes transmitted diseases, the problem of determination of thresholds for invasion in this equation is of tremendous importance.

    The issue of survival probability of invading species has attracted a lot of attention by many researchers. Among such we may cite [6] and [31]. We stress, however, that this is not the direction followed in this paper. In the cited articles indeed, the basic underlying model is either a stochastic PDE or its discretization, and the uncertainty concerning the initial state is not considered.

    In other words, although in a deterministic model as ours one can in principle numerically check for a specific initial configuration whether the invasion by the Wolbachia-infected mosquitoes will be successful or not, in practice such a specific initial condition is subject to uncertainty, and therefore the uncertainty quantification of the success probability is a natural question.

    Our modeling goes as follows: We consider on a domain $\Omega \subseteq \mathbb{R}^d$ (usually $d \in \{1, 2\}$ and $\Omega = \mathbb{R}^d$), a frequency $p : \Omega \to [0, 1]$ that models the prevalence of the Wolbachia infection trait. More specifically, in the case of cytoplasmic incompatibility caused in Aedes mosquitoes by the endo-symbiotic bacterium Wolbachia, $p$ is the proportion of mosquitoes infected by the bacterium (e.g. $p=1$ means that the whole population is infected). Then, this frequency obeys a bistable reaction-diffusion equation. We aim at estimating the invasion success probability with respect to the initial data (= release profile).

    In [4, 32] it was obtained an expression for the reaction term $f$ in the limit Allen-Cahn equation

    $tpσΔp=f(p)$ (1)

    in terms of the following biological parameters: $\sigma$ diffusivity (in square-meters per day, for example), $s_f$ (effect of Wolbachia on fecundity, $=0$ if it has no effect); $s_h$ (strength of the cytoplasmic incompatibility, $=1$ if it is perfect); $\delta$ (effect on death rate, $d_i = \delta d_s$ where $d_s$ is the regular death rate without Wolbachia) and $\mu$ (imperfection of vertical transmission, expected to be small). It reads as follows:

    $f(p)=δdspshp2+(1+sh(1sf)(1μδ+μ))p+(1sf)1μδ1shp2(sf+sh)p+1.$ (2)

    Bistable reaction terms are such that $f < 0$ on $(0, \theta)$ and $f > 0$ on $(\theta, \theta_+)$. Usually, we consider $\theta_+ = 1$. This is the case if $\mu = 0$.

    The outline of the paper is the following. In the next section, we explain how to use a threshold property for bistable reaction-diffusion equation in order to obtain explicit sufficient conditions for invasion success of a release protocol (Theorem 2.1). In a relevant stochastic framework, we show in Section 2.2 how these conditions provide uncertainty quantification for invasion success when release locations are random. Thanks to this, we prove in Section 2.3 that if the release domain is wide enough (with an explicit bound), the success probability goes to $1$ as the number of releases goes to $+\infty$. Our main tool is the construction of compactly supported radially symmetric functions (in Section 2.4 for any dimension, and in Section 3 for the $1$-dimensional case) such that if the initial data is above one of such functions, then invasion occurs. Section 3 and the following are devoted to the one dimensional case. We prove in Section 4.1 that the sufficient conditions for invasion are very hard to meet with a single release point (Proposition 5), and this leads to considering multiple release locations. For a deterministic (Section 4.2, Lemma 4.1) and a stochastic (Sections 4.3 and 4.4, Proposition 7) set of release profiles, we give analytical formulae for uncertainty quantification. Numerical simulations illustrate these results in dimension $1$ in Section 5. We conclude in Section 6. Finally an appendix is devoted to the study of the minimization of the perimeter of release in one dimension.


    2. Setting the problem: How to use a threshold property to design a release protocol?


    2.1. The threshold phenomenon for bistable equations

    In Equation (1), we assume that

    ${ θ(0,1), f(0)=f(θ)=f(1)=0,f<0 on (0,θ),f>0 on (θ,1),10f(x)dx>0.$ (3)

    A consequence of this hypothesis is the existence of invading traveling waves. From now on, we denote $F$ the anti-derivative of $f$ which vanishes at $0$,

    $F(x):=x0f(y) dy.$ (4)

    Since we have assumed $F(1)>0$, by the bistability of the function $f$, there exists a unique $\theta_c\in(0,1)$ such that

    $F(θc)=θc0f(x)dx=0.$

    In all numerical simulations we use the following values taken from [20, 12, 27] for the Wolbachia and mosquito parameters:

    $ds=0.27day1, sf=0.1, μ=0, sh=0.8,δ=0.3/0.27=10/9 and σ=877m2.day1.$ (5)

    In particular we obtain the profiles for $f$ and its anti-derivative in Figure 1. In [20], the authors used the notations $\phi = 1 - s_f$, $\delta = d_i / d_s$, $u = 1 - s_h$ and $v= 1 - \mu$. They gave a range of values of these parameters for three Wolbachia strains, namely wAlbB, which has no impact on death ($\delta = 1$) but reduces fecundity, wMelPop which highly increases death rate but isn't detrimental to fecundity, and wMel which has a moderate impact on both. Values are given in Table 3 of the cited article (which contains also a parameter $r$, standing for differential vector competence of Wolbachia-infected mosquitoes for dengue, a feature we do not include in our modelling since we focus on the mosquito population dynamics), see the references therein for more details. According to the aforementioned references, the authors always assumed perfect CI and maternal transmission, that is, with our notations $s_h=1$ and $\mu=0$. Our notations mimic those of [4, 14], where they did not give as detailed tables for the parameters as in [20], although we refer the reader to the references they gave, which contain some quantitative estimations of these parameters. Our choices in (5) for $d_s, s_f$ and $\delta$ reflect the field data exposed in [12], for the (life-shortening) wMel strain in the context of the city of Rio de Janeiro, in Brazil.

    Figure 1. Profile of $f$ defined in (2) (left) and of its anti-derivative $F$ (right) with parameters given by (5).

    We will always assume $\mu=0$ (perfect vertical transmission) in the following. Complex dynamical behaviors can arise in the case when $\mu$ exceeds some threshold, as was proved in [39] for a system of two ordinary differential equations. For such values of $\mu$, in particular, population replacement may not be guaranteed by invasion success. Note however that our results apply when $\mu > 0$ is small. In this case the "invasion" state is not exactly $p=1$, but $p=p_+ (\mu) < 1$, because of the flaw in Wolbachia vertical (=maternal) transmission.

    Moreover, following estimates from [12, 35] for Aedes aegypti in Rio de Janeiro (Brazil), and general literature review and discussion in Section 3 of [27] we consider that mosquitoes spread at around $\sigma = 830 \mathrm{m}^2 / \text{day}$ (see the references given in [27] for more details). With these estimations of the parameters, the quantitative results we get are satisfactory because they appear to be relevant for practical purposes. For example, in order to get a significant probability of success, the release perimeter we find is around $595 \mathrm{m}$ wide (in one dimension). In the example from Figure 1, $\theta_c \simeq 0.36$.

    We say that a radially symmetric function $\phi$ on $\mathbb{R}^d$ is non-increasing if $\phi(x)=g(|x|)$ for some $g$ that is non-increasing on $\mathbb{R}^+$.

    The following result gives a criterion on the initial data to guarantee invasion.

    Theorem 2.1.  Let us assume that $f$ is bistable in the sense of (3). Then, for all $\alpha \in (\theta_c, 1]$ there exists a compactly supported, radially symmetric non-increasing function$v_\alpha(|x|)$, with $v_\alpha:\mathbb{R}_+ \to \mathbb{R}_+$ non-increasing, $v_\alpha (0) = \alpha$ (called "$\alpha$-bubble"), such that if $p$ is a solution of

    $tpσΔp=f(p),$ (6)
    $p(t=0,x)=p0(x)vα(|x|),$

    then $p \xrightarrow[t \to \infty]{} 1$ locally uniformly. Moreover, we can take Supp$(v_\alpha)=B_{R_\alpha}$ with

    $Rα=σinfρΓ1ρd(1ρ)21(α0(11ραx)df(x)dx)+,$ (7)

    where $\Gamma=\{\rho\in(0,1),\ \int_0^\alpha (1-\frac{1-\rho}{\alpha}x)^d f(x)dx >0\}$.

    In one dimension, we have the sharper estimate Supp$(v_\alpha)=[-L_\alpha,L_\alpha]$ with

    $Lα=σ2α0dvF(α)F(v).$ (8)

    Remark 1. Clearly, the set $\Gamma$ is nonempty. Indeed for $\rho \sim 1$,

    $α0(11ραx)df(x)dx>0,$

    since $F(\alpha) > 0$. However, it is hard to say more unless we consider a specific function $f$.

    (Sharp) threshold phenomena are well-known for bistable reaction-diffusion equations (see [11, 24, 25, 30, 40]). In Theorem 2.1, we use this property to derive the new formula (7), and (8), which are very useful to quantify invasion success uncertainty. We postpone to Section 2.4 the proof of this result for dimension $d \geq 1$, which is based upon an energy method developed in [25]. When $d=1$, we give an alternative proof using sharp critical bubbles and a result of [11] in Section 3.1. To the best of our knowledge, we give in Section 3.2 the first comparison between the two approaches.

    We recall the definition of a "ground state" as a positive stationary solution $v$ of (1), i.e.:

    $Δv=f(v)$

    that decays to $0$ at infinity. In dimension $d=1$ (and in some special cases in higher dimensions, see [25]), such a ground state is unique up to translations. When $d=1$ we denote $v_{\theta_c}$ the ground state which is maximal at $x=0$. It is symmetric decreasing and $v_{\theta_c} (0) = \theta_c$, which is consistent with the notation $v_{\alpha}$ in Theorem 2.1. Although we won't use it in the rest of the paper, we note that with a similar argument, we have a sufficient condition for the extinction:

    Proposition 1. In dimension $d=1$, let $p$ be a solution of equation (1), associated with the initial value $p_0$. If $p_0 < \theta$ or $p_0 < v_{\theta_c}(\cdot - \zeta)$ for some $\zeta \in \mathbb{R}$, then $p$ goes extinct: $p \xrightarrow[t \to \infty]{} 0$ uniformly on $\mathbb{R}$.


    2.2. The stochastic framework for release profiles

    When mosquitoes are released in the field, the actual profile of Wolbachia infection in the days right after the release is very uncertain. In order to quantify this uncertainty, we define in this section an adequate space of release profiles. The preexisting mosquito population is assumed to be homogeneously dense, at a level $N_0 \in \mathbb{R}_+$.

    From now on, we assume that we have fixed a space unit, so that we may talk of numbers or densities of mosquitoes without any trouble.

    We define a spatial process $X_{\cdot} (\omega)=X(\cdot,\omega) : \mathbb{R}^d \to \mathbb{R}_+$ as the introduced mosquitoes profile.

    We expect that the time dynamics of the infection frequency will be given by

    ${tpσΔp=f(p),p(t=0,τ;ω)=Xτ(ω)Xτ(ω)+N0.$ (9)

    We want to measure the probability of establishment associated with this set of initial profiles.

    Making use of Theorem 2.1, we want to give a lower bound for the probability of non-extinction (which is equivalent to the probability of invasion, by the sharpness of threshold solutions, as described in [24, 25]).

    An initial condition $X_{\tau}$ ensures non-extinction if

    $α(θc,1], τ0R, τRd, XτXτ+N0vα(τ+τ0),$

    where $v_{\alpha}$ is the "$\alpha$-bubble" used in Theorem 2.1.

    Now, we assume that we have a fixed number of mosquitoes to release, say $N$. When we release mosquitoes in the field (out of boxes), they will spread out to find vertebrates to feed on (if not fed in the lab prior to the release), to mate or to rest. Many environmental factors may influence their spread (see [23]). As a very rough estimate we consider that the distribution of the released mosquitoes can be described by a Gaussian. A Gaussian profile is typically the result of a diffusion process. However, we shall not use very fine properties of these profiles, and mainly focus on a "significant spread radius", so that this assumption is not too restrictive.

    Due to the above simplification, the set of releases profiles ("RP") for a total of $N$ mosquitoes at $k$ locations in a domain $[-L, L]^d$ is defined as

    $RPdk(N):={τNkki=1e(ττi)22σi(2πσi)d/2, with τi[L,L]d, σi[σ0ϵ,σ0+ϵ]},$ (10)

    where $\sigma_0$ is an estimated diffusion coefficient and $\epsilon > 0$ represents the uncertainty on this parameter ($\sigma_i$ is the "significant spread radius"). In other words, for any $i$ between $1$ and $k$, the release profile is locally at the $i$-th release point a centered Gaussian with fixed amplitude $N/k$ and variance $\sigma_i$.

    The basic requirement for a release profile is that $\int_{\mathbb{R}^d} X_{\tau} d \tau = N$. It is obviously satisfied for the elements in $RP_k^d(N)$.

    We use uniform measure on $\big( [-L, L]^d \times [\sigma_0 - \epsilon, \sigma_0 + \epsilon] \big)^k$ to equip $RP_k^d (N)$ with a probability measure, denoted by $\mathcal{M}$ in the following.

    According to our estimate, the success probability satisfies

    $P[Non-extinction after releasing N mosquitoes at k locations ]P[Xτ(ω) satisfies (NEC)],$

    where $X_{\tau} (\omega)$ is taken in $RP_k^d (N)$ according to the uniform probability measure.


    2.3. First result: Relevance of under-estimating success

    Though it may seem naive, our under-estimation by radii given in Theorem 2.1 is rather good, and this can be quantified in any dimension $d$. Indeed, in any dimension we can prove convergence of our under-estimation in (SP) to $1$ as the number of releases goes to infinity, if we fix the number of mosquitoes per release.

    More precisely, we define for a domain $\Omega\subset\mathbb{R}^d$,

    $Pdk(N,Ω):=M{(xi)1ik,α(θc,1),x0Ω,x0+BRαΩ and xx0+BRα, Nkki=1Gσ,d(xxi)α},$ (11)

    where $G_{\sigma,d} (y) = \frac{1}{(2 \pi \sigma)^{d/2}} e^{-|y|^2/ 2 \sigma}$ and $B_{R_{\alpha}}=B_{R_{\alpha}} (0)$ is the ball of radius $R_{\alpha}$, centered at $0$. Then, the probability of success of a random (in the sense of Section 2.2) $k$-release of $N$ mosquitoes in the $d$-dimensional domain $\Omega$ is bigger than $P^d_k (N, \Omega)$, because of Theorem 2.1.

    Fixing the number of mosquitoes per release and letting the number of releases go to $\infty$ yield:

    Proposition 2.  Let $1 > \alpha > \theta_c$, $N \geq N^* := (2 \pi \sigma)^{d/2} \frac{\alpha}{1 - \alpha} N_0$ and $\Omega \subset \mathbb{R}^d$ be a compact set containing a ball of radius $R_{\alpha}$. Then,

    $Pdk(kN,Ω)k1.$ (12)

    Proof. There are two ingredients for the proof: First, we minimize a Gaussian at $x$ on a ball centered at $x$ by its value on the border of the ball. Second, if we pick uniformly an increasing number of balls with fixed radius and center in a compact domain, then their union covers almost-surely any given subset (this second ingredient is connected with the well-known coupon collector's problem). Namely,

    $y2σlog(2)ey2/2σ1/2.$

    Now, when we pick uniformly in a compact set the centers of balls of fixed radius $\alpha$, the probability of covering a given subset $\Omega_c \subset \Omega$ increases with the number $k$ of balls. Therefore it has a limit as $k \to + \infty$. In fact, this limit is equal to $1$.

    One can prove this claim using the coupon collector problem (see the classical work [13] for the main results on this problem), after selecting a mesh for the compact domain $\Omega_c$. We take this mesh such that each cell has diameter less than $\sqrt{2 \sigma \log(2)}/2$, and positive measure. The domain $\Omega$ is compact, hence finitely many cells is enough. Picking the center of a random ball in a given cell of the mesh has probability $>0$, and we simply need to have picked one center in each element to be done. It remains to choose the (compact) set $\Omega_c = B_{R_{\alpha}} + x_0 \subset \Omega$ to conclude the proof.

    Remark 2. We could have been a little more precise, and get an estimate for the expected value of the number $k$ of small balls required to cover the domain. According to classical results [13] on the coupon collector problem, it typically grows as $N_c \log(N_c)$, where $N_c$ is the number of cells. If the domain $\Omega$ has diameter $R$, $N_c$ is typically $(2 R/\sqrt{2 \sigma \log(2)})^{d}$, in dimension $d$.

    Therefore we should expect $\mathbb{E} [k] \sim d \big( \frac{2 R}{\sqrt{2 \sigma \log(2)} }\big)^d \log( \frac{2 R}{\sqrt{2 \sigma \log(2)}})$, and for a typical release area $R$ should be of the same order as $R_{\alpha}$.

    In fact, any $N > 0$ enjoys the same property, but then we need to assume that each cell contains a large enough number of release points.

    Corollary 1. For any $N > 0$ and $\alpha \in (\theta_c, 1)$, for $\Omega \subset \mathbb{R}^d$ a compact set containing a ball of radius $R_{\alpha}$, then for any compact subset $\Omega_c \subset \Omega$ containing a ball of radius $R_{\alpha}$ we have

    $Pdk(kN,Ωc)k1.$

    Proof. Let $\iota = \lceil\!\!\lceil \frac{N^*}{N} \rceil\!\!\rceil$. With the same technique as for proving Proposition 2, we get a coupon collector problem where $\iota$ coupons of each kind must be collected, whence the result.


    2.4. Proof of invasiveness in Theorem 2.1 in any dimension

    We consider in this section the proof of Theorem 2.1 in any dimension. The case $d=1$ is postponed to the next section.

    We use an approach based on the energy as proposed by [25]. In the present context, the energy is defined by

    $E[u]=Rd(σ2|u|2F(u(x)))dx.$ (13)

    It is straightforward to see that if $p$ is a solution to (6), then the energy is non-increasing along a solution, i.e.,

    $ddtE[p]=Rd(σΔp+f(p))2 dx0.$

    Thus, $E[p](t)\leq E[p^0]$ for all nonnegative $t$ and for $p$ solution to (6). Moreover, Theorem 2 of [25] states that if $\lim_{t\to +\infty} E[p(t,\cdot)] <0$, then $p(t,\cdot)\to 1$ locally uniformly in $\mathbb{R}^d$ as $t\to +\infty$. Thus, since $t\mapsto E[p(t,\cdot)]$ is non-increasing, it is enough to choose $p^0$ such that $E[p^0]<0$ to conclude the proof of Theorem 2.1.

    For any $\alpha > \theta_c$, we construct $p^0(x)=v_\alpha(|x|)$ as defined in the statements of Theorem 2.1. To do so, consider the family of non-increasing radially symmetric functions, compactly supported in $B_{R_0}$ with $R_0>0$, indexed by a small radius $0 < r_0 < R_0$, defined by $\phi(r) = 1$ if $r \leq r_0$, $\phi(r) = \frac{R_0 - r}{R_0 - r_0}$ if $r_0 < r < R_0$, and $\phi(r) \equiv 0$ if $r > R_0$.

    For any $0 < r_0 < R_0$, $\phi$ is continuous and piecewise linear. We define $v_\alpha(r)=\alpha \phi(r)$, for $r\geq 0$. By the comparison principle, it suffices to find $(r_0, R_0)$ such that $E[\alpha \phi] < 0$ to ensure that $R_{\alpha} = R_0$ is suitable in Equation (7) of Theorem 2.1. To do so, we introduce

    $Jd(r0,R0,α,ϕ):=E[αϕ]|Sd1|=α2σ0rd1|ϕ(r)|2dr(rd0dF(α)+R0r0rd1αϕ(r)0f(s)dsdr).$ (14)

    Now, we use our specific choice of non-increasing radially symmetric function $\phi$. Introducing $\rho := r_0 / R_0$, and with obvious abuses of notation, $J_d$ stands again for

    $Jd(ρ,R0,α):=Rd0(σdR201ρd(1ρ)2F(α)ρdd1ραα0(11ραx)d1F(x)dx),$ (15)

    where $F$ is the antiderivative of $f$ (as introduced in (4)). After an integration by parts, we have

    $Jd(ρ,R0,α)=Rd0(σdR201ρd(1ρ)2α0(11ραx)df(x)dx).$

    We choose $\rho\in (0,1)$ such that

    $α0(11ραx)df(x)dx>0$ (16)

    Then the energy $J_d(\rho,R_0,\alpha)$ decreases to $- \infty$ with $R_0$ and is positive for $R_0 \to 0$, so the minimal scaling ensuring negative energy is obtained for some known value of $R_0 =: R^{(d)}_{\alpha} (\rho)$, such that $J_d (\rho, R^{(d)}_{\alpha} (\rho), \alpha) = 0$. Namely,

    $(R(d)α(ρ))2=σ1ρd(1ρ)21α0(11ραx)df(x)dx,$ (17)

    which is a rational fraction in $\rho$. Thus we recover formula (7) in Theorem 2.1 by minimizing with respect to those $\rho$ satisfying constraint (16).

    We examine in particular formula (17) in the case $d=1$. To do so, we introduce

    $U(α):=F(α)1αα0F(x)dx,V(α):=1αα0F(x)dx.$ (18)

    Since $F(x) \leq F(\alpha)$ for $x \leq \alpha$, we know that $U$ is positive and $V$ is increasing with respect to $\alpha$ ($V'(\alpha) = \frac{1}{\alpha} U(\alpha)$). Moreover, $V(\theta_c) < 0$. We get

    $R(1)α(ρ)=ασ(1ρ)(V(α)+ρU(α)),$ (19)

    under the constraint $V(\alpha) + \rho U(\alpha) > 0$. The optimal choice for $\rho$ is then $\rho^*_1 (\alpha) := \frac{1}{2} - \frac{1}{2} \frac{V(\alpha)}{U(\alpha)}$. It satisfies $V(\alpha) + \rho^*_1 (\alpha) U(\alpha) > 0$ since $U(\alpha) = F(\alpha) - V(\alpha) > 0$ and $F(\alpha) > 0$.

    Finally, $\rho^*_1$ corresponds to a minimal radius

    $R(1),α:=R(1)α(ρ1(α))=2σαU(α)F(α),$ (20)

    with $U(\alpha)$ as in (18).

    Remark 3. We emphasize that $R_\alpha$ quantifies the minimal radius which ensures invasion from level $\alpha$, in the sense that it provides an upper bound for it. However, we were not able to perform an analytical computation of the actual optimal radius (=support size) of a critical bubble.

    Remark 4. We note in passing that the same energy (13) appears for instance in the review paper [5] and in associated literature, but is used in a different spirit (stemming from statistical physics).

    Before restricting to dimension $1$ in the sequel, we end the general exposition in this section with a numerical illustration. In order to help the reader getting a clearer picture of the invasion problem we investigate in the present paper, Figure 2 displays the time dynamics of equation (1) in two spatial dimensions, with three different initial conditions. In this simulation we use the function $f$ defined in (2) with parameter values given in (5). It illustrates the fact that with a fixed number of release points taken uniformly in a rectangle, invasion typically appears only if the size of the rectangle is well chosen.

    Figure 2. Time dynamics with three different initial releases belonging to the set $RP_{50}^2(N)$ of (10), with $N/(N+N_0) = 0.75$. Integration is performed on the domain $[-L, L]$ with $L = 50 \textrm{km}$. The release box is plotted in dashed red on the first picture of each configuration. Left: Release box $[-2 L/3, 2 L/3]^2$. Center: Release box $[-L/2, L/2]^2$. Right: Release box $[-L/12.5, L/12.5]^2$. From top to bottom: increasing time $t \in \{0, 1, 25, 50, 75\}$, in days. The color indicates the value of $p$ (with the scale on the right).

    If it is too small (Figure 2-Right) the pressure of the surrounding Wolbachia-free environment is too strong for the infection to propagate. If it is too large (Figure 2-Left), the release points are likely to be too scattered and never reach and invasion threshold. Whereas in Figure 2-Center, the release area and the number of releases is sufficient to generate a wide enough domain of Wolbachia-infected mosquitoes which spreads for larger times.


    3. Critical bubbles of non-extinction in dimension 1


    3.1. Construction

    In this section, we consider the particular one dimensional case for which we can construct a sharp critical bubble. To do so, we consider the following differential system:

    $σuα"+f(uα)=0 in R+,uα(0)=α, uα(0)=0.$ (21)

    Proposition 3.  System (21) admits a unique maximal solution $u_{\alpha}$; it is global and can be extended by symmetry on $\mathbb{R}$ as a function of class $\mathcal{C}^2$. Moreover, if $\alpha > \theta_c$, then $L_\alpha$ defined in (8) is finite and$u_\alpha$ is monotonically decreasing on $\mathbb{R}_+$ and vanishes at $L_\alpha$.

    Definition 3.1.  For $\alpha\in (\theta_c,1]$, we denote by an $\alpha$-bubble in one dimension the function $v_\alpha$ defined by

    $vα(x)=uα(|x|)+:=max{0,uα(|x|)} .$

    From Proposition 3 and Definition 3.1 we have that $v_\alpha$ is compactly supported with supp$(v_\alpha)=[-L_\alpha,L_\alpha]$.

    Proof. Local existence is granted by Cauchy-Lipschitz theorem. Then, we multiply Equation (21) by $u'_\alpha$,

    $σ2((uα)2)+(F(uα))=0,$

    which implies (since $u_{\alpha}' (0) = 0, u_{\alpha}(0) = \alpha$ and the domain is connected) that:

    $σ2(uα)2=F(α)F(uα).$

    Recall that $F(x) = \int_0^x f(y)dy$ is positive and increasing from $\theta_c$. Hence, for $\alpha > \theta_c$, $u_{\alpha}$ stays strictly below $\alpha$ except at $0$; $u_{\alpha}'$ cannot vanish unless $u_{\alpha}= \alpha$. Hence, $u_{\alpha}$ is decreasing on $\mathbb{R}_+$.

    Because $u_{\alpha}$ is decreasing, its derivative is negative and thus:

    $σduαdx=2(F(α)F(uα)).$ (22)

    Then, $u_{\alpha}$, being monotonic, is invertible on its range. Let us define $\chi_{\alpha}(u_{\alpha}(x))=x$, so that $u_{\alpha}(\chi_{\alpha}(\omega))=\omega$. By the chain rule, we have

    $dχαdω=σ2(F(α)F(ω)),$

    so that,

    $χα(ω)=αωσ2(F(α)F(v))dv.$ (23)

    Thus the function $\chi_{\alpha}$ evaluated at $\omega$ is equal to the unique radius at which the solution of (21) takes the value $\omega$. It remains to prove that $L_{\alpha} = \chi_{\alpha} (0)$ is finite, i.e. that $v \mapsto \frac{1}{\sqrt{ F (\alpha) - F(v)} }$ is integrable on $(0, \alpha)$. This function has the following equivalents at $\alpha$ and $0$:

    $1F(α)F(v)vα1f(α)1αv,1F(α)F(v)v0+{1v2f(0) if α=θc,1F(α) if α>θc.$

    Therefore $L_{\alpha}$ is finite if and only if $\alpha > \theta_c$.

    Proposition 4.  The limit bubble $u_{\theta_c}$ (also known as the "ground state") is exponentially decaying at infinity.

    Proof. The function $u_{\theta_c}$ satisfies the following equation:

    $σ2(uθc)2=F(θc)F(uθc)=F(uθc).$

    Hence,

    $σuθc=2F(uθc) on R+.$

    Moreover, for small $\epsilon$, $\sqrt{-2 F(\epsilon)} = \epsilon \sqrt{-f'(0)} + o (\epsilon)$.

    As a consequence, as $u_{\theta_c}$ gets small (at infinity), it is equivalent to the solution of

    $y=f(0)y,$

    that is $x \mapsto e^{-\sqrt{-f'(0)} x}$.

    Proof of Theorem 2.1 in dimension d=1. Let $\alpha\in(\theta_c,1]$, and let us assume that the initial data for system (1) satisfies $p(0,\cdot)\geq v_\alpha$ where $v_\alpha$ is the $\alpha$-bubble defined in Definition 3.1. From Proposition 3, it suffices to prove that $p(t,\cdot)\to 1$ locally uniformly on $\mathbb{R}$ as $t\to +\infty$.

    We first notice that the $\alpha$-bubble $v_\alpha$ is a sub-solution for (1). Indeed it is the minimum between the two sub-solutions $0$ and $u_\alpha$. Therefore, by the comparison principle, if $p(0,\cdot)\geq v_\alpha$, then for all $t>0$, $p(t,\cdot)\geq v_\alpha$.

    Then, the proof follows from the "sharp threshold phenomenon" for bistable equations, as exposed for example in [11, Theorem 1.3], which we recall below:

    Theorem 3.2 [11, Theorem 1.3]. Let $\phi_\lambda$, $\lambda>0$ be a family of $L^\infty(\mathbb{R})$ nonnegative, compactly supported initial data such that

    (ⅰ) $\lambda\mapsto \phi_\lambda$ is continuous from $\mathbb{R}^+$ to $L^1(\mathbb{R})$;

    (ⅱ) if $0<\lambda_1<\lambda_2$ then $\phi_{\lambda_1}\leq \phi_{\lambda_2}$ and$\phi_{\lambda_1}\neq \phi_{\lambda_2}$;

    (ⅲ) $\lim_{\lambda\to 0} \phi_\lambda(x)=0$ a.e. in $\mathbb{R}$.

    Let $p_\lambda$ be the solution to (1) with initial data $p_\lambda(0,\cdot)=\phi_\lambda$. Then, one of the following alternative holds:

    (a) $\lim_{t\to \infty} p_\lambda(t,x) = 0$ uniformly in $\mathbb{R}$ for every $\lambda>0$;

    (b) there exists $\lambda^* \geq 0$ and $x_0 \in \mathbb{R}$ such that

    $limtpλ(t,x)={0 uniformly in R(0λ<λ),uθc(xx0) uniformly in R(λ=λ),1 locally uniformly in R(λ>λ).$

    In our case, we define $\phi_\lambda(x)=v_\alpha(\frac{x}{\lambda})$ for $\lambda>0$. We have $\phi_1=v_\alpha$. Since $v_\alpha$ is a sub-solution to (1), the solution to this equation with initial data $\phi_1$ stays above $v_\alpha$ for all positive time. From the alternative in the above Theorem, we deduce that the solution to (1) with initial data $v_\alpha$ converges to $1$ as time goes to $+\infty$ locally uniformly on $\mathbb{R}$. (Indeed, the ground state $u_{\theta_c}$ is bounded from above by $\theta_c < \alpha$.) By the comparison principle, we conclude that if $p(0,\cdot) \geq v_\alpha$, then $\lim_{t\to+\infty} p(t,\cdot) = 1$ locally uniformly as $t\to+\infty$.


    3.2. Comparison of the energy and critical bubble methods

    Our construction of a critical $\alpha$-bubble, inspired by [4], holds in dimension $1$. In this context we may compare the "minimal invasion radius" at level $\alpha$ for initial data, given by the two sufficient conditions: being above an $\alpha$-bubble (which is the maximum of two stationary solutions), or being above an initial condition with negative energy.

    We first compute the energy of the critical $\alpha$-bubble $v_\alpha$ of Definition 3.1,

    $E[vα]=R(σ2|vα|2F(vα)) dx.$

    From Equation (21), we have

    $E[vα]=LαLα(σ|vα|2F(α))dx=2Lα0σ|vα|2 dx2LαF(α).$

    Performing the change of variable $v=v_\alpha(x)$ we have

    $Lα0|vα|2dx=α0vα(v1α(v)) dv=1σα02(F(α)F(v))dv,$

    where we use Equation (22) for the last equality. Finally, using the expression of $L_\alpha$ in (8) we arrive at

    $E[vα]=2σα0F(α)2F(v)2(F(α)F(v)) dv.$

    To emphasize the difference between the two sufficient conditions, we observe that when $\alpha\to \theta_c$, since $F(\theta_c)=0$, we obtain

    $E[vθc]=2σθc02F(v) dv>0.$

    By continuity of $\alpha \mapsto E[v_{\alpha}]$ we deduce

    Lemma 3.3. The $\alpha$-bubbles $v_{\alpha}$ have positive energy if $\alpha$ is close to $\theta_c$.

    Remark 5. In particular, the energy estimate alone does not imply invasiveness of the $\alpha$-bubbles, which justifies the interest of our particular approach in one dimension. We do not claim that the "energy" or the "bubble" method is better, but we highlight the fact that they do not perfectly overlap.

    Figure 3 gives a numerical illustration of the fact that $\alpha$-bubbles give smaller radii at level $\alpha$, except for $\alpha \sim 1$, and at any rate provide a smaller minimal radius for invasion when the same parameters as in Figure 1 are used.

    Figure 3. Comparison of minimal invasion radii $R_{\alpha}$ (obtained by energy) in dashed line and $L_{\alpha}$ (obtained by critical bubbles) in solid line, varying with the maximal infection frequency level $\alpha$. The scale is such that $\sigma=1$.

    4. Specific study of a relevant set of release profiles

    In this section we discuss a specific release protocol, with a total of $N$ mosquitoes divided equally into $k$ locations, in a space of dimension $1$. It yields a release profile in the set $RP_k^d (N)$ we defined in (10).


    4.1. Analytical study of the case of a single release

    In the case of a single release ($k=1$), we can easily describe the relationship between the mosquito diffusivity $\sigma$ and the total number of mosquitoes to release. Morally, as long as the mosquitoes diffuse they could theoretically invade (in dimension $1$) by a single release, by introducing a sufficiently large amount of mosquitoes. This is the object of the next proposition:

    Proposition 5.  Let $G_\sigma(\tau):=G_{\sigma,1}(\tau) = \frac{1}{\sqrt{2\pi \sigma}} e^{-\tau^2/2\sigma}$. The following equivalent properties hold:

    (ⅰ)  There exists $\sigma_+ : \mathbb{R}_+^* \to \mathbb{R}_+^*$ such that $N G_{\sigma}$ satisfies (NEC) with $\tau_0 = 0$ if and only if $\sigma \in (0, \sigma_+(N)]$.

    (ⅱ)  There exists $N_m : \mathbb{R}_+^* \to \mathbb{R}_+^*$ such that $N G_{\sigma}$ satisfies (NEC) with $\tau_0 = 0$ if and only if $N \geq N_m (\sigma)$.

    Moreover, $\sigma_+$ and $N_m = \sigma_+^{-1}$ are increasing and in both cases, evolution in (1) with initial data $p_{\sigma, N} := \frac{N G_{\sigma}}{N G_{\sigma} + N_0}$ yields invasion by the introduced population.

    Part () of Proposition 5 asserts that if we fix the total number $N$ of mosquitoes to introduce, single introduction is a failure if diffusivity is too large. Part () is just the converse viewpoint: if we know estimates on the diffusivity (thanks to field experiments like mark-release-recapture for example [35]), then we can define a minimal number $N_m$ of mosquitoes to introduce at a single location to succeed.

    Remark 6. If $\alpha \in (\theta_c, 1)$ makes $N G_{\sigma}$ satisfy (NEC) ("be above the $\alpha$-bubble"), then necessarily (evaluating at $0$ to take the maximum of $G_{\sigma}$), $\alpha \leq \frac{N}{N + \sqrt{2 \pi \sigma} N_0}$. In particular, our under-estimation of the probability is equal to $0$ as soon as

    $N<2πσN0θc1θc.$

    Equivalently, the density of mosquitoes at the center of the single release location $\frac{N}{\sqrt{2 \pi \sigma}}$ should exceed $\frac{\theta_c}{1 - \theta_c} N_0$ for our estimate to prove useful. (If $\theta_c = 0.8$, this is already $4$ times the existing mosquito density. If $\theta_c = \frac{2}{3}$, then it is only $2$ times; in the case of Figure 1, $\theta_c=0.36$ and then the ratio is only $0.56$).

    Proof of Proposition 5. Both the introduction profile given by the fraction $\displaystyle\frac{N G_{\sigma} (\tau)}{N G_{\sigma} (\tau) + N_0}$ and non-extinction bubbles from Theorem 2.1 built by (21) $(u_{\alpha} (\tau))_{\alpha}$ are symmetric, radial-decreasing functions. Instead of comparing them, we compare their reciprocals. We define $T_{\sigma, N}$ such that for all $p \in [0, \alpha]$,

    $NGσ(Tσ,N(p))NGσ(Tσ,N(p))+N0=p,$

    and $\chi_{\alpha}$ such that $u_{\alpha} (\chi_{\alpha} (p)) = p$. Respectively, they read

    ${Tσ,N(p)=2σlog(NN02πσ1pp),χα(p)=σ2αpdvF(α)F(v).$ (24)

    By construction, the following equivalence holds

    $τR+, NGσ(τ)NGσ(τ)+N0uα(τ)p s.t. 0pα, χα(p)Tσ,N(p).$

    Using (24) this rewrites as

    $log(NN02πσ)(αpdv2F(α)F(v))2log(1pp),p[0,α].$ (25)

    From (25), we define

    $Jα(p):=log(p)log(1p)+(αpdv2F(α)F(v))2,$ (26)
    $I(σ,N):=log(N2πσN0).$ (27)

    For any given $N$, the problem we want to solve amounts at finding couples $(\alpha, \sigma)$ such that

    $p[0,α], Jα(p)I(σ,N).$ (28)

    We study the function $J_{\alpha}$. First, we note that $J_{\alpha} (p) \xrightarrow[p \to 0]{} - \infty$, $J_{\alpha} (\alpha) = \log \big( \frac{\alpha}{1 -\alpha} \big)$ and it is continuous. Moreover,

    $Jα(p)=1p(1p)1F(α)F(p)αpdv2F(α)F(v),$

    and we may compute $\lim_{p \to \alpha} J_{\alpha}' (p) = \frac{1}{\alpha (1-\alpha)} - \frac{1}{f(\alpha)}$. Then, we can define

    $jα:=maxp[0,α]Jα(p),j:=minα(θc,1]jα.$

    Thus there exists $\alpha \in (\theta_c, 1]$ such that (25) holds if and only if $N \geq N_0 \sqrt{2 \pi \sigma} e^{j^*}$. This gives Proposition 5 () with $\sigma_+ (N) = \frac{e^{-2 j^*}}{2 \pi} \big( \frac{N}{N_0} \big)^2$ and Proposition 5 () with $N_m = N_0 \sqrt{2 \pi \sigma_+} e^{j^*}$.

    Remark 7. With parameter values from (5), the expected number of mosquitoes to release is huge, since we need to have $\frac{N_m}{N_0 \sqrt{2 \pi \sigma}} = e^{j^*} \simeq 7 \cdot 10^{10}$ with $j^* \simeq 25$, where $\frac{N_m}{N_0 \sqrt{2 \pi \sigma}}$ is the quotient between total mosquitoes to release and wild initial population in an area of typical size $\sqrt{2 \pi \sigma}$. (This is approximately the distance diffused in $1$ day, equal to $72 \textrm{m}$ in this case). To obtain $j^*$ numerically, we used MATLAB function fminbnd. Here, the model has a clear and crucial conclusion: it is very hard to invade an area with a single, localized release. Therefore, we must model several releases (whether in time or in space). In the rest of the paper we discuss the case of releases at multiple locations at same time $t=0$.


    4.2. Equally spaced releases

    If we space the $k$ release points regularly in the interval $[-L_{\alpha}, L_{\alpha}]$, we want to check that (NEC) holds for

    $Xτ=Nkk1i=0Gσ(τ+Lα(1+2ik1)).$

    Within a fairly good approximation, this is the case if

    $τ[Lα,Lα],XτXτ+N0α.$

    This holds in particular if

    $N˜N(k,α,σ)=N02πσ2α1αkeL2α2σ(k1)2.$

    If we fix $\sigma$ then we may try to find optimal $k$ and $\alpha$ in order to minimize $\widetilde{N}$. Alternatively, we can do the same, fixing $N$ or $N/k$ (number of mosquitoes per release), and find the optimal number of releases $k$.

    It is straightforward, keeping in mind that $L_{\alpha}$ is proportional to $\sqrt{\sigma}$, that the optimal $\alpha$ here merely depends on $k$, not on $\sigma$. We may introduce

    $j(k):=minα(θc,1)α1αeL2α/(2σ(k1)2).$

    and find the minimal (in view of our sufficient criterion) value $\widetilde{N}^*$ for $\widetilde{N}$:

    Lemma 4.1.  For $k$ equally spaced releases on the line, there exists an invading release profile with $L^1$ norm:

    $˜N(k,σ)=N02πσk2j(k).$ (29)

    However, we want to take into account the uncertainties and variability in the release protocol and population fixation. Namely, the release points might not be exactly equally spaced, so that introducing $\widetilde{N}^*$ mosquitoes would only give some probability of success. This is what we want to quantify now and shall be addressed in Section 4.3.


    4.3. Multiple release locations: Towards a geometric problem

    When we sum several Gaussians, the profile is neither symmetric (in general), nor monotone. Therefore the previous analytical argument does not apply. However, at the cost of fixing $\sigma$ we are left with a simple geometric problem.

    First step: fixing $\sigma$ and bounding by level rather than profile. We assume first that there is no uncertainty on $\sigma$, which is taken as $\sigma_0$ ($\epsilon = 0$ in (10)). As a further simplification, we shall not compare the introduction frequency profile to some $\alpha$-bubble (because it is too hard), but rather to the very simple upper bound of an $\alpha$-bubble: the characteristic function $\tau \mapsto \alpha 1_{- L_{\alpha} \leq \tau \leq L_{\alpha}}$.

    Moreover, we assume that our $k$ release locations $(x_i)_{1 \leq i \leq k}$ are within the compact set $[-L, L]$, for some $L>0$. As above, we write

    $Gσ(y):=12πσey2/2σ$

    and

    $G=Nkki=1Gσ(xi).$

    We define

    $P(σ,Nk,(xi)1ik,L0,α):=min[Lα+L0,Lα+L0]G$ (30)

    Then, the probability of success for the release of $N$ mosquitoes in a total of $k$ different sites in $[-L, L]^k$, when they all spread according to $\sigma$ diffusivity, and the initial population density was $N_0$, is given by:

    $Pk(N,L)=P[L0R, α(θc,1), P(σ,Nk,(xi)1ik,L0,α)α1αN0].$ (31)

    Here, the probability $\mathbb{P}$ is taken over all the real $k$-uples $(x_l)_{1 \leq l \leq k}$ such that $-L < x_1 \leq \dots \leq x_k < L$, and $[-L, L]^k$ is equipped with the uniform measure.

    Second step: transformation into a geometric problem. In order to get a more tractable bound, we make use of the following property:

    Proposition 6.  Let $(x_i)_{i} \in [-L, L]^k$ with $x_1 \leq \dots \leq x_k$. We define$\mathcal{G} = \frac{N}{k}\sum_{i=1}^k G_{\sigma} (\cdot - x_i)$.

    If there is $\alpha \in (\theta_c, 1)$ such that

    $Nk12πσα1αN0$

    and $1 \leq l < m \leq k$ such that

    (ⅰ)  $\forall l \leq j \leq m-1, ~ x_{j+1} - x_j \leq 2 \sqrt{2 \log(2)} \sqrt{\sigma}$,

    (ⅱ)  $ x_m - x_l \geq 2 L_{\alpha}$,

    then

    $GG+N0vα(xm+xl2) .$

    We notice that the constant $2 \sqrt{2 \log (2)} \simeq 2.35$ is optimal with this property: if two translated Gaussians centered at $x_0, x_1$ are at a distance $x_1 - x_0 = \lambda \sqrt{\sigma}$, with $\lambda > 2 \sqrt{2 \log (2)}$, then their sum is smaller at $\frac{x_0 + x_1}{2}$ than at $x_0$.

    Proof. This property relies on the simple computation of the sum of two $G_{\sigma}$s, centered at $-h$ and $h$ ($h>0$), is greater than $G_{\sigma} (0)$ on $[-h, h]$ as soon as $h \leq \sqrt{2 \log(2)}\sqrt{\sigma}$. Figure 4 illustrates this property.

    Figure 4. Two $G_{\sigma}$ profiles and their sum (in thick line). The level $G_{\sigma} (0)$ is the dashed line. On the left, $h=\sqrt{2\log(2)\sigma}$. On the right, $h>\sqrt{2\log(2)\sigma}$.

    Indeed, considering the sum of two Gaussian $G_\sigma$,

    $ξ(x)=12πσ(e(x+h)22σ+e(xh)22σ)=2eh22σGσ(x)cosh(xhσ).$

    Then, recalling that $\sigma G_{\sigma}'(z) = - z G_{\sigma} (z) $, we compute

    $12eh22σσξ(x)=xGσ(x)cosh(xhσ)+hGσ(x)sinh(xhσ)12eh22σσ2ξ"(x)=(h2+x21σ)Gσ(x)cosh(xhσ)2hxGσ(x)sinh(xhσ).$

    As a consequence, the sign of $\xi" (x)$ is that of

    $γ(x):=h2+x22hxtanh(xhσ)1σ.$

    We notice that $\gamma(0) = h^2 - \frac{1}{\sigma}$. Hence, $\xi$ has a local maximum (resp. a local minimum) at $x=0$ if $h < \sqrt{\sigma}$ (resp. $h > \sqrt{\sigma}$). Since $\xi(0) = 2 e^{-\frac{h^2}{2 \sigma}} G_{\sigma} (0)$, the maximal $h > 0$ that ensures $\xi(0) \geq G_{\sigma} (0)$ is $h=h_0 :=\sqrt{2 \log(2) \sigma}$.

    Now, we examine the necessary condition $\xi'(x) = 0$ for a local extremum on $(-h, h)$. It implies

    $x=htanh(xhσ).$

    This is true for $x=0$ (and we have seen the condition on $h - \sqrt{\sigma}$ to have a local extremum indeed). Then, there is a solution $x_+ > 0$ if, and only if, $\frac{h^2}{\sigma} > 1$, i.e. $h > \sqrt{\sigma}$. In this case, $x_+$ is unique (and $x_- := - x_+$ is a solution as well).

    So, for $h=h_0>\sqrt{\sigma}$, we know that $\xi$ has a local minimum at $x=0$, is smooth, has at most one local extremum on $(0, + \infty)$, and goes to $0$ at $+\infty$. Hence, this local extremum exists and is a maximum. Therefore (and by symmetry), the minimum of $\xi$ on $(-h, h)$ is attained at $x=0$ or $x=h$. Since $h=h_0$, $\xi(h) > \xi(0) = G_{\sigma} (0)$. We deduce that $\xi > G_{\sigma} (0)$ on $(-h, h)$.

    We may use this property to prove Proposition 6. By condition (ⅰ) the above lower-bound holds between $x_l$ and $x_m$, and not only between two adjacent locations $x_j, x_{j+1}$. Now, the first condition implies that $G_{\sigma} (0) \geq \frac{\alpha}{1-\alpha} N_0$. Combining these two facts with $x_m - x_l \geq 2 L_{\alpha}$ implies that

    $GG+N0α,$

    on $[x_l,x_m]$ which is an interval of length at least $2 L_{\alpha}$. Precisely, for all $x \in \mathbb{R}$,

    $G(xxm+xl2)G(xxm+xl2)+N0αvα(xxm+xl2).$

    As a consequence, we may translate the generic inequality (SP) into:

    $P1k(N,(L,L))=Pk(N,L)P[α(θc,11+N0Nk2πσ),1l<mk,xmxl2Lα and ljm1,xj+1xj22log(2)σ]$ (32)

    Then, we define

    $L:=minθc<α11+N0Nk2πσLα,$

    and equivalently estimate (32) reads

    $Pk(N,L)P[1l<mk,xmxl2L and maxljm1(xj+1xj)22log(2)σ].$ (33)

    The study of the minimization of $L_\alpha$ with respect to $\alpha$ is discussed further in Appendix.

    Remark 8. Note that for this estimate, we only consider initial data that are above a characteristic function at level $\alpha$ on an interval of length $2 L_{\alpha}$. This is far from being the optimal way to be above the $\alpha$-bubble $v_{\alpha}$.

    Remark 9. It is easy to check that our estimate yields $0$ (no information) as long as $k$ is too small, namely $k \sqrt{2 \log(2)} \sqrt{\sigma} \leq L^*$. A necessary condition for our estimate not to yield $0$ may read:

    $k12log(2)minθc<α1α0dv2(F(α)F(v)).$

    Specific discussion for $\alpha = \theta_c$. By Proposition 4, $u_{\theta_c}$ decays exponentially. As a consequence, no sum of $G_{\sigma}$s may be above it. This is why this profile cannot be used in our approach (because we consider that introduction profiles should be Gaussian).


    4.4. Analytical computations of the probability of success: Recursive formulae

    In order to compute analytically the right-hand-side in (33), we may introduce the following notations:

    ●  $\mathcal{T}_k (u, v)$ is the set of ordered $k$-uples between $u$ and $v$ ($u<v \in \mathbb{R}$), the measure of which is

    $τk(u,v)=(vu)k+k!.$

    ●  $\mathcal{C}_k^{\lambda} (u,v) \subseteq \mathcal{T}_k (u,v)$ is the subset of $k$-uples such that $y_1= u$, $y_k = v$ and for all $l \in [\![ 1, k-1 ]\!]$, $y_{l+1} - y_l \leq \lambda$. Its measure is denoted $\gamma_k^\lambda(u,v)$.

    ●  $\mathcal{B}^{\lambda, R^*}_k (u,v) \subseteq \mathcal{T}_k (u,v)$ is the subset of $k$-uples such that $\exists 1 \leq l < m \leq k, y_m - y_l \geq R^*$ and $\max_{l \leq j \leq m-1} (y_{j+1} - y_j) \leq \lambda$. We denote $\beta_k^{L,R^*}(u,v)$ its measure.

    Remark 10. Back to problem (33), we recover the problem of estimating $\beta$ with the notations of Proposition 7 through a simple change of variables. We divide all positions ($x_1, \dots, x_k$) by $\sqrt{2 \sigma}$. Then in the right-hand side of (33) we replace $2 L^*$ by

    $R:=minαα0dvF(α)F(v) ,$

    and $2 \sqrt{2 \log(2) \sigma}$ by $\lambda := 2 \sqrt{\log(2)}$. This was done in order to simplify computations. Moreover, it shows that the success probabilities do not depend on diffusivity. In fact, scaling in $\sigma$ as we did merely amounts at choosing a space scale such that $\sigma = 1$. Even though probabilities themselves do not make $\sigma$ appear, one must keep in mind that the corresponding release protocols (including the space between release points or the size of the release box) are proportional to $\sqrt{\sigma}$.

    We want to under-estimate the probability of success with $k$ releases in the box $[-L, L]$. In view of (SP), it amounts to computing $\frac{\beta_k^{\lambda, R^*} (-L, L)}{\tau_k (-L, L)}$. In fact, we get a general recursive formula for $\beta$ in the following proposition.

    Proposition 7.  Let $k_0 := \lceil\!\!\lceil \frac{R^*}{\lambda} \rceil\!\!\rceil + 1$. Then,

    $βλ,Rk(L,χ)=ki=k0ki+1j=1χRLmin(χ,u+(k1)λ)u+Rγλi(u,v)$
    $(τj1(L,uλ)βλ,Rj1(L,uλ))τk(i+j1)(v+λ,χ)dvdu.$ (34)

    Proof. The idea is simple: we count each "positive initial data", that is an ordered $k$-uple $(y_i)_i$ such that a subfamily satisfies $y_m - y_l \geq R^*$ and $y_{i+1} - y_i \leq \lambda$ between $l$ and $m$, according to its leftmost "positive sub-family", which is then taken of maximal length.

    We shall use the index $i$ to denote the length of this maximal family (between $k_0$ and $k$), and $j$ its first rank ($1 \leq j \leq k-i+1$). Then,

    $βλ,Rk(L,χ)=[L,χ]k1{y1y2yk}1{(y1,,yk)Bλ,Rk(L,χ)}dy1dyk.$ (35)

    Now, we split:

    $1{(y1,,yk)Bλ,Rk(L,χ)}=ki=k0ki+1j=11{yi+j1yjR}j+i2l=j1{yl+1ylλ}1{(y1,,yj1)Bλ,Rj1(L,χ)}1{yjyj1>λ}1{yi+jyi+j1>λ}.$ (36)

    This identity requires some explanations. It comes from the partition of $\mathcal{B}$ using maximal leftmost positive sub-family, as described above. Then, the term $1_{\{ (y_1, \dots, y_{j-1}) \not\in \mathcal{B}^{\lambda, R^*}_{j-1} (-L, \chi) \}}$ simply comes from the definition of $\mathcal{B}$. Since we consider the leftmost positive subfamily, no family on its left should be positive. Moreover no element on its left can be added, which justifies the $1_{\{y_j - y_{j-1} > \lambda \}}$. Then, we have in addition that for $j > 1$ and $y_j \leq \chi$,

    $1{(y1,,yj1)Bλ,Rj1(L,χ)}1{yj1yj}1{yjyj1>λ}=1{(y1,,yj1)Bλ,Rj1(L,yjλ)},$

    with the obvious convention that $\mathcal{B} (u, v) = \emptyset$ if $v < u$.

    In addition, for $i + j - 1 < k$

    $[L,χ]k(i+j1)1{yi+j1yk}1{yi+jyi+j1>λ}dyi+jdyk=τk(i+j1)(yi+j1+λ,χ)=(χyi+j1λ)k(i+j1)+(k(i+j1))!.$

    Combining these results, and using (35) and (36) yields

    $βλ,Rk(L,χ)=ki=k0ki+1j=1χLχxi+j21{yj+i1yjR}j+i2l=j1{0yl+1ylλ}(τj1(L,yjλ)βλ,Rj1(L,yjλ))τk(i+j1)(yi+j1+λ,χ)dyjdyi+j1,$ (37)

    with conventions $\tau_0 = 1$ and $\beta_0 = 0$, regardless of their arguments.

    We assume $\chi \geq - L + R^*$ (otherwise $\beta_k^{\lambda,R^*} (-L, \chi) = 0$). Using the notation $\gamma$ we introduced, Equation (37) is simplified again into:

    $βλ,Rk(L,χ)=ki=k0ki+1j=1χRLmin(χ,u+(k1)λ)u+Rγλi(u,v)(τj1(L,uλ)βλ,Rj1(L,uλ))τk(i+j1)(v+λ,χ)dvdu,$

    where $u$ stands for $y_j$ and $v$ for $y_{i+j-1}$. This is our recursive formula (34).

    Now, we may give an explicit formula for $\gamma^{\lambda}_i (u, v)$. We should notice that by definition,

    $γλi+2(u,v)=u+λuu1+λu1ui1+λui11vuivλduidu1,$

    that is

    $γλi+2(u,v)=u+λuγλi+1(u1,v)du1.$ (38)

    Hence, we deduce the recursive formula,

    Lemma 4.2. For all $i, \lambda, u, v$ as above,

    $γλi+2(u,v)=λi+i+1k=1(1)ki!((ik1)(vukλ)i++(1)i+1(i1k1)(kλ(vu))i+).$ (39)

    Proof. Obviously, $\gamma_2^\lambda (u, v)= 1_{v \geq u \geq v-\lambda}$ and we deduce from (38)

    $γλ3(u,v)=λ+(vu2λ)+(λ(vu))+(vuλ)+$

    Then, using (38) again proves (39) by induction.

    Remark 11.  For $k < 2 k_0$, formula (34) simplifies a lot for it is no longer recursive. It enables us to compute $\beta_{k_0}^{\lambda, R^*} (-L, L)$.

    $βλ,Rk0(L,L)=LRLmin(L,u+(k01)λ)u+Rγλk0(u,v)dvdu.$ (40)

    Then by (39) we know $\gamma_{k_0}^{\lambda} (u, v)$. With the change of variables $w = v + u$, when $L > - L + (k_0 - 1) \lambda$, equation (40) becomes

    $βλ,Rk0(L,L)=L(k01)λL(k01)λR(λk02+k01k=1(1)k(k02)!((k02k1)(wkλ)k02++(1)k01(k03k1)(kλw)k02+))dwdu+LRL(k01)λLuRγλk0(u,u+w)dwdu.$ (41)

    Clearly, the first integral in the right-hand side of (41) may be written as

    $(2L(k01)λ)f1(λ,R),$

    where $f_1$ does not depend on $L$. With the change of variables $z = L - u$, the second term in the right-hand side of (41) becomes

    $f2(λ,R):=(k01)λRzR(λk02+k01k=1(1)k(k02)!((k01k1)(wkλ)k02++(1)k01(k03k1)(kλw)k02+))dwdz.$

    In particular, it appears that it does not depend on $L$. (Recall that by definition, $k_0~=~\lceil\!\!\lceil \frac{R^*}{\lambda} \rceil\!\!\rceil~+~1$).

    For $\chi \in (-L + R^*, -L + (k_0 - 1) \lambda)$, we can compute similarly

    $βλ,Rk0(L,χ)=χ(L)RzRγλk0(0,w)dwdz,$

    and notice that our expressions are consistent since

    $βλ,Rk0(L,L+(k01)λ)=L+(k01)λ(L)RzRγλk0(0,w)dwdz=f2(λ,R).$

    All in all, $\beta_{k_0}$ is expressed as follows:

    $βλ,Rk0(L,χ)={0 if χ+LRχ(L)RzRγλk0(0,w)dwdz if χ+L(R,(k01)λ),(χ+L(k01)λ))f1+f2 if χ+L>(k01)λ$ (42)

    (This is an affine function for $\chi + L > (k_0 - 1 ) \lambda$, with pent $f_1 (\lambda, R^*)$).

    Then, we obtain a bound on the probability of success with $k_0$ (the minimal number of) releases after dividing by $\tau_{k_0} (-L, L)$ :

    $Pk0(L)βλ,Rk0τk0(L,L)=k0!(2L)k0((2L(k01)λ)f1(λ,R)+f2(λ,R)).$

    In particular, we see that this underestimation of the success probability is increasing and then decreasing, and thus reaches a unique maximum at $L = \widehat{L}$.

    We find

    $2ˆL=λ(Rλ+1)k0k01f2(λ,R)f1(λ,R).$

    We may note that introducing the non-negative and non-decreasing function

    $Γλ,Rk(z):=zRγλk(0,w)dw$

    we get

    $f1(λ,R)=Γλ,Rk0((k01)λ),f2(λ,R)=(k01)λRΓλ,Rk0(z)dz.$

    As a consequence, $f_2 \leq \big( (k_0-1)\lambda - R^* \big) f_1$ and thus

    $2ˆLk0k01R.$

    5. Numerical results

    Now, we present some numerical results we obtained on this set of release profiles. Numerical simulations confirm the intuition of Proposition 2. Our under-estimation is not very bad. Indeed, as one increases the number of release points ($k$) in a fixed perimeter, with a fixed number of mosquitoes per release, then our under-estimation of the probability of success converges to $1$.

    Figure 5 shows the probability profile as a function of the size $L$ of the release box, for $20$, $40$ and $80$ release points. With parameter values from (5), $R^* = 10.981$, $\lambda = 1.665$ and thus $k_0 = 8$. The curves are obtained by a simple Monte-Carlo method. They lead to the appearance of an optimal size for the release box ($6.3$ in this example), that does not seem to depend on the number of release points between $20$ and $80$.

    Figure 5. Under-estimation $\beta^{\lambda, R^*} (-L, L)$ of introduction success probability for $L$ ranging from $R^*/2 = 5.49$ to $3 R^*/2 = 16.47$. The seven curves correspond to increasing number of release points. (From bottom to top: $20$ to $80$ release points).

    However, for small (relatively to $k_0$) numbers of releases, the probabilities are very small. In the case of $10$ release points, the maximal probability we find is about $1. 10^{-5}$.

    Our numerical values are somehow consistent with field experiments: typically, the space between release points is less than $\lambda \sqrt{2 \sigma}$, which is about $68 \mathrm{m}$, and the optimal box size is approximately equal to $6.3 \times \sqrt{2 \sigma} \simeq 257 \mathrm{m}$, with the values from (5).

    The factor $2 \sqrt{2 \log(2)}$ is crucial with this respect. Losing it changes $\lambda$ from $2\sqrt{ \log(2)} \simeq 1.665$ to $1/\sqrt{2} \simeq 0.707$ and makes $k_0$ ("the minimal theoretical number of releases to make our under-estimation of the probability of success positive") increase from $8$ to $17$. We show in Figure 6 the probability profile for $80$ releases in this case, to illustrate the loss with this "worse" geometric estimation. It culminates at around $50 \%$ only and is comparable with the green curve (for $40$ release points) of Figure 5.

    Figure 6. Effect of losing the constant $2 \sqrt{2 \log(2)}$ in Proposition 6: under-estimation $\beta^{\lambda, R^*} (-L, L)$ of introduction success probability for $L$ ranging from $R^*/2 = 5.49$ to $3 R^*/2 = 16.47$, with $80$ release points.

    6. Conclusion and perspectives

    We considered spatial aspects of a biological invasion mechanism associated to release programs and their uncertainty. We validated the framework in the one-dimensional case, and the two-dimensional case is the natural extension.

    Two difficulties must be tackled in higher dimensions. First, the radial-symmetric "$\alpha$-bubbles" may still exist, but we no longer have an exact formula like (8) for their support. Second, the geometric problem underlying our estimation gets harder, but not impossible to manage. To deal with it, we need an analogue of Proposition 6 in order to get a lower bound for a sum of Gaussians in two dimensions.

    An interesting feature of the approach we introduced is that it can be extended to cases when neither sub-solutions nor geometric properties are available. Heuristically, we need first a criterion to tell us if a given initial data belongs to a "set of interest". Second, we need to put a probability measure on the set of "feasible initial data". Combining these, we compute the probability that the criterion is satisfied. This probability gives an insight into the role any given aspect of the release protocol plays.

    We used a sufficient condition for invasion, the criterion from Theorem 2.1. However, we proved that our under-estimation of probability is rather good: in particular, it converges to $1$ when the number $k$ of releases goes to $\infty$. This fact is the object of Proposition 2, holds true in any dimension, and is supported by numerical simulations in dimension $1$.

    We have always considered a homogeneous "context of introduction", so that the stochasticity would only affect the release process itself. Another natural continuation of this work, trying to go further into spatial stochasticity for release protocols, is the use of other stochastic parameters, such as the diffusion process (here it is given by a deterministic diffusivity $\sigma$), or the local carrying capacity. We let this open for further research.

    Some other questions remain open. For instance: in one dimension, we considered releases in $[-L, L]$. We know that if $2 L < L^*$ then our condition in the right-hand side of (33) is zero. On the other hand, this right-hand side goes to $0$ as $L \to + \infty$. This suggests that there exists a (non-necessarily unique) size $\widehat{L}$ that maximizes this right-hand side. Back to (40), we obtained in Remark 11 a lower bound for $\widehat{L}$ in this case:

    $ˆLR1+RλRλ.$ (43)

    It is a numerical conjecture that the optimal value of $L$ is close to $\frac{1}{2} ( \lambda + R^*)$ for any $k$. For this particular protocol feature (the optimal size of the release area), our approach already provides an interesting indication which - to the best of our knowledge - has not been used in previous release experiments.

    As a possible follow-up to this work, one can set up several optimization problems. First, on a purely theoretical side, how to optimize the threshold functions in Theorem 2.1 with respect to a cost functional such as the $L^1$ norm (for the total number of released mosquitoes)? Then, if we fix a cost, how to maximize the under-estimated probability of success with respect to the size of the release area? Ultimately, how to optimize a release protocol (playing on the probability law of the release profiles space)?


    Appendix: Uniqueness of the minimal radius

    In this appendix we investigate sufficient conditions for the uniqueness of a minimal radius among the $\alpha$ bubbles we constructed in Section 3. More precisely, we establish the number of bubbles of a given radius (which is typically $2$). General results in any dimension on the exact multiplicity of solutions for such problems (semilinear elliptic Dirichlet problems) have been obtained in [28] and [29], so in essence the results below are not new and are even contained in the cited articles. However we emphasize that our proof, limited to dimension $1$, uses very simple arguments and even provides an equivalent formulation of the problem in terms of a single real function $h$ built from $f$ and $F$, see formula (45) below.

    Let $f \in \mathcal{C}^2 ([0, 1], \mathbb{R})$ be a bistable function in the sense of (3) and $F (x) = \int_0^x f(y) dy$ its antiderivative as introduced in (4).

    We make the following assumptions:

    $f(0)<0,f(θ)>0,f(1)<0,F(1)>0,x[0,1],(f(x)+xf"(x))f(x)x(f(x))2.$

    Under assumption (B1), there exists a unique $\theta_c \in (\theta, 1)$ such that $F(\theta_c) = 0$. We introduce

    $g(x):=xf(x)/f(x).$ (44)

    By simple computation we have

    Lemma .1.  Under assumption (B0), (B2),$g$ is decreasing on $[0, \theta)$ and on $(\theta, 1]$. In addition, $g(0) = 1$, $g(\theta_-) = - \infty$, $g(\theta_+) = + \infty$ and $g(1) = - \infty$. As a consequence, there exists a unique $\alpha_1 \in (\theta, 1)$ such that

    $g(α1)=1.$

    We add the following assumption:

    $α>max(θc,α1),F(α)(f(α)+αf(α))α(f(α))2.$

    Now, we recall the $\alpha$-bubble radius, as introduced before, for $\alpha \in (\theta_c, 1]$:

    $Lα=σα0dv2(F(α)F(v)).$

    Proposition 8.  Under conditions (B0), (B1), the bistable (in the sense of (3)) function $f$ is such that$L_{\alpha}$ reaches its minimum on $(\theta_c, 1]$ (which is well-defined) at points in $(\theta_c, 1)$.

    If in addition (B2), (B3) hold, then there exists a unique $\alpha_0 \in (\theta_c, 1)$ such that

    $Lα0=minαLα,$

    and for all $L > L_{\alpha_0}$, there exists unique $\alpha_{\pm} (L)$ with $\alpha_- (L) \in (\theta_c, \alpha_0)$ and $\alpha_+ (L) \in (\alpha_0, 1)$ such that $L_{\alpha_{\pm} (L)} = L$.

    Remark 12. Although assumptions (B0) and (B1) are very general, (B2) and (B3) are debatable. They yield a simple sufficient condition for uniqueness of minimum (which is the object of Proposition 8), but are by no means necessary to get it. We expect that they can be refined and improved in order to get uniqueness for a wider class of bistable functions.

    Using $f$ defined by (2) with values from (5), we verified numerically that (B2)-(B3) are satisfied. Indeed, using MATLAB we found that $x(f'(x))^2 - f(x) (f'(x) + xf"(x))$ and $x (f(x))^2 - F(x) (f(x)+xf'(x))$ are increasing on $[0, 1]$ and $[\alpha_1, 1]$ (with $\max(\theta_c, \alpha_1) = \alpha_1$), respectively. The former is equal to $0$ at $0$, and the latter is approximately equal to $2 \cdot 10^{-4} > 0$ at $\alpha_1$ in this case, hence the two assumptions hold.

    Generally, we can check that (B2)-(B3) hold for the classical bistable function $f(x) = x (1 - x) (x - \theta)$ with $\theta \in (0, 1/2)$. We first compute

    $f(x)+xf"(x)=9x2+4(1+θ)xθ.$

    Then (B2) is equivalent to

    $(9x24(1+θ)x+θ)x(x2(1+θ)x+θ)x(3x22(1+θ)x+θ)213(1+θ)x2+10θx5θ(1+θ)12(1+θ)x2+6θx4θ(1+θ)0(1+θ)x24θx+θ(1+θ).$

    The discriminant of this second-order polynomial is $-4 \theta (1 - \theta)^2 < 0$, so this inequality holds for any $\theta \in (0, 1)$. Then a straightforward computation shows that $\alpha_1 = \frac{1 + \theta}{2}$.

    Now, we want to check (B3). To do so we compute

    $F(x)=14x4+1+θ3x3θ2x2.$

    Then (B3) is equivalent to

    $x2(14x21+θ3x+θ2)(4x23(1+θ)x+2θ)x3(x2(1+θ)x+θ)2x2(1+θ)(23443)+θ2x+θ(1+θ)(23223)0.$

    Then we recall that $2 - \frac{3}{4} - \frac{4}{3} = -\frac{1}{12} < 0$, so we just need to show that the discriminant is negative. This discriminant is equal to

    $θ24θ(1+θ)29=θ4(θ49(1+θ)2)<0.$

    Hence simplest bistable functions of the form $f(x) = x (1 - x) (x - \theta)$ satisfy our assumptions (B2) and (B3), and in particular the set of such functions is non-empty.

    Proof. Without loss of generality we assume $\sqrt{\sigma} = \sqrt{2}$ to get rid of the constant. From (8), we deduce the equivalent expression:

    $Lα=α0(1F(α)F(v)1f(α)(αv))dv+α0dvf(α)(αv)=1f(α)(α0(f(α)F(α)F(v)1αv)dv+2α)$

    Hence

    $ddαLα=1αf(α)+12f(α)α0(1(αv)3/2(f(α)F(α)F(v))3/2)dv,$

    which is a continuous function from $(\theta_c, 1)$ to $\mathbb{R}$. It is easily seen that $\frac{d}{d\alpha} L_{\alpha}$ goes to $- \infty$ as $\alpha \to \theta_c^+$, and to $+ \infty$ as $\alpha \to 1^-$ (recalling $f(1) = 0$). Therefore, we know that $L_{\alpha}$ reaches its minimum (which is well-defined) at points strictly in the interior of $(\theta_c, 1)$. This is the first point of Proposition 8.

    Then, $\frac{d}{d \alpha} L_{\alpha} = 0$ if and only if

    $1α+12α0(1(αv)3/2(f(α)F(α)F(v))3/2)dv=0.$

    For $\alpha \in (\theta_c, 1)$, we introduce

    $h(α):=10(1(1w)3/2(αf(α)F(α)F(αw))3/2)dw.$ (45)

    Then $\frac{d}{d \alpha} L_{\alpha} = 0$ if and only if $h(\alpha) = -2$. In addition, $h(\theta_c) = - \infty$ and $h(1)= +\infty$ are well-defined by continuity.

    We compute

    $h(α)=3210(αf(α))1/2(F(α)F(αw))5/2((f(α)+αf(α))(F(α)F(αw))αf(α)(f(α)wf(αw)))dw,$ (46)

    and introduce

    $z(α,w):=(f(α)+αf(α))(F(α)F(αw))αf(α)(f(α)wf(αw)).$

    Now, we are going to prove that under conditions (B2), (B3), for all $\alpha \in (\theta_c, 1]$, $w \in [0, 1]$,

    $z(α,w)0,$

    with strict inequality almost everywhere. First, we notice that $z(\alpha, 1) = 0$ and

    $z(α,0)=F(α)(f(α)+αf(α))αf(α)2.$

    Then we compute

    $wz=αf(αw)(f(α)+αf(α))+αf(α)f(αw)+α2wf(α)f(αw)=α2wf(α)f(αw)α2f(αw)f(α).$

    Now, denoting $g(x) = x f'(x) / f(x)$, we get

    $wz=αf(αw)f(α)(g(αw)g(α)).$ (47)

    We are going to make use of the assumptions on $f$ and equation (47) to prove that $z \leq 0$.

    Recall that there exists a unique $\alpha_1 \in (\theta, 1)$ such that $g(\alpha_1) = 1$. If $\alpha \leq \alpha_1$, then for all $w \in [0, \alpha/\theta)$, $g(\alpha w) \leq g(\alpha)$ while for all $w \in (\alpha / \theta, 1]$, $g(\alpha w) \geq g(\alpha)$ (these facts are stated in Lemma .1).

    Hence $w \mapsto z (\alpha, w)$ is increasing on $[0, 1]$. Since $z(\alpha, 1) = 0$, it implies that $z \leq 0$.

    Now, if $\alpha > \alpha_1$, there exists a unique $\beta(\alpha) \in (0, \theta)$ such that $g(\beta(\alpha)) = g(\alpha)$. In this case, if $w \in [0, \alpha / \beta(\alpha)] \cup (\theta, 1]$, $g(\alpha w) \geq g(\alpha)$. If $w \in (\alpha / \beta(\alpha), \theta)$, then $g(\alpha w) < g(\alpha)$. Hence, $\partial_w z \leq 0$ on $[0, \beta(\alpha)/\alpha]$ and $\partial_w z \geq 0$ on $[\beta(\alpha) / \alpha, 1]$. It implies that $z \leq 0$ if, and only if, $z(\alpha, 0) \leq 0$ for all $\alpha > \alpha_1$. This is assumption (B3).

    All in all, we proved that $z \leq 0$ for all $\alpha, w$. Hence $h'(\alpha) > 0$, and there exists a unique $\alpha_0 \in (\theta_c, 1)$ such that $h(\alpha_0) = -2$.

    We conclude that $L_{\alpha}$ is decreasing on $(\theta_c, \alpha_0)$ and increasing on $(\alpha_0, 1]$. Hence $\alpha_0$ is the unique minimum point of $L_{\alpha}$, and the uniqueness of $\alpha_{\pm} (L)$ follows.


    Acknowledgments

    The authors acknowledge partial support from the Programme Convergence Sorbonne Universités / FAPERJ Control and identification for mathematical models of Dengue epidemics and from the Capes-Cofecub project Ma-833 15 Modeling innovative control method for Dengue fever. MS and NV acknowledge partial funding from the ANR blanche project Kibord: ANR-13-BS01-0004 funded by the French Ministry of Research, from the Emergence project from Mairie de Paris, Analysis and simulation of optimal shapes - application to lifesciences and from Inria, France and CAPES, Brazil (processo 99999.007551/2015-00), in the framework of the STIC AmSud project MOSTICAW. JPZ was supported by CNPq grants 302161/2003-1 and 474085/2003-1, by FAPERJ through the programs Cientistas do Nosso Estado, and by the Brazilian-French Network in Mathematics.



    Acknowledgments



    The authors would like to thank Dee Mclean for providing the artwork for Figure 3.

    Conflict of interest



    All authors declare no conflicts of interest in this paper.

    [1] Donaldson L (2009) 150 years of the Annual Report of the Chief Medical Officer: On the state of public health 2008. London: Department of Health.
    [2] Ahangari A (2014) Prevalence of chronic pelvic pain among women: An updated review. Pain Physician 17: E141–147.
    [3] van Wilgen CP, Keizer D (2012) The sensitization model to explain how chronic pain exists without tissue damage. Pain Manag Nurs 13: 60–65. doi: 10.1016/j.pmn.2010.03.001
    [4] Lamvu G (2011) Role of hysterectomy in the treatment of chronic pelvic pain. Obstet Gynecol 117: 1175–1178. doi: 10.1097/AOG.0b013e31821646e1
    [5] Mathias SD, Kuppermann M, Liberman RF, et al. (1996) Chronic pelvic pain: prevalence, health-related quality of life, and economic correlates. Obstet Gynecol 87: 321–327. doi: 10.1016/0029-7844(95)00458-0
    [6] Royal Collge of Obstetricians & Gynaecologists (RCOG) (2015) Scientific Impact Paper No. 46. Therapies Targeting the Nervous System for Chronic Pelvic Pain Relief. RCOG: London
    [7] Lee TT, Yang LC (2008) Pelvic denervation procedures: A current reappraisal. Int J Gynaecol Obstet 101: 304–308. doi: 10.1016/j.ijgo.2008.02.010
    [8] Huber SA, Northington GM, Karp DR (2015) Bowel and bladder dysfunction following surgery within the presacral space: an overview of neuroanatomy, function, and dysfunction. Int Urogynecol J 26: 941–946. doi: 10.1007/s00192-014-2572-x
    [9] Chen FP, Soong YK (1997) The efficacy and complications of laparoscopic presacral neurectomy in pelvic pain. Obstet Gynecol 90: 974–977. doi: 10.1016/S0029-7844(97)00484-5
    [10] Lichten EM, Bombard J (1987) Surgical treatment of primary dysmenorrhea with laparoscopic uterine nerve ablation. J Reprod Med 32: 37–41.
    [11] Daniels JP, Middleton L, Xiong T, et al. (2010) International LUNA IPD Meta-analysis Collaborative Group. Individual patient data meta-analysis of randomized evidence to assess the effectiveness of laparoscopic uterosacral nerve ablation in chronic pelvic pain. Hum Reprod Update 16: 568–576.
    [12] El-Din Shawki H (2011) The efficacy of laparoscopic uterosacral nerve ablation (LUNA) in the treatment of unexplained chronic pelvic pain: a randomized controlled trial. Gynecol Surg 8: 31–39. doi: 10.1007/s10397-010-0612-1
    [13] Daniels J, Gray R, Hills RK, et al. (2009) LUNA Trial Collaboration. Laparoscopic uterosacral nerve ablation for alleviating chronic pelvic pain: A randomized controlled trial. JAMA 302: 955–961.
    [14] Jedrzejczak P, Sokalska A, Spaczynski RZ, et al. (2009) Effects of presacral neurectomy on pelvic pain in women with and without endometriosis. Ginekol Pol 80: 172–178.
    [15] Rouholamin S, Jabalameli M, Mostafa A (2015) The effect of preemptive pudendal nerve block on pain after anterior and posterior vaginal repair. Adv Biomed Res 27: 153. doi: 10.4103/2277-9175.161580
    [16] Chanrachakul B, Likittanasombut P, O-Prasertsawat P, et al. (2001) Lidocaine versus plain saline for pain relief in fractional curettage: A randomized controlled trial. Obstet Gynecol 98: 592–595.
    [17] Naghshineh E, Shiari S, Jabalameli M (2015) Preventive effect of ilioinguinal nerve block on postoperative pain after cesarean section. Adv Biomed Res 4: 229. doi: 10.4103/2277-9175.166652
    [18] Binkert CA, Hirzel FC, Gutzeit A, et al. (2015) Superior hypogastric nerve block to reduce pain after uterine artery embolization: Advanced technique and comparison to epidural anesthesia. Cardiovasc Intervent Radiol 38: 1157–1161. doi: 10.1007/s00270-015-1118-z
    [19] Rapp H, Ledin Eriksson S, Smith P (2017) Superior hypogastric plexus block as a new method of pain relief after abdominal hysterectomy: Double-blind, randomised clinical trial of efficacy. BJOG 124: 270–276. doi: 10.1111/1471-0528.14119
    [20] Fujii M, Sagae S, Sato T, et al. (2002) Investigation of the localization of nerves in the uterosacral ligament: Determination of the optimal site for uterosacral nerve ablation. Gynecol Obstet Invest 54: discussion 16–7. doi: 10.1159/000066289
    [21] Matalliotakis IM, Katsikis IK, Panidis DK (2005) Adenomyosis: What is the impact on fertility? Curr Opin Obstet Gynecol 17: 261–264. doi: 10.1097/01.gco.0000169103.85128.c0
    [22] Desrosiers JA, Faucher GL (1964) Uterosacral block: A new diagnostic procedure. Obstet Gynecol 23: 671–677.
    [23] Rana MV, Candido KD, Raja O, et al. (2014) Celiac plexus block in the management of chronic abdominal pain. Curr Pain Headache Rep 18: 394. doi: 10.1007/s11916-013-0394-z
    [24] Soysal ME, Soysal S, Gurses E, et al. (2003) Laparoscopic presacral neurolysis for endometriosis-related pelvic pain. Hum Reprod 18: 588–592. doi: 10.1093/humrep/deg127
    [25] Byrd D, Mackey S (2008) Pulsed radiofrequency for chronic pain. Curr Pain Headache Rep 12: 37–41. doi: 10.1007/s11916-008-0008-3
  • This article has been cited by:

    1. Camille Pouchol, Emmanuel Trélat, Enrique Zuazua, Phase portrait control for 1D monostable and bistable reaction–diffusion equations, 2019, 32, 0951-7715, 884, 10.1088/1361-6544/aaf07e
    2. Pierre-Alexandre Bliman, Daiver Cardona-Salgado, Yves Dumont, Olga Vasilieva, Implementation of control strategies for sterile insect techniques, 2019, 314, 00255564, 43, 10.1016/j.mbs.2019.06.002
    3. Martin Strugarek, Hervé Bossin, Yves Dumont, On the use of the sterile insect release technique to reduce or eliminate mosquito populations, 2019, 68, 0307904X, 443, 10.1016/j.apm.2018.11.026
    4. Benoît Perthame, Martin Strugarek, Cécile Taing, Selection–mutation dynamics with asymmetrical reproduction kernels, 2022, 222, 0362546X, 112947, 10.1016/j.na.2022.112947
    5. Michel Duprez, Romane Hélie, Yannick Privat, Nicolas Vauchelet, G. Buttazzo, E. Casas, L. de Teresa, R. Glowinski, G. Leugering, E. Trélat, X. Zhang, Optimization of spatial control strategies for population replacement, application toWolbachia, 2021, 27, 1292-8119, 74, 10.1051/cocv/2021070
    6. Yantao Shi, Bo Zheng, Modeling Wolbachia infection frequency in mosquito populations via a continuous periodic switching model, 2023, 12, 2191-950X, 10.1515/anona-2022-0297
    7. L. Roques, T. Boivin, J. Papaïx, S. Soubeyrand, O. Bonnefon, Dynamics of Aedes albopictus invasion insights from a spatio-temporal model, 2023, 1387-3547, 10.1007/s10530-023-03062-y
    8. Samia Ben Ali, Mohamed Lazhar Tayeb, Nicolas Vauchelet, Influence of the competition in the spatial dynamics of a population of Aedes mosquitoes, 2025, 421, 00220396, 208, 10.1016/j.jde.2024.12.002
    9. Nicolas Vauchelet, On a reaction–diffusion system modeling strong competition between two mosquito populations, 2025, 1972-6724, 10.1007/s40574-025-00492-5
  • Reader Comments
  • © 2019 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(4834) PDF downloads(614) Cited by(0)

Article outline

Figures and Tables

Figures(3)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog