Modeling eating disorders in young people

Nonlinear Analysis: Real World Applications 53 (2020) 103064

Contents lists available at ScienceDirect

Nonlinear Analysis: Real World Applications www.elsevier.com/locate/nonrwa

Modeling eating disorders in young people Andrea Giacobbe a , Giuseppe Mulone a ,∗,1 , Wendi Wang b ,2 a

Dipartimento di Matematica e Informatica, Città Universitaria, Viale A. Doria, 6, 95125, Catania, Italy b Key Laboratory of Eco-environments in Three Gorges Reservoir Region, School of Mathematics and Statistics, Southwest University, Chongqing 400715, PR China

article

info

Article history: Received 9 July 2019 Received in revised form 7 October 2019 Accepted 24 October 2019 Available online xxxx Keywords: Nonlocal Reproduction number Body mass index Stability Distribution Eating disorders

abstract A mathematical model is proposed to simulate the eating disorders of bulimia or anorexia. Earlier models are extended to incorporate the body mass index, which plays a key role in the eating attitude of self thinners. The global existence and ultimate boundedness of solutions of the nonlocal model are proved by using estimates of solutions. The basic reproduction number of eating disorder contagion is shown to be the invasion threshold. The testable linear and nonlinear stability conditions are established by Lyapunov functions. Further numerical simulations are given to reveal how self-forces and peer pressures to be thinner affect the emergence and distributions of eating disorders. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction Eating disorders of anorexia nervosa and bulimia nervosa have become serious problems for young people both in developed countries and developing countries such as China. It was estimated that about 1%–3% of adolescent and young adult women in the United States are diagnosed with an eating disorder [1]. It was reported in a review article on eating disorders in Europe that the prevalence of anorexia nervosa can reach an incidence up to 1%–4% and bulimia nervosa could be up to 1%–2% [2]. China has a similar proportion to that of Western societies [3]. These social diseases impact many aspects of life and result in significant functional impairment, emotional distress, and medical problems for the infected people [4]. Mathematical models have been proposed to study the dynamics of anorexia and bulimia in young populations [5,6]. Basically, these models split the population into some compartments: susceptible, infected ∗ Corresponding author. E-mail addresses: [email protected] (A. Giacobbe), [email protected] (G. Mulone), [email protected] (W. Wang). 1 The research that led to the present paper was partially supported by the group GNFM of Istituto Nazionale di Alta Matematica Francesco Severi, by the University of Catania, Grant No. PTRDMI-53722122113, and by Grant 2017YBKNCE of national project PRIN of Italian Ministry for University and Research. 2 The research is supported by the NSF of China (11571284).

https://doi.org/10.1016/j.nonrwa.2019.103064 1468-1218/© 2019 Elsevier Ltd. All rights reserved.

2

A. Giacobbe, G. Mulone and W. Wang / Nonlinear Analysis: Real World Applications 53 (2020) 103064

(with or without the progression stages), and recovered. Note that the body mass index (BMI) is the key quantity to characterize body image and body dissatisfaction, and is one of the major driving forces to reduce body’s weights [7–10]. Furthermore, the influence of body image varies with the body mass index. Therefore, it is important to incorporate this variable into the mathematical models to understand how the distributions of BMI affect the propagation of eating disorder behaviors, and how intervention strategies can be designed to prevent the disease spread on the basis of BMI distributions. The organization of this paper is as follows. In the next section, we present the formulation of mathematical model. Section 3 includes the analysis of basic properties of the system. In Section 4, we derive the basic reproduction number of the disease and show that it is the threshold value of disease invasion. In Section 5, testable conditions for the stability of disease-free steady state and its attractive basin are established by Lyapunov method. Section 6 gives the numerical simulation to demonstrate how the distributions of BMI affect the disease spread. The paper ends with brief discussions. 2. Model formulations We consider a population of young people in which some individuals have eating disorder behaviours (anorexia or bulimia). For simplicity, we use B to denote the collection of individuals with eating disorders. It is assumed that every individual in B group pursues the goal of becoming as thin as possible. The group of young people with the risk to develop into the class of eating disorders is denoted by S, in which all individuals admit normal eating behaviours. In practice, the illness of eating disorders is treated by medication or psychological educations. It is estimated that about half of the individuals can be cured within ten years of their diagnosis and about 30%–50% of them experience relapse [11]. Thus, we let T denote the group of individuals in treatment and R denote the group of cured individuals. It is assumed that an individual in class T stops the eating disorder behaviours, but can relapse or progress into the cured class. We denote by x the body mass index and t the time. Let S(x, t), B(x, t), T (x, t) and R(x, t) be respectively the sizes of susceptible individuals, anorexic or bulimic individuals, individuals in treatment and cured individuals with body mass index x at time t. Assume that the maximal BMI and minimal BMI are x1 and x0 , respectively, and that a person can only live if its BMI is in the interval [x0 , x1 ] outside which life is forbidden. As a result, we impose the boundary conditions: ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ S⏐ = S⏐ = B⏐ = B⏐ = T⏐ = T⏐ = R⏐ = R⏐ = 0. x=x0

x=x1

x=x0

x=x1

x=x0

x=x1

x=x0

x=x1

We now extrapolate how susceptible individuals are affected to enter the group of eating disorders. In the present paper, we focus on the case that susceptible individuals are mainly tempted to become anorexic or bulimic by the images from classmates, and assume that the social influences in body images are neutral to susceptible individuals under positive educations and negative medium forces. Let us consider a weight function g1 (a) which is positive, continuous and satisfies ∫ x1 g1 (a)da = 1. x0

This function describes how infectious individuals with different BMI affect susceptible individuals. The force of infection is modelled by the law ∫ x 1

β

g1 (a)B(a, t)da. x0

The variation of BMI is affected by two sources: random effect and directed effect. The random effect results from ordinary diet fluctuations, whereas the directed transition is induced by the intention to get an ideal body image. The BMI of an individual in B class varies not only from random fluctuations, but also

A. Giacobbe, G. Mulone and W. Wang / Nonlinear Analysis: Real World Applications 53 (2020) 103064

3

from the directed eating behaviours towards the ideal BMI x from the point of view of bulimic or anorexic people. In the present paper, we consider only the case where x = x0 for the convenience of mathematical analysis. 2 2 2 Let d be the random diffusion coefficient of body mass index for each individual. Then d ∂∂xS2 , d ∂∂xB2 , d ∂∂xT2 , 2 and d ∂∂xR2 denote the diffusions of BMI for individuals from the four classes, respectively. Notice that the eating behaviours of thinness are mainly driven by self-forces and peer pressures. In the case where the eating behaviours are only affected by self-forces, we may assume that the directed BMI transition flux is given by Jds := −η0 B,

(2.1)

where η0 > 0 is the self-force diffusion coefficient. For the case where peer pressure is the major force for an individual in group B to reduce weights, we suppose that the directed BMI transition flux is described by Jdp := −η1 B

∫

x1

g2 (a)B(a, t)da,

(2.2)

x0

where η1 > 0 is the peer-pressure diffusion coefficient and g2 is a non-negative, continuous weight function that satisfies ∫ x1 g2 (a)da = 1. x0

