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
[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 |
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
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
In [4, 32] it was obtained
an expression for the reaction term
$∂tp−σΔp=f(p)$ | (1) |
in terms of the following biological parameters:
$f(p)=δdsp−shp2+(1+sh−(1−sf)(1−μδ+μ))p+(1−sf)1−μδ−1shp2−(sf+sh)p+1.$ | (2) |
Bistable reaction terms are such that
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
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(x):=∫x0f(y) dy.$ | (4) |
Since we have assumed
$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.27day−1, sf=0.1, μ=0, sh=0.8,δ=0.3/0.27=10/9 and σ=877m2.day−1.$ | (5) |
In particular we obtain the profiles for
We will always assume
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
We say that a radially symmetric function
The following result gives a criterion on the initial data to guarantee invasion.
Theorem 2.1.
Let us assume that
$∂tp−σΔp=f(p),$ | (6) |
$p(t=0,x)=p0(x)≥vα(|x|),$ |
then
$Rα=√σinfρ∈Γ√1−ρd(1−ρ)21(∫α0(1−1−ραx)df(x)dx)+,$ | (7) |
where
In one dimension, we have the sharper estimate Supp
$Lα=√σ2∫α0dv√F(α)−F(v).$ | (8) |
Remark 1. Clearly, the set
$∫α0(1−1−ραx)df(x)dx>0,$ |
since
(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
We recall the definition of a "ground state" as a positive stationary solution
$−Δv=f(v)$ |
that decays to
Proposition 1. In dimension
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
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
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
$∃α∈(θc,1], ∃τ0∈R, ∀τ∈Rd, XτXτ+N0≥vα(τ+τ0),$ |
where
Now, we assume that we have a fixed number of mosquitoes to release, say
Due to the above simplification, the set of releases profiles ("RP") for
a total of
$RPdk(N):={τ↦Nkk∑i=1e−(τ−τi)22σi(2πσi)d/2, with τi∈[−L,L]d, σi∈[σ0−ϵ,σ0+ϵ]},$ | (10) |
where
The basic requirement for a release profile is that
We use uniform measure on
According to our estimate, the success probability satisfies
$P[Non-extinction after releasing N mosquitoes at k locations ]≥P[Xτ(ω) satisfies (NEC)],$ |
where
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
More precisely, we define for a domain
$Pdk(N,Ω):=M{(xi)1≤i≤k,∃α∈(θc,1),∃x0∈Ω,x0+BRα⊂Ω and ∀x∈x0+BRα, Nkk∑i=1Gσ,d(x−xi)≥α},$ | (11) |
where
Fixing the number of mosquitoes per release and letting the number of releases go to
Proposition 2.
Let
$Pdk(kN,Ω)→k→∞1.$ | (12) |
Proof. There are two ingredients for the proof: First, we minimize a Gaussian at
$‖y‖≤√2σlog(2)⟹e−‖y‖2/2σ≥1/2.$ |
Now, when we pick uniformly in a compact set the centers of balls of fixed radius
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
Remark 2. We could have been a little more precise, and get an estimate for the expected value of the number
Therefore we should expect
In fact, any
Corollary 1. For any
$Pdk(kN,Ωc)→k→∞1.$ |
Proof. Let
We consider in this section the proof of Theorem 2.1 in any dimension.
The case
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|2−F(u(x)))dx.$ | (13) |
It is straightforward to see that if
$ddtE[p]=−∫Rd(σΔp+f(p))2 dx≤0.$ |
Thus,
For any
For any
$Jd(r0,R0,α,ϕ):=E[αϕ]|Sd−1|=α2σ∫∞0rd−1|∇ϕ(r)|2dr−(rd0dF(α)+∫R0r0rd−1∫αϕ(r)0f(s)dsdr).$ | (14) |
Now, we use our specific choice of non-increasing radially symmetric function
$Jd(ρ,R0,α):=Rd0(σdR201−ρd(1−ρ)2−F(α)ρdd−1−ρα∫α0(1−1−ραx)d−1F(x)dx),$ | (15) |
where
$Jd(ρ,R0,α)=Rd0(σdR201−ρd(1−ρ)2−∫α0(1−1−ραx)df(x)dx).$ |
We choose
$∫α0(1−1−ραx)df(x)dx>0$ | (16) |
Then the energy
$(R(d)α(ρ))2=σ1−ρd(1−ρ)21∫α0(1−1−ραx)df(x)dx,$ | (17) |
which is a rational fraction in
We examine in particular formula (17) in the case
$U(α):=F(α)−1α∫α0F(x)dx,V(α):=1α∫α0F(x)dx.$ | (18) |
Since
$R(1)α(ρ)=α√σ√(1−ρ)(V(α)+ρU(α)),$ | (19) |
under the constraint
Finally,
$R(1),∗α:=R(1)α(ρ∗1(α))=2√σα√U(α)F(α),$ | (20) |
with
Remark 3. We emphasize that
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
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.
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
Definition 3.1.
For
$vα(x)=uα(|x|)+:=max{0,uα(|x|)} .$ |
From Proposition 3 and Definition 3.1 we have that
Proof. Local existence is granted by Cauchy-Lipschitz theorem.
Then, we multiply Equation (21) by
$σ2((u′α)2)′+(F(uα))′=0,$ |
which implies (since
$σ2(u′α)2=F(α)−F(uα).$ |
Recall that
Because
$√σduαdx=−√2(F(α)−F(uα)).$ | (22) |
Then,
$dχαdω=−√σ2(F(α)−F(ω)),$ |
so that,
$χα(ω)=∫αω√σ2(F(α)−F(v))dv.$ | (23) |
Thus the function
$1√F(α)−F(v)∼v→α1√f(α)1√α−v,1√F(α)−F(v)∼v→0+{1v√−2f′(0) if α=θc,1√F(α) if α>θc.$ |
Therefore
Proposition 4.
The limit bubble
Proof. The function
$σ2(u′θc)2=F(θc)−F(uθc)=−F(uθc).$ |
Hence,
$√σu′θc=−√−2F(uθc) on R+.$ |
Moreover, for small
As a consequence, as
$y′=−√−f′(0)y,$ |
that is
Proof of Theorem 2.1 in dimension d=1.
Let
We first notice that the
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
(ⅰ)
(ⅱ) if
(ⅲ)
Let
(a)
(b) there exists
$limt→∞pλ(t,x)={0 uniformly in R(0≤λ<λ∗),uθc(x−x0) uniformly in R(λ=λ∗),1 locally uniformly in R(λ>λ∗).$ |
In our case, we define
Our construction of a critical
We first compute the energy of the critical
$E[vα]=∫R(σ2|v′α|2−F(vα)) dx.$ |
From Equation (21), we have
$E[vα]=∫Lα−Lα(σ|v′α|2−F(α))dx=2∫Lα0σ|v′α|2 dx−2LαF(α).$ |
Performing the change of variable
$∫Lα0|v′α|2dx=∫α0v′α(v−1α(v)) dv=1√σ∫α0√2(F(α)−F(v))dv,$ |
where we use Equation (22) for the last equality.
Finally, using the expression of
$E[vα]=2√σ∫α0F(α)−2F(v)√2(F(α)−F(v)) dv.$ |
To emphasize the difference between the two sufficient conditions, we observe that
when
$E[vθc]=2√σ∫θc0√−2F(v) dv>0.$ |
By continuity of
Lemma 3.3. The
Remark 5. In particular, the energy estimate alone does not imply invasiveness of the
Figure 3 gives a numerical illustration of the fact that
In this section we discuss a specific release protocol, with a total of
In the case of a single release (
Proposition 5.
Let
(ⅰ) There exists
(ⅱ) There exists
Moreover,
Part (ⅰ) of Proposition 5 asserts that if we fix the total number
Remark 6. If
$N<√2πσN0θc1−θc.$ |
Equivalently, the density of mosquitoes at the center of the single release location
Proof of Proposition 5.
Both the introduction profile given by the fraction
$NGσ(Tσ,N(p))NGσ(Tσ,N(p))+N0=p,$ |
and
${Tσ,N(p)=√2σ√log(NN0√2πσ1−pp),χα(p)=√σ2∫αpdv√F(α)−F(v).$ | (24) |
By construction, the following equivalence holds
$∀τ∈R+, NGσ(τ)NGσ(τ)+N0≥uα(τ)⟺∀p s.t. 0≤p≤α, χα(p)≤Tσ,N(p).$ |
Using (24) this rewrites as
$log(NN0√2πσ)≥(∫αpdv2√F(α)−F(v))2−log(1−pp),∀p∈[0,α].$ | (25) |
From (25), we define
$Jα(p):=log(p)−log(1−p)+(∫αpdv2√F(α)−F(v))2,$ | (26) |
$I(σ,N):=log(N√2πσN0).$ | (27) |
For any given
$∀p∈[0,α], Jα(p)≤I(σ,N).$ | (28) |
We study the function
$J′α(p)=1p(1−p)−1√F(α)−F(p)∫αpdv2√F(α)−F(v),$ |
and we may compute
$jα:=maxp∈[0,α]Jα(p),j∗:=minα∈(θc,1]jα.$ |
Thus there exists
Remark 7. With parameter values from (5),
the expected number of mosquitoes to release is huge, since we need to have
If we space the
$Xτ=Nkk−1∑i=0Gσ(τ+Lα(−1+2ik−1)).$ |
Within a fairly good approximation, this is the case if
$∀τ∈[−Lα,Lα],XτXτ+N0≥α.$ |
This holds in particular if
$N≥˜N(k,α,σ)=N0√2πσ2α1−αkeL2α2σ(k−1)2.$ |
If we fix
It is straightforward, keeping in mind that
$j∗(k):=minα∈(θc,1)α1−αeL2α/(2σ(k−1)2).$ |
and find the minimal (in view of our sufficient criterion) value
Lemma 4.1.
For
$˜N∗(k,σ)=N0√2πσ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
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
First step: fixing
Moreover, we assume that our
$Gσ(y):=1√2πσe−y2/2σ, $ |
and
$G=Nkk∑i=1Gσ(⋅−xi).$ |
We define
$P(σ,Nk,(xi)1≤i≤k,L0,α):=min[−Lα+L0,Lα+L0]G$ | (30) |
Then, the probability of success for the release of
$Pk(N,L)=P[∃L0∈R, ∃α∈(θc,1), P(σ,Nk,(xi)1≤i≤k,L0,α)≥α1−αN0].$ | (31) |
Here, the probability
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
If there is
$Nk1√2πσ≥α1−αN0$ |
and
(ⅰ)
(ⅱ)
then
$GG+N0≥vα(⋅−xm+xl2) .$ |
We notice that the constant
Proof. This property relies on the simple computation of the sum of two
Indeed, considering the sum of two Gaussian
$ξ(x)=1√2πσ(e−(x+h)22σ+e−(x−h)22σ)=2e−h22σGσ(x)cosh(xhσ).$ |
Then, recalling that
$12eh22σσξ′(x)=−xGσ(x)cosh(xhσ)+hGσ(x)sinh(xhσ)12eh22σσ2ξ"(x)=(h2+x2−1σ)Gσ(x)cosh(xhσ)−2hxGσ(x)sinh(xhσ).$ |
As a consequence, the sign of
$γ(x):=h2+x2−2hxtanh(xhσ)−1σ.$ |
We notice that
Now, we examine the necessary condition
$x=htanh(xhσ).$ |
This is true for
So, for
We may use this property to prove Proposition 6.
By condition (ⅰ) the above lower-bound
holds between
$GG+N0≥α,$ |
on
$G(x−xm+xl2)G(x−xm+xl2)+N0≥α≥vα(x−xm+xl2).$ |
As a consequence, we may translate the generic inequality (SP) into:
$P1k(N,(−L,L))=Pk(N,L)≥P[∃α∈(θc,11+N0Nk√2πσ),∃1≤l<m≤k,xm−xl≥2Lα and ∀l≤j≤m−1,xj+1−xj≤2√2log(2)√σ]$ | (32) |
Then, we define
$L∗:=minθc<α≤11+N0Nk√2πσLα,$ |
and equivalently estimate (32) reads
$Pk(N,L)≥P[∃1≤l<m≤k,xm−xl≥2L∗ and maxl≤j≤m−1(xj+1−xj)≤2√2log(2)√σ].$ | (33) |
The study of the minimization of
Remark 8. Note that for this estimate, we only consider initial data that are above a characteristic function at level
Remark 9. It is easy to check that our estimate yields
$k≥1√2log(2)minθc<α≤1∫α0dv√2(F(α)−F(v)).$ |
Specific discussion for
In order to compute analytically the right-hand-side in (33), we may introduce the following notations:
●
$τk(u,v)=(v−u)k+k!.$ |
●
●
Remark 10. Back to problem (33), we recover the problem of estimating
$R∗:=minα∫α0dv√F(α)−F(v) ,$ |
and
We want to under-estimate the probability of success with
Proposition 7.
Let
$βλ,R∗k(−L,χ)=k∑i=k0k−i+1∑j=1∫χ−R∗−L∫min(χ,u+(k−1)λ)u+R∗γλi(u,v)$ |
$(τj−1(−L,u−λ)−βλ,R∗j−1(−L,u−λ))τk−(i+j−1)(v+λ,χ)dvdu.$ | (34) |
Proof. The idea is simple: we count each "positive initial data", that is an ordered
We shall use the index
$βλ,R∗k(−L,χ)=∫[−L,χ]k1{y1≤y2≤⋯≤yk}1{(y1,…,yk)∈Bλ,R∗k(−L,χ)}dy1…dyk.$ | (35) |
Now, we split:
$1{(y1,…,yk)∈Bλ,R∗k(−L,χ)}=k∑i=k0k−i+1∑j=11{yi+j−1−yj≥R∗}j+i−2∏l=j1{yl+1−yl≤λ}1{(y1,…,yj−1)∉Bλ,R∗j−1(−L,χ)}1{yj−yj−1>λ}1{yi+j−yi+j−1>λ}.$ | (36) |
This identity requires some explanations.
It comes from the partition of
$1{(y1,…,yj−1)∉Bλ,R∗j−1(−L,χ)}1{yj−1≤yj}1{yj−yj−1>λ}=1{(y1,…,yj−1)∉Bλ,R∗j−1(−L,yj−λ)},$ |
with the obvious convention that
In addition, for
$∫[−L,χ]k−(i+j−1)1{yi+j−1≤⋯≤yk}1{yi+j−yi+j−1>λ}dyi+j…dyk=τk−(i+j−1)(yi+j−1+λ,χ)=(χ−yi+j−1−λ)k−(i+j−1)+(k−(i+j−1))!.$ |
Combining these results, and using (35) and (36) yields
$βλ,R∗k(−L,χ)=k∑i=k0k−i+1∑j=1∫χ−L…∫χxi+j−21{yj+i−1−yj≥R∗}j+i−2∏l=j1{0≤yl+1−yl≤λ}(τj−1(−L,yj−λ)−βλ,R∗j−1(−L,yj−λ))τk−(i+j−1)(yi+j−1+λ,χ)dyj…dyi+j−1,$ | (37) |
with conventions
We assume
$βλ,R∗k(−L,χ)=k∑i=k0k−i+1∑j=1∫χ−R∗−L∫min(χ,u+(k−1)λ)u+R∗γλi(u,v)(τj−1(−L,u−λ)−βλ,R∗j−1(−L,u−λ))τk−(i+j−1)(v+λ,χ)dvdu,$ |
where
Now, we may give an explicit formula for
$γλi+2(u,v)=∫u+λu∫u1+λu1…∫ui−1+λui−11v≥ui≥v−λdui…du1,$ |
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+2(u,v)=λi+i+1∑k=1(−1)ki!((ik−1)(v−u−kλ)i++(−1)i+1(i−1k−1)(kλ−(v−u))i+).$ | (39) |
Proof. Obviously,
$γλ3(u,v)=λ+(v−u−2λ)+−(λ−(v−u))+−(v−u−λ)+$ |
Then, using (38) again proves (39) by induction.
Remark 11.
For
$βλ,R∗k0(−L,L)=∫L−R∗−L∫min(L,u+(k0−1)λ)u+R∗γλk0(u,v)dvdu.$ | (40) |
Then by (39) we know
$βλ,R∗k0(−L,L)=∫L−(k0−1)λ−L∫(k0−1)λR∗(λk0−2+k0−1∑k=1(−1)k(k0−2)!((k0−2k−1)(w−kλ)k0−2++(−1)k0−1(k0−3k−1)(kλ−w)k0−2+))dwdu+∫L−R∗L−(k0−1)λ∫L−uR∗γλk0(u,u+w)dwdu.$ | (41) |
Clearly, the first integral in the right-hand side of (41) may be written as
$(2L−(k0−1)λ)f1(λ,R∗),$ |
where
$f2(λ,R∗):=∫(k0−1)λR∗∫zR∗(λk0−2+k0−1∑k=1(−1)k(k0−2)!((k0−1k−1)(w−kλ)k0−2++(−1)k0−1(k0−3k−1)(kλ−w)k0−2+))dwdz.$ |
In particular, it appears that it does not depend on
For
$βλ,R∗k0(−L,χ)=∫χ−(−L)R∗∫zR∗γλk0(0,w)dwdz,$ |
and notice that our expressions are consistent since
$βλ,R∗k0(−L,−L+(k0−1)λ)=∫−L+(k0−1)λ−(−L)R∗∫zR∗γλk0(0,w)dwdz=f2(λ,R∗).$ |
All in all,
$βλ,R∗k0(−L,χ)={0 if χ+L≤R∗∫χ−(−L)R∗∫zR∗γλk0(0,w)dwdz if χ+L∈(R∗,(k0−1)λ),(χ+L−(k0−1)λ))f1+f2 if χ+L>(k0−1)λ$ | (42) |
(This is an affine function for
Then, we obtain a bound on the probability of success with
$Pk0(L)≥βλ,R∗k0τk0(−L,L)=k0!(2L)k0((2L−(k0−1)λ)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
We find
$2ˆL=λ(⌈⌈R∗λ⌉⌉+1)−k0k0−1f2(λ,R∗)f1(λ,R∗).$ |
We may note that introducing the non-negative and non-decreasing function
$Γλ,R∗k(z):=∫zR∗γλk(0,w)dw$ |
we get
$f1(λ,R∗)=Γλ,R∗k0((k0−1)λ),f2(λ,R∗)=∫(k0−1)λR∗Γλ,R∗k0(z)dz.$ |
As a consequence,
$2ˆL≥k0k0−1R∗.$ |
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 (
Figure 5 shows the probability profile as a function of the size
However, for small (relatively to
Our numerical values are somehow consistent with field experiments: typically, the space between release points is less than
The factor
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 "
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
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
Some other questions remain open. For instance:
in one dimension, we considered releases in
$ˆL≥R∗1+⌈⌈R∗λ⌉⌉⌈⌈R∗λ⌉⌉.$ | (43) |
It is a numerical conjecture that the optimal value of
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
In this appendix we investigate sufficient conditions for the uniqueness of a minimal radius among the
Let
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
$g(x):=xf′(x)/f(x).$ | (44) |
By simple computation we have
Lemma .1.
Under assumption (B0), (B2),
$g(α1)=1.$ |
We add the following assumption:
$∀α>max(θc,α1),F(α)(f(α)+αf′(α))≤α(f(α))2.$ |
Now, we recall the
$Lα=√σ∫α0dv√2(F(α)−F(v)).$ |
Proposition 8.
Under conditions (B0), (B1), the bistable (in the sense of (3)) function
If in addition (B2), (B3) hold,
then there exists a unique
$Lα0=minαLα,$ |
and for all
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
Generally, we can check that (B2)-(B3) hold for the classical bistable function
$f′(x)+xf"(x)=−9x2+4(1+θ)x−θ.$ |
Then (B2) is equivalent to
$(9x2−4(1+θ)x+θ)x(x2−(1+θ)x+θ)≤x(3x2−2(1+θ)x+θ)2⟺−13(1+θ)x2+10θx−5θ(1+θ)≤−12(1+θ)x2+6θx−4θ(1+θ)⟺0≤(1+θ)x2−4θx+θ(1+θ).$ |
The discriminant of this second-order polynomial is
Now, we want to check (B3). To do so we compute
$F(x)=−14x4+1+θ3x3−θ2x2.$ |
Then (B3) is equivalent to
$x2(14x2−1+θ3x+θ2)(4x2−3(1+θ)x+2θ)≤x3(x2−(1+θ)x+θ)2⟺x2(1+θ)(2−34−43)+θ2x+θ(1+θ)(2−32−23)≤0.$ |
Then we recall that
$θ24−θ(1+θ)29=θ4(θ−49(1+θ)2)<0.$ |
Hence simplest bistable functions of the form
Proof. Without loss of generality we assume
$Lα=∫α0(1√F(α)−F(v)−1√f(α)(α−v))dv+∫α0dv√f(α)(α−v)=1√f(α)(∫α0(√f(α)√F(α)−F(v)−1√α−v)dv+2√α)$ |
Hence
$ddαLα=1√αf(α)+12√f(α)∫α0(1(α−v)3/2−(f(α)F(α)−F(v))3/2)dv,$ |
which is a continuous function from
Then,
$1√α+12∫α0(1(α−v)3/2−(f(α)F(α)−F(v))3/2)dv=0.$ |
For
$h(α):=∫10(1(1−w)3/2−(αf(α)F(α)−F(αw))3/2)dw.$ | (45) |
Then
We compute
$h′(α)=−32∫10(α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
$z(α,w)≤0,$ |
with strict inequality almost everywhere.
First, we notice that
$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
$∂wz=αf(αw)f(α)(g(αw)−g(α)).$ | (47) |
We are going to make use of the assumptions on
Recall that there exists a unique
Hence
Now, if
All in all, we proved that
We conclude that
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.
[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
![]() |
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 |