Citation: Micaela Morettini, Christian Göbl, Alexandra Kautzky-Willer, Giovanni Pacini, Andrea Tura, Laura Burattini. Former gestational diabetes: Mathematical modeling of intravenous glucose tolerance test for the assessment of insulin clearance and its determinants[J]. Mathematical Biosciences and Engineering, 2020, 17(2): 1604-1615. doi: 10.3934/mbe.2020084
[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] | American Diabetes Association, 2. Classification and diagnosis of diabetes: Standards of medical care in diabetes-2019, Diabetes Care, 42 (2019), S13-S28. |
[2] | C. Kim, Maternal outcomes and follow-up after gestational diabetes mellitus, Diabetic Med., 31 (2014), 292-301. |
[3] | B. E. Metzger, T. A. Buchanan, D. R. Coustan, A. de Leiva, D. B. Dunger, D. R. Hadden, et al., Summary and recommendations of the Fifth International Workshop-Conference on Gestational Diabetes Mellitus, Diabetes Care, 30 (2007), S251-S260. |
[4] | A. Kautzky-Willer, R. Prager, W. Waldhäusl, G. Pacini, K. Thomaseth, O. F Wagner, et al., Pronounced insulin resistance and inadequate β-cell secretion characterize lean gestational diabetes during and after pregnancy, Diabetes Care, 20 (1997), 1717-1723. |
[5] | C. S. Göbl, L. Bozkurt, T. Prikoszovich, C. Winzer, G. Pacini, A. Kautzky-Willer, Early possible risk factors for overt diabetes after gestational diabetes mellitus, Obstet. Gynecol., 118 (2011), 71-78. |
[6] | A. Tura, A. Grassi, Y. Winhofer, A. Guolo, G. Pacini, A. Mari, et al., Progression to type 2 diabetes in women with former gestational diabetes: Time trajectories of metabolic parameters, PLoS One, 7 (2012), e50419. |
[7] | L. Bozkurt, C. S. Göbl, A. Tura, M. Chmelik, T. Prikoszovich, L. Kosi, et al., Fatty liver index predicts further metabolic deteriorations in women with previous gestational diabetes, PLoS One, 7 (2012), e32710. |
[8] | M. Morettini, C. Castriota, C. Göbl, A. Kautzky-Willer, G. Pacini, L. Burattini, et al., Glucose effectiveness from short insulin-modified IVGTT and its application to the study of women with previous gestational diabetes mellitus, Diabetes Metab. J., in press. |
[9] | D. C. Polidori, R. N. Bergman, S. T. Chung, A. E. Sumner, Hepatic and extrahepatic insulin clearance are differentially regulated: Results from a novel model-based analysis of intravenous glucose tolerance data, Diabetes, 65 (2016), 1556-1564. |
[10] | E. Van Cauter, F. Mestrez, J. Sturis, K. S. Polonsky, Estimation of insulin secretion rates from C-peptide levels: Comparison of individual and standard kinetic parameters for C-peptide clearance, Diabetes, 41 (1992), 368-377. |
[11] | G. Pacini, G. Tonolo, M. Sambataro, M. Maioli, M. Ciccarese, E. Brocco, et al., Insulin sensitivity and glucose effectiveness: Minimal model analysis of regular and insulin-modified FSIGT, Am. J. Physiol. Metab., 274 (2017), E592-E599. |
[12] | S. E. Kahn, R. L. Prigeon, D. K. McCulloch, E. J. Boyko, R. N. Bergman, M. W. Schwartz, et al., Quantification of the relationship between insulin sensitivity and β-cell function in human subjects: Evidence for a hyperbolic function, Diabetes, 42 (1993), 1663-1672. |
[13] | A. Tura, A. Mari, T. Prikoszovich, G. Pacini, A. Kautzky-Willer, Value of the intravenous and oral glucose tolerance tests for detecting subtle impairments in insulin sensitivity and beta-cell function in former gestational diabetes, Clin. Endocrinol., 69 (2008), 237-243. |
[14] | A. Mari, A. Tura, G. Pacini, A. Kautzky-Willer, E. Ferrannini, Relationships between insulinsecretion after intravenous and oral glucose administration in subjects with glucose toleranceranging from normal to overt diabetes, Diabet. Med., 25 (2008), 671-677. |
[15] | A. Caumo, L. Luzi, First-phase insulin secretion: Does it exist in real life? Considerations on shape and function, Am. J. Physiol. Endocrinol. Metab., 287 (2004), E371-E385. |
[16] | I. Blumer, E. Hadar, D. R. Hadden, L. Jovanovič, J. H. Mestman, M. H. Murad, et al., Diabetes and pregnancy: An endocrine society clinical practice guideline, J. Clin. Endocrinol. Metab., 98 (2013), 4227-4249. |
[17] | R. Retnakaran, Y. Qi, C. Ye, A. J. G. Hanley, P. W. Connelly, M. Sermer, et al., Hepatic insulin resistance is an early determinant of declining β-cell function in the first year postpartum after glucose intolerance in pregnancy, Diabetes Care, 34 (2011), 2431-2434. |
[18] |
S. D. Mittelman, G. W. Van Citters, S. P. Kim, D. A. Davis, M. K. Dea, M. Hamilton-Wessler, et al., Longitudinal compensation for fat-induced insulin resistance includes reduced insulin clearance and enhanced β-cell response, Diabetes, 49 (2000), 2116-2125. doi: 10.2337/diabetes.49.12.2116
![]() |
[19] | M. Ader, D. Stefanovski, S. P. Kim, J. M. Richey, V. Ionut, K. J. Catalano, et al., Hepatic insulin clearance is the primary determinant of insulin sensitivity in the normal dog, Obesity, 22 (2014), 1238-1245. |
[20] | M. O. Goodarzi, C. D. Langefeld, A. H. Xiang, Y. I. Chen, X. Guo, A. J. G. Hanley, et al., Insulin sensitivity and insulin clearance are heritable and have strong genetic correlation in mexican americans, Obesity, 22 (2014), 1157-1164. |
[21] | K. Ohashi, M. Fujii, S. Uda, H. Kubota, H. Komada, K. Sakaguchi, et al., Increase in hepatic and decrease in peripheral insulin clearance characterize abnormal temporal patterns of serum insulin in diabetic subjects, NPJ Syst. Biol. Appl., 4 (2018), 14. |
[22] | J. Ling, L. Ge, D. H. Zhang, Y. Wang, Z. Xie, J. Tian, et al., DPP-4 inhibitors for the treatment of type 2 diabetes: A methodology overview of systematic reviews, Acta Diabetol., 56 (2019), 7-27. |
[23] | R. M. Goldenberg, L. Berard, Adding prandial GLP-1 receptor agonists to basal insulin: A promising option for type 2 diabetes therapy, Curr. Med. Res. Opin., 34 (2018), 1-10. |
[24] | A. Tura, G. Pacini, Y. Yamada, Y. Seino, B. Ahrén, Glucagon and insulin secretion, insulin clearance, and fasting glucose in GIP receptor and GLP-1 receptor knockout mice, Am. J. Physiol. Integr. Comp. Physiol., 316 (2019), R27-R37. |
[25] | A. Tura, R. Bizzotto, Y. Yamada, Y. Seino, G. Pacini, B. Ahrén, Increased insulin clearance in mice with double deletion of glucagon-like peptide-1 and glucose-dependent insulinotropic polypeptide receptors, Am. J. Physiol. Integr. Comp. Physiol., 314 (2018), R639-R646. |
[26] | A. Shah, M. M. Holter, F. Rimawi, V. Mark, R. Dutia, J. McGinty, et al., Insulin clearance after oral and intravenous glucose following gastric bypass and gastric banding weight loss, Diabetes Care, 42 (2019), 311-317. |
[27] | M. Morettini, F. Di Nardo, L. Ingrillini, S. Fioretti, C. Göbl, A. Kautzky-Willer, et al., Glucose effectiveness and its components in relation to body mass index, Eur. J. Clin. Invest., 49 (2019), e13099. |
[28] | G. Toffolo, R. N. Bergman, D. T. Finegood, C. R. Bowden, C. Cobelli, Quantitative estimation of beta cell sensitivity to glucose in the intact organism. A minimal model of insulin kinetics in the dog, Diabetes, 29 (1980), 979-990. |
[29] | F. Di Nardo, M. Mengoni, M. Morettini, MATLAB-implemented estimation procedure for model-based assessment of hepatic insulin degradation from standard intravenous glucose tolerance test data, Comput. Methods Programs Biomed., 110 (2012), 215-225. |
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 |