This function describes how infectious individuals with different BMI affect individuals in B class to reduce their weights. In the present paper, we assume that the directed BMI transition flux Jd is the additive effect of the two sources, that is, ∫ x1 Jd = Jds + Jdp = −η0 B − η1 B g2 (a)B(a, t)da, (2.3) x0

which captures the main forces to reduce BMI for individuals with eating disorders, and is similar to those in the construction of ecological models [12]. Similarly, the directed BMI transition flux of treated individuals is assumed to be η2 T where η2 > 0 is the coefficient of BMI recovery. Let Λ be the recruitment rate of the population, µ be the removing rate of the population, ϵ be the disease-induced death rate, α1 be the treatment rate, α0 be the relapse rate of treated individuals, and γ be the transition rate of individuals with treatment to the recovered class. Summarising the above hypotheses, we get the mathematical model: ⎧ ∫x ∂2S ∂S ⎪ ⎪ = d + Λ − µS − βS x 1 g1 (a)B(a, t)da, ⎪ ⎪ 2 0 ⎪ ∂t ∂x ⎪ ⎪ ( ) ∂B 2 ⎪ ∫ ∂ B ∂B ⎪ x1 ⎪ ⎪ = d + η + η g (a)B(a, t)da + 2 0 1 x ⎪ 0 ⎪ ∂t ∂x2 ∂x ⎨ ∫ x1 + βS x g1 (a)B(a, t)da + α0 T − (µ + ϵ + α1 )B, 0 ⎪ ⎪ ⎪ 2 ⎪ ∂T ∂ T ∂T ⎪ ⎪ =d − η2 + α1 B − (µ + γ + α0 )T, ⎪ 2 ⎪ ∂t ∂x ∂x ⎪ ⎪ ⎪ 2 ⎪ ⎪ ⎩ ∂R = d ∂ R + γT − µR. ∂t ∂x2

(2.4)

We will focus our analysis on the first three equations that are decoupled from the last equation. Let us scale the variable x by x ˆ = (x − x0 )/(x1 − x0 ) and set dˆ = d(x1 − x0 )−2 , βˆ = β(x1 − x0 ), ηˆ0 = η0 (x1 − x0 )−1 ,

4

A. Giacobbe, G. Mulone and W. Wang / Nonlinear Analysis: Real World Applications 53 (2020) 103064

ηˆ1 = η1 , ηˆ2 = η2 (x1 − x0 )−1 . Dropping the hats we have ⎧ ∫1 ∂2S ∂S ⎪ ⎪ = d 2 + Λ − µS − βS 0 g1 (a)B(a, t)da, ⎪ ⎪ ∂t ∂x ⎪ ⎪ ( ) 2 ⎪ ⎪ ⎨ ∂B = d ∂ B + η + η ∫ 1 g (a)B(a, t)da ∂B + 0 1 0 2 2 ∂t ∂x ∂x ∫1 ⎪ ⎪ ⎪ +βS g (a)B(a, t)da + α0 T − (µ + ϵ + α1 )B, ⎪ 0 1 ⎪ ⎪ ⎪ 2 ⎪ ⎩ ∂T = d ∂ T − η ∂T + α B − (µ + γ + α )T, 2 1 0 ∂t ∂x2 ∂x which are subject to the Dirichlet conditions ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ = B⏐ = S⏐ S⏐ x=1

x=0

x=0

⏐ ⏐ = B⏐

x=1

⏐ ⏐ = T⏐

x=0

⏐ ⏐ = T⏐

= 0,

(2.5)

(2.6)

x=1

and initial conditions S(x, 0) = S0 (x),

B(x, 0) = B0 (x),

T (x, 0) = T0 (x).

(2.7)

