Abstract
We study an epidemiological mathematical model formulated in terms of an ODE system taking into account both social and nonsocial contagion risks of obesity. Analyzing first the case in which the model presents only the effect due to social contagion and using qualitative methods of the stability analysis, we prove that such system has at the most three equilibrium points, one disease-free equilibrium and two endemic equilibria, and also that it has no periodic orbits. Particularly, we found that when considering (the basic reproductive number) as a parameter, the system exhibits a backward bifurcation: the disease-free equilibrium is stable when and unstable when , whereas the two endemic equilibria appear from (a specific positive value reached by and less than unity), one being asymptotically stable and the other unstable, but for values, only the former remains inside the feasible region. On the other hand, considering social and nonsocial contagion and following the same methodology, we found that the dynamic of the model is simpler than that described above: it has a unique endemic equilibrium point that is globally asymptotically stable.
1. Introduction
Obesity has gone from being an isolated health problem, related to some people, to a global problem. Considered as “The Pandemia of the 21st Century” [1], it is present in both developed and underdeveloped countries, eclipsing in the latter the problem of malnutrition to become today one of its main priorities [2]. Around the world countries spend a huge amount of their year budgets as well as qualified human resources to fight this disease, which is frequently associated with serious health pathologies such as diabetes mellitus, high blood pressure, and lung and heart diseases and, moreover, it is the cause of several kinds of cancer [3, 4]. Also it affects the psychological condition of the individuals, because it can damage their self-esteem and social relationships [1, 2, 5, 6].
Clearly, obesity is a huge and difficult problem, with the aggravating factor that it is present in all sectors of society regardless of the income, ethnicity, age, gender, or another socioeconomic status of their individuals. Besides, it is often associated with the wrong diet, sedentary lifestyle, or genetic predisposition of individuals, and if it was not enough, it has also been found that obesity can be produced by a large variety of causes that are linked to cultural, social, and economic conditions of the environment in which people develop. Today, the way in which the latter influence the origin of obesity is far from being understood. In summary, the problem of pandemic obesity is complex and multifactorial and it has increased in recent years over all the world. Thus, obesity has become a relevant current research topic in which different fields of human knowledge converge for the purpose of understanding its causes, knowing its consequences, and, as far as possible, keeping it under control or eradicating it.
In this way, mathematical modelling is a means to provide a general insight for the dynamics of obesity and, as such, could hopefully become a useful device to develop control strategies. With regard to its causes of social origin, the dynamics of obesity can be well modelled by epidemic-type models as a process of social contagion, as was evidenced by Christakis and Fowler who studied the spread of obesity in a large social network over 32 years and established that obesity can spread through social ties [7]. This approach has resulted in a wide range of papers of mathematical modelling in which obesity is studied as a social epidemic [8–15]. Social obesity epidemic models typically divide the population into two or several classes or subpopulations. In [12] the classical model is extended, where infection occurs by nonsocial mechanisms as well as through social transmission. There are models in which it has been considered a bilinear incidence rate [14] (for subpopulations of normal weight, overweight, and obese individuals), obtaining as a result a unique stable equilibrium point; in [9–11, 13] this effect was considered, but for six subpopulations: normal weight, latent, overweight, obese, becoming overweight, on diet, and obese on diet individuals. Other models have incorporated the effects of the time delay [15] and have also formulated nonautonomous obesity epidemic models [9], in which periodic positive solutions were found under some sufficient conditions using a continuation theorem based on coincidence degree theory.
In this paper we analyze the model proposed by Ejima et al. [8]: a variant of the SIRI model in which the individuals who recover temporarily may get recurrence to infectious state and is formulated on the premise that obesity is caused by both social and nonsocial contagion routes [16]. The objective of this work is to analyze this system using the methods of the qualitative theory of ordinary differential equations.
The rest of this paper is organized as follows: in Section 2, we present the Ejima et al. model [8] and we reduce it in a two-dimensional system. Section 3 focuses on the case in which only the risk of social contagion of obesity is considered: we perform a local analysis in order to establish the equilibrium points and their corresponding local stabilities as well as bifurcations; we also obtain a global analysis by means of an appropriate Lyapunov function to establish the global stability of the endemic equilibrium point and, by using Dulac’s criterion, the nonexistence of periodic orbits. In Section 4 we will study the case with risk of both social and nonsocial contagion of obesity: we prove the existence of a unique endemic equilibrium point that is globally asymptotically stable, and also through a suitable Dulac function, we determine the nonexistence of periodic orbits. Section 5 investigates several important aspects of our model from a numerical point of view. Finally, in Section 6 we collect some observations and conclusions.
2. Obesity Mathematical Model
The model proposed by Ejima et al. [8] for the dynamics of obesity is given by the following set of three differential equations:wherein , , and , respectively, denote the susceptible (never-obese), infectious (obese), and recovered (ex-obese) individuals in a population. In (1) the natural death and birth rates are assumed to be equal and denoted by ; thus, we have for all time (the population size is constant). Also, the parameter is the transmission rate due to social contagion risk of obesity, describes the rate at which the infectious individuals become recovered individuals, is the hazard of obesity due to nonsocial contagion reasons, and is the relative risk of weight regain among ex-obese individuals which typically takes a value greater than unity () due to high risk of coming back to the obese state. All the involved parameters are positive.
In Figure 1 the flow diagram is shown in which the model assumptions are based and system (1) is deduced. There may be seen how it is a variant of a SIRI model, since the term , given in the second equation (1), represents the effect of relapse; that is, the recovered people (ex-obese), after some time, become infected (obese) again. It should be noted that, in total absence of the relapse term, (1) is reduced to the corresponding well-known SIR model.
Adding the three equations (1) we obtain the differential equation for the total population of this systemIt may be shown easily that the region of biological sense is given bythat is, all solutions starting in remain there for all . Clearly, the set is positively invariant with respect to (1).
2.1. Reduction to a Two-Dimensional IR System
We can eliminate from the equations of system (1), using the identity , in order to obtain the following two-dimensional system written in terms of the variables :It is straightforward to prove that the regionis also a positively invariant set for the reduced model (4).
With the aim to carry out a qualitative analysis of the dynamics for system (4), we study it considering firstly that transmission only is produced by social contagion route () and, afterwards, taking into account both social and nonsocial contagion hazards of the obesity (). In each case we performed a local and global analysis of the behaviour of solutions of the system.
3. Case with Only Social Contagion Risk of Obesity
3.1. Equilibria, Local Stability, and Backward Bifurcation
If we impose the condition in system (4) it takes the following form:By equalizing to zero the right members of system (6), as is commonly done to find equilibrium points, we obtain the cubic polynomial equationwhere the coefficients have been definedIn (8) is defined aswhich is named as the basic reproductive number of the system and, as is well known, is a dimensionless quantity that represents the average number of secondary infections caused by an infective individual introduced into a group of susceptible individuals. From the second equation (6) is obtained in terms of by means of
According to (7) and (10), we have several points of equilibrium for model (6): one disease-free equilibrium, , and at least two endemic equilibria, . Besides, the solutions of (7), given by andwherein is a discriminant, can be graphed as functions of in the first quadrant of the plane where they represent different families of equilibrium points: one formed by a line of disease-free equilibrium along the horizontal axis and the other two of endemic equilibrium constituted by quadratic branches. In this situation plays the role of a bifurcation parameter. See Figure 2.
In order to determine the stability of these families of equilibrium points, we need to calculate the Jacobian matrix of system (6) and evaluate it on them. This matrix overall is given asNow we perform the local stability analysis considering separately the disease-free equilibrium and the endemic equilibria cases.
(1) The Disease-Free Equilibrium, . The Jacobian matrix (12) evaluated in the disease-free equilibrium takes the simple form whose eigenvalues are and . Note that if , then and therefore is locally asymptotically stable; while when , and is unstable. In both cases the equilibrium point is hyperbolic. This situation is depicted in Figure 2 by the horizontal line on the axis, wherein for we have stability (indicated by a solid line) and for instability (showed by a dashed line). It should be noted that if , then ; as a consequence, the equilibrium point is nonhyperbolic and could be the value in which a bifurcation is produced. Hence, we have the following result.
Theorem 1. If , then is an equilibrium of system (6) and it is locally asymptotically stable. Otherwise, if , then is unstable.
(2) The Endemic Equilibria, . By substituting the endemic equilibria in system (6) the identities are obtainedwhich allows writing the Jacobian matrix (12), evaluated in such points, asThe eigenvalues of (15) are given bywhereare the trace and determinant of (15), respectively. Because in the region both coordinates of are positive, then always is negative, not so with , which could be positive or negative. According to the discriminant of (11) in the corresponding analysis of the equilibrium points, in order to ensure that , the following three cases may be identified:(i)If , then considering in this condition expressions , , and given by (8), we obtain , whereConsequently, in accordance with (11), in this interval we have that there is not endemic equilibrium points. Since and are two positive real numbers, the relation between its arithmetic and geometric means is always given bytherefore, from (19) and (20), it may be concluded that .(ii)If , then substituting , , and given by (8) in this equality, we obtain . Thus, according to (11), when reaches this value we have only one endemic equilibrium point and , whereis the corresponding infected population of Moreover, if (19) and (21) are substituted in determinant (18) and trace (17), these are simplified asrespectively. Consequently, taking into account (22) and (23), eigenvalues (16) are reduced to and , one of which is negative and the other is zero. Therefore, this equilibrium point is nonhyperbolic, and could also be the value in which a bifurcation is produced. See Figure 2.(iii)If , then similarly using , , and given by (8), this inequality leads to . As a consequence, from (11) it follows that in this interval we have families of endemic equilibrium points determined by two quadratic branches. It may be shown that, for all endemic equilibrium points that form the upper curve, since , from (18) it results in the fact that ; thereby, in (16) it follows that both eigenvalues and are negative, and therefore, these equilibrium points are asymptotically stable. On the other hand, it can also be shown that for endemic equilibrium points that are part of the lower curve, since , from (18) it occurs that ; in this way, again from (16), it is obtained that and , and consequently such equilibrium points are unstable. Note that in the latter case (in order that ), quadratic branch is delimited by the point in which, by the way, and therefore and ; that is, the said point is nonhyperbolic and in it a bifurcation could also occur. These results are illustrated in Figure 2.It is necessary to point out that from the above arguments it follows that if , then the disease-free family of equilibrium points is the only one in the feasible region . Based on the previous analysis, we have the following theorem which summarizes the local stability of all endemic equilibria points of system (6) contained in .
Theorem 2. If , then there are not endemic equilibrium points in the feasible region . Also, If , then in such region there is a unique , which is not hyperbolic. Moreover, if , then there are families of endemic equilibrium points determined by two quadratic branches: an upper, in which all points are asymptotically stable, and the other lower, formed by unstable points and with epidemiological significance only when .
In Figure 2 the results indicated by Theorems 1 and 2 are displayed schematically. There it may be appreciated that the trajectories of solutions of system (4) change abruptly in the nonhyperbolic points and ; that is, in these points really bifurcations occur. In that sense, Figure 2 can be considered as the bifurcation diagram in the feasible region of said system, and the type of bifurcation shown in it corresponds to a backward bifurcation.
3.2. Nonexistence of Periodic Orbits and Global Stability
Global analysis of system (6) includes the study of periodic orbits. The following result shows that this system does not have periodic orbits.
Theorem 3. System (6) does not have periodic orbits in the interior of .
Proof. Consider system (6) for and . We consider the Dulac function Given that that is, in the interior of , it follows from the Dulac criterion that system (6) has no periodic orbits in that region.
On the other hand, global analysis of system (6) also involves the study of global stability of its disease-free equilibrium and its endemic equilibria. The former may be proved by arguments of local stability, as can be seen in the following theorem.
Theorem 4. If , then the unique disease-free equilibrium of system (6) is globally asymptotically stable in .
Proof. According to the discussion given previously during the formulation of Theorems 1 and 2, it was found that if , then the disease-free equilibrium point is the only one, and it also is asymptotically stable in the feasible region . Besides, according to Theorem 3, there are not periodic orbits here. As a result, all initial condition in , which satisfies this interval, originates paths tending to ; consequently, such point must be globally asymptotically stable.
Finally, in order to prove the global stability of a unique endemic equilibrium in the interior of when , we use the method of Lyapunov functions. This is accomplished in the next theorem.
Theorem 5. Assume that . Then a unique endemic equilibrium of system (6) is globally asymptotically stable in the interior of .
Proof. To prove the global asymptotic stability of the unique endemic equilibrium , we define a Lyapunov function given byLet . By using , we haveLet . By using , we haveThe derivative of (26) along solution of (6) is given byBy means of (27) and (28), we obtainUsing we haveTherefore, if , then is negative definite. By the Lyapunov asymptotic stability theorem [17], this result implies that is globally asymptotically stable in the interior of .
Remark 6. Recently, Vargas-De-León [18] used the Lyapunov function to prove global stability of the coexistence equilibrium of two-species mutualism models.
4. Case with Social and Nonsocial Contagion Risks of Obesity ()
4.1. Equilibria and Nonexistence of Periodic Orbits
We might follow the usual method of setting the right-hand side of (4) equal to zero in order to obtain their equilibria. This procedure is not the most appropriate, because it leads to a cubic equation whose analytical solutions are quite complicated. Instead, an alternative way to find the equilibrium points is to perform a geometric analysis based on the intersection of the nullclines of system (4). The nullclines (vertical directions) and (horizontal directions) are given by the functionsrespectively.
Note that and , respectively, are discontinuous in and . Besides, it is clear that (resp., ) as and that the function (resp., ) is increasing and concave down. Given that and , then the functions and intersect at a single point (the equilibrium point) in the first quadrant (see Figure 3). This result indicates that in an endemic equilibrium point exists and is unique.
Regarding whether or not system (4) has periodic orbits, we have precisely the following result.
Theorem 7. System (4) does not have periodic orbits in the interior of .
Proof. Consider system (4) for and . From the Dulac function it follows thatSince in , taking into account the Dulac criterion, we conclude that system (4) has no periodic orbits.
4.2. Local and Global Stability
In order to determine the local and global stability of the single equilibrium point in , we consider that it is inside of the triangular region bounded by the dashed line shown in Figure 4. This is a trapping region: the vector field of system (4) on the boundary point into the box. For the points on the horizontal side, and , where and are the positive and negative roots of , respectively. Thus, these inequalities imply that the vector field points up and right on the interval , while it goes up and left on . Moreover, in the vertical side, and . This implies that the vector field is directed downward and right; except at the origin, where it goes to the right.
Finally, on the diagonal side of slope extending from the point to , it can be shown thatwhere and are the positive and negative roots of , respectively. Hence, the vector field, for large values of , is almost parallel to the diagonal line. In a more precise analysis, always and if or if . According to (37), the first case implies that the vector field points inward (to the right) on the diagonal line, because is more negative than , while the latter indicates that the vector field is directed inside (to the left) the diagonal side, because is less negative than . (See Figure 4.) Therefore, we conclude that the triangular region is effectively a trapping region.
On the other hand, Theorem 7 prohibits the existence of closed orbits inside the triangular region. Consequently, the trajectories entering the triangular region, according to the Poincaré-Bendixson theorem [19], must converge to the single equilibrium point located inside it. Therefore, the equilibrium point inside the triangular region must be asymptotically stable. Moreover, since this point is the only point of equilibrium in , it is also globally asymptotically stable.
The results of this section may be summarized in the following theorem.
Theorem 8. System (4) has a unique equilibrium point in the triangular region , which is globally asymptotically stable.
5. Numerical Results
In this section we perform a series of numerical simulations and graphs of system (4) to illustrate the different results obtained for each of the two cases of interest and previously analyzed. We will use the parameter values reported in the work developed by Ejima et al. [8]; their epidemiological meanings and magnitudes are indicated in Table 1.
5.1. Case
According to the magnitudes of the parameters reported in Table 1, (19) takes the value , whereas (21) acquires the value . Furthermore, taking into consideration (10) and (11), the resulting two endemic equilibria are which, as discussed previously, are asymptotically stable and unstable, respectively.
In Figures 5–9 we illustrated representative phase portraits of the dynamic system (6), for some different values of between 0 and just over 1, according to the bifurcation diagram shown in Figure 2. These values are obtained by maintaining fixed in (9) , , and , but varying the parameter . Thus, if takes values between 0 and , say , we have only one equilibrium point (the disease-free equilibrium) which is globally asymptotically stable (see Figure 5). If , we have two equilibrium points: one asymptotically stable disease-free equilibrium and the other endemic equilibria point which is nonhyperbolic (see Figure 6). When takes values between and 1, say , we have three equilibria: one disease-free equilibrium (asymptotically stable) and two endemic equilibria (asymptotically stable and unstable), as shown in Figure 7. If , there are two equilibrium points: one disease-free equilibrium which is nonhyperbolic and the other asymptotically stable endemic equilibrium, as shown in Figure 8. Finally, if takes a slightly greater value to 1, for example, , we have again two equilibrium points: one disease-free equilibrium (unstable) and the other endemic equilibrium (asymptotically stable), but both are hyperbolic (see Figure 9).
On the other hand, in Figures 10–14 examples of the evolution in time of and are presented, for the values , , , , and , which illustrate, respectively, some of the dynamics shown in the phase portraits in Figures 5–9.
5.2. Case
Taking into account the value given in Table 1, in addition to the others there reported, nullclines (33) and (34) intersect at the only endemic equilibrium pointwhich is globally asymptotically stable. In Figure 15 the corresponding phase portrait with several representative trajectories around this equilibrium point is shown. Also, performing another series of numerical simulations of system (4), we show for this case in Figure 16 an example of the evolution in time of and .
6. Discussion and Conclusions
According to the results of our analysis of the Ejima et al. model, we arrive at the following conclusions. It has been shown that the Ejima et al. model [8], if only the social contagion () is considered, presents the two bifurcation values and which results in what is known as backward bifurcation. From a mathematical point of view, the backward bifurcation is apparently constituted by the combination of two different kinds of bifurcations: a saddle-node bifurcation that occurs at , whose upper (asymptotically stable) and lower (unstable) quadratic branches are defined for ; and a transcritical bifurcation that takes place in , whose horizontal branches on the axis are asymptotically stable for and unstable for . For values of , only the upper quadratic branch asymptotically stable of the saddle-node bifurcation and the unstable horizontal branch of the transcritical bifurcation have a meaning in the feasible region . The diagram shown in Figure 2 suggests the presence of these two bifurcations in the vicinities of the points and . Of course, we should go beyond the graphic evidence and employ the qualitative analysis of local bifurcations in those points to verify our statement; even without the realization of this depth analysis, our discussions and conclusions can be supported firmly from the results shown in this work of investigation.
Indeed, the presence of the different branches of families of equilibrium points, shown in Figure 2, causes the system to display, as the parameter is varied, a bistable region and, as a consequence, the hysteresis phenomenon. In order to show in more detail the presence of hysteresis in the bifurcation diagram mentioned before, suppose we start the system in the state (see Figure 17) and then slowly increase (which is indicated by an arrow to the right just beneath the axis of Figure 17). We remain at the origin until , when the origin loses its stability. At this point, at the slightest disturbance, the system will jump to the state located in the stable upper branch. Insofar as increases its value (, the equilibrium points of the system will move along this branch away from the state . If now continuously reduced its value, equilibrium states will return to the state , cross it, and reach the state (in which ). At this point, at the slightest change again, equilibrium states will jump back to the origin (in ) and if continues to decrease we will return to the state . Consequently, the system exhibits a lack of reversibility or it is said that it has memory (or that it presents hysteresis), because the forward and backward paths are not identical.
From an epidemiological point of view, in this backward bifurcation, if , no matter where the initial conditions are taken, solutions always tend to the straight line of stable disease-free points (that correspond to a never-obese population). When , if the initial conditions are located below the unstable quadratic branch of endemic equilibrium points (obese people), then solutions continue going to the horizontal line of stable disease-free points; otherwise, they will be directed to the stable curve of endemic points located at the upper quadratic branch. Finally, if , no matter again where the initial conditions are situated, solutions will continue heading for the stable endemic points of such upper branch. Should be noted that, unlike what happens in systems that exhibit the phenomenon of forward bifurcation, in which the endemic equilibrium exists only for [20], in systems exhibiting a backward bifurcation, under certain initial conditions it is possible to have endemic states when . In our case, the endemic equilibrium that exists for just above one has a large obese population, so the result of rising above one would be a sudden and dramatic jump in the number of infective. Moreover, reducing back below one would not eradicate the obesity if the size of the infected population is greater or equal than , the corresponding value of . In order to eradicate the disease, must be further reduced, so it just needs to be slightly less than , to suddenly enter the region where only stable disease-free equilibrium (never-obese population) exists. Both such different abrupt behaviours are the result of hysteresis phenomenon.
It can be seen from the bifurcation diagram shown in Figure 2 that the way in which obesity infected population can be lowered (or raised) is through the reduction (or increase) of parameter, which, in accordance with (9), depends on the four quantities , , , and . Since we have assumed in (1) that the natural death and birth rates are equal and denoted by , the population size is constant. Therefore, can only change according to the following forms. Firstly, diminishes if the natural rate (of death or birth) increases, the transmission rate decreases, or the rate grows at which the infected individuals become recovered. On the other hand, increment its value if the opposite occurs: decreases, goes up, or declines. It should be noted that changes in could be slower and achieved in the long term, probably in decades, which is not a viable option. Conversely, changes in or might be faster and perhaps occur in fewer years. Since the focus is on reducing obesity levels and therefore the values of , it is of interest to seek mechanisms to decrease and/or increase . The former could be achieved by sensitizing at the population in general about the benefits associated with having a healthy and balanced diet and also by designing effective health and informative public prevention policies. The latter could be done by providing better medical assistance and making it available to the obese population. It should be emphasized that this last strategy is not feasible in countries with high levels of obesity as Mexico, because at the present time it is very expensive and will become even more expensive as the affected population increases. In this regard, we believe that it is much more convenient that the governments of the affected countries implement broader and more effective preventive policies.
It is worth mentioning that in Mexico there have been studies to quantify more precisely the prevalence of people suffering from overweight and obesity. It is estimated that 71.3% of Mexican adults aged 20 or more are found in this condition [21]. While trends show slowing down of the increase in the obesity prevalence, there is no evidence to infer that prevalence will decrease in the coming years. In this sense, studies like this allow us to analyze the trends of the population with obesity and study the different possible scenarios. Particularly, with the model studied this could be done for Mexico if we knew accurately the values of its parameters and . To our knowledge, there is not such information. Nevertheless, these parameters could be estimated in a short or medium term, based on statistical information collected on obesity in our country. The knowledge of how these parameters behave over time would contribute in quantifying the effectiveness of the Mexican public policies against obesity.
With respect to the model studied, it is pertinent to comment about the term of relapse of the second equation (1). In the first place, in this work it has been considered that the parameter , the relative risk of weight regain among ex-obese individuals, is greater than unity (). This condition represents the high risk that recovered (ex-obese) individuals will become obese again. However, it is also plausible to consider that , which represents that ex-obese population is more resistant to being obese again; that is, they are more aware of the health risks that this would represent. The consideration that does not change the results of the stability analysis and periodic orbits performed in Sections 3 and 4. With respect to the latter, the nullcline (given by (33)) changes to , where , while (34) remains the same. The intersection of the nullclines and again occurs at only one point in the region , and its stability is the same. Regarding this aspect, we could conclude that, regardless of the value of the parameter , there is always the possibility that the ex-obese population is at risk of relapse and becoming obese again; of course, this possibility is less if and greater if . In this sense, public health programs against obesity designed by governments may be more effective if they take into account the awareness of the ex-obese population.
Finally, with regard to the mentioned relapse term , if we consider only the effects of social contagion (which corresponds to ), it is reduced to the nonlinear term which causes system (4) to exhibit a backward bifurcation. If the effects of nonsocial contagion (that is, ) are also incorporated, the presence of the linear term causes the backward bifurcation to disappear in system (4) and in its place it only presents an equilibrium point global asymptotically stable. The latter case suggests that the way in which the effects of nonsocial contagion are introduced into the model by means of the linear part of the relapse term should be modified in order to maintain much of the interesting dynamic that is present in the backward bifurcation. Certainly, an alternative way is considering, instead of a linear, a nonlinear contribution. Perhaps it is possible to build a more general nonlinear model, where in the case in which the contribution due to nonsocial effects gradually fades, the backward bifurcation of the Ejima et al. model could be obtained as a limit case.
Competing Interests
The authors declare that they have no conflict of interests.
Acknowledgments
This research was supported by the UACM-SECITI Research Grant PI-2013-34. Cruz Vargas-De-León would like to express his indebtedness to Dr. Ramón Reyes Carreto for his hospitality during his stay as visiting professor at the UAGro in Chilpancingo.