- Email: [email protected]
Modeling the social obesity epidemic with stochastic networks
Physica A 389 (2010) 3692–3701
Contents lists available at ScienceDirect
Physica A journal homepage: www.elsevier.com/locate/physa
Modeling the social obesity epidemic with stochastic networks Gilberto González-Parra a,b,∗ , L. Acedo b , Rafael-J. Villanueva Micó b , Abraham J. Arenas b,c a
Dpto. de Cálculo, Fac. Ingeniería, Universidad de los Andes, Mérida, Venezuela
b
Instituto de Matemática Multidisciplinar, Universidad Politécnica de Valencia, Edif. 8G, 2o , 46022 Valencia, Spain
c
Dpto. de Matemáticas y Estadística, Universidad de Córdoba, Montería, Colombia
article
info
Article history: Received 27 January 2010 Received in revised form 18 April 2010 Available online 29 April 2010 Keywords: Networks dynamics Populations Mean-field Gamma distribution Obesity epidemic model
abstract In this paper we extend a compartmental model to the case of a homogenous network epidemic model for a study of the dynamics of obese populations. The social epidemic network-based approach developed here uses different algorithms and points of views regarding the simulation of the dynamics of the network. First, Monte Carlo simulations for homogeneous networks using a traditional constant probability transition rates and a mean-field-like approach are presented. We show that these networks evolve towards an approximately stationary state, which coincides with the one obtained by the underlying classical compartmental continuous model. A mean-field-like approach is applied in order to reduce the large computation time required when dealing with large contact networks. We also investigate, using homogenous contact networks, the effect of the realistic assumption that the waiting times between subpopulations follow a gamma distribution instead of the traditional exponential distribution. It is concluded that careful attention must be paid to the distributions assumed for the state periods. © 2010 Elsevier B.V. All rights reserved.
1. Introduction Classical models of disease dynamics rely on systems of differential equations that represent the number of individuals in various categories through continuous variables, allowing for infinitesimal population densities. The origin of these models is commonly traced back to the well-known pioneering work of Kermack and McKendrick. Systems of ordinary differential equations (ODEs) are well-known tools that have been used to model different types of disease [1–3]. The most discussed type of infection spread is the SIR system, in which individuals are susceptible (S), infective (I) or removed/immune (R). In addition, several differential equation models have been proposed for modeling social behavior, such as rumors, social behavior and ideologies [1,3–5]. Understanding the nature and origin process of epidemic behavior has been a subject of considerable interest to mathematical epidemiologists. In order to model these epidemics it is important to choose a convenient framework. In this context we must acknowledge that the controversy between discrete individual models and continuous population models is an old one, and now a revival is taking place. The classical continuous approach to the modeling of epidemic spread suffers from some drawbacks; for example, it neglects the local character of epidemic spread. There are several important studies that deal with epidemics using networks [6–10]. In addition, different types of network such as small world and scale-free ones have been studied in order to understand epidemics [11–15]. In [16], a survey of the most important measurements of networks is presented. These measures give an insight into the understanding of the evolution and topological features of different network models.
∗
Corresponding author at: Dpto. de Cálculo, Fac. Ingeniería, Universidad de los Andes, Mérida, Venezuela. Tel.: +34 963879144; fax: +34 963879144. E-mail addresses: [email protected], [email protected] (G. González-Parra), [email protected] (L. Acedo), [email protected] (R.-J. Villanueva Micó). 0378-4371/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.physa.2010.04.024
G. González-Parra et al. / Physica A 389 (2010) 3692–3701
3693
In [17], there are shown some advantages of modeling epidemics with networks and as well with the continuous model. It is important to remark that in some cases continuous models and discrete network results may differ. For instance, the Halloran et al. results [18] do not agree with the recent findings of [17,19]. In [20], the authors develop an interesting study of an SIR model using networks under different algorithms and conclude that simulated epidemics generated under some model architectures are insensitive to the average degree of contact amongst nodes. The efficiency of the epidemic spread related to community structure and the high clustering coefficient has been also studied [21]. Continuous population models have the advantage of being easier to analyze theoretically than discrete individual models. In addition, simple deterministic models are still at the core of theoretical epidemiology despite the increasing evidence for the importance of contact networks underlying transmission at the individual level. These compartmental mathematical models based on homogeneous mixing have made, and are still making, important contributions to the epidemiology and the ecology of infectious diseases, but they fail to reproduce many of the features observed for the spreading of diseases in contact networks [9]. In this paper, we introduce several stochastic networks describing the obesity population dynamics in order to incorporate different assumptions and compare the results with the classical mathematical model. It is important to remark that most of the previous studies developed using networks focused on the general classical SI, SIS and SIR epidemic models. Here, we use a more complex network regarding classes or states, but heterogenous networks such as random, scale-free or small world are not considered in this paper. An interesting study is presented in [22] about the progression of dengue disease in Singapore, where the authors show that the dengue epidemic organized itself into a scale-free network. Other interesting applications of scale-free networks have been applied in different areas [23–25]. The study of graphs, which originated with the abstract theory of Erdös and Rényi of the 1950s and 1960s, has become important more recently in many areas, including social contacts and computer networks, as well as the spread of communicable diseases [7]. In this paper we describe the stochastic contact network by a complete graph with members of the population represented by vertices and with contacts between individuals represented by edges. An edge is a contact between vertices that can transmit infection. The number of edges of a graph at a vertex is called the degree of the vertex. A well-known definition is the degree distribution of a graph pk , where pk is the fraction of vertices having degree k. The degree distribution is fundamental in the description of the spread of the model using heterogenous networks. However, in this study we assumed a homogenous network, where all vertices have contact (connecting edges) with other vertices. Here, we assume that all contacts between an infective and a susceptible transmit the infection with the same probability. In [26], a compartmental deterministic mathematical model to study the dynamics of overweight and obesity in a childhood population was introduced. The authors consider obesity as a social transmission epidemic concern produced by unhealthy habits. Thus, our first aim is to simulate in detail the dynamics of the overweight and obese population on a homogenous network through numerical studies, and investigate the impact of the interaction rules on the efficiency and reliability of the simulation process. Therefore, different algorithms to simulate the network are performed in order to characterize the system dynamics. Simulation of large network models requires heavy computational cost; therefore, when large populations are simulated it is necessary to consider more efficient algorithms, for instance a mean-field approach or discrete events oriented simulation. One main assumption commonly employed in the classical continuous differential equation models is that the time for which individuals remain in one state can be described by an exponential distribution. This assumption avoids more complex equations such integro-differential equations. For instance, if a fixed lifetime is assumed in the deterministic model, then integro-differential equations arise in the model. However, real-world waiting times may have different probability distribution functions. Lloyd states that this distribution is biologically unrealistic, because it corresponds to the assumption that the chance of changing the state at a given time interval is independent of the time spend on the actual state. In addition, this assumption translates in that the distribution of remaining times in a state are too dispersed. In reality, the state periods are fairly closely centered about the mean duration of the period to state change [27]. Thus, another aim of this study is to investigate the effect of introducing the assumption that the waiting times between subpopulations (states) of the model follow a gamma distribution [27,28]. Using the simulation of the networks we investigate the effect that this has on the population dynamics. It is important to remark that simple modeling assumptions, such as the aforementioned, may be a driving force on the dynamics of the model; therefore their study is well justified. The layout of this paper is as follows. In Section 2, the underlying classical compartmental mathematical model is introduced. Section 3 is devoted to a brief introduction to networks modeling. Simulations of the proposed homogenous contact networks are presented in Section 4. A discussion and conclusions are presented in Section 5. 2. Classical mathematical obesity population model In several previous works, it is considered that the obesity problem is a result of social and cultural forces, including social networks [26,29–31]. In [26], a deterministic mathematical model to investigate the dynamics of overweight and obesity in a childhood population was introduced. In this model, the childhood population was divided into six disjoint classes or subpopulations: children with normal weight N (t ), latent children (children with unhealthy habits but who are still normal weight) L(t ), overweight children S (t ), obese children O(t ), overweight children on a diet DS (t ) and obese children on a diet DO (t ). The proposed model considers the 3–5-year-old population in order to study a homogenous group that is commonly studied. The model was normalized for the sake of clarity and since the qualitative behavior of the different classes is more
3694
G. González-Parra et al. / Physica A 389 (2010) 3692–3701
Table 1 Parameters of model (1), for the 3–5-year-old child population in the Spanish region of Valencia. Parameter
Value (weeks)
Parameter
Value (weeks)
µ γS α σ γD
0.0064 0.003085 4.068 × 10−3 4.4379 × 10−3 4.5045 × 10−4
γL ε ϕ δ β
0.0089 2.7768 × 10−3 0.12735 0.15974 0.0222949
easily observed using proportions. Therefore, without loss of generality and for the sake of clarity, one gets for all time t, N (t ) + L(t ) + S (t ) + O(t ) + DS (t ) + DO (t ) = 1. The compartmental deterministic mathematical model can be represented analytically by the following nonlinear system of ordinary differential equations: N˙ (t ) = µ + ε DS (t ) − µN (t ) − β N (t ) [L(t ) + S (t ) + O(t )] , L˙ (t ) = β N (t ) [L(t ) + S (t ) + O(t )] − [µ + γL ] L(t ), S˙ (t ) = γL L(t ) + ϕ DS (t ) − [µ + γS + α ] S (t ),
(1)
˙ (t ) = γS S (t ) + δ DO (t ) − [µ + σ ] O(t ), O D˙S (t ) = γD DO (t ) + α S (t ) − [µ + ε + ϕ ] DS (t ), D˙O (t ) = σ O(t ) − [µ + γD + δ]DO (t ) where the parameters of the model are
• • • • • • • • • •
β , social transmission rate to unhealthy eating habits due to social environment pressure, µ, inflow rate population for the system, inversely proportional to the time spent by a child in the system, γL , rate at which a latent child becomes an overweight child, γS , rate at which an overweight child becomes an obese child by unhealthy habits, ε , rate at which an overweight child with healthy habits becomes a normal weighted child, α , rate at which an overweight child stops or reduces unhealthy habits, ϕ , rate at which an overweight child relapses to unhealthy habits, σ , rate at which an obese child stops or reduces unhealthy habits, δ , rate at which an obese child relapses to unhealthy habits, γD , rate at which an obese child with healthy habits become overweight with healthy habits.
The parameter values of the model (1) are presented in Table 1, and based on survey data, the initial conditions are assumed as N (0) = 0.462, L(0) = 0.194, S (0) = 0.2176, O(0) = 0.09, DS (0) = 0.0249 and DO (0) = 0.0115 [32]. Using these values, the model shows an increasing trend of obese and overweight populations. A linear stability analysis of system (1) was developed in [33]. The authors show that the system (1) has two steady states: a trivial steady state called the obesity-free equilibrium (OFE) and a non-trivial steady state called the obesity-endemic equilibrium (OEE). The nonlinear model (1) with the parameter values shown in Table 1 has the following two equilibrium points: OFE = (1, 0, 0, 0, 0, 0),
and
OEE ∼ = (0.293, 0.297, 0.271, 0.126, 0.008, 0.003).
Throughout this paper we denote by N ∗ , L∗ , S ∗ , O∗ , D∗S , D∗O the equilibrium points of system (1).
3. Network models The class of epidemiological models that is most widely used is the so-called homogenous models [1,2]. A homogenous model assumes that every individual has contacts to everyone else in the population, and that the rate of infection is largely determined by the density of the infected population. These models provide a good approximation of the disease spread in networks where the contact among individuals is sufficiently homogenous. Networks are perhaps the most natural abstract substrates to implement models of dynamical processes in many complex systems in which the relevant concept is the connectivity of the nodes and not their Euclidean distance. Here, we describe the stochastic contact network by a graph with N members of the population represented by nodes and with transits between states represented by N (N − 1)/2 edges. Each node is connected to all the other nodes; hence we have a complete graph [34]. An edge can be seen as a transit between states; that is, two vertices are connected by an edge if they have a flux of individuals in a regular way, and there exist a potential to transmit the disease if one of the nodes is infectious. The number of edges of a graph at a vertex is called the degree of the vertex. A well-known definition is the degree distribution of a graph pk , where pk is the fraction of vertices having degree k. The degree distribution is fundamental in the description of the spread of the model using heterogenous networks. However, in this study we assumed a homogenous network, where all vertices have contact with other vertices. Here, we assume that all contacts between an infective and a susceptible transmit infection with the same probability.
G. González-Parra et al. / Physica A 389 (2010) 3692–3701
3695
Simulating networks have the disadvantage of being computationally intensive: an epidemic in a population of just a few hundred thousand nodes can take hours of computation to obtain averaged results over a range of model parameters. However, that is a one-off cost of doing away with the assumption of a fully mixed population, and the range of possibilities it provides makes it a necessary trade-off. When large populations are simulated, a viable option is to consider a mean-field approach. A challenge in mathematical biology is to develop models that incorporate parameters with biological meaning. The exponential distribution plays a main role in the traditional classical mathematical epidemic models. Consider a death process and let the constant µ denote the per capita death rate of a population. Then its dynamics are modeled by the equation N˙ (t ) = −µN (t ),
0 ≤ µ < ∞,
N (0) = N0 ,
(2)
where N (t ) denotes the population size at time t. Hence, N (t ) = N0 e−µt , −µt
that is, e
for t ≥ 0;
(3)
denotes the probability of a individual being alive at time t ≥ 0 given that he was alive at time t = 0. Thus,
F (t ) =
1 − e−µt , 0,
for t ≥ 0, otherwise,
where F (t ) is a probability distribution that denotes the probability of dying in the time interval [0, t ). In this way, modeling the time to death X with an exponential probability distribution is equivalent to the following probability statement: Prob [X ≤ t] = F (t ).
(4)
The average time before death, or the life expectancy, is given by the mean or expected value of X ; that is,
Z E [X ] ≡
tf (t ) dt ,
(5)
where f (t ) = dF /dt denotes the probability density function of X . In order to simulate our models on the networks, we rely on Bayes’ theorem to obtain the following useful probabilistic relation [1]: Prob [X ≤ t + 1t | X > t] Prob [X > t] = Prob [t ≤ X ≤ t + 1t]
(6)
or Prob [X ≤ t + 1t | X > t] =
Prob [t ≤ X ≤ t + 1t] Prob [X > t]
.
(7)
In the particular case of the exponential distribution,
µe−µt 1t
= µ1t . (8) e−µt It is important to remark that in the case of the exponential distribution, Prob [X ≤ t + 1t | X > t] is independent of time. However, in the case of other distribution functions, such as the gamma function, this is not necessarily true. Exponential distributions are not only used to describe population dynamics, but also other biological processes such the dynamics of infectious diseases in populations. Intuitively, this does not seem correct because from the same time that individuals get infected there is a non-zero probability of recovery. Gamma distributions have been proposed in a more realistic assumption that suppresses early recovery by a particular function proportional to a power of time [28,27]: Prob [X ≤ t + 1t | X > t] ≈
f (x; k, θ ) = xk−1
e−x/θ
θ k Γ (k)
for x > 0 and k, θ > 0,
(9)
where the gamma distribution is parameterized in terms of a shape parameter k and scale parameter θ . The average waiting time is then given by ht i = kθ . In our case, the metabolic transitions by which the people change from latent to overweight and obese states should be better described by Eq. (9). Thus, if a gamma distribution is assumed for the waiting times between states, one gets Prob [X ≤ t + 1t | X > t] =
t k−1 e−x/θ
θ
kΓ
(k, x/θ )
,
(10)
where Γ (k, x/θ ) is the Plica function:
Γ (a, z ) =
∞
Z
t a−1 e−t dt
for x > 0 and k, θ > 0.
(11)
z
In order to simulate our models on the networks, we obtain different useful probabilistic relations for the particular cases k = 2, 3, 4. Thus, one gets for k = 2, Prob [X ≤ t + 1t | X > t] =
t
θ (θ + t )
,
(12)
3696
G. González-Parra et al. / Physica A 389 (2010) 3692–3701
where θ is the scale parameter of the gamma distribution and t is the time. Notice that, if t → ∞, the probability becomes 1/θ , and unlike the exponential distribution, the probability of death or of changing state at t = 0 is null. For k = 3, one gets Prob [X ≤ t + 1t | X > t] =
t 2 e−t /θ
θ 3 e−t /θ (2 + 2t /θ + (t /θ )2 )
=
t2
θ (2θ 2 + 2t θ + t 2 )
.
(13)
For k = 4, we obtain Prob [X ≤ t + 1t | X > t] =
t 3 e−t /θ
θ 4 Γ (4, t /θ )
=
t3
θ (6θ 3 + 6t θ 2 + 3t 2 θ + t 3 )
.
(14)
In order to allow that the mean duration of time periods on a state will be equal to the exponential distribution mean, the value of θ needs to be modified for each value of k, since kθ is the mean of the gamma distribution. 4. Simulations of models This section is devoted to presenting the simulation results from the application of constant probability and mean-field approach algorithms to a network representing the different classes of a child population using an epidemiological point of view. In addition, simulations are made assuming exponential and gamma probability distributions for the transition times between states (classes). Comparisons between the algorithms and different networks regarding probability distributions are performed in order to investigate the differences and advantages of each one. The algorithms are coded in the computer using Fortran language. Moreover, the calculations are based upon different numbers of runs and population size. However, only some representative cases are shown. 4.1. Simulation process In the simulation process, one means of avoiding passing through the entire network and reducing the simulation time is to perform the simulation according to the next event, that is, store in a separate list the events ordered by time. By running through this list, only the first event of the list is performed and then a new event is generated and substituted orderly on the list. In this way, one avoids having to check each node. However, since our graph size is not large, we rely on checking each node for each time step 1t. In order to investigate the dynamics of the population network it is necessary to reproduce the initial conditions of the population. Therefore, initially there is allocated in each epidemic class the corresponding nodes that mimic these conditions. In the simulation process we perform several time steps until the simulation time is achieved. At each successive time step 1t, each node is checked to see if a state transition occurs. Here 1t = 1 time unit because more than one transition state is able to take place in each step. We now explain simulation operations in greater detail. At each time step, children who have attained the maximum possible age of our 3–5-year-old population are removed. Since we assume that there is a constant population, a new child with normal weight is introduced into the network. The simulation process is performed in the following way at each time step.
• Every node of each one of the classes is visited in order to check if his state (class) changes. • A normal N (t ) node (normal weight child) becomes infected with unhealthy eating habits with probability 0 < β < 1 and • • • • • • • • •
becomes a latent node (latent individual). In order to check this social contagion process, we visit every latent, overweight and obese node in the population and select a pseudo-random number ξ . If ξ < β the normal weight node is infected. An L(t ) node (child) becomes an overweight S (t ) node with probability γL . An S (t ) node (child) becomes an obese O (t ) node with probability γS . An S (t ) node (child) moves to diet class DS (t ) node with probability α . A O(t ), node (child) moves to a diet class DO (t ) node with probability σ . An overweight on diet DS (t ) node (child) return to be an overweight S (t ) node with probability ϕ . An overweight on diet DS (t ) node (child) becomes a normal N (t ) node with probability ε . An obese on diet DO (t ) node (child) returns to being an obese O(t ) node with probability δ . An obese on diet DO (t ) node (child) becomes an overweight on diet DS (t ) node with probability γD . Every node (child) is removed with probability µ and replaced by another normal individual, in order to represent when a child becomes six years old (going out of the network).
We must also take into account that in our compartmental model several probabilistic transitions are possible for the next states. Here, the programming tree has been developed to randomly chose among all possible sequences of transition arrows a particular sequence for a given individual at a given time. Then, we check if the transition really takes place; if so, the remaining transitions are not checked. Otherwise, we follow the sequence until a transition occurs or every transition arrow is checked (by always selecting the same transition arrow for a given compartmental, every time step and every individual we could bias the simulation results). An important issue that needs to be considered in the simulation of the networks is that in the contagion process the value of the probability of contagion per unit time β needs to be modified for each network size N in order to enable comparison with the classical deterministic continuous model [35].
G. González-Parra et al. / Physica A 389 (2010) 3692–3701 0.5
0.5 Normal Weight Latent Overweight Obese Overweight Diet Obese Diet
0.45 0.4
0.4 0.35
0.3 0.25 0.2
0.3 0.25 0.2
0.15
0.15
0.1
0.1
0.05
0.05 0
100
200
300
400
500
Normal Weight Latent Overweight Obese Overweight Diet Obese Diet
0.45
Subpopulations
Subpopulations
0.35
0
3697
600
Time t (Weeks)
0
0
100
200
300
400
500
600
Time t (Weeks)
Fig. 1. Comparison of the numerical approximations solutions between the classical deterministic model (left) and the evolution of the constant rate homogenous network with N = 103 and 15 Monte Carlo simulations (right). Notice that the network evolves in a noisy way to the same steady-state OEE N ∗ = 0.293, L∗ = 0.2971, S ∗ = 0.271, O∗ , 0.126, D∗S = 0.008, D∗O = 0.003 of the deterministic model (1).
4.2. Networks with constant probability (CP) In this first simulation we consider a network with exponential transitions times in order to mimic the counterpart deterministic classical differential equation model. This simulation is applied using the conditional probability of transition from Eq. (8). In Fig. 1 there can be seen the comparison of the numerical approximations solutions between the classical deterministic model and the evolution of the constant probability homogenous network with N = 103 and 15 Monte Carlo simulations. Notice that the network evolves in a noisy way to the same steady-state OEE of the deterministic model. However, the inherent noisy evolution can be smoothed using more Monte Carlo simulations. When the number of nodes (children) is increased the computation time increases exponentially; therefore one needs to use another approach for larger networks such as the mean-field approach of the next section. As an example of computation time, a network with N = 104 and 15 Monte Carlo simulations can take several days and the mean-field approach takes only minutes. 4.3. Network model with mean-field approach (MF) In this simulation we also consider a network with exponential transition times in order to mimic the counterpart deterministic classical differential equation model. However, in this case we apply the mean-field algorithm to simulate the dynamics network in order to reduce the simulation time. The evolution of the epidemic network using the mean-field approach proceeds similarly to the constant probability network, but in the mean-field approach every node is checked only once each time step 1t. Instead, in the constant probability network each node is visited N 2 times each time step. Thus, the computation time of the mean-field approach is reduced considerably with respect to the previous network. Nevertheless, the mean-field approach can be applied only when the contagion probability β per unit time is small. Taking into account the simulation rules we can easily write the following mean-field evolution equation for L(t ): L(t + 1) = L(t ) + 1 − (1 − β)(L(t )+S (t )+O(t )) N (t ) − γL L(t ).
(15)
Here we have assumed the statistical independence of the social contagion process, i.e., the probability for a given normal node of not becoming contagious is given by (1 − β)(L(t )+S (t )+O(t )) [36]. Fig. 2 shows the mean-field approach simulation result. This result shows very good agreement with the classical differential equation model (1) and with the constant probability case, as was expected since the value of β is sufficiently small to be a good approximation. In addition, the profiles of the evolution of the different populations are smoother than for the constant probability network. This result is reasonable since in the mean-field approach the populations that are contagious are grouped as a mean value. However, averaging over a large number of independent realizations reduces the fluctuations, even if there are a finite number of nodes (individuals). 4.4. Gamma distribution network (GDN) Here, we investigate the effect of introducing the assumption that the waiting times between subpopulations (states) of the model follow a gamma distribution. This simple modeling assumption may be a driving force to change the dynamics of the classical exponential distribution model. In the differential equation model (1), one can assume a gamma distribution on the different transitions between states (classes). However, in order to study the dynamical changes which result from the inclusion of gamma distributions in the waiting time between states, we only consider the gamma distribution on the metabolic transitions by which the children change from latent to overweight and obese states.
3698
G. González-Parra et al. / Physica A 389 (2010) 3692–3701 0.5 Normal Weight Latent Overweight Obese Overweight Diet Obese Diet
0.45 0.4
Subpopulations
0.35 0.3 0.25 0.2 0.15 0.1 0.05 0
0
200
400
600
800
1000
Time t (Weeks)
Fig. 2. Evolution of the network model with the mean-field approach with N = 103 nodes and 15 Monte Carlo runs. Notice that the network evolves in a small noisy way to the same steady-state OEE N ∗ = 0.293, L∗ = 0.2971, S ∗ = 0.271, O∗ , 0.126, D∗S = 0.008, D∗O = 0.003 of the deterministic model (1). However, the noise is reduced with respect to the previous method (see Fig. 1). 0.5 0.45 0.4
Subpopulations
0.35 0.3 0.25 Normal Weight Latent Overweight Obese Overweight Diet Obese Diet
0.2 0.15 0.1 0.05 0
0
200
400
600
800
1000
Time t (Weeks)
Fig. 3. Evolution of the GDNLS (k = 2) network model with the CP algorithm, N = 103 nodes and 100 runs. Notice that the steady state for the latent subpopulation L∗ ≈ 0.32 has increased and the overweight subpopulation S ∗ ≈ 0.25 has decreased with respect to the deterministic model (1).
In addition to performing different simulations, some common settings are assumed. Throughout what follows, to enable comparisons between the previous networks and the gamma distribution network, the mean duration of time periods on a state ht i = kθ will always be taken to be equal to the corresponding mean of the exponential distribution networks. Therefore, we rely on the shape parameter k and scale parameter θ of the gamma function to accomplish this comparison requirement. The simulations are performed for different values of the parameters and with the values k = 2, 3, 4 and 5. Different initial conditions were also used but they showed no significant change in the results. Fig. 3 shows the network simulation result when a gamma distribution only on the transition time from the latent state L to the overweight state S (GDNLS) is assumed. This result shows that the steady state for the latent subpopulation L∗ ≈ 0.32 has been increased and the overweight subpopulation S ∗ ≈ 0.25 has decreased with respect to the deterministic model (1). This fact can be explained based on that fast transitions from a state have less probability with the gamma distribution. Nevertheless, the mean expected time for the transition is the same as that in the exponential case. Fig. 4 shows the network simulation result when a gamma distribution on the transition time from the overweight state S to the obese state O (GDNSO) is assumed. This result shows now that the steady state for the overweight subpopulation S ∗ ≈ 0.32 has increased and the obese subpopulation O∗ ≈ 0.08 has decreased with respect to the deterministic model (1). Fig. 5 shows the network simulation result when a gamma distribution on the two previously studied transition times is assumed. This result shows now that the steady state for the latent subpopulation L∗ ≈ 0.32 has increased and the obese subpopulation O∗ ≈ 0.08 has decreased with respect to the deterministic model (1). This particular result is the combination of the two previous cases and the results were expected. The last simulation results are obtained using the GDNLS network
G. González-Parra et al. / Physica A 389 (2010) 3692–3701
3699
0.5 0.45 0.4
Subpopulations
0.35 0.3 Normal Weight Latent Overweight Obese Overweight Diet Obese Diet
0.25 0.2 0.15 0.1 0.05 0
0
200
400
600
800
1000
Time t (Weeks)
Fig. 4. Evolution of the GDNLS (k = 2) network model with the CP algorithm, N = 103 nodes and 20 runs. Notice that the steady state for the overweight subpopulation S ∗ ≈ 0.32 has increased and the obese subpopulation O∗ ≈ 0.08 has decreased with respect to the deterministic model (1).
0.5 0.45 0.4
Subpopulations
0.35 0.3 Normal Weight
0.25
Latent 0.2
Overweight Obese
0.15
Overweight Diet 0.1
Obese Diet
0.05 0
0
200
400
600
800
1000
Time t (Weeks)
Fig. 5. Evolution of the combination of GDNLS and GDNSO previous networks with the CP algorithm, N = 103 nodes and 20 runs. Notice that the steady state for the latent subpopulation L∗ ≈ 0.32 has increased and the obese subpopulation O∗ ≈ 0.08 has decreased with respect to the deterministic model (1). Table 2 Steady states for the different simulated obesity dynamics networks using 103 nodes and 103 runs. Notice that the gamma distribution decreases the steady state for the overweight population when the parameter k is increased. Network
N∗
L∗
S∗
O∗
D∗S
D∗O
Classical Constant probability Mean field Gamma on L → S (k = 2) Gamma on L → S (k = 3) Gamma on L → S (k = 4) Gamma on L → S (k = 5) Gamma on S → O (k = 2)
0.293 0.293 0.293 0.2958 0.295 0.295 0.295 0.29
0.2971 0.2971 0.2971 0.3242 0.336 0.342 0.344 0.29
0.271 0.271 0.271 0.254 0.245 0.243 0.242 0.32
0.126 0.126 0.126 0.115 0.112 0.108 0.1076 0.08
0.008 0.008 0.008 0.007 0.007 0.007 0.007 0.008
0.003 0.003 0.003 0.003 0.003 0.003 0.003 0.001
with different values for the parameter k (k = 2, 3, 4 and 5). Fig. 6 shows the dynamics of the latent L(t ) and overweight S (t ) populations under a gamma distribution with k = 2 and k = 5. It can be observed that k = 5 leads to more rapid changes due to the less dispersed distribution. These results are in good agreement with the works of Lloyd [27,28]. The steady states for the network simulation are shown in Table 2. Notice the effect of introducing the assumption that the waiting times between subpopulations (states) of the model follow a gamma distribution. It can be observed that the increase of the values of k did not change the steady states significantly.
3700
G. González-Parra et al. / Physica A 389 (2010) 3692–3701 0.4
Subpopulations
0.35
0.3
0.25
0.2
Latent(k=2) Overweight(k=2) Latent(k=5) Overweight(k=5) 0
100
200
300
400
500
Time t (Weeks)
Fig. 6. Dynamics of latent L(t ) and overweight S (t ) populations under the gamma distribution with k = 2 and k = 5 with the CP algorithm, N = 103 nodes and 500 runs. It can be observed that k = 5 leads to more rapid changes due to less dispersed distributions, but the steady states are similar.
5. Discussion and conclusions In the context of epidemic diseases, several mathematical models have been used to study different demographic and epidemiological mechanisms that characterize disease dynamics. In this paper we extend a compartmental model to the case of a homogenous network epidemic model for the dynamics of an obese population. The social epidemic networks developed here use different approaches and points of view regarding the simulation of the dynamics of the network. First, Monte Carlo simulations for homogeneous networks using the traditional constant probability function and mean-field-like rate equations were performed. We show that these networks evolve towards an approximately stationary state, which is exactly equal to the one obtained by the underlying classical mathematical continuous model. A mean-field-like approach is applied in order to reduce the large computation time required when dealing with a large complete contact network. Another issue observed through this paper is that the profiles obtained for the populations by the simulation of the obesity epidemics on the homogenous network have a noisy signal. Hence, the onset of the realistic epidemic would be characterized by noisy dynamics as shown in Section 4, and this could be a driving force on the expansion of the epidemic for other models. However, averaging over a large number of realizations reduces the noise but increases the computation time. In addition, using homogenous contact networks, we investigated the effect of introducing the assumption that the waiting times between subpopulations (states) of the model follow a gamma distribution. Based on numerical results, we show that different distributions for the state periods lead to different results. However, the steady states of the subpopulations are relatively insensitive to the different values of k in the gamma distribution and this results agree with the result obtained by Lloyd [27]. This issue is of paramount importance in mathematical modeling since a model’s predictions need to be independent of the detailed assumptions of the underlying model. In this way, it is concluded that careful attention must be paid to the distributions assumed for the state periods in the mathematical models in order to obtain robust results. Finally, it can be stressed out that the approach presented here can provide a guide to examining alternative modeling of other epidemic-like models and that it provides new scenarios of epidemic propagation. Nevertheless, there are still many problems under debate, such as the role of network topology (random or scale-free) on epidemic dynamics. Acknowledgements We acknowledge the collaboration of the Oficina de Plan de Salud of the Conselleria de Sanitat of the Comunidad Valenciana. References [1] [2] [3] [4]
F. Brauer, C. Castillo-Chavez, Mathematical Models in Population Biology and Epidemiology, Springer Verlag, 2001. Herbert W. Hethcote, The mathematics of infectious diseases, SIAM Review 42 (4) (2000) 599–653. J.D. Murray, Mathematical Biology: I. An Introduction, Springer, Berlin, 2002. A. Noymer, The transmission and persistence of ‘urban legends’: sociological application of age-structured epidemic models, Journal of Mathematical Sociology 25 (3) (2001) 299–323. [5] Francisco J. Santonja, Ana C. Tarazona, Rafael J. Villanueva, A mathematical model of the pressure of an extreme ideology on a society, Computer & Mathematics with Applications 56 (3) (2008) 836–846. [6] Matt J. Keeling, Ken T.D. Eames, Networks and epidemic models, Journal of The Royal Society Interface 2 (4) (2005) 295–307.
G. González-Parra et al. / Physica A 389 (2010) 3692–3701
3701
[7] Fred Brauer, An Introduction to Networks in Epidemic Modeling, 2008. [8] Chengbin Peng, Xiaogang Jin, Meixia Shi, Epidemic threshold and immunization on generalized networks, Physica A: Statistical Mechanics and its Applications 389 (3) (2010) 549–560. [9] Juan Pablo Aparicio, Mercedes Pascual, Building epidemiological models from R0 : an implicit treatment of transmission in networks, Proceedings of the Royal Society B: Biological Sciences 274 (1609) (2007) 505–512. [10] Pieter Trapman, On analytical approaches to epidemics on networks, Theoretical Population Biology 71 (2) (2007) 160–173. [11] J. Verdasca, M.M. Telo da Gama, A. Nunes, N.R. Bernardino, J.M. Pacheco, M.C. Gomes, Recurrent epidemics in small world networks, Journal of Theoretical Biology 233 (4) (2005) 553–561. [12] F.C. Santos, J.F. Rodrigues, J.M. Pacheco, Epidemic spreading and cooperation dynamics on homogeneous small-world networks, Physical Review E 72 (5) (2005) 056128. [13] Thilo Gross, Carlos J. Dommar D’Lima, Bernd Blasius, Epidemic dynamics on an adaptive network, Phys. Rev. Lett. 96 (20) (2006) 208701. [14] Manojit Roy, Mercedes Pascual, On representing network heterogeneities in the incidence rate of simple epidemic models, Ecological Complexity 3 (1) (2006) 80–90. [15] Tao Zhou, Zhong-Qian Fu, Bing-Hong Wang, Epidemic dynamics on complex networks, Progress in Natural Science 16 (452) (2006). [16] Luciano da F. Costa, Francisco A. Rodrigues, Gonzalo Travieso, P.R. Villas Boas, Characterization of complex networks: a survey of measurements, Advances in Physics 56 (2007) 167–242. [17] J. Koopman, Controlling smallpox, Science 298 (2002) 1342–1344. [18] M.E. Halloran, I.M. Longini Jr., A. Nizam, Y. Yang, Containing bioterrorist smallpox, Science 298 (2002) 1428–1432. [19] E.H. Kaplan, D.L. Craft, L.M. Wein, Emergency response to a smallpox attack: the case for mass vaccination, Proceedings of the National Academy of Science 99 (2002) 10935–10940. [20] Darren M. Green, Istvan Z. Kiss, Rowland R. Kao, Parameterization of individual-based models: comparisons with deterministic mean-field models, Journal of Theoretical Biology 239 (3) (2006) 289–297. [21] Xiaoyan Wu, Zonghua Liu, How community structure influences epidemic spread in social networks, Physica A: Statistical Mechanics and its Applications 387 (2–3) (2008) 623–630. [22] Eduardo Massad, Stefen Ma, Mark Chen, Cláudio José Struchiner, Nico Stollenwerk, Maíra Aguiar, Scale-free network of a dengue epidemic, Applied Mathematics and Computation 195 (2) (2008) 376–381. [23] Romualdo Pastor-Satorras, Alessandro Vespignani, Epidemic spreading in scale-free networks, Physics Review Letters 86 (14) (2001) 3200–3203. [24] C.P. Warren, L.M. Sander, I.M. Sokolov, Geography in a scale-free network model, Physical Review E 66 (5) (2002) 056105. [25] J. Wu, Z. Gao, H. Sun, H. Huang, Urban transit system as a scale-free network, Modern Physics Letters B 18 (2004) 1043–1049. [26] L. Jódar, F. Santonja, G. González-Parra, Modeling dynamics of infant obesity in the region of Valencia, Spain, Computer & Mathematics with Applications 56 (3) (2008) 679–689. [27] Alun L. Lloyd, Destabilization of epidemic models with the inclusion of realistic distributions of infectious periods, Proceedings of the Royal Society B: Biological Sciences 268 (1470) (2001) 985–993. [28] Alun L. Lloyd, Realistic distributions of infectious periods in epidemic models: changing patterns of persistence and dynamics, Theoretical Population Biology 60 (1) (2001) 59–71. [29] Cara B. Ebbeling, Dorota B. Pawlak, David S. Ludwig, Childhood obesity: public-health crisis, common sense cure, The Lancet 360 (9331) (2002) 473–482. [30] Nicholas A. Christakis, James H. Fowler, The spread of obesity in a large social network over 32 years, N. Engl. J. Med. 357 (4) (2007) 370–379. [31] Maura MacPhee, Global childhood obesity: how to curb an epidemic, Journal of Pediatric Nursing 23 (1) (2008) 1–4. [32] A. Fullana, P. Momparler, J. Quiles, M.J. Redondo, Present situation of infant and adolescent obesity and prevention strategies 2005–2009, www.sp.san.gva.es/DgspPortal/docs/Informe_Obesidad.pdf. [33] G. González, L. Jódar, R. Villanueva, F. Santonja, Random modeling of population dynamics with uncertainty, Wseas Trans. Biol. Biomed. 5 (2) (2008) 34–45. [34] Reinhard Diestel, Graph Theory, Springer, 2000. [35] Yang Wang, Deepayan Chakrabarti, Christos Faloutsos, Chenxi Wang, Chenxi Wang, Epidemic spreading in real networks: an eigenvalue viewpoint, In SRDS (2003) 25–34. [36] L. Acedo, A second-order phase transition in the complete graph stochastic epidemic model, Physica A: Statistical Mechanics and its Applications 370 (2) (2006) 613–624.
Contents lists available at ScienceDirect
Physica A journal homepage: www.elsevier.com/locate/physa
Modeling the social obesity epidemic with stochastic networks Gilberto González-Parra a,b,∗ , L. Acedo b , Rafael-J. Villanueva Micó b , Abraham J. Arenas b,c a
Dpto. de Cálculo, Fac. Ingeniería, Universidad de los Andes, Mérida, Venezuela
b
Instituto de Matemática Multidisciplinar, Universidad Politécnica de Valencia, Edif. 8G, 2o , 46022 Valencia, Spain
c
Dpto. de Matemáticas y Estadística, Universidad de Córdoba, Montería, Colombia
article
info
Article history: Received 27 January 2010 Received in revised form 18 April 2010 Available online 29 April 2010 Keywords: Networks dynamics Populations Mean-field Gamma distribution Obesity epidemic model
abstract In this paper we extend a compartmental model to the case of a homogenous network epidemic model for a study of the dynamics of obese populations. The social epidemic network-based approach developed here uses different algorithms and points of views regarding the simulation of the dynamics of the network. First, Monte Carlo simulations for homogeneous networks using a traditional constant probability transition rates and a mean-field-like approach are presented. We show that these networks evolve towards an approximately stationary state, which coincides with the one obtained by the underlying classical compartmental continuous model. A mean-field-like approach is applied in order to reduce the large computation time required when dealing with large contact networks. We also investigate, using homogenous contact networks, the effect of the realistic assumption that the waiting times between subpopulations follow a gamma distribution instead of the traditional exponential distribution. It is concluded that careful attention must be paid to the distributions assumed for the state periods. © 2010 Elsevier B.V. All rights reserved.
1. Introduction Classical models of disease dynamics rely on systems of differential equations that represent the number of individuals in various categories through continuous variables, allowing for infinitesimal population densities. The origin of these models is commonly traced back to the well-known pioneering work of Kermack and McKendrick. Systems of ordinary differential equations (ODEs) are well-known tools that have been used to model different types of disease [1–3]. The most discussed type of infection spread is the SIR system, in which individuals are susceptible (S), infective (I) or removed/immune (R). In addition, several differential equation models have been proposed for modeling social behavior, such as rumors, social behavior and ideologies [1,3–5]. Understanding the nature and origin process of epidemic behavior has been a subject of considerable interest to mathematical epidemiologists. In order to model these epidemics it is important to choose a convenient framework. In this context we must acknowledge that the controversy between discrete individual models and continuous population models is an old one, and now a revival is taking place. The classical continuous approach to the modeling of epidemic spread suffers from some drawbacks; for example, it neglects the local character of epidemic spread. There are several important studies that deal with epidemics using networks [6–10]. In addition, different types of network such as small world and scale-free ones have been studied in order to understand epidemics [11–15]. In [16], a survey of the most important measurements of networks is presented. These measures give an insight into the understanding of the evolution and topological features of different network models.
∗
Corresponding author at: Dpto. de Cálculo, Fac. Ingeniería, Universidad de los Andes, Mérida, Venezuela. Tel.: +34 963879144; fax: +34 963879144. E-mail addresses: [email protected], [email protected] (G. González-Parra), [email protected] (L. Acedo), [email protected] (R.-J. Villanueva Micó). 0378-4371/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.physa.2010.04.024
G. González-Parra et al. / Physica A 389 (2010) 3692–3701
3693
In [17], there are shown some advantages of modeling epidemics with networks and as well with the continuous model. It is important to remark that in some cases continuous models and discrete network results may differ. For instance, the Halloran et al. results [18] do not agree with the recent findings of [17,19]. In [20], the authors develop an interesting study of an SIR model using networks under different algorithms and conclude that simulated epidemics generated under some model architectures are insensitive to the average degree of contact amongst nodes. The efficiency of the epidemic spread related to community structure and the high clustering coefficient has been also studied [21]. Continuous population models have the advantage of being easier to analyze theoretically than discrete individual models. In addition, simple deterministic models are still at the core of theoretical epidemiology despite the increasing evidence for the importance of contact networks underlying transmission at the individual level. These compartmental mathematical models based on homogeneous mixing have made, and are still making, important contributions to the epidemiology and the ecology of infectious diseases, but they fail to reproduce many of the features observed for the spreading of diseases in contact networks [9]. In this paper, we introduce several stochastic networks describing the obesity population dynamics in order to incorporate different assumptions and compare the results with the classical mathematical model. It is important to remark that most of the previous studies developed using networks focused on the general classical SI, SIS and SIR epidemic models. Here, we use a more complex network regarding classes or states, but heterogenous networks such as random, scale-free or small world are not considered in this paper. An interesting study is presented in [22] about the progression of dengue disease in Singapore, where the authors show that the dengue epidemic organized itself into a scale-free network. Other interesting applications of scale-free networks have been applied in different areas [23–25]. The study of graphs, which originated with the abstract theory of Erdös and Rényi of the 1950s and 1960s, has become important more recently in many areas, including social contacts and computer networks, as well as the spread of communicable diseases [7]. In this paper we describe the stochastic contact network by a complete graph with members of the population represented by vertices and with contacts between individuals represented by edges. An edge is a contact between vertices that can transmit infection. The number of edges of a graph at a vertex is called the degree of the vertex. A well-known definition is the degree distribution of a graph pk , where pk is the fraction of vertices having degree k. The degree distribution is fundamental in the description of the spread of the model using heterogenous networks. However, in this study we assumed a homogenous network, where all vertices have contact (connecting edges) with other vertices. Here, we assume that all contacts between an infective and a susceptible transmit the infection with the same probability. In [26], a compartmental deterministic mathematical model to study the dynamics of overweight and obesity in a childhood population was introduced. The authors consider obesity as a social transmission epidemic concern produced by unhealthy habits. Thus, our first aim is to simulate in detail the dynamics of the overweight and obese population on a homogenous network through numerical studies, and investigate the impact of the interaction rules on the efficiency and reliability of the simulation process. Therefore, different algorithms to simulate the network are performed in order to characterize the system dynamics. Simulation of large network models requires heavy computational cost; therefore, when large populations are simulated it is necessary to consider more efficient algorithms, for instance a mean-field approach or discrete events oriented simulation. One main assumption commonly employed in the classical continuous differential equation models is that the time for which individuals remain in one state can be described by an exponential distribution. This assumption avoids more complex equations such integro-differential equations. For instance, if a fixed lifetime is assumed in the deterministic model, then integro-differential equations arise in the model. However, real-world waiting times may have different probability distribution functions. Lloyd states that this distribution is biologically unrealistic, because it corresponds to the assumption that the chance of changing the state at a given time interval is independent of the time spend on the actual state. In addition, this assumption translates in that the distribution of remaining times in a state are too dispersed. In reality, the state periods are fairly closely centered about the mean duration of the period to state change [27]. Thus, another aim of this study is to investigate the effect of introducing the assumption that the waiting times between subpopulations (states) of the model follow a gamma distribution [27,28]. Using the simulation of the networks we investigate the effect that this has on the population dynamics. It is important to remark that simple modeling assumptions, such as the aforementioned, may be a driving force on the dynamics of the model; therefore their study is well justified. The layout of this paper is as follows. In Section 2, the underlying classical compartmental mathematical model is introduced. Section 3 is devoted to a brief introduction to networks modeling. Simulations of the proposed homogenous contact networks are presented in Section 4. A discussion and conclusions are presented in Section 5. 2. Classical mathematical obesity population model In several previous works, it is considered that the obesity problem is a result of social and cultural forces, including social networks [26,29–31]. In [26], a deterministic mathematical model to investigate the dynamics of overweight and obesity in a childhood population was introduced. In this model, the childhood population was divided into six disjoint classes or subpopulations: children with normal weight N (t ), latent children (children with unhealthy habits but who are still normal weight) L(t ), overweight children S (t ), obese children O(t ), overweight children on a diet DS (t ) and obese children on a diet DO (t ). The proposed model considers the 3–5-year-old population in order to study a homogenous group that is commonly studied. The model was normalized for the sake of clarity and since the qualitative behavior of the different classes is more
3694
G. González-Parra et al. / Physica A 389 (2010) 3692–3701
Table 1 Parameters of model (1), for the 3–5-year-old child population in the Spanish region of Valencia. Parameter
Value (weeks)
Parameter
Value (weeks)
µ γS α σ γD
0.0064 0.003085 4.068 × 10−3 4.4379 × 10−3 4.5045 × 10−4
γL ε ϕ δ β
0.0089 2.7768 × 10−3 0.12735 0.15974 0.0222949
easily observed using proportions. Therefore, without loss of generality and for the sake of clarity, one gets for all time t, N (t ) + L(t ) + S (t ) + O(t ) + DS (t ) + DO (t ) = 1. The compartmental deterministic mathematical model can be represented analytically by the following nonlinear system of ordinary differential equations: N˙ (t ) = µ + ε DS (t ) − µN (t ) − β N (t ) [L(t ) + S (t ) + O(t )] , L˙ (t ) = β N (t ) [L(t ) + S (t ) + O(t )] − [µ + γL ] L(t ), S˙ (t ) = γL L(t ) + ϕ DS (t ) − [µ + γS + α ] S (t ),
(1)
˙ (t ) = γS S (t ) + δ DO (t ) − [µ + σ ] O(t ), O D˙S (t ) = γD DO (t ) + α S (t ) − [µ + ε + ϕ ] DS (t ), D˙O (t ) = σ O(t ) − [µ + γD + δ]DO (t ) where the parameters of the model are
• • • • • • • • • •
β , social transmission rate to unhealthy eating habits due to social environment pressure, µ, inflow rate population for the system, inversely proportional to the time spent by a child in the system, γL , rate at which a latent child becomes an overweight child, γS , rate at which an overweight child becomes an obese child by unhealthy habits, ε , rate at which an overweight child with healthy habits becomes a normal weighted child, α , rate at which an overweight child stops or reduces unhealthy habits, ϕ , rate at which an overweight child relapses to unhealthy habits, σ , rate at which an obese child stops or reduces unhealthy habits, δ , rate at which an obese child relapses to unhealthy habits, γD , rate at which an obese child with healthy habits become overweight with healthy habits.
The parameter values of the model (1) are presented in Table 1, and based on survey data, the initial conditions are assumed as N (0) = 0.462, L(0) = 0.194, S (0) = 0.2176, O(0) = 0.09, DS (0) = 0.0249 and DO (0) = 0.0115 [32]. Using these values, the model shows an increasing trend of obese and overweight populations. A linear stability analysis of system (1) was developed in [33]. The authors show that the system (1) has two steady states: a trivial steady state called the obesity-free equilibrium (OFE) and a non-trivial steady state called the obesity-endemic equilibrium (OEE). The nonlinear model (1) with the parameter values shown in Table 1 has the following two equilibrium points: OFE = (1, 0, 0, 0, 0, 0),
and
OEE ∼ = (0.293, 0.297, 0.271, 0.126, 0.008, 0.003).
Throughout this paper we denote by N ∗ , L∗ , S ∗ , O∗ , D∗S , D∗O the equilibrium points of system (1).
3. Network models The class of epidemiological models that is most widely used is the so-called homogenous models [1,2]. A homogenous model assumes that every individual has contacts to everyone else in the population, and that the rate of infection is largely determined by the density of the infected population. These models provide a good approximation of the disease spread in networks where the contact among individuals is sufficiently homogenous. Networks are perhaps the most natural abstract substrates to implement models of dynamical processes in many complex systems in which the relevant concept is the connectivity of the nodes and not their Euclidean distance. Here, we describe the stochastic contact network by a graph with N members of the population represented by nodes and with transits between states represented by N (N − 1)/2 edges. Each node is connected to all the other nodes; hence we have a complete graph [34]. An edge can be seen as a transit between states; that is, two vertices are connected by an edge if they have a flux of individuals in a regular way, and there exist a potential to transmit the disease if one of the nodes is infectious. The number of edges of a graph at a vertex is called the degree of the vertex. A well-known definition is the degree distribution of a graph pk , where pk is the fraction of vertices having degree k. The degree distribution is fundamental in the description of the spread of the model using heterogenous networks. However, in this study we assumed a homogenous network, where all vertices have contact with other vertices. Here, we assume that all contacts between an infective and a susceptible transmit infection with the same probability.
G. González-Parra et al. / Physica A 389 (2010) 3692–3701
3695
Simulating networks have the disadvantage of being computationally intensive: an epidemic in a population of just a few hundred thousand nodes can take hours of computation to obtain averaged results over a range of model parameters. However, that is a one-off cost of doing away with the assumption of a fully mixed population, and the range of possibilities it provides makes it a necessary trade-off. When large populations are simulated, a viable option is to consider a mean-field approach. A challenge in mathematical biology is to develop models that incorporate parameters with biological meaning. The exponential distribution plays a main role in the traditional classical mathematical epidemic models. Consider a death process and let the constant µ denote the per capita death rate of a population. Then its dynamics are modeled by the equation N˙ (t ) = −µN (t ),
0 ≤ µ < ∞,
N (0) = N0 ,
(2)
where N (t ) denotes the population size at time t. Hence, N (t ) = N0 e−µt , −µt
that is, e
for t ≥ 0;
(3)
denotes the probability of a individual being alive at time t ≥ 0 given that he was alive at time t = 0. Thus,
F (t ) =
1 − e−µt , 0,
for t ≥ 0, otherwise,
where F (t ) is a probability distribution that denotes the probability of dying in the time interval [0, t ). In this way, modeling the time to death X with an exponential probability distribution is equivalent to the following probability statement: Prob [X ≤ t] = F (t ).
(4)
The average time before death, or the life expectancy, is given by the mean or expected value of X ; that is,
Z E [X ] ≡
tf (t ) dt ,
(5)
where f (t ) = dF /dt denotes the probability density function of X . In order to simulate our models on the networks, we rely on Bayes’ theorem to obtain the following useful probabilistic relation [1]: Prob [X ≤ t + 1t | X > t] Prob [X > t] = Prob [t ≤ X ≤ t + 1t]
(6)
or Prob [X ≤ t + 1t | X > t] =
Prob [t ≤ X ≤ t + 1t] Prob [X > t]
.
(7)
In the particular case of the exponential distribution,
µe−µt 1t
= µ1t . (8) e−µt It is important to remark that in the case of the exponential distribution, Prob [X ≤ t + 1t | X > t] is independent of time. However, in the case of other distribution functions, such as the gamma function, this is not necessarily true. Exponential distributions are not only used to describe population dynamics, but also other biological processes such the dynamics of infectious diseases in populations. Intuitively, this does not seem correct because from the same time that individuals get infected there is a non-zero probability of recovery. Gamma distributions have been proposed in a more realistic assumption that suppresses early recovery by a particular function proportional to a power of time [28,27]: Prob [X ≤ t + 1t | X > t] ≈
f (x; k, θ ) = xk−1
e−x/θ
θ k Γ (k)
for x > 0 and k, θ > 0,
(9)
where the gamma distribution is parameterized in terms of a shape parameter k and scale parameter θ . The average waiting time is then given by ht i = kθ . In our case, the metabolic transitions by which the people change from latent to overweight and obese states should be better described by Eq. (9). Thus, if a gamma distribution is assumed for the waiting times between states, one gets Prob [X ≤ t + 1t | X > t] =
t k−1 e−x/θ
θ
kΓ
(k, x/θ )
,
(10)
where Γ (k, x/θ ) is the Plica function:
Γ (a, z ) =
∞
Z
t a−1 e−t dt
for x > 0 and k, θ > 0.
(11)
z
In order to simulate our models on the networks, we obtain different useful probabilistic relations for the particular cases k = 2, 3, 4. Thus, one gets for k = 2, Prob [X ≤ t + 1t | X > t] =
t
θ (θ + t )
,
(12)
3696
G. González-Parra et al. / Physica A 389 (2010) 3692–3701
where θ is the scale parameter of the gamma distribution and t is the time. Notice that, if t → ∞, the probability becomes 1/θ , and unlike the exponential distribution, the probability of death or of changing state at t = 0 is null. For k = 3, one gets Prob [X ≤ t + 1t | X > t] =
t 2 e−t /θ
θ 3 e−t /θ (2 + 2t /θ + (t /θ )2 )
=
t2
θ (2θ 2 + 2t θ + t 2 )
.
(13)
For k = 4, we obtain Prob [X ≤ t + 1t | X > t] =
t 3 e−t /θ
θ 4 Γ (4, t /θ )
=
t3
θ (6θ 3 + 6t θ 2 + 3t 2 θ + t 3 )
.
(14)
In order to allow that the mean duration of time periods on a state will be equal to the exponential distribution mean, the value of θ needs to be modified for each value of k, since kθ is the mean of the gamma distribution. 4. Simulations of models This section is devoted to presenting the simulation results from the application of constant probability and mean-field approach algorithms to a network representing the different classes of a child population using an epidemiological point of view. In addition, simulations are made assuming exponential and gamma probability distributions for the transition times between states (classes). Comparisons between the algorithms and different networks regarding probability distributions are performed in order to investigate the differences and advantages of each one. The algorithms are coded in the computer using Fortran language. Moreover, the calculations are based upon different numbers of runs and population size. However, only some representative cases are shown. 4.1. Simulation process In the simulation process, one means of avoiding passing through the entire network and reducing the simulation time is to perform the simulation according to the next event, that is, store in a separate list the events ordered by time. By running through this list, only the first event of the list is performed and then a new event is generated and substituted orderly on the list. In this way, one avoids having to check each node. However, since our graph size is not large, we rely on checking each node for each time step 1t. In order to investigate the dynamics of the population network it is necessary to reproduce the initial conditions of the population. Therefore, initially there is allocated in each epidemic class the corresponding nodes that mimic these conditions. In the simulation process we perform several time steps until the simulation time is achieved. At each successive time step 1t, each node is checked to see if a state transition occurs. Here 1t = 1 time unit because more than one transition state is able to take place in each step. We now explain simulation operations in greater detail. At each time step, children who have attained the maximum possible age of our 3–5-year-old population are removed. Since we assume that there is a constant population, a new child with normal weight is introduced into the network. The simulation process is performed in the following way at each time step.
• Every node of each one of the classes is visited in order to check if his state (class) changes. • A normal N (t ) node (normal weight child) becomes infected with unhealthy eating habits with probability 0 < β < 1 and • • • • • • • • •
becomes a latent node (latent individual). In order to check this social contagion process, we visit every latent, overweight and obese node in the population and select a pseudo-random number ξ . If ξ < β the normal weight node is infected. An L(t ) node (child) becomes an overweight S (t ) node with probability γL . An S (t ) node (child) becomes an obese O (t ) node with probability γS . An S (t ) node (child) moves to diet class DS (t ) node with probability α . A O(t ), node (child) moves to a diet class DO (t ) node with probability σ . An overweight on diet DS (t ) node (child) return to be an overweight S (t ) node with probability ϕ . An overweight on diet DS (t ) node (child) becomes a normal N (t ) node with probability ε . An obese on diet DO (t ) node (child) returns to being an obese O(t ) node with probability δ . An obese on diet DO (t ) node (child) becomes an overweight on diet DS (t ) node with probability γD . Every node (child) is removed with probability µ and replaced by another normal individual, in order to represent when a child becomes six years old (going out of the network).
We must also take into account that in our compartmental model several probabilistic transitions are possible for the next states. Here, the programming tree has been developed to randomly chose among all possible sequences of transition arrows a particular sequence for a given individual at a given time. Then, we check if the transition really takes place; if so, the remaining transitions are not checked. Otherwise, we follow the sequence until a transition occurs or every transition arrow is checked (by always selecting the same transition arrow for a given compartmental, every time step and every individual we could bias the simulation results). An important issue that needs to be considered in the simulation of the networks is that in the contagion process the value of the probability of contagion per unit time β needs to be modified for each network size N in order to enable comparison with the classical deterministic continuous model [35].
G. González-Parra et al. / Physica A 389 (2010) 3692–3701 0.5
0.5 Normal Weight Latent Overweight Obese Overweight Diet Obese Diet
0.45 0.4
0.4 0.35
0.3 0.25 0.2
0.3 0.25 0.2
0.15
0.15
0.1
0.1
0.05
0.05 0
100
200
300
400
500
Normal Weight Latent Overweight Obese Overweight Diet Obese Diet
0.45
Subpopulations
Subpopulations
0.35
0
3697
600
Time t (Weeks)
0
0
100
200
300
400
500
600
Time t (Weeks)
Fig. 1. Comparison of the numerical approximations solutions between the classical deterministic model (left) and the evolution of the constant rate homogenous network with N = 103 and 15 Monte Carlo simulations (right). Notice that the network evolves in a noisy way to the same steady-state OEE N ∗ = 0.293, L∗ = 0.2971, S ∗ = 0.271, O∗ , 0.126, D∗S = 0.008, D∗O = 0.003 of the deterministic model (1).
4.2. Networks with constant probability (CP) In this first simulation we consider a network with exponential transitions times in order to mimic the counterpart deterministic classical differential equation model. This simulation is applied using the conditional probability of transition from Eq. (8). In Fig. 1 there can be seen the comparison of the numerical approximations solutions between the classical deterministic model and the evolution of the constant probability homogenous network with N = 103 and 15 Monte Carlo simulations. Notice that the network evolves in a noisy way to the same steady-state OEE of the deterministic model. However, the inherent noisy evolution can be smoothed using more Monte Carlo simulations. When the number of nodes (children) is increased the computation time increases exponentially; therefore one needs to use another approach for larger networks such as the mean-field approach of the next section. As an example of computation time, a network with N = 104 and 15 Monte Carlo simulations can take several days and the mean-field approach takes only minutes. 4.3. Network model with mean-field approach (MF) In this simulation we also consider a network with exponential transition times in order to mimic the counterpart deterministic classical differential equation model. However, in this case we apply the mean-field algorithm to simulate the dynamics network in order to reduce the simulation time. The evolution of the epidemic network using the mean-field approach proceeds similarly to the constant probability network, but in the mean-field approach every node is checked only once each time step 1t. Instead, in the constant probability network each node is visited N 2 times each time step. Thus, the computation time of the mean-field approach is reduced considerably with respect to the previous network. Nevertheless, the mean-field approach can be applied only when the contagion probability β per unit time is small. Taking into account the simulation rules we can easily write the following mean-field evolution equation for L(t ): L(t + 1) = L(t ) + 1 − (1 − β)(L(t )+S (t )+O(t )) N (t ) − γL L(t ).
(15)
Here we have assumed the statistical independence of the social contagion process, i.e., the probability for a given normal node of not becoming contagious is given by (1 − β)(L(t )+S (t )+O(t )) [36]. Fig. 2 shows the mean-field approach simulation result. This result shows very good agreement with the classical differential equation model (1) and with the constant probability case, as was expected since the value of β is sufficiently small to be a good approximation. In addition, the profiles of the evolution of the different populations are smoother than for the constant probability network. This result is reasonable since in the mean-field approach the populations that are contagious are grouped as a mean value. However, averaging over a large number of independent realizations reduces the fluctuations, even if there are a finite number of nodes (individuals). 4.4. Gamma distribution network (GDN) Here, we investigate the effect of introducing the assumption that the waiting times between subpopulations (states) of the model follow a gamma distribution. This simple modeling assumption may be a driving force to change the dynamics of the classical exponential distribution model. In the differential equation model (1), one can assume a gamma distribution on the different transitions between states (classes). However, in order to study the dynamical changes which result from the inclusion of gamma distributions in the waiting time between states, we only consider the gamma distribution on the metabolic transitions by which the children change from latent to overweight and obese states.
3698
G. González-Parra et al. / Physica A 389 (2010) 3692–3701 0.5 Normal Weight Latent Overweight Obese Overweight Diet Obese Diet
0.45 0.4
Subpopulations
0.35 0.3 0.25 0.2 0.15 0.1 0.05 0
0
200
400
600
800
1000
Time t (Weeks)
Fig. 2. Evolution of the network model with the mean-field approach with N = 103 nodes and 15 Monte Carlo runs. Notice that the network evolves in a small noisy way to the same steady-state OEE N ∗ = 0.293, L∗ = 0.2971, S ∗ = 0.271, O∗ , 0.126, D∗S = 0.008, D∗O = 0.003 of the deterministic model (1). However, the noise is reduced with respect to the previous method (see Fig. 1). 0.5 0.45 0.4
Subpopulations
0.35 0.3 0.25 Normal Weight Latent Overweight Obese Overweight Diet Obese Diet
0.2 0.15 0.1 0.05 0
0
200
400
600
800
1000
Time t (Weeks)
Fig. 3. Evolution of the GDNLS (k = 2) network model with the CP algorithm, N = 103 nodes and 100 runs. Notice that the steady state for the latent subpopulation L∗ ≈ 0.32 has increased and the overweight subpopulation S ∗ ≈ 0.25 has decreased with respect to the deterministic model (1).
In addition to performing different simulations, some common settings are assumed. Throughout what follows, to enable comparisons between the previous networks and the gamma distribution network, the mean duration of time periods on a state ht i = kθ will always be taken to be equal to the corresponding mean of the exponential distribution networks. Therefore, we rely on the shape parameter k and scale parameter θ of the gamma function to accomplish this comparison requirement. The simulations are performed for different values of the parameters and with the values k = 2, 3, 4 and 5. Different initial conditions were also used but they showed no significant change in the results. Fig. 3 shows the network simulation result when a gamma distribution only on the transition time from the latent state L to the overweight state S (GDNLS) is assumed. This result shows that the steady state for the latent subpopulation L∗ ≈ 0.32 has been increased and the overweight subpopulation S ∗ ≈ 0.25 has decreased with respect to the deterministic model (1). This fact can be explained based on that fast transitions from a state have less probability with the gamma distribution. Nevertheless, the mean expected time for the transition is the same as that in the exponential case. Fig. 4 shows the network simulation result when a gamma distribution on the transition time from the overweight state S to the obese state O (GDNSO) is assumed. This result shows now that the steady state for the overweight subpopulation S ∗ ≈ 0.32 has increased and the obese subpopulation O∗ ≈ 0.08 has decreased with respect to the deterministic model (1). Fig. 5 shows the network simulation result when a gamma distribution on the two previously studied transition times is assumed. This result shows now that the steady state for the latent subpopulation L∗ ≈ 0.32 has increased and the obese subpopulation O∗ ≈ 0.08 has decreased with respect to the deterministic model (1). This particular result is the combination of the two previous cases and the results were expected. The last simulation results are obtained using the GDNLS network
G. González-Parra et al. / Physica A 389 (2010) 3692–3701
3699
0.5 0.45 0.4
Subpopulations
0.35 0.3 Normal Weight Latent Overweight Obese Overweight Diet Obese Diet
0.25 0.2 0.15 0.1 0.05 0
0
200
400
600
800
1000
Time t (Weeks)
Fig. 4. Evolution of the GDNLS (k = 2) network model with the CP algorithm, N = 103 nodes and 20 runs. Notice that the steady state for the overweight subpopulation S ∗ ≈ 0.32 has increased and the obese subpopulation O∗ ≈ 0.08 has decreased with respect to the deterministic model (1).
0.5 0.45 0.4
Subpopulations
0.35 0.3 Normal Weight
0.25
Latent 0.2
Overweight Obese
0.15
Overweight Diet 0.1
Obese Diet
0.05 0
0
200
400
600
800
1000
Time t (Weeks)
Fig. 5. Evolution of the combination of GDNLS and GDNSO previous networks with the CP algorithm, N = 103 nodes and 20 runs. Notice that the steady state for the latent subpopulation L∗ ≈ 0.32 has increased and the obese subpopulation O∗ ≈ 0.08 has decreased with respect to the deterministic model (1). Table 2 Steady states for the different simulated obesity dynamics networks using 103 nodes and 103 runs. Notice that the gamma distribution decreases the steady state for the overweight population when the parameter k is increased. Network
N∗
L∗
S∗
O∗
D∗S
D∗O
Classical Constant probability Mean field Gamma on L → S (k = 2) Gamma on L → S (k = 3) Gamma on L → S (k = 4) Gamma on L → S (k = 5) Gamma on S → O (k = 2)
0.293 0.293 0.293 0.2958 0.295 0.295 0.295 0.29
0.2971 0.2971 0.2971 0.3242 0.336 0.342 0.344 0.29
0.271 0.271 0.271 0.254 0.245 0.243 0.242 0.32
0.126 0.126 0.126 0.115 0.112 0.108 0.1076 0.08
0.008 0.008 0.008 0.007 0.007 0.007 0.007 0.008
0.003 0.003 0.003 0.003 0.003 0.003 0.003 0.001
with different values for the parameter k (k = 2, 3, 4 and 5). Fig. 6 shows the dynamics of the latent L(t ) and overweight S (t ) populations under a gamma distribution with k = 2 and k = 5. It can be observed that k = 5 leads to more rapid changes due to the less dispersed distribution. These results are in good agreement with the works of Lloyd [27,28]. The steady states for the network simulation are shown in Table 2. Notice the effect of introducing the assumption that the waiting times between subpopulations (states) of the model follow a gamma distribution. It can be observed that the increase of the values of k did not change the steady states significantly.
3700
G. González-Parra et al. / Physica A 389 (2010) 3692–3701 0.4
Subpopulations
0.35
0.3
0.25
0.2
Latent(k=2) Overweight(k=2) Latent(k=5) Overweight(k=5) 0
100
200
300
400
500
Time t (Weeks)
Fig. 6. Dynamics of latent L(t ) and overweight S (t ) populations under the gamma distribution with k = 2 and k = 5 with the CP algorithm, N = 103 nodes and 500 runs. It can be observed that k = 5 leads to more rapid changes due to less dispersed distributions, but the steady states are similar.
5. Discussion and conclusions In the context of epidemic diseases, several mathematical models have been used to study different demographic and epidemiological mechanisms that characterize disease dynamics. In this paper we extend a compartmental model to the case of a homogenous network epidemic model for the dynamics of an obese population. The social epidemic networks developed here use different approaches and points of view regarding the simulation of the dynamics of the network. First, Monte Carlo simulations for homogeneous networks using the traditional constant probability function and mean-field-like rate equations were performed. We show that these networks evolve towards an approximately stationary state, which is exactly equal to the one obtained by the underlying classical mathematical continuous model. A mean-field-like approach is applied in order to reduce the large computation time required when dealing with a large complete contact network. Another issue observed through this paper is that the profiles obtained for the populations by the simulation of the obesity epidemics on the homogenous network have a noisy signal. Hence, the onset of the realistic epidemic would be characterized by noisy dynamics as shown in Section 4, and this could be a driving force on the expansion of the epidemic for other models. However, averaging over a large number of realizations reduces the noise but increases the computation time. In addition, using homogenous contact networks, we investigated the effect of introducing the assumption that the waiting times between subpopulations (states) of the model follow a gamma distribution. Based on numerical results, we show that different distributions for the state periods lead to different results. However, the steady states of the subpopulations are relatively insensitive to the different values of k in the gamma distribution and this results agree with the result obtained by Lloyd [27]. This issue is of paramount importance in mathematical modeling since a model’s predictions need to be independent of the detailed assumptions of the underlying model. In this way, it is concluded that careful attention must be paid to the distributions assumed for the state periods in the mathematical models in order to obtain robust results. Finally, it can be stressed out that the approach presented here can provide a guide to examining alternative modeling of other epidemic-like models and that it provides new scenarios of epidemic propagation. Nevertheless, there are still many problems under debate, such as the role of network topology (random or scale-free) on epidemic dynamics. Acknowledgements We acknowledge the collaboration of the Oficina de Plan de Salud of the Conselleria de Sanitat of the Comunidad Valenciana. References [1] [2] [3] [4]
F. Brauer, C. Castillo-Chavez, Mathematical Models in Population Biology and Epidemiology, Springer Verlag, 2001. Herbert W. Hethcote, The mathematics of infectious diseases, SIAM Review 42 (4) (2000) 599–653. J.D. Murray, Mathematical Biology: I. An Introduction, Springer, Berlin, 2002. A. Noymer, The transmission and persistence of ‘urban legends’: sociological application of age-structured epidemic models, Journal of Mathematical Sociology 25 (3) (2001) 299–323. [5] Francisco J. Santonja, Ana C. Tarazona, Rafael J. Villanueva, A mathematical model of the pressure of an extreme ideology on a society, Computer & Mathematics with Applications 56 (3) (2008) 836–846. [6] Matt J. Keeling, Ken T.D. Eames, Networks and epidemic models, Journal of The Royal Society Interface 2 (4) (2005) 295–307.
G. González-Parra et al. / Physica A 389 (2010) 3692–3701
3701
[7] Fred Brauer, An Introduction to Networks in Epidemic Modeling, 2008. [8] Chengbin Peng, Xiaogang Jin, Meixia Shi, Epidemic threshold and immunization on generalized networks, Physica A: Statistical Mechanics and its Applications 389 (3) (2010) 549–560. [9] Juan Pablo Aparicio, Mercedes Pascual, Building epidemiological models from R0 : an implicit treatment of transmission in networks, Proceedings of the Royal Society B: Biological Sciences 274 (1609) (2007) 505–512. [10] Pieter Trapman, On analytical approaches to epidemics on networks, Theoretical Population Biology 71 (2) (2007) 160–173. [11] J. Verdasca, M.M. Telo da Gama, A. Nunes, N.R. Bernardino, J.M. Pacheco, M.C. Gomes, Recurrent epidemics in small world networks, Journal of Theoretical Biology 233 (4) (2005) 553–561. [12] F.C. Santos, J.F. Rodrigues, J.M. Pacheco, Epidemic spreading and cooperation dynamics on homogeneous small-world networks, Physical Review E 72 (5) (2005) 056128. [13] Thilo Gross, Carlos J. Dommar D’Lima, Bernd Blasius, Epidemic dynamics on an adaptive network, Phys. Rev. Lett. 96 (20) (2006) 208701. [14] Manojit Roy, Mercedes Pascual, On representing network heterogeneities in the incidence rate of simple epidemic models, Ecological Complexity 3 (1) (2006) 80–90. [15] Tao Zhou, Zhong-Qian Fu, Bing-Hong Wang, Epidemic dynamics on complex networks, Progress in Natural Science 16 (452) (2006). [16] Luciano da F. Costa, Francisco A. Rodrigues, Gonzalo Travieso, P.R. Villas Boas, Characterization of complex networks: a survey of measurements, Advances in Physics 56 (2007) 167–242. [17] J. Koopman, Controlling smallpox, Science 298 (2002) 1342–1344. [18] M.E. Halloran, I.M. Longini Jr., A. Nizam, Y. Yang, Containing bioterrorist smallpox, Science 298 (2002) 1428–1432. [19] E.H. Kaplan, D.L. Craft, L.M. Wein, Emergency response to a smallpox attack: the case for mass vaccination, Proceedings of the National Academy of Science 99 (2002) 10935–10940. [20] Darren M. Green, Istvan Z. Kiss, Rowland R. Kao, Parameterization of individual-based models: comparisons with deterministic mean-field models, Journal of Theoretical Biology 239 (3) (2006) 289–297. [21] Xiaoyan Wu, Zonghua Liu, How community structure influences epidemic spread in social networks, Physica A: Statistical Mechanics and its Applications 387 (2–3) (2008) 623–630. [22] Eduardo Massad, Stefen Ma, Mark Chen, Cláudio José Struchiner, Nico Stollenwerk, Maíra Aguiar, Scale-free network of a dengue epidemic, Applied Mathematics and Computation 195 (2) (2008) 376–381. [23] Romualdo Pastor-Satorras, Alessandro Vespignani, Epidemic spreading in scale-free networks, Physics Review Letters 86 (14) (2001) 3200–3203. [24] C.P. Warren, L.M. Sander, I.M. Sokolov, Geography in a scale-free network model, Physical Review E 66 (5) (2002) 056105. [25] J. Wu, Z. Gao, H. Sun, H. Huang, Urban transit system as a scale-free network, Modern Physics Letters B 18 (2004) 1043–1049. [26] L. Jódar, F. Santonja, G. González-Parra, Modeling dynamics of infant obesity in the region of Valencia, Spain, Computer & Mathematics with Applications 56 (3) (2008) 679–689. [27] Alun L. Lloyd, Destabilization of epidemic models with the inclusion of realistic distributions of infectious periods, Proceedings of the Royal Society B: Biological Sciences 268 (1470) (2001) 985–993. [28] Alun L. Lloyd, Realistic distributions of infectious periods in epidemic models: changing patterns of persistence and dynamics, Theoretical Population Biology 60 (1) (2001) 59–71. [29] Cara B. Ebbeling, Dorota B. Pawlak, David S. Ludwig, Childhood obesity: public-health crisis, common sense cure, The Lancet 360 (9331) (2002) 473–482. [30] Nicholas A. Christakis, James H. Fowler, The spread of obesity in a large social network over 32 years, N. Engl. J. Med. 357 (4) (2007) 370–379. [31] Maura MacPhee, Global childhood obesity: how to curb an epidemic, Journal of Pediatric Nursing 23 (1) (2008) 1–4. [32] A. Fullana, P. Momparler, J. Quiles, M.J. Redondo, Present situation of infant and adolescent obesity and prevention strategies 2005–2009, www.sp.san.gva.es/DgspPortal/docs/Informe_Obesidad.pdf. [33] G. González, L. Jódar, R. Villanueva, F. Santonja, Random modeling of population dynamics with uncertainty, Wseas Trans. Biol. Biomed. 5 (2) (2008) 34–45. [34] Reinhard Diestel, Graph Theory, Springer, 2000. [35] Yang Wang, Deepayan Chakrabarti, Christos Faloutsos, Chenxi Wang, Chenxi Wang, Epidemic spreading in real networks: an eigenvalue viewpoint, In SRDS (2003) 25–34. [36] L. Acedo, A second-order phase transition in the complete graph stochastic epidemic model, Physica A: Statistical Mechanics and its Applications 370 (2) (2006) 613–624.