3. Preliminary analysis In this section, we present the mathematical analysis for the existence and uniqueness of solution of (2.5) with b. c. (2.6). Let X = C02+α ([0, 1], R3 ) with 0 < α < 1 be the space of twice continuously differentiable functions on [0, 1] with values in R3 that satisfy the Dirichlet conditions (2.6) and whose second partial derivatives are H¨ older continuous of exponent α. In view of the biological context, we consider the dynamics 0 of (2.5)–(2.6) in X+ which is the nonnegative cone of space X. Let X+ be the interior of X+ . Set Ω = (0, 1) and Qτ = Ω × (0, τ ) where τ > 0 is a fixed time. Let m be a nonnegative integer, 0 < α < 1 and 0 < δ < 1. Assume that function u(x, t) and its derivatives up to order m in x are bounded on the closure of Qτ , that the m derivatives with respect to x are H¨older continuous with exponent α and they are m+α,δ continuous in t and H¨ older continuous with exponent δ. Following [13,14], we denote by Cx,t (Qτ ) the space of function u(x, t) with the norm ∥u∥C m+α,δ x,t

] m [ ∑ l m (δ) sup |Dx u| + ⟨Dxm u⟩(α) = x + ⟨Dx u⟩t , l=0

Qτ

where ⟨w⟩(α) x =

sup (x,t),(y,t)∈Qτ

(δ)

⟨w⟩t

=

sup (x,t),(x,t1 )∈Qτ

|w(x, t) − w(y, t)| , α |x − y| |w(x, t) − w(x, t1 )| δ

|t − t1 |

.

2+α,1+δ Moreover, we let Cx,t (Qτ ) be the space of function u(x, t) with the finite norm:

∥u∥C 2+α,δ + ∥ut ∥C α,δ . x,t

x,t

In what follows, the subscripts x and t in the spaces and norms will be suppressed when there is no confusion. Lemma 3.1. Let (S(x, t), B(x, t), T (x, t)) be a solution of system (2.5) with the boundary condition (2.6) and (S(·, 0), B(·, 0), T (·, 0)) = (S0 (·), B0 (·), T0 (·)) ∈ X+ . Then each component of (S(x, t), B(x, t), T (x, t)) remains nonnegative for t > 0 in the existence interval of solution. Furthermore, if there is an 0 < x0 < 1 0 such that S0 (x0 ) > 0 and B0 (x0 ) > 0, then (S(·, t), B(·, t), T (·, t)) stays in X+ for t > 0 in the existence interval of solution.

A. Giacobbe, G. Mulone and W. Wang / Nonlinear Analysis: Real World Applications 53 (2020) 103064

5

Proof . We present only the discussions to the second case where S0 and B0 are positive at x0 . Note that the first equation and third equation of (2.5) are standard to use the maximum principle of parabolic equation [15]. It is easily seen that S(x, t) > 0 for t > 0 and x ∈ (0, 1). We claim that B(x, t) > 0 for t > 0 and x ∈ (0, 1). Indeed, the second equation of (2.5) can be rewritten as ∫ 1 ∂B + LB = βS g1 (a)B(a, t)da + α0 T, ∂t 0 where

( ) ∫ 1 ∂2B ∂B LB = −d 2 − η0 + η1 g2 (a)B(a, t)da + (µ + ϵ + α1 )B. ∂x ∂x 0

Thus we have

∂B + LB ≥ 0 ∂t for small t > 0. It follows from the maximum principle [15] that B(x, t) > 0 for small t > 0 and x ∈ (0, 1). Then it is easy to conclude from the maximum principle that B(x, t) > 0 for x ∈ (0, 1) and those t > 0 in the existence interval of the solution. The positivity of T follows easily because Eq. (2.5)3 is a linear equation. □ Lemma 3.2. For any initial function in X+ , there exists a unique solution of system (2.5), which satisfies the boundary condition (2.6) and the initial data (2.7). Furthermore, each component of the solution is in C 2+α,1+α/2 (Qτ ) for small τ > 0. Proof . We start with the nonlinear and nonlocal problem: ⎧ ∫ ∫ 1 ∂B ∂B 1 ⎪ ⎪ ⎪ + LB = η g (a)B(a, t)da + βS g1 (a)B(a, t)da + α0 T, 1 2 ⎪ ⎨ ∂t ∂x 0 0 B(x, 0) = B0 (x), ⎪ ⎪ ⏐ ⏐ ⎪ ⎪ ⏐ ⎩ B ⏐⏐ = B⏐ = 0, x=0

(3.1)

x=1

where

∂B ∂2B − η0 + (µ + ϵ + α1 )B. ∂x ∂x2 We claim that for each nonnegative S ∈ C α,α/2 (Qτ ) and nonnegative T ∈ C α,α/2 (Qτ ), there exists a unique nonnegative solution B ∈ C 2+α,1+α/2 (Qτ ) to (3.1). To this end, we consider a closed and convex set in C 2+α,1+α/2 (Qτ ): { } CM = B ∈ C 2+α,1+α/2 (Qτ ) : B(x, 0) = B0 (x), ∥B∥C 1+α,α/2 ≤ M , LB = −d

˜ be the unique solution of the following linear problem: where M = ∥B0 ∥X + 1. For any B ∈ CM , we let B ⎧ ∫ 1 ∫ 1 ˜ ∂B ⎪ ⎪ ˜ = η1 ∂B + L B g (a)B(a, t)da + βS g1 (a)B(a, t)da + α0 T, ⎪ 2 ⎪ ⎨ ∂t ∂x 0 0 (3.2) ˜ 0) = B0 (x), B(x, ⎪ ⎪ ⏐ ⏐ ⎪ ⎪ ⎩B ˜ ⏐⏐ ˜ ⏐⏐ =B = 0. x=0

x=1

˜ = G(B). We show below that G(CM ) ⊂ CM and G is contractive for a Then we define a mapping G by B small τ > 0. First, by the parabolic Schauder theory [16], we have ˜ 2+α,1+α/2 ≤ M1 (∥B0 ∥X + M0 ), ∥B∥ C

(3.3)

A. Giacobbe, G. Mulone and W. Wang / Nonlinear Analysis: Real World Applications 53 (2020) 103064

6

where

1

∫ 1  ∂B  + βS g1 (a)B(a, t)da + α0 T  α,α/2 ≤ M0 , ∂x C 0 0 and M0 , M1 are positive numbers. It is easy to see that there exists a positive constant N1 such that  ∫  η1

g2 (a)B(a, t)da

˜ t) − B(x, ˜ 0)∥ 1+α,α/2 ≤ N1 τ ∥B∥ ˜ 2+α,1+α/2 . ∥B(x, C C As a result, it follows from (3.3) that ˜ t)∥ 1+α,α/2 ≤ ∥B0 ∥ 1+α,α/2 + N1 τ ∥B∥ ˜ 2+α,1+α/2 ∥B(x, C C C ≤ ∥B0 ∥X + N1 T M1 (∥B0 ∥X + M0 ) < M for small τ > 0. This verifies that G maps CM into CM . ˜1 = G(B1 ), B ˜2 = G(B2 ), Next, we show that the map G is contractive in CM . For B1 , B2 ∈ CM , we let B ˜1 − B ˜2 . Then we have u = B1 − B2 and v = B ∂v + Lv = f, ∂t (3.4) v(x, 0) = 0, x ∈ Ω , ⏐ ⏐ ⏐ ⏐ v⏐ = v⏐ = 0, x=0

where

x=1

∫ ∫ ∂B1 1 ∂u 1 f = βS g1 (a)u(a, t)da + η1 g2 (a)u(a, t)da + η1 g2 (a)B2 (a, t)da. ∂x 0 ∂x 0 0 Evidently, there is a constant M2 > 0 such that (  ∂u  )   ∥f ∥L∞ ≤ M2 ∥u∥C 0 +   = M2 ∥u∥C 1,0 , for 0 ≤ t ≤ τ. ∂x C 0 ∫

1

(3.5)

As a result, by the parabolic Lp theory, there exists an M3 > 0 such that ∥v∥W 2,1 ≤ M3 M2 ∥u∥C 1,0

(3.6)

p

for any p > 1. Taking large p and using the Sobolev embedding, we see that there exists an M4 > 0 such that ∥v∥C 1+α,α/2 ≤ M4 ∥v∥W 2,1 ≤ M2 M3 M4 ∥u∥C 1,0 . p

Consequently, we have ∥v∥C 1+α,α/2 ≤ M2 M3 M4 τ 1−α/2 ∥u∥C 1,α/2 ≤ M2 M3 M4 τ 1−α/2 ∥u∥C 1+α,α/2 . Hence, G is contractive in CM if τ > 0 is small. As a result, we conclude that (3.1) admits a unique solution for small τ > 0 and the claim is examined. Now, we let Y = C α,α/2 (Qτ ) × C α,α/2 (Qτ ) × C α,α/2 (Qτ ) with small τ > 0. Then we define a closed and convex set in Y : DM ={(S, B, T ) ∈ Y : S(x, 0) = S0 (x), B(x, 0) = B0 (x), T (x, 0) = T0 (x), (S(·, t), B(·, t), T (·, t)) ∈ X+ , ∥S∥C α,α/2 + ∥B∥C α,α/2 + ∥T ∥C α,α/2 ≤ M }, where M is a positive constant to be determined later. For any fixed (u, v, w) ∈ DM , we consider the following problem in X+ : ∫ 1 ∂S ∂2S = d 2 + Λ − µS − βS g1 (a)v(t, a)da + γB, ∂t ∂x 0 ∫ 1 ∫ 1 ∂B ∂B + LB = η1 g2 (a)B(a, t)da + βu g1 (a)B(t, a)da, ∂t ∂x (3.7) 0 0 ∂T ∂2T ∂T =d − η2 + α1 B − (µ + γ + α0 )T, ∂t ∂x⏐ 2 ∂x ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ S⏐ = S⏐ = 0, B ⏐ = B⏐ = 0, T ⏐ = T⏐ = 0. x=0

x=1

x=0

x=1

x=0

x=1

A. Giacobbe, G. Mulone and W. Wang / Nonlinear Analysis: Real World Applications 53 (2020) 103064

7

Since the equations in S and T are linear, it follows from the arguments in the first step that (3.7) admits a unique solution (S, B, T ) with S, B, T ∈ C 2+α,1+α/2 (Qτ ) for small τ > 0. Then we define a mapping F by (S, B, T ) = F (u, v, w). By similar discussions to those in the first step, we see that there exists a M > 0 such that the mapping F maps DM into itself and is contractive in DM . Therefore, we conclude that F has a unique fixed point in DM , which is the unique solution of system (2.5) with the boundary condition (2.6). Clearly, each component of the solution is in C 2+α,1+α/2 (Qτ ) from the parabolic Schauder estimates. □ Our next objective is to examine that a nonnegative solution (S, B, T ) of (2.5) with the boundary condition (2.6) exists on [0, ∞). To this end, we let N = S + B + T . Then (N, B, T ) solves the problem: ∂N ∂B ∂2N ∂T = d 2 + Q(t) − η2 + Λ − µN − ϵB − γT, ∂t ∂x ∂x ∂x ∫ 1 ∂2B ∂B ∂B g1 (a)B(a, t)da + α0 T − (µ + ϵ + α1 )B, = d 2 + Q(t) + βS ∂t ∂x ∂x 0 ∂T ∂2T ∂T =d − η2 + α1 B − (µ + γ + α0 )T, ∂t⏐ ∂x2⏐ ∂x⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ N⏐ = N⏐ = B⏐ = B⏐ = T⏐ = T⏐ = 0, x=0

x=1

x=0

where

x=1

( Q(t) =

x=0

∫ η0 + η1

1

(3.8)

x=1

) g2 (a)B(a, t)da .

(3.9)

0

Since S, B and T are nonnegative, the solution (N, B, T ) of (3.8) satisfies N ≥ S ≥ 0,

N ≥ B ≥ 0,

N ≥ T ≥ 0.

(3.10)

Henceforth, we show that a solution of (3.8) satisfying (3.10) exists on [0, ∞). First, we state that the population size N is bounded by a positive constant. Lemma 3.3. Let (S, B, T ) be a nonnegative solution of (2.5) with (2.6) which exists on [0, τ0 ). Then there exists a positive constant C0 , which is independent of τ0 , such that ∫ 1 N (x, t)dx < C0 , for 0 ≤ t < τ0 . (3.11) 0

Proof . Let N0 (x) = N (x, 0). Integrating the first equation of (3.8) with respect to x over [0, 1] and using the strong parabolic maximum principle, we get ∫ ∫ 1 d 1 ∂N ⏐⏐1 N (x, t)dx ≤ d N (x, t)dx ⏐ +Λ−µ dt 0 ∂x 0 0 ∫ 1 ≤Λ−µ N (x, t)dx. 0

It follows that

∫ 0

1

N (x, t)dx ≤ e−µt

∫

1

N0 (x)dx + 0

Λ , µ

0 ≤ t < τ0 .

(3.12)

Consequently, there exists a positive constant C0 such that (3.11) holds. □ Lemma 3.4. Let (S, B, T ) be a nonnegative solution of (2.5) with (2.6) which exists on [0, τ0 ). Then there exists a positive constant C1 , which is independent of τ0 , such that ∥N (·, t)∥L∞ (Ω) < C1 , for 0 ≤ t < τ0 .

(3.13)

8

A. Giacobbe, G. Mulone and W. Wang / Nonlinear Analysis: Real World Applications 53 (2020) 103064

Proof . Note that (2.5) can be rewritten as ⎧ ∫1 ∂2S ∂S ⎪ ⎪ = d 2 + Λ − µS − βS 0 g1 (a)B(a, t)da, ⎪ ⎪ ∂t ∂x ⎪ ⎨ ∂N = d ∂ 2 N + Q(t) ∂B − η ∂T + Λ − µN − ϵB − γT, 2 ∂t ∂x ⎪ ∂x2 ∂x ⎪ ⎪ 2 ⎪ ⎪ ∂T ∂T ∂ T ⎩ − η2 =d + α1 B − (µ + γ + α0 )T, ∂t ∂x2 ∂x

(3.14)

where B = N − S − T and Q(t) is defined by (3.9). Multiplying the second equation of (3.14) by N p−1 with p ≥ 2 and then integrating about x over [0, 1], we obtain )2 p ∫ ∫ ( ∫ 1 ( ) 4d(p − 1) 1 ∂N 2 1 d 1 p N dx ≤ − dx + ΛN p−1 − µN p dx p dt 0 p2 ∂x 0 0 (3.15) ∫ 1 p−2 ∂N − (Q(t) − η2 )(p − 1) (B + T )N dx. ∂x 0 By Lemma 3.3, there exists a positive constant C21 such that (Q(t) − η2 )(p − 1) ≤ C21 p,

t ∈ [0, τ0 ).

By computing the maximum of the function ΛN p−1 − µN p , it is easy to see that there exists a positive constant C22 such that ( )p ∫ 1 ( ) Λ , t ∈ [0, τ0 ). ΛN p − µN p+1 dx ≤ C22 µ 0 Hence, from (3.15) we obtain ( )p ∫ 1 d 1 p Λ (3.16) . N dx ≤K + C22 p dt 0 µ where )2 p p ∫ ( ∫ 1 p ∂N 2 4d(p − 1) 1 ∂N 2 2 K=− dx + 2C | |dx. N 21 p2 ∂x ∂x 0 0 In order to estimate K, we need Young’s inequality ab ≤ ϵam + c(ϵ)bn , a ≥ 0, b ≥ 0, m > 0, n > 0, ϵ > 0

(3.17)

where c(ϵ) = (ϵm)−n/m /n and 1/m+1/n = 1. We also need the Gagliardo–Nirenberg interpolation inequality for one-dimensional space [17,18] 1,r ∥w∥Ls ≤ C∥w∥δW 1,r ∥w∥1−δ ∩ Lq Lq , ∀w ∈ W

(3.18)

where C is a constant, δ ∈ (0, 1), 1 ≤ q, r ≤ ∞ and ( ) 1 1 1 =δ − 1 + (1 − δ) . s r q Using Young’s inequality with m = n = 2 and ϵ = d(p − 1)/(p2 C21 ), we get )2 ( )p p ∫ ∫ ( ∫ 1 2 1 d 1 p 2d(p − 1) 1 ∂N 2 p2 C21 Λ p N dx ≤ − dx + N dx + C22 . 2 p dt 0 p ∂x 2d(p − 1) 0 µ 0

(3.19)

Using the Gagliardo–Nirenberg interpolation inequality with s = r = 2, q = 1, δ = 1/3 and the Poincar´e inequality, we obtain ⎛ )2 ⎞ 13 (∫ ) 43 p ∫ 1 ∫ 1( 1 2 p ∂N p ⎝ ⎠ 2 N dx ≤ C3 N dx dx , (3.20) ∂x 0 0 0

A. Giacobbe, G. Mulone and W. Wang / Nonlinear Analysis: Real World Applications 53 (2020) 103064

9

where C3 is a positive constant. With the help of Young’s inequality where m = 3, n = 3/2, we get from (3.20) that )2 (∫ 1 )2 p ∫ 1( ∫ 1 p ∂N 2 2 1 p N dx ≤ ϵ N 2 dx , dx + √ (3.21) C3 0 ∂x 3 3ϵ 0 0 where ϵ > 0. As a result, it follows from (3.19) that ( )p (∫ 1 )2 ∫ p Λ 4d(p − 1) 1 d 1 p 2 dx N dx ≤C22 + N 3 p dt 0 µ 0 p2 (3ϵ) 2 )∫ 1 ( 2 2 p C21 2d(p − 1) N p dx. + − 2d(p − 1) ϵC3 p2 0 Let us select

2d2 (p − 1)2 2 C3 p4 C21

ϵ= and define

1

∫

N p dx.

N (p) = sup 0≤t<τ0

0

It follows from (3.22) that ( )p ∫ ∫ 1 ( ( p ))2 2 1 d 1 p p2 C21 Λ N dx ≤C22 + C4 N − N p dx, p dt 0 µ 2 2d(p − 1) 0 where 3 4p4 C21 C4 = 2 d (p − 1)2

As a result, we obtain

(

C3 6

where C5 = 1 + 0

Set

1

N0p (x)dx

2d(p − 1)C22 + 2 p2 C21

(3.23)

) 32 .

( ( p ))2 N (p) ≤ C5 + C6 N , 2 ∫

(3.22)

(3.24)

( )p ( )3 Λ 8p2 C21 C3 2 , C6 = . µ d(p − 1) 6

1 ( )1 M (p) = max{C52p , N (p) p }.

Then it is easy to examine

1

M (p) ≤ (1 + C6 ) 2p M

(p) 2

.

Clearly, there is a positive constant C7 > 1 such that (p) (p) 1 1 M (p) ≤ (1 + C7 p) 2p M ≤ (2C7 p) p M . 2 2 Taking p = 2k where k is a positive integer, we get 1

k

1 +

1 +···+ 1 2 2k−1

M (2k ) ≤ (2C7 ) 2k 2 2k M (2k−1 ) ≤ (2C7 ) 2k

k + k−1 +···+ 1 2 2k−1

2 2k

M (1).

Consequently, there is a positive constant C1 such that M (2k ) ≤ C1 for any positive integer. Hence, we conclude ∥N (·, t)∥∞ ≤ lim M (2k ) ≤ C1 , 0 ≤ t < τ0 . k→∞

This completes the proof.

□

(3.25)

10

A. Giacobbe, G. Mulone and W. Wang / Nonlinear Analysis: Real World Applications 53 (2020) 103064

A consequence of Lemma 3.4 is that a positive solution of (2.5) with (2.6) exists on [0, ∞). Notice that (3.12) implies that there exists a positive constant C0∗ such that any nonnegative solution of (2.5) with (2.6) satisfies ∫ 1 N (x, t)dx ≤ C0∗ .

lim sup t→∞

(3.26)

0

Then by similar arguments to those in the proof of Lemma 3.4, we see that there exists a positive constant C0 such that any nonnegative solution of (2.5) with (2.6) satisfies lim sup ∥N (·, t)∥∞ ≤ C0 .

(3.27)

t→∞

Hence, we can state the main results of this section. Theorem 3.1. The nonnegative solutions of (2.5) with (2.6) exist on [0, ∞) and satisfy (3.27). 4. Basic reproduction number In this section, we formulate the basic reproduction number of the model and show that it gives the invasion threshold of the disease. Let S0 (x) be the unique positive solution of the equation: ∂2S + Λ − µS ∂x2 with zero Dirichlet boundary conditions on [0, 1]. Then direct calculations lead to ( (√ ) ( √ )) µ (2x − 1) µ Λ √ √ S0 (x) = 1 − sech cosh . µ 2 d 2 d 0=d

(4.1)

(4.2)

Clearly, E0 = (S0 , 0, 0) is the disease-free steady state. Let us introduce a few infected individuals with the distribution ϕ(x) near the disease-free state E0 . As time evolves, the distribution of the infective individuals under the synthetical influences of mortality, mobility, and recovery is described by ∂B ∂2B ∂B = d 2 + η0 + α0 T − (µ + α1 + γ)B, ∂t ∂x ∂x 2 ∂T ∂ T ∂T =d − η2 + α1 B − (µ + γ + α0 )T, ∂t⏐ ∂x⏐2 ∂x⏐ ⏐ ⏐ ⏐ ⏐ ⏐ B⏐ = B⏐ = T⏐ = T⏐ = 0. x=0

x=1

x=0

(4.3)

x=1

Let H(t) be the solution map of (4.3) in X2 = C02 ([0, 1], R2 ), which is the space of continuous functions, together to its first and second derivatives, on [0, 1] with values in R2 and satisfying the Dirichlet conditions. If ψ = (ψ1 , ψ2 ) is the distribution of infectious individuals, then the new infection rate is described by ( ) ∫1 βS0 (x) 0 g1 (a)ψ1 (a)da . F ◦ψ = 0 Note that φ = (ϕ, 0) ∈ X2 is the initial distribution of infectious individuals and H(t)φ gives the distribution of the infective individuals at time t through the mortality, recovery and mobility. It follows that the new infection rate at time t is F ◦ H(t)φ and the total new infections resulting from the initial introduction ϕ of infected individuals are ∫ ∞ L(ϕ) := F ◦ H(t)ϕdt.

(4.4)

0

Evidently, L is a positive and compact operator which maps the initial infection distribution ϕ to the distribution of the total infective members produced during the infection period. Following the idea of next generation operators (see, e.g., [19–24]), we give the following definition.

A. Giacobbe, G. Mulone and W. Wang / Nonlinear Analysis: Real World Applications 53 (2020) 103064

11

Definition 4.1. The spectral radius of L is defined as the basic reproduction number of model (2.5), which is denoted by R0 . Theorem 4.1. The disease-free equilibrium E0 is asymptotically stable if R0 < 1, and is unstable if R0 > 1. Proof . Let us define an operator B on X2 by ⎞ ⎛ ∂ψ1 ∂ 2 ψ1 + η d + α ψ − (µ + α + γ)ψ 0 0 2 1 1 ∂x ⎟ ⎜ ∂x2 Bψ := ⎝ ⎠ ∂ψ2 ∂ 2 ψ2 − η d + α ψ − (µ + γ + α )ψ 2 1 1 0 2 ∂x2 ∂x + where ψ = (ψ1 , ψ2 ), and let X2 denote the nonnegative subset of X2 . Clearly, B is the generator of H(t) on X2 , which is positive in the sense that H(t)X2+ ⊂ X2+ . Furthermore, the spectral bound s(B) of B is negative. It follows from [25] that B is resolvent-positive and ∫ ∞ −1 −B ψ = H(t)ψ dt, ∀ψ ∈ X2 . (4.5) 0

It follows from (4.4) and (4.5) that L = −F ◦ B −1 . Let us define an operator A by Aψ = Bψ + F ◦ ψ. It is evident that A is resolvent-positive. As a consequence, we conclude from [25] that the spectral bound s(A) has the same sign as R0 − 1. If R0 < 1, then s(A) < 0. As a result, it follows from [25] that the exponential growth bound of H(t) is negative. Then it is easy to see that E0 is asymptotically stable. On the other hand, E0 becomes unstable when R0 > 1 since s(A) > 0. □ 5. Stability of disease-free steady state The purpose of this section is to provide powerful testable conditions for the stability of the disease-free equilibrium. We consider here, for the sake of simplicity, the case α0 = 0, the more general case can be done by using another appropriate Lyapunov function. 5.1. Linear Lyapunov stability Substituting S = S − S0 (x), B = B − 0, T = T − 0 into (2.5), we get the nonlinear perturbation equations: ∫ 1 ∫ 1 ∂2S ∂S =d 2 − µS − βS0 (x) g1 (a)B(a, t)da − βS g1 (a)B(a, t)da, ∂t ∂x 0 0 ∫ ∫ 1 ∂2B ∂B 1 ∂B ∂B =d 2 + η0 + η1 g2 (a)B(a, t)da + βS0 (x) g1 (a)B(a, t)da ∂t ∂x ∂x 0 ∂x 0 (5.1) ∫ 1

g1 (a)B(a, t)da − (µ + ϵ + α1 )B.

+ βS 0

∂T ∂2T ∂T =d − η2 + α1 B − (µ + γ)T . ∂t ∂x2 ∂x The linearized system about E0 is given by: ∫ 1 ∂2S ∂S =d 2 − µS − βS0 (x) g1 (a)B(a, t)da, ∂t ∂x 0 ∫ 1 ∂B ∂2B ∂B =d 2 + η0 + βS0 (x) g1 (a)B(a, t)da − (µ + ϵ + α1 )B ∂t ∂x ∂x 0 ∂T ∂2T ∂T =d − η2 + α1 B − (µ + γ)T . 2 ∂t ∂x ∂x

(5.2)

A. Giacobbe, G. Mulone and W. Wang / Nonlinear Analysis: Real World Applications 53 (2020) 103064

12

To this system we add zero boundary conditions for x = 0 and x = 1 and initial conditions S(x.0) = S0 (x), B(x, 0) = B0 (x), T (x, 0) = T0 (x), where S0 (x), B0 (x), T0 (x) are sufficiently regular functions and vanish at x = 0 and x = 1. Multiplying the second equation of (5.2) by B, integrating in the interval [0, 1], and taking into account the boundary conditions, we have: [ ] ( ) ∫ 1 ∂B 2 d ∥B∥2 = −d∥ ∥ + β S0 (x) g1 (a)B(a, t) da, B − (µ + ϵ + α1 )∥B∥2 , dt 2 ∂x 0 ∫1 ∫1 where (f, g) = 0 f (x)g(x) dx and ∥f ∥2 = 0 f 2 (x) dx. By using the Cauchy–Schwarz inequality, we have (∫ 1 )1/2 ∫ 1 g1 (a)B(a, t) da ≤ g12 (a) da ∥B∥ = g0 ∥B∥ , 0

where g0 =

0

)1/2 1 2 . g (a) da 1 0

