A negative answer to a conjecture arising in the study of selection-migration models in population genetics
Elisa Sovrano




[t1]This work was supported by the auspices of
Grup­po Na­zio­na­le per l'Anali­si Ma­te­ma­ti­ca, la Pro­ba­bi­li­tà e le lo­ro
Appli­ca­zio­ni (GNAMPA) of the Isti­tu­to Na­zio­na­le di Al­ta Ma­te­ma­ti­ca (INdAM).



Department of Mathematics, Computer Science and Physics, University of Udine, 
via delle Scienze 206, 33100 Udine, Italy 


We deal with the study of the evolution of the allelic frequencies, at a single locus, 
for a population distributed continuously over a bounded habitat. 
We consider evolution which occurs under the joint action of selection and arbitrary migration, 
that is independent of genotype, in absence of mutation and random drift. 
The focus is on a conjecture, that was raised up in literature of population genetics, 
about the possible uniqueness of polymorphic equilibria, which are known as clines, 
under particular circumstances. 
We study the number of these equilibria, making use of topological tools, 
and we give a negative answer to that question by means of two examples. 
Indeed, we provide numerical evidence of multiplicity of positive solutions
for two different Neumann problems satisfying the requests of the conjecture.


 
Migration Selection Cline Polymorphism Indefinite weight Neumann problem  
[2010] 92D25 35K57 34B18.








§ INTRODUCTION

Population genetics is a field of biology concerning the genetic structure inside the populations. 
Its main interest is the understanding of evolutionary processes that make the complexity of the Nature so intriguing. 
One of the main causes of the diversity among organisms are the changes in the genetic sequence. 
The genome evolution is influenced by selection, recombination, harmful and beneficial mutations, among others. 
This way, population genetics becomes helpful in order to tackle a very broad class of issues, 
from epidemiology, animal or plant breeding, demography and ecology.

The birth of the “modern population genetics” can be traced back at the need to interlace the Darwin's evolution theory with the Mendelian laws of inheritance. 
This has taken place in the 1920s and early 1930s, when Fisher, Haldane and Wright 
have developed mathematical models in order to analyze how the natural selection, along with other factors, 
would modify the genetic composition of a population over time. 
Accordingly, an impressive moment on the history of this field of the genetics is the “Sixth International Congress of Genetics” in Ithaca, 
where all the three fathers of the genetical theory of evolution have given a presentation of their pioneering works, see <cit.>.

Mathematical models of population genetics can be described by relative genotypic frequencies or 
relative allelic frequencies, that may depend on both space and time. 
A common assumption is that individuals mate at random in a habitat 
(which can be bounded or not) with respect to the locus under consideration. 
Furthermore, the population is usually considered large enough so that frequencies can be treated as deterministic. 
This way, a probability is associated to the relative frequencies of genotypes/alleles.
The dynamics of gene frequencies are the result of some genetic principles along with several environmental influences, 
such as selection, segregation, migration, mutation, recombination and mating, that lead to different evolutionary processes 
like adaptation and speciation, see <cit.>. 

Amongst these influences, by natural selection we mean that some genotypes enjoy a survival or reproduction advantage over other ones. 
This way, the genotypic and allelic frequencies change in accord to the proportion of progeny to the next generation of 
the various genotypes which is named fitness. 
Thinking to model real-life populations, we have to take into account which is unusual that the selection factor acts alone. 
Since every organism lives in environments that are heterogeneous, 
another considerable factor is the natural subdivision of the population that mate at random only locally. 
Thus, migration is often considered as a factor that affects the amount of genetic change. 
There are two different ways, in order to model the dispersion or the migration of organisms: 
one is of discrete type and the other one of a continuous nature. 
If the population size is sufficiently large and the selection is restricted to a single locus with two alleles, 
then deterministic models continuous in time and space lead to mathematical problems 
which involve a single nonlinear partial differential equation of reaction-diffusion type.

In this direction, a seminal paper was given by <cit.>. 
In that work, it was studied the frequency of an advantageous gene for a uniformly distributed population in a one-dimensional habitat 
which spreads through an intensity constant selection term. 
Accordingly, a mathematical model of a cline was built up as a non-constant stationary solution 
of the nonlinear diffusion equation in question. 
The term cline was coined by J. Huxley in <cit.>: 
“Some special term seems desirable to direct attention to variation within groups, and I propose the term cline, 
meaning a gradation in measurable characters.”
One of the major causes of cline's occurrence is the migration or the selection which favors an allele in a region of the habitat 
and, a different one in another region. 
The steepness of a cline is considered as an indicating character of the level of the geographical variation.
Another contribution comes from J.B. Haldane in <cit.>, who has studied the cline's stability by considering 
as a selection term a stepwise function which depends on the space and changes its sign. 

Some meaningful generalizations of these models have been performed, for example, 
in <cit.> by introducing a linear spatial dependence in the selection term; in
<cit.> by considering a different diffusion term that can model barriers and in 
<cit.> by taking into account population not necessarily uniformly distributed and 
terms of migration-selection that depend on both space and time.
During the past decade, these mathematical treatments have opened the door to a great amount of works that 
investigated the existence, uniqueness and stability of clines. 
Since a complete list of references of further analysis on clines is out of the scope of this work, 
we limit ourselves to cite some of the earliest contributions in the literature that have inspired the succeeding ones, 
see for instance <cit.>.

Understanding the processes that act in order to have non-constant genetic polymorphisms 
(i.e., loci that occur in more than one allelic form) is an important challenge in population genetics. 
In the present work, we deal with a class of diallelic migration-selection models in continuous space and time 
introduced by W.H. Fleming in <cit.> and D. Henry in <cit.>.
We focus on a conjecture stated in <cit.>,
that, for such a kind of reaction diffusion equations, guess the uniqueness of a cline (instead of the existence of multiple ones). 
In a one-dimensional setting, we will give a negative answer to that conjecture, 
by providing two examples with multiplicity of non-constant steady states. 
This type of treatment is inspired by the result, about multiplicity of positive solutions for indefinite weight problems 
with Dirichlet boundary conditions, performed in <cit.>.
Although the problem approach has a topological feature, numerical simulations are given in order to support it.

The plan of the paper is the following. 
In Section <ref>, we present the class of migration-selection models 
considered and the state of the art which has lead to the formulation of the conjecture of Lou and Nagylaki,
with reference to the genetical and mathematical literature.
In Section <ref>, we build up two examples giving a negative answer to this conjecture.
In Section <ref>, we conclude with a discussion.




§ MIGRATION-SELECTION MODEL: THE CONJECTURE OF LOU AND NAGYLAKI


2.equation To ease understanding the conjecture raised up in <cit.>, 
we introduce some notations. 
We also provide an overview of the classical migration-selection model, continuous in space and in time,
of a population in which the genetic diversity occurs in one locus with two alleles, A_1 and A_2. 

Let us consider a population continuously distributed in a bounded habitat, say Ω. 
In our context, genetic diversity is the result only of the joint action of dispersal within Ω 
and selective advantage  for some genotypes, 
so that, no mutation nor genetic drift will be considered. 
This way, the gene frequencies, after random mating, are given by the Hardy-Weinberg relation. 
The genetic structure of the population is measured by the frequencies p(x,t) and 
q(x,t):=(1-p(x,t)) at time t and location x∈Ω of A_1 and A_2, respectively. 
Thus, by the assumptions made, the mathematical formulation of this migration-selection model 
leads to the following semilinear parabolic PDE:

    ∂ p/∂ t = Δ p + λ w(x) f(p)    in Ω×]0,∞[,

where Δ denotes the Laplace operator and Ω⊆ℝ^N is a bounded open connected set, with N≥1,
whose boundary ∂Ω is C^2.
The term λ w(x) f(u) models the effect of the natural selection. 
More in detail, the real parameter λ > 0 plays the role of the ratio of the selection intensity and
the function w∈ L^∞(Ω) represents the local selective advantage (if w(x) > 0), 
or disadvantage (if w(x) < 0), of the gene at the position x∈Ω. 
Moreover, following <cit.> and <cit.>, the nonlinear term we treat is a general function f: [0,1]→ℝ of class C^2 satisfying 

    f(0) = f(1) =0,    f(s) > 0   ∀ s∈ ]0,1[ ,    f'(0) > 0 > f'(1).
    (f_*)

We also impose that there is no-flux of genes into or out of the habitat Ω, 
namely we assume that

    ∂ p/∂ν = 0    on ∂Ω×]0,∞[,

where ν is the outward unit normal vector on ∂Ω.
Since p(t,x) is a frequency, then we are interested only in positive solutions of (<ref>)–(<ref>) such that 0≤ p≤ 1.

By the analysis developed in <cit.>, 
we know that, if the conditions in (f_*) hold and 0≤ p(·,0) ≤ 1 in Ω, then 0≤ p(x,t) ≤ 1 
for all (x,t)∈Ω×]0,∞[ and equation (<ref>) defines a dynamical system in 

    X:={p∈ H^1(Ω): 0 ≤ p(x) ≤ 1,  a.e. in Ω},

where H^1(Ω) is the standard Sobolev space of integrable functions whose first derivative is also square integrable.
Moreover, the stability of the solutions is determined by the equilibrium solutions in the space X. 
Clearly, a stationary solution of the problem  is a function p(·) satisfying 0 ≤ p≤ 1, 

    -Δ p = λ w(x) f(p)    in Ω

and the Neumann boundary condition 

    ∂ p/∂ν = 0    on ∂Ω.


Notice that p≡ 0 and p≡ 1 are constant trivial solutions 
to problem , that correspond to monomorphic equilibria, 
namely when, in the population, the allele A_2 or A_1, respectively, is gone to fixation. 
So, one is interested in finding non-constant stationary solutions or, in other words, polymorphic equilibria. 
Indeed, our main interest is the existence of clines for system (<ref>)–(<ref>). 

The maintenance of genetic diversity is examined by seeking for the existence 
of polymorphic stationary solutions/clines, 
that are solutions p(·) to system 

    -Δ p = λ w(x) f(p)     in Ω,
    ∂ p/∂ν = 0     on ∂Ω,
    (𝒩_λ)

with 0 < p(x) < 1 for all x∈Ω.

In this respect, the assumption f(s) > 0 for every s > 0 implies that a necessary condition 
for positive solutions of problem (𝒩_λ) is that the function w changes its sign. 
In fact, by integrating (<ref>) over Ω, we obtain

    0=∫_ΩΔ p +λ w(x) f(p)  dx=
    ∫_∂Ω∂ p/∂ν dx+λ∫_Ω w(x) f(p) dx =
    λ∫_Ω w(x) f(p) dx.

Notice that we can see the function w in (𝒩_λ) 
as a weight term which attains both positive and negative values, 
so that such a kind of system is usually known as problem with indefinite weight. 

It is a well-known fact that the existence of positive solutions of (𝒩_λ) 
depends on the sign of 

    w̅:=∫_Ωw(x) dx.

Indeed, for the linear eigenvalue problem -Δ p(x)=λ w(x)p(x), 
under Neumann boundary condition on Ω, the following facts hold:
if w̅<0, then there exists a unique positive eigenvalue having an associated eigenfunction which does not change sign; 
on the contrary, if w̅≥0 such an eigenvalue does not exist and 0 is the only non-negative eigenvalue 
for which the corresponding eigenfunction does not vanish, see <cit.>.

Furthermore, under the additional assumption of concavity for the nonlinearity:

    f”(s)≤0,   ∀ s>0,

it follows that, if w̅<0,
then there exists λ_0>0 such that for each λ>λ_0 
problem  has at most one nonconstant stationary positive solution 
(i.e., cline) which is asymptotically stable, see <cit.>.

After these works a great deal of contributions appeared in order 
to complement these results of existence and uniqueness on population genetics, see for instance <cit.>; 
or to consider also unbounded habitats as done in <cit.>; 
or even to treat more general uniformly elliptic operators, as in <cit.>.
Taking into account these works, in <cit.> the migration-selection model 
with an isotropic dispersion, that is identified with the Laplacian operator,
was generalized to an arbitrary migration, which involves a strongly uniformly elliptic differential operator of second order 
(see also <cit.> for the derivation of this model as a continuous approximation of the discrete one).

By modeling single locus diallelic populations, there is an interesting family of nonlinearities 
which satisfies the conditions in (f_*) and 
allows to consider different phenotypes of alleles, A_1 and A_2. 
This family can be obtained by considering
the map f_k:ℝ^+→ℝ^+ such that

    f_k(s):= s(1-s)(1+k-2ks),

where -1≤ k≤1 represents the degree of dominance of the alleles independently of the space variable, see <cit.>.
In this special case, if k=0 then the model does not present any kind of dominance, 
instead, if k=1 or k=-1 then the allelic dominance is relative to A_1, in the first case, and to A_2 in the second one 
(the last is also equivalent to said that A_1 is recessive).
In view of this, we can make mainly the following two observations.


In the case of no dominance, i.e. k=0, from (<ref>) we have f_0(s)=s(1-s) which is a concave function. 
Therefore, we can enter in the settings considered by <cit.>.
So if w(x)>0 on a set of positive measure in Ω and w̅<0, 
then for λ sufficiently large there exists a unique positive non trivial equilibrium of the equation
∂ p/ ∂ t = Δ u + λ w(x) p(1-p) for every (x,t)∈Ω×]0,∞[ under 
the boundary condition (<ref>).



In the case of completely dominance of allele A_2, i.e. k=-1, from (<ref>)
we have f_-1(s)=2s^2(1-s) which is not a concave function. 
Thanks to the results in <cit.>, if w(x)>0 on a set of positive measure in Ω and w̅<0, 
then for λ sufficiently large there exist at least two positive non trivial equilibrium of the equation
∂ p/ ∂ t = Δ p + λ w(x) 2 p^2(1-p) for every (x,t)∈Ω×]0,∞[ under 
the boundary condition (<ref>).


We observe that the map s↦ f_0(s)/s is strictly decreasing with f_0(s) concave. 
On the contrary, the map s↦ f_-1(s)/s is not strictly decreasing with f_-1(s) not concave. 
Thus, from Remark <ref> and Remark <ref>, it arises a natural question which involves the possibility 
to weaken the concavity assumption (<ref>) further to the monotonicity of the map s↦ f(s)/s, in order to get 
uniqueness results of nontrivial equilibria for problem . 

This is still an open question, firstly appeared in <cit.>, known as the “conjecture of Lou and Nagylaki”. 




Conjecture “Suppose that w(x)>0 on a set of positive measure in Ω and such that w̅=∫_Ω w  dx<0. 
If the map s↦ f(s)/s is monotone decreasing in ]0,1[, 
then  has at most one nontrivial equilibrium p(t,x) with 0<p(0,x)<1 for every x∈Ω̅,
which, if it exists, is globally asymptotically stable.” <cit.>.



The study of existence, uniqueness and multiplicity of positive solutions for nonlinear indefinite weight problems 
is a very active area of research, starting from the Seventies, 
and several types of boundary conditions along with a wide variety of nonlinear functions, classified according to growth conditions, 
were taken into account.
Several authors have addressed this topic, see <cit.>, 
just to recall the first main papers dedicated.
Instead, the recent literature about multiplicity results for positive solutions of indefinite weight problems with 
Dirichlet or Neumann boundary conditions is really very rich. 
In order to cover most of the results achieved with different techniques so far, 
we give reference of the following bibliography <cit.>. 

Nevertheless, as far as we known, there is no answer about the conjecture of Lou and Nagylaki.
It is interesting to notice that the study of the concavity of f(s) versus the monotonicity of f(s)/s has significance also 
in the investigation on the uniqueness of positive solutions for a particular class of indefinite weight problems with Dirichlet boundary conditions.
More in detail, these problems involve positive nonlinearities which have linear growth at zero and sublinear growth at infinity,
namely

    -Δ p = λ w(x) g(p)     in Ω,
    
    p = 0     on ∂Ω,
    (𝒟_λ)

where g:ℝ^+→ℝ^+ is a continuous function satisfying 

    g(0)=0 ,   g(s) > 0   ∀ s>0 ,   lim_s→0^+g(s)/s>0 ,  lim_s→+∞g(s)/s=0.
    (g_*)


The state of the art on this topic refers mainly on two papers.
From the results achieved in <cit.>, it follows that, 
if (g_*) holds and, moreover, the map s→ g(s)/s is strictly decreasing, 
then there exist at most one positive solution of (𝒟_λ) only 
if the weight function w(x) > 0 for a.e. x∈Ω.
On the other hand, from <cit.>, if the conditions in (g_*) are satisfied for a smooth concave nonlinear term g  
and the weight w is a smooth and changing sign function, then there exists at most one positive solutions of (𝒟_λ).
Therefore, if the weight function is positive, then the hypothesis of Brezis-Oswald, concerning the monotonicity of g(s)/s, 
is more general than the requirement of Brown-Hess about the concavity of g(s).

At this point one could query whether something similar to the conjecture of Lou and Nagylaki happens 
also for this family of Dirichlet problems.
This was done in <cit.>, where it was shown that the monotonicity of the map s↦ g(s)/s is not enough 
to guarantee the uniqueness of positive solutions for problems as in (𝒟_λ) with an indefinite weight. 
Through numerical evidence, more than one positive solution has been detected for an exemplary 
two-point boundary value problem (𝒟_λ).




§ MULTIPLICITY OF CLINES: THE CONJECTURE HAS NEGATIVE ANSWER


3.equation In this section we look at the framework of the conjecture of Lou and Nagylaki. 
So, from now on we tacitly consider a nonlinear function f: [0,1]→ℝ of class C^2 
which is not concave, satisfies (f_*) and is such that the map s↦ f(s)/s is strictly decreasing. 

We concentrate on the one-dimensional case N=1 
and we take as a habitat an open interval Ω:=]ω_1,ω_2[ with ω_1,ω_2∈ℝ 
such that ω_1<0<ω_2.
This type of habitats, confined to one-dimensional spaces, have an intrinsic interest in modeling phenomena which occur, 
for example, in neighborhoods of rivers, sea shores or hills, see <cit.>.
As in  <cit.>, we assume that the weight term w is step-wise.
Hence, let us consider the following class of indefinite weight functions 

    w(x):=
    
    -α    x∈[ω_1,0[,
    
    1	    x∈]0,ω_2],

such that 

    w̅ = -ω_1α +ω_2<0,
 
with w̅ defined as in (<ref>).
In these settings, the indefinite Neumann problem (𝒩_λ) reads as follows

    p” + λ w(x) f(p) = 0,
    
    p'(ω_1)=0=p'(ω_2),

with 0 < p(x) < 1 for all x∈[ω_1,ω_2].

Inspired by the results in <cit.>, we will consider two particular functions f 
in order to provide a negative reply to the conjecture under examination.
In both cases, we are going to use a topological argument, that is called shooting method, 
and, with the aid of some numerical computations, we give evidence
of multiplicity of positive solutions for the corresponding problems in (<ref>).
The shooting method relies on the study of the deformation of planar continua under the action of the vector field associated to 
the second order scalar differential equation in (<ref>), whose formulation, in the phase-plane (u,v),
is equivalent to the first order planar system

    u '=v, 	
    
    v '= -λ w(x) f(u).

Solutions p(·) of problem (<ref>) we are looking for are also solutions (u(·),v(·)) of system (<ref>),
such that v(ω_1)=0=v(ω_2).

We set the interval [0,1] contained in the u-axis as follows

    ℒ_{v=0}:={(u,v)∈ℝ^2: 0≤ u≤ 1,  v=0}.

This way, as a real parameter r ranges between 0 and 1, we are interested in the solution, 
(u(·  ;ω_1,(r,0)),v(·  ;ω_1,(r,0))), of the Cauchy problem
with initial conditions

    u(ω_1)=r,
    
    v(ω_1)=0,

such that (u(ω_2 ;ω_1,(r,0)),v(ω_2 ;ω_1,(r,0)))∈ℒ_{v=0}.
Hence, let us consider the planar continuum Γ obtained by shooting ℒ_{v=0} forward
from ω_1 to ω_2, namely

    Γ:={(u(ω_2;r),v(ω_2;r))∈ℝ^2 : r∈[0,1]	}.

We define the set of the intersection points between this continuum and the segment [0,1] contained in the u-axis, as

    𝒮:=	Γ∩ℒ_{v=0}.

Then, there exists an injection form the set of the solutions p(·) of (<ref>) such that 0 < p(x) < 1 for all x∈[ω_1,ω_2]
and the set 𝒮∖( {(0,0)}∪{(1,0)}).

More formally, we denote 
by ζ(· ; ω_0,z_0)=(u(· ; ω_0,z_0),v(· ; ω_0,z_0)) 
the solution of (<ref>) with ω_0∈[ω_1,ω_2] and initial condition
ζ(ω_0 ; ω_0,z_0)=z_0=(u_0,v_0)∈ℝ^2.
The uniqueness of the solutions for the associated initial value problems guarantee that 
the Poincaré map associated to system (<ref>) is well defined.
Recall that, for any τ_1,τ_2∈[ω_1,ω_2], the Poincaré map for system (<ref>), 
denoted by Φ_τ_1^τ_2,
is the planar map which at any point z_0=(u_0,v_0)∈ℝ^2 associates the point (u(τ_2),v(τ_2))
where (u(·),v(·)) is the solution of (<ref>) with (u(τ_1),v(τ_1))=z_0.
Notice that Φ^τ_2_τ_1 is a global diffeomorphism of the plane onto itself.

Under these notations, the recipe of the shooting method is the following.
A solution p(·) of (<ref>) such that 0<p(x)<1 for all x∈[ω_1,ω_2]
is identified by a point (c,0)∈ℒ_{v=0} whose image through the action of the Poincaré map,
say C:=Φ_ω_1^ω_2((c,0))∈Γ, belongs to ℒ_{v=0}.
This way, the solution p(·) of the Neumann problem with p(ω_1)=c is found 
looking at the first component of the map 

    x↦Φ_ω_1^x((c,0))=(u(x),v(x)),

since, by construction, p'(ω_1)=v(ω_1)=0 and p'(ω_2)=v(ω_2)=0.
This means that the set 𝒮 is made by points such that, each of them determines univocally 
an initial condition, of the form (<ref>), for which the solution (u(·),v(·)) of the Cauchy problem 
associated to (<ref>) verifies v(ω_1)=0=v(ω_2).

The study of the uniqueness of the clines is based on the study, in the phase plane (u,v), of the qualitative properties 
of the shape of the continuum Γ which is the image of ℒ_{v=0} 
under the action of the Poincaré map Φ_ω_1^ω_2. 
More in detail, we are interested in find real values c∈]0,1[ such that 

    Φ^ω_2_ω_1((c,0))∈Φ_ω_1^ω_2(ℒ_{v=0})∩ℒ_{v=0}.

Indeed, our aim is looking for values c∈]0,1[ such that the point C=Φ^ω_2_ω_1((c,0)) belongs to 𝒮.

So, let us show now that there exist Neumann problems as in (<ref>) that admit more than one positive solution.
Namely, there exist more than one polymorphic stationary solution for the equation:

    ∂ p/∂ t=p”+λ w(x)f(p).

Roughly speaking, if Γ crosses the u-axis more than one time, out of the points (0,0) and (1,0),
then #(𝒮∖( {(0,0)}∪{(1,0)}))>1 and so, we expect a result of non-uniqueness of clines
for equation (<ref>).




 §.§ First example

Taking into account the definition of the functions in (<ref>), given a real parameter h>0, let us consider the family of maps 
f̂_h:[0,1]→ℝ of class C^2 such that 

    f̂_h(s):=s(1-s)(1-h s+h s^2).
 
By definition f̂(0)=0=f̂(1). Moreover, to have s↦f̂_h(s)/s monotone decreasing in ]0,1[ it is sufficient to assume 0<h≤ 3. 
If the parameter h ranges in ]0,3], then it is straightforward to check that f̂_h is not concave and 
f̂_h(s)>0 for every s∈]0,1[.

Let us fix h=3. Then, in this case, all the conditions in (f_*) are verified and 
f̂_3(s)=s(1-s)(1-3s+3s^2). 
As a consequence, we point out the following result of multiplicity.


Let f:[0,1]→ℝ be such that 

    f(s):=s(1-s)(1-3s+3s^2).

Assume w:[ω_1,ω_2]→ℝ be defined as in (<ref>) 
with α=1, ω_1=-0.21 and ω_2=0.2. 
Then, for λ=45 the problem (<ref>) 
has at least 3 solutions such that 0 < p(x) < 1 for all x∈[ω_1,ω_2].


Notice that w̅=-0.01<0, so we are in the hypotheses of the conjecture.
Now we follow the scheme of the shooting method, in order to detect three polymorphic stationary solutions
for the equation (<ref>). 
This approach, with the help of numerical estimates, will enable us to prove Proposition <ref>.


In the phase-plane (u,v), Figure <ref> shows the existence of at least four points 
(r_i,0)∈ℒ_{v=0} with i=1,…,4 such that, 
by defining their images through the Poincaré map Φ^ω_2_ω_1 as
R_i:=(R_i^u,R_i^v)=Φ^ω_2_ω_1((r_i,0))∈Γ
for every i∈{1,…,4},
the following conditions  

    R_i^v<0	 for i=1,3,  R_i^v>0	 for i=2,4,

are satisfied. 
This is done, for example, with the choice of the values r_1=0.1, r_2=0.4, r_3=0.65 and r_4=0.75.
The solutions of the Cauchy problems associated to system (<ref>), with initial conditions (r_i,0) for i=1,…,4, 
assume at x=ω_2 the values R_1=(0.230, -0.066), R_2= (0.922, 0.165), R_3=(0.790,  0.036) 
and R_4=(0.533, 0.055), truncated at the third significant digit. 
Therefore, we have R_1^v<0<R_2^v, R_2^v>0>R_3^v and R_3^v<0<R_4^v.
Then, by a continuity argument (that means an application of the Mean Value Theorem), 
there exist at least three real values c_1,c_2 and c_3 such that

    r_j<c_j<r_j+1  and  
    C_j:=Φ^ω_2_ω_1((c_j,0))∈𝒮∖( {(0,0)}∪{(1,0)}),

for every j∈{1,…,i-1}. So, let us see how to find such values.



The curve Γ is the result of the integration of several system of differential equation (<ref>),
with initial conditions taken within a uniform discretization of the interval [0,1], 
followed by the interpolation of the approximated values of each solution ζ(x ;ω_0,z_0) at x=ω_2.
Hence, Γ represents the approximation of the image of the interval [0,1] 
under the action of the Poincaré map Φ^ω_2_ω_1.

As the Figure <ref> suggests, the projection of Γ on its first component 
is not necessarily contained in the interval [0,1], which includes the only values of biological pertinence. 
Nonetheless, this does not avoid the existence of solutions p(·) of the problem (<ref>) such that 0 <p(x)<1 
for all x∈[ω_1,ω_2].
This way, by means of a fine discretization of ℒ_{v=0}, 
we have found the approximate values of the intersection points C_j∈Γ∩ℒ_{v=0}, 
with j=1,2,3. In this case they are:
C_1=(0.273,0), C_2=(0.601,0) and C_3=(0.833,0), truncated at the third significant digit (see Figure <ref>).
The intersection points between ℒ_{v=0} and its image Γ 
through the Poincaré map Φ^ω_2_ω_1, namely C_j with j=1,2,3, are in agreement with the previous predictions.

At last, we computed the values c_1=0.125, c_2=0.479 and c_3=0.683, which verify the required conditions (<ref>). 
For j=1,2,3, in Figure <ref> are represented the trajectories of the solutions of the initial value problem

    p”+λ w(x)f(p)=0,
    
    p(ω_1)=c_j,
    
    p'(ω_1)=0,

that, by construction, satisfy p'(ω_2)=0. 
We observe also that the values of each solution p(·) of the three different
initial value problems ranges in ]0,1[ as desired.

Once found the values c_j with j=1,2,3, a numerically result of multiplicity of clines is achieved.
Indeed, in Figure <ref>, we display the approximation of the three nontrivial stationary solution p(·)
of equation (<ref>) that are identified by the points
C_j∈(𝒮∖( {(0,0)}∪{(1,0)})),
with j=1,2,3. 







 §.§ Second example

We refer now to the application given in <cit.> and we adapt it to our purposes. 
So we consider, the nonlinear term f̃:ℝ^+→ℝ^+ defined by

    f̃(s):=(10 s e^-25s^2+s/|s|+1).

It is straightforward to check that f̃ is not concave and the map s→f̃(s)/s is strictly decreasing. 
Moreover, f̃(0)=0 and f̃(s)>0 for every s>0, but f̃ does not take value zero in s=1, 
since f̃(1)=10 e^-25+1≠0.
To satisfy all the conditions in (f_*), it is sufficient to multiply f̃ by the term arctan(m(1-x)) with m>0.
This way, the following result holds.


Let f:[0,1]→ℝ be such that 

    f(s):=(10 s e^-25s^2+s/|s|+1) arctan(10-10s).

Assume w:[ω_1,ω_2]→ℝ be defined as in (<ref>) 
with α=2.4, ω_1=-0.255 and ω_2=0.6. 
Then, for λ=3 the problem (<ref>) 
has at least 3 solutions such that 0 < p(x) < 1 for all x∈[ω_1,ω_2].


Notice that, under the assumptions of Proposition <ref>, the hypotheses of the conjecture are now all satisfied since w̅=-0.012<0. 
To prove the existence of at least three clines for the equation (<ref>),
we exploit again the shooting method.

So, our main interest is in finding real values r_i∈]0,1[ with i∈ℕ such that,
given R_i:=(R_i^u,R_i^v)=Φ^ω_2_ω_1((r_i,0)), it follows

    R_i^v<0	 for i=2ℓ +1,  R_i^v>0	 for i=2ℓ,  with ℓ∈ℕ.


In this case, the features of the nonlinear term along with the joint action of the indefinite weight 
give rise to an involved deformation of the segment ℒ_{v=0}. 
Nevertheless, looking at Figure <ref>, we can see that there exist more than one intersection point between the continuum
Γ and the u-axis such that their abscissa is contained in the interval open ]0,1[.

This way, the previous observation suggests us the following analysis. By choosing the values
 r_1=0.01, r_2=0.1, r_3=0.45 and r_4=0.9 we compute the points R_i for i=1,…,4. 
All the results achieved are truncated at the third significant digit and so we obtain 
R_1^v=-0.639<0, R_2^v=2.160>0, R_3^v<-0.036 and R_4^v=1.392>0.
The numerical details are thus represented in Figure <ref>.




At this point, an application of the Intermediate Value Theorem guarantees the existence of at least three initial conditions 
(c_j,0) with j=1,2,3, such that each respective solution of the initial value problem (<ref>) is also a positive solution 
of the Neumann problem (<ref>) we are looking for. 
Indeed, the values c_1=0.436, c_2=0.776 and c_3=0.854 satisfy the conditions in (<ref>).
Finally, we display the approximation of the three nontrivial stationary solution p(·)
of equation (<ref>) in Figure <ref>. 




§ DISCUSSION

5.equation 
When the selection gradient is described by a piecewise constant coefficient,
we have studied the conjecture proposed in <cit.> within a finite and one-dimensional environment. 

Summing up: we have found multiplicity of positive solutions
for two different indefinite Neumann problems defined as in (<ref>) where the nonlinear term is an application
f:[0,1]→ℝ which assumes two particular forms, the one in (<ref>) or that in (<ref>).
In our examples, the nonlinearity f is a function of class C^2 such that

    f(0) = f(1) =0,    f(s) > 0   ∀ s∈ ]0,1[ ,    f'(0) > 0 > f'(1),
    (f_*)

and

    f is not concave, s↦ f(s)/s is strictly decreasing.(H)

Hence, uniqueness of positive solutions in general is not guaranteed for indefinite Neumann problems 
whose nonlinear term f is a function satisfying (f_*) and (H), and
the indefinite weight w is defined on a bounded domain Ω with ∫_Ωw(x) dx<0.
Nevertheless, in the case of the family of functions f_k(s)=s(1-s)(1+k-2ks)
with k∈[-1,1], there are issues in this direction that have not yet been answered. 

So, a question still open is the following: 
under the action of gene flow, what is the minimal set of assumptions 
under which a selection gradient will maintain a unique gene frequency cline?

The delicate matter of the comparison between the concavity 
versus a condition about monotonicity arises also in other context than the Neumann one.
With this respect, indefinite weight problems under Dirichlet boundary conditions 
have been considered in <cit.> where an example of multiplicity of positive solutions was given.
As far as we know, the mathematical literature lacks of a rigorous multiplicity result in both the two cases.
This way, we see how these issues deserve to be studied in deep for both a natural mathematical 
and genetical interest.



§ ACKNOWLEDGMENTS

I thank prof. Reinhard Bürger for inspiring me to work on this problem.  
I am also very grateful to profs. Fabio Zanolin, Carlota Rebelo and Alessandro Margheri 
for providing useful comments 
that greatly improved the manuscript.



*
elsart-num-sort
   