(∫

Denoting by S1 = ∥S0 (x)∥ and using the Poincar´e inequality, we have ] d ∥B∥2 ≤ − dπ 2 ∥B∥2 + (βg0 S1 − µ − ϵ − α1 )∥B∥2 dt 2 ( ) βg0 S1 2 =(dπ + µ + ϵ + α1 ) − 1 ∥B∥2 dπ 2 + µ + ϵ + α1 =(dπ 2 + µ + ϵ + α1 )(R0 − 1)∥B∥2 , [

(5.3)

where

βg0 S1 . dπ 2 + µ + ϵ + α1 If R0 < 1, it follows from (5.3) that ∥B∥ → 0 according to the following inequality R0 =

∥B∥ ≤ ∥B(x, 0)∥e(dπ

2

+µ+ϵ+α1 )(R0 −1)t

.

Now we prove that ∥S∥ and ∥T ∥ decay also exponentially to zero if R0 < 1. To this end, multiplying the first equation of (5.2) by S and the third equation of (5.2) by T and integrating over [0, 1], we obtain ( [ ] ) ∫ 1 d ∥S∥2 ∂S 2 ∥ − µ∥S∥2 − β S0 (x) = − d∥ g1 (a)B(a, t) da, S dt 2 ∂x 0 (5.4) ≤ − (dπ 2 + µ)∥S∥2 + βg0 S1 ∥S∥∥B∥ [ ] 1 ϵ ≤ − (dπ 2 + µ)∥S∥2 + βg0 S1 ∥S∥2 + ∥B∥2 . 2 2ϵ Choosing ϵ =

dπ 2 + µ , we have βg0 S1 [ ] d ∥S∥2 dπ 2 + µ ≤− ∥S∥2 + dt 2 2 dπ 2 + µ ≤− ∥S∥2 + 2

(βg0 S1 )2 ∥B∥2 2(dπ 2 + µ) 2 (βg0 S1 )2 ∥B(x, 0)∥2 e2(dπ +µ+ϵ+α1 )(R0 −1)t . 2 2(dπ + µ)

Moreover, we easily have [ ] d ∥T ∥2 dπ 2 + µ + γ α12 ≤− ∥T ∥2 + ∥B∥2 2 dt 2 2 2(dπ + µ + γ) 2 dπ 2 + µ + γ α12 ≤− ∥T ∥2 + ∥B(x, 0)∥2 e2(dπ +µ+ϵ+α1 )(R0 −1)t . 2 2 2(dπ + µ + γ) Integrating the two last inequalities with respect to t, we have

(5.5)

(5.6)

A. Giacobbe, G. Mulone and W. Wang / Nonlinear Analysis: Real World Applications 53 (2020) 103064

13

Theorem 5.1. The disease-free equilibrium is linearly asymptotically stable in the energy norm E(t) =

∥S∥2 + ∥B∥2 + ∥T ∥2 2

if R0 < 1. Remark 5.1. We note that if R0 < 1 also ∥R∥2 → 0 asymptotically, where R is the perturbation to R. 5.2. Nonlinear stability Here we study the nonlinear stability of E0 . For this purpose, we consider the nonlinear system of the perturbations: ∫ 1 ∫ 1 ∂2S ∂S =d 2 − µS − βS0 (x) g1 (a)B(a, t)da − βS g1 (a)B(a, t)da, ∂t ∂x 0 0 ∫ ∫ 1 ∂2B ∂B ∂B ∂B 1 =d 2 + η0 + η1 g2 (a)B(a, t)da + βS0 (x) g1 (a)B(a, t)da ∂t ∂x ∂x 0 ∂x 0 (5.7) ∫ 1 + βS g1 (a)B(a, t)da − (µ + ϵ + α1 )B. 0

∂2T ∂T ∂T =d − η2 + α1 B − (µ + γ)T . 2 ∂t ∂x ∂x Recalling that S = S − S0 (x), B = B, T = T , and observing that B and T are nonnegative, multiplying (5.7)1 and (5.7)2 , (5.7)3 by S, B and T , respectively, integrating over the interval [0, 1] and taking into account of the boundary conditions, we have [ ] ( ) ∫ 1 d ∥S∥2 ∂S 2 2 = − d∥ ∥ − µ∥S∥ − β S0 (x) g1 (a)B(a, t) da, S dt 2 ∂x 0 ( ∫ 1 ) −β S g(a)B(a, t) da, S 0 [ ] ( ) ∫ 1 ∂B 2 d ∥B∥2 (5.8) = − d∥ ∥ + β S0 (x) g1 (a)B(a, t) da, B dt 2 ∂x 0 ( ∫ 1 ) +β S g1 (a)B(a, t) da, B − (µ + ϵ + α1 )∥B∥2 0 [ ] ∂T 2 d ∥T ∥2 = − d∥ ∥ + α1 (B, T ) − (µ + γ)∥T ∥2 . dt 2 ∂x Since

( ∫ −β S

1

) g(a)B(a, t) da, S

≤ 0,

0

we have

[ ] ( ) ∫ 1 ∂S 2 d ∥S∥2 2 ≤ − d∥ ∥ − µ∥S∥ − β S0 (x) g(a)B(a, t) da, S dt 2 ∂x 0 [ ] d ∥B∥2 ∂B 2 ≤ − d∥ ∥ − (µ + ϵ + α1 − βS1 g0 )∥B∥2 + βg0 ∥S∥∥B∥2 dt 2 ∂x

Multiplying (5.9)1 by a positive number λ to be chosen, and using the Poicar´e inequality, we have [ ] d ∥S∥2 λ ≤ − λ(dπ 2 + µ)∥S∥2 + λβg0 S1 ∥B∥∥S∥ dt 2 [ ] d ∥B∥2 ≤ − (dπ 2 + µ + ϵ + α1 − βS1 g0 )∥B∥2 + βg0 ∥S∥∥B∥2 dt 2

(5.9)

(5.10)

A. Giacobbe, G. Mulone and W. Wang / Nonlinear Analysis: Real World Applications 53 (2020) 103064

14

Now, we introduce the Energy norm E(t) =

λ∥S∥2 ∥B∥2 + . 2 2

By summing (5.10)1 and (5.10)2 , we have E˙ ≤ − λ(dπ 2 + µ)∥S∥2 + λβg0 S1 ∥B∥∥S∥ − (dπ 2 + µ + ϵ + α1 − βS1 g0 )∥B∥2 + βg0 ∥S∥∥B∥2 ) ( ∥B∥2 ϵ¯∥S∥2 + ≤ −λ(dπ 2 + µ)∥S∥2 + λβg0 S1 2 2¯ ϵ 2 2 − (dπ + µ + ϵ + α1 − βS1 g0 )∥B∥ + βg0 ∥S∥∥B∥2 .

(5.11)

We now suppose that R0 =

βS1 g0 <1 dπ 2 + µ + ϵ + α1

and choose ϵ¯ =

dπ 2 + µ , βg0 S1

λ=

(dπ 2 + µ + ϵ + α1 )(dπ 2 + µ)(1 − R0 ) , (βg0 S1 )2

(5.12)

we have 1 λ(dπ 2 + µ) E˙ ≤ − (dπ 2 + µ + ϵ + α1 )(1 − R0 )∥B∥2 − ∥S∥2 + βg0 ∥S∥∥B∥2 . 2 2 The last inequality may be written in this way √ 1 2βg0 λ(dπ 2 + µ) 2 2 2 ˙ E ≤ − (dπ + µ + ϵ + α1 )(1 − R0 )∥B∥ − ∥S∥ + 2 √ E 3/2 2 2 λ √ ( ) 2 βg0 E 3/2 ≤ max −(dπ 2 + µ + ϵ + α1 )(1 − R0 ), −(dπ 2 + µ) E + 2 λ = AE + CE 3/2 ,

(5.13)

where A and C are the coefficients of E and E 3/2 in (5.13), respectively. By assuming R0 < 1 and E(0) < (−

A 2 ) C

we have ˙ E(0) < 0.

(5.14)

By a recursive argument (see, e.g., [26]), we have E˙ ≤ (C

√

E(0) + A)E(t), ∀t

(5.15)

A. Giacobbe, G. Mulone and W. Wang / Nonlinear Analysis: Real World Applications 53 (2020) 103064

15

Integrating the last inequality we obtain the main theorem: Theorem 5.2. Assuming R0 = and

dπ 2

βS1 g0 <1 + µ + ϵ + α1

)]2 [ ( max h(1 − R0 ), dπ 2 + µ h(dπ 2 + µ)(1 − R0 ) E(0) < , 8β 2 g02

where h = (dπ 2 + µ + ϵ + α1 ), the disease-free equilibrium is asymptotically exponentially stable in the energy norm λ∥S∥2 + ∥B∥2 E(t) = , 2 where λ is given by (5.12), according to the inequality √ E(t) ≤ E(0) exp[(C E(0) + A)t]. As for the linear Lyapunov stability, also here it can be proved that the perturbations to the individuals ¯ 0 < 1. in treatment and to those cured go to zero exponentially fast if R 6. Numerical simulations In this section, we use numerical simulations to reveal how self-forces and peer pressures to be thinner affect the emergence and distributions of eating disorders. First, the distribution of susceptible individuals in disease-free steady state is shown in Fig. 1. An endemic solution, if it exists, is the steady states solution in which the disease exists. The determination of such a solution is, in this case, non-trivial. In fact, one has to determine the stationary solution of the non-local equations (2.5) ⎧ ∂2S ⎪ ⎪ 0 = d 2 + Λ − µS − β I S, ⎪ ⎪ ⎪ ∂x ⎪ ⎪ 2 ⎪ ∂ B ∂B ⎪ ⎪ ⎨0 = d + (η0 + η1 J ) + β I S + α0 T − (µ + ϵ + α1 )B, 2 ∂x ∂x (6.1) 2 ⎪ ∂ T ∂T ⎪ ⎪0 = d − η + α B − (µ + γ + α )T, 2 1 0 ⎪ ⎪ ∂x2 ∂x ⎪ ⎪ ⎪ 2 ⎪ ∂ R ⎪ ⎩0 = d + γT − µR ∂x2 ∫1 ∫1 where I = 0 g1 (x)B(x)dx and J = 0 g2 (x)B(x)dx. These equations can be solved numerically by substituting I, J with (positive) real numbers i, j. This gives solutions Si,j , Bi,j , Ti,j , Ri,j . These functions are steady state solutions if and only if the two equations ∫ 1 ∫ 1 i= g1 (x)Bi,j (x)dx, j= g2 (x)Bi,j (x)dx, 0

0

hold. In our simulations we choose g1 = g2 = g so that I = J , and we will indicate with Si , Bi , Ti the corresponding solutions. We first choose g = 1. Fixing all parameters except β we can define a function ∫1 fβ (i) that, to every i, associates the number 0 Bi dx. The graph of fβ for increasing values of β is shown in Fig. 2. These graphs imply that there is a critical value of β, called βc above which an endemic solution exists. Below such a threshold there is only the disease-free steady state. With the same choice of parameters, and when β = 4.3 > βc , the shape of body mass index distribution of the different classes for the endemic solution is given in Fig. 3. We also perform the same computations as

16

A. Giacobbe, G. Mulone and W. Wang / Nonlinear Analysis: Real World Applications 53 (2020) 103064

Fig. 1. The BMI distribution of susceptibles for the disease-free steady state SDF , where d = 0.25, Λ = 4.1, µ = 0.4.

Fig. 2. The dashed line is the diagonal of I,III quadrant, the solid lines are the graph of fβ for choices of β = 1.5 (left) β = 2.8 (centre) and β = 4.5 (right). These values are below, on, and above the critical value βc . The other parameters are chosen to be d = 0.25, Λ = 4.1, µ = 0.4, β = 4, η0 = 0.12, η1 = 0.11, η2 = 0.08, α0 = 0.3, α1 = 0.2, ε = 0.023, γ = 1.0.

Fig. 3. The plot in the right pane represents the BMI distribution of susceptible (solid) and bulimic (dashed) for the endemic distribution. The right pane represents the treated (solid) and the removed (dashed). The values of the parameters are the same as above, the values of β is 4.3 (above threshold).

A. Giacobbe, G. Mulone and W. Wang / Nonlinear Analysis: Real World Applications 53 (2020) 103064

17

Fig. 4. In this plot the BMI distribution of bulimic are plotted for two different choices of the function g1 = g2 . In the small pane is the graph of the two choices. The function gl is plotted as a solid line, gr is plotted dashed. The corresponding BMI distributions

(

)2

for B have the same convention. The explicit form of the function g1 = g2 is exp[−1/ ((x − m)/σ)2 − 1

] prolonged to zero outside

the interval x ∈ [−σ + m, σ + m]. We have chosen ml = 1/5, mr = 4/5 and σ = 1/5.

Fig. 5. The three panes represent susceptible and bulimic for increasing values of η1 = 0.5 (left) η1 = 1.5 (centre) η1 = 2.5 (right). All other parameters are kept fixed d = 0.1, Λ = 4.1, µ = 0.4, β = 2, η0 = 0.12, η2 = 0.08, α0 = 0.3, α1 = 0.2, γ = 1, ϵ = 0.023.

above with two different choices of g, that we denote gl and gr . The function gl is a function with compact support around 1/5, plotted as a solid curve, the function gr is the same curve but translated so to have 4/5 as mean (see Fig. 4). Observe that a population that gives more importance to a low body mass index has much higher number of bulimic individuals. We finally perform the computations with fixed β above threshold βc and with increasing η1 (Fig. 5). Observe that high values of η1 correspond to lower impact of bulimia, and that the class of bulimic is more and more shifted towards low levels of BMI. 7. Conclusions In this article, for the first time, we investigate a behavioural-epidemiological model, related in particular to eating behaviours, in which a continuum parameter, the body mass index, plays a role. The existence, uniqueness and the stability of the disease-free equilibrium is analysed, and the emergence of an endemic equilibrium is proven numerically. What is more important for applications, some characteristic features of the endemic equilibrium have been exposed, and their dependence on some environmental controls have been discussed.

18

A. Giacobbe, G. Mulone and W. Wang / Nonlinear Analysis: Real World Applications 53 (2020) 103064

Despite the fact that this model needs some reliable accommodation dictated by the analysis of datasets from the medical community (for example on the form of the functions g1 , g2 , or on the treatment and relapse parameters αi ), the results could indicate some strategies for treating the disease, and provide valuable information that could allow an inverse-problem analysis to model the disease from, for example, the real BMI distribution of healthy and affected individuals. References [1] F.L. Denmark, K. Schaffer, E.M. Baron, et al., Women in the United States, in: C.M. Brown, et al. (Eds.), Women’s Evolving Lives: Global and Psychosocial Perspectives, Springer, Cham, 2017, pp. 257–276. [2] A. Keski-Rahkonen, L. Mustelin, Epidemiology of eating disorders in Europe: Prevalence, incidence, comorbidity, course, consequences, and risk factors, Curr. Opin. Psychiatry 29 (2016) 340–345. [3] T. Jackson, H. Chen, Sociocultural experiences of bulimic and non-bulimic adolescents in a school-based Chinese sample, J. Abnorm. Child Psychol. 38 (2010) 69–76. [4] P. Rohde, E. Stice, C.N. Marti, Development and predictive effects of eating disorder risk factors during adolescence: Implications for prevention efforts, Int. J. Eat. Disord. 48 (2015) 187–198. [5] C. Ciarci` a, P. Falsaperla, A. Giacobbe, G. Mulone, A mathematical model of anorexia and bulimia, Math. Methods Appl. Sci. 38 (2015) 2937–2952. [6] B. Gonzalez, E. Huerta-Sanchez, A. Ortiz-Nieves, et al., Am I too fat? Bulimia as an epidemic, J. Math. Psych. 47 (2003) 515–526. [7] Y. Fan, Y. Li, A. Liu, et al., Associations between body mass index, weight control concerns and behaviors, and eating disorder symptoms among non-clinical Chinese adolescents, BMC Public Health 10 (2010) 314. [8] M.L. Hill, A. Masuda, R.D. Latzman, Body image flexibility as a protective factor against disordered eating behavior for women with lower body mass index, Eat. Behav. 14 (2013) 336–341. [9] L.M. Thornton, J.E. Dellava, T.L. Root, et al., Anorexia nervosa and generalized anxiety disorder: Further explorations of the relation between anxiety and body mass index, J. Anxiety Disord. 25 (2011) 727–730. [10] M.J. Tov´ ee, J.L. Emery, E.M. Cohen-Tov´ ee, The estimation of body mass index and physical attractiveness is dependent on the observer’s own body mass index, Proc. R. Soc. Lond. B 267 (2000) 1987–1997. [11] P.K. Keel, J.E. Mitchell, Outcome in bulimia nervosa, Am. J. Psychiatry 154 (1997) 313–321. [12] C. Cosner, Y. Lou, Does movement toward better environments always benefit a population? J. Math. Anal. Appl. 277 (2003) 489–503. [13] A. Friedman, G. Lolas, Analysis of a mathematical model of tumor lymphangiogenesis, Math. Models Methods Appl. Sci. 1 (2005) 95–107. [14] Y. Tao, Global existence of classical solutions to a predator–prey model with nonlinear prey-taxis, Nonlinear Anal. RWA 11 (2010) 2056–2064. [15] M.H. Protter, H.F. Weinberger, Maximum Principles in Differential Equations, second ed., Springer-Verlag, Berlin, 1984. [16] O.A. Ladyzenskaja, V.A. Solonnikov, N.N. Uralceva, Linear and Quasi-Linear Equations of Parabolic Type, in: Amer. Math. Soc. Transl., vol. 23, Providence, RI, 1968. [17] L. Nirenberg, An extended interpolation inequality, Ann. Sc. Norm. Super. Pisa 20 (1966) 733–737. [18] A. Friedman, Partial Differential Equations, Holt, Rinehart and Winston, New York, Montreal, Que., London, 1969. [19] O. Diekmann, J.A.P. Heesterbeek, J.A.J. Metz, On the definition and the computation of the basic reproduction ratio R0 in the models for infectious disease in heterogeneous populations, J. Math. Biol. 28 (1990) 365–382. [20] P. van den Driessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci. 180 (2002) 29–48. [21] W. Wang, X.Q. Zhao, Threshold dynamics for compartmental epidemic models in periodic environments, J. Dynam. Differential Equations 20 (2008) 699–717. [22] W. Wang, X.-Q. Zhao, A nonlocal and time-delayed reaction–diffusion model of Dengue transmission, SIAM J. Appl. Math. 71 (2011) 147–168. [23] W. Wang, X.-Q. Zhao, Basic reproduction numbers for reaction–diffusion epidemic models, SIAM J. Appl. Dyn. Syst. 11 (2012) 1652–1673. [24] W. Wang, X.-Q. Zhao, Spatial invasion threshold of lyme disease, SIAM J. Appl. Math. 75 (2015) 1142–1170. [25] H.R. Thieme, Spectral bound and reproduction number for infinite-dimensional population structure and time heterogeneity, SIAM J. Appl. Math. 70 (2009) 188–211. [26] G. Mulone, S. Rionero, Necessary and sufficient conditions in the magnetic B´ enard problem, Arch. Ration. Mech. Anal. 166 (3) (2003) 197–218.