The analysis of the SIRS alcoholism models with relapse on weighted networks.

Two SIRS alcoholism models with relapse on networks with fixed and adaptive weight are introduced. The spread of alcoholism threshold [Formula: see text] is calculated by the next generation matrix method. For the model with fixed weight, we prove that when [Formula: see text] the alcohol free equilibrium is globally asymptotically stable, then the drinking crowd gradually disappear. When [Formula: see text], the alcoholism equilibrium is global attractivity, then the density of alcoholics will remain in a stable value. For the model with adaptive weight, we only make some numerical simulations. We also give two effective strategies. Our results show that the treatment of recuperator for stopping relapsing and preventing the susceptible people to drink are two effective measures to eliminate alcoholism problem, and preventing the susceptible people to drink is more effective when the proportion of recuperator to accept treatment is equal to the proportion of susceptible people to refuse drinking alcohol.

Background

Alcohol use and abuse have been part of human society for centuries. The College Alcohol Study defines alcoholism as male students who had five or more and female students who had four or more drinks in a row at least once in a 2-week period (Wechsler 2000). Similar to other drug addictions, alcoholism can cause to a series of serious consequences. The World Health Organization estimates that about 140 million people around the world under the influence of alcohol-related problems, such as being sick, losing a job and so on (Saunders et al. 1993). What is more serious is that chronic alcohol consumption damages almost all parts of the body and contribute to a number of human diseases including but not limited to liver cirrhosis, pancreatitis, heart disease, and sexual dysfunction and eventually be deadly (Glavas and Weinberg 2005). Early recognized since the 1800s that alcoholism produced not only impairment of the senses but also higher predisposition for tuberculosis. William Osler reported in 1905 that patients who misused alcohol had higher predisposition to pneumonia (Giraldina et al. 2015). And studies over the last 30 years have also demonstrated that chronic alcohol consumption impairs the functions of both T cells and B cells (Sumana et al. 2015). Between 2006 and 2010 in the United States, alcohol abuse resulted in approximately 88,000 deaths, and the average death rate associated with alcoholism was 28.5 per 100,000 population (Gonzales et al. 2014). Importantly, the majority of alcohol-related deaths were among adults aged 20–64 years old (Giraldina et al. 2015). In view of the above situation, alcoholism has become a issue that need to be solved urgently.

Since mathematical model can mimic the process of alcoholism and provide useful methods to control the spread of drinking behavior. Several different mathematical models for alcoholism have been formulated and studied recently. Sanchez et al. (2007) presented a simple model for alcohol treatment. Their model were based on studying binge drinking in a college system and assumed the same “leaving rate”. Manthey et al. (2008) built a model to capture the dynamics of campus drinking and to study the spread of drinking on campus. Benedict (2007) proposed an SIR model and used standard contact rate between susceptible and alcoholics. Furthermore they obtained the alcoholism reproductive number and discussed the existence and stability of all the equilibria. Huo and Song (2012) introduced a two-stage model for binge drinking problem, which the youths with alcohol problems were divided into those who admit the problem and those who do not admit it. Wang et al. (2014) presented a deterministic SATQ-type mathematical model for the spread of alcoholism with two control strategies and analyzed some properties of the solutions including positivity, existence and stability. Huo and Wang (2014) developed a nonlinear mathematical model with the effect of awareness programs on the binge drinking. Their results showed that awareness programs is an effective measure in reducing alcohol problems.

Quit drinking is usually temporary. Some drinking people may relapse since contacting with alcoholics or weak self-control ability. Sharma and Samanta (2015) developed an alcohol abuse model by introducing a treatment programm in the population and considered all possible relapses. They assumed that the drinkers in treatment most commonly relapse due to contact with heavy drinkers who are not in treatment. For the other mathematical models for alcoholism or smoking, please see Wechsler and Nelson (2008), Room et al. (2005), Mushayabasa and Bhunu (2011), Huo and Zhu (2013) and references cited therein. They commonly assume that communities are homogeneous, that is, communities are made up of individuals who mix uniformly and randomly with each other in the above models. These assumptions make the analysis tractable but not realistic (Bansal et al. 2007).

In contrast to classical compartment models, a lot of studies on complex networks have been investigated during the past years. Liu et al. (2013) presented an SIR model with individual’s birth and death on scale-free networks and analyzed the stability of three equilibria. They also gave out two immunization schemes. But they didn’t consider recuperator’s temporarily immune, that is to say, recuperator is likely to become infected or susceptible because of the loss of immune. Zhu et al. (2012a) investigated a new epidemic SIS model with nonlinear infectivity on heterogeneous networks. The global behavior of the model is studied. Wang et al. (2012) proposed a modified SIS model with an infective vector on complex networks. They treated direct human contacts as a social network and assumed spatially homogeneous mixing between vector and human populations. Huo and Liu (2016) proposed an alcoholism model on complex heterogeneous networks and proved stability of all the equilibria.

Nodes usually represent individuals and links represent potential contacts among those individuals on the complex network (Zhu et al. 2012a). The connectivity of a node is defined as the number of the links connected to the node, represent by k. The degree distribution of a network is defined as the probability of a randomly chosen node to have a degree k, represent by P(k). Many networks (Zhu et al. 2012b; Liu and Zhang 2011; Zhang and Jin 2011; Liu et al. 2004) have been found to be scale-free networks, that is to say the degree distribution follows a power law distribution , where c is any constant satisfy the equation .

Notice that in many real contact networks, there are groups of nodes with a high density of edges within them and a lower density of links between groups. The differences between links within a contact network can be described by link weights, which can represent the amount of time two individuals interact or the intimacy between individuals. The larger the weight is, the more the two nodes communicate, while, the more possible a susceptible individual will be infected through the edge (Chu et al. 2009). The usual assumption is that weights are constant and driven by the network connectivity, which is fixed as time goes on. For example, the weight between two nodes with degrees i and j are represented by a function of their degrees (Barrat et al. 2004a, b, c). However, as the disease becomes severe, individuals tend to be more cautious in social contacts and make some reflection such as decreasing the out going visits, cutting down the meeting time and reducing the intimacy. Such behaviors will change the strengths of nodes and the weights of links, which can be seen as an adaptive weight network. Further, it was found that the infectivity exponent has a stronger effect on the epidemic threshold and the epidemic prevalence than the weight exponent (Chu et al. 2011). For the other mathematical models on network with weights, please see Macdonald et al. (2005), Zhu et al. (2013) and references cited therein.

Motivated by the Liu et al. (2013), Zhu et al. (2012b, 2013), we introduce the individuals’ birth and death in our model, and set up an SIRS alcoholism model in complex network. Furthermore, we study the impact of the fixed weight and adaptive weight on the spread of alcoholism. We not only introduce general forms of the weight function to account for different cases of transmission but also to analyze the influence of weights on alcoholism spreading. In addition, we add the group of recuperator and study the relapse of the recuperator. We also give some control strategies against drinking, our results show that the treatment of recuperator for stopping relapsing and preventing the susceptible people to drink are two effective control strategy, and the latter has more effective than the former when the proportion of recuperator to accept treatment is equal to the proportion of susceptible people to refuse drinking alcohol.

The paper is organized as follows: in “Model formulation” section, we set up the model via differential equations. Then we present a global analysis of the model in “Global dynamics of the model” section. In “Control strategy” section, we perform two control strategies. In “Sensitivity analysis and numerical simulations” section, we perform sensitivity analysis and numerical simulations. We finally conclude the paper and give some measures to control alcoholism in “Conclusion and discussions” section.

Model formulation

In epidemic model, the total population N generally is divided into susceptible, represented by S, infections, represented by I and recovery, represented by R. An SIRS model, susceptible people usually infected by infections and become infected individual, after infection acquired immunity, it will be recovery. It allows members of the recovered class to be free of infection and rejoin the susceptible class. Zhu et al. (2012b) introduced an SEIRS epidemic model with the incubation period on complex network. Zhu et al. (2013) proposed a modified epidemic SIS model on an adaptive and weighted contact network, they introduced the general forms of the weight function and presented a new weight called “adaptive weight”. Liu et al. (2013) presented an SIR model with individual’s birth and death on scale-free networks and analyzed the stability of three equilibria. Motivated by these work, we set up a new SIRS alcoholism model on complex network. First, we introduce the individual’s birth and death in our model, Second, we take into account adaptive weight in our model. At last, we studied the effect of alcohol relapse on the spread of the alcoholism. In our model, we divide the whole population into three compartments: the susceptible S(t), denote the people who do not drink or drink limited; the problem alcoholic I(t),  denote the people who drink more than daily and weekly limit; the recuperator R(t),  denote the people who recover from alcoholism after treatment. In order to reflect the heterogeneity of the contact network, it is necessary to consider the node with different degree. Let , and denote the densities of susceptible individuals, problem alcoholics and recuperator individuals with degree k at time t respectively. Then , and are the average densities of susceptible individuals, problem alcoholics and recuperator individuals respectively, where P(k) is the probability that a randomly chosen node has degree k.

On complex networks, as alcoholism spread in a crowd, every site of N is empty or occupied by only one individual. Just as Liu et al. (2013), we give each site a number: 0, 1, 2, 3. We interpret the four states as: state 0: vacant; state 1: a susceptible individual occupied; state 2: a problem alcoholic occupied; state 3: a recuperator individual occupied. The states of the system at time t can be described by a set of numbers Each site can change its state at a certain rate. We assume that a birth event occurs at a vacant node at rate b. A susceptible individual can be infected through contact with a problem alcoholic. While a problem alcoholic can be cured at rate or lead to relapse at rate through contact with problem alcoholics or other reasons. All individuals’ death rate is and we assume that alcoholism is not fatal. If a person is dead, the corresponding side becomes vacant. Therefore, the dynamics of , and are described by the following differential equationswith initial conditionswhereand the parameters are all positive constants. The model structure is shown in Fig. 1.

Fig. 1

Transfer diagram for alcoholism model

The meanings of the parameters and variables in model (1) and (3) are as follows

represents the new born susceptible individuals per unit time, which is proportional to vacant nodes’ birth rate b and the density of vacant nodes .

represents the new problem alcoholic individuals per unit time, which is proportional to the degree k, the density of susceptible individuals and the probability that alcoholism transmits through a link , while is the transmission rate from nodes with degree i to nodes with degree k, is the infectivity of the problem alcoholic nodes with degree i. So is the link’s average infectivity of the problem alcoholic nodes with degree i. P( i|k ) is the probability that a node of degree k connected to a node of degree i. In this paper, we focuses on degree uncorrelated networks. Hence, , where is the average degree of the network.

is the natural death rate. Since the disease of alcoholism is assumed not fatal, so there is no disease related death. We assumed that if an individual dies, the corresponding side will become empty. represents the recovery rate of the problem alcoholics. Some recuperators are likely to recur drinking. The density of relapse alcoholics is , where means the recurrence rate. represents the transfer rate from recuperator to susceptible people.

There is little literature about the network model with links’ or nodes’ weights, but the weighted patterns on complex networks have various formats. Weighted patterns are used to represent the different intensities of infection by contact. Usually, the weight between two nodes with degree i and j are measured by a function of their degrees (Barrat et al. 2004a, b, c; Macdonald et al. 2005), where and m depend on the specific network. In the Escherichia coli metabolic network in the US airport network in the scientist collaboration network Here, we use a different expression for the weight function (Zhu et al. 2013), where g(k) is an increasing function of k, because the nodes with more connections will be more influential and gain larger weights. Since estimates the links’ weight, the weight of each node can also be measured by summing up the weights of links connected to it. Thus, . On uncorrelated networks, . We assume that the node with degree i has a fixed transmission rate given by , and the transmission by the link from the i-degree node to a k-degree node is measured by the proportion of this link’s weight accounting for the k-degree nodes’ weight (Chu et al. 2011). So, we have

In this paper, we also consider people’s health-conscious behavior, so the value of weight function will become less and less as the alcoholism progresses. In particular, if a person has more neighbors, it will be more cautious, therefore the weight will decrease more obviously. Thus, the weight function can be expressed as , where h(k) is an increase function of k. The corresponding becomeswhere .

Substituting (4) and (5) into (3) respectively, we get two respectivelyandit is clear that when , (7) reduces to (6).

Then, substituting (6) into (1), we obtain the fixed weight systemand substituting (7) into (1), we obtain the adaptive weight systemwhere in (8) and (9).

Let be the density of the whole individuals with degree k, . Then adding the three equations in (8) or (9) gives

By Eq. (10), we get that , where represents the initial density of whole population with degree k. Hence, , then for all .

Due to the limit system and original system have the same dynamic behaviors for a long time. And at steady-state, it is sufficient to study the limiting systemsand

It is easy to obtain that and for . So the region is the positive invariant for both (11) and (12).

Global dynamics of the model

The basic reproduction number R0

Here, we first calculate the fix weight model’s basic reproduction number. Using the next generation method in Driessche and Watmough (2002), it is clear that model (11) has an alcohol free equilibrium . System (11) can be written asandwhere the rate of appearance of new infections isand the rate of transfer of individuals out of compartments is

The Jacobian matrices of and at the alcohol free equilibrium arewhereand E is identity matrix, 0 is zero matrix. It is clear that V is a nonsingular M-matrix and F is a nonnegative matrix. According to the concept of next generation matrix and reproduction number given in Driessche and Watmough (2002), the reproduction number of (11) equals towhere .

Coincidentally, we get that matrices F and V in model (12) are the same as that in model (11). Therefore, the reproduction number of model (12) is also given by (15), which implies that the adaptive weights cannot change the propagation threshold.

According to the above process and Theorem 2 in Driessche and Watmough (2002), we obtain the following results.

Theorem 1

For the two alcoholism models (11) and (12), we have

Both of their basic reproductive number are equal toin (15).

If, the alcohol free equilibriumof (11) and (12) is locally asymptotically stable, but unstable if, whereis defined by (15).

Next, we will investigate the global stability of the alcohol free equilibrium and the globally attractive of the alcoholism equilibrium of model (11). Since the stability of the equilibria in model (12) are difficult to demonstrate, so we only give some numerical simulations to discuss it at “Sensitivity analysis and numerical simulations” section.

Uniqueness of the alcoholism equilibrium

We first give the following Lemma which guarantee that the density of population in each compartment cannot become negative or greater than .

Let , we study the system for .

Lemma 1

The setis the positively invariant for system (11).

Proof

We will show that if , then for all . Denotewhere . Let the ‘outer normals’ be denoted by and . We use the Nagumo’s result in Yorke (1967). Since is a 2n-dimensional rectangle. From (11), it is easy to obtain that for .

So, any solution that starts in stays inside .

Furthermore, we will ascertain the uniqueness of the alcoholism equilibrium. We give the following theorem.

Theorem 2

There exists a unique alcoholism equilibriumin model (11) when.

The alcoholism equilibrium of system (11) is determined by equationsa direct calculation yields

Then we get a self-consistency equation of as follows

Obviously, satisfies above equation, then , which is the alcohol free equilibrium of (11). We transform the self-consistency equation form as , wheresinceandthe equation has a unique non-trivial solution if and only if then , through which the alcoholism equilibrium is admitted. The proof is completed.

Global stability of the alcohol free equilibrium

Here, we use the method in Lajmanovich and Yorke (1976), d’Onofrio (2008) to demonstrate the global behavior of the system (11). By letting . Then, equations in (11) can be rewritten as a formwhere Ay is the linear part, H( y) is the nonlinear part andwhere E is unit matrix and

Lemma 2

Letbe anmatrix, and assumewhenever. Then there exists an eigenvectorofAsuch that, and the corresponding eigenvalue isS(A) .

The stability modulusS(A) is defined by, , whereare the eigenvalues ofA.

Choose , such that , for . Then is an nonnegative matrix. Therefore, by Theorem 2.20 from Varga (2000), there exists a nonnegative eigenvector with nonnegative real eigenvalue equal to its spectral radius . So we have , where E is the unit matrix. Then, , so is also an eigenvector of A, and the corresponding eigenvalue is . If is any eigenvalue of A, then is an eigenvalue of , so , then we have , and , therefore , that is to say is the maximum real part of all eigenvalues of A, so . The proof is completed.

Lemma 3

(Lajmanovich and Yorke 1976) Consider the systemwhereAis anmatrix andN(y) is continuously differentiable in a region. AssumeThen eitheris globally asymptotically stable inC, or for anythe solutionof (17) satisfies, independent of. Moreover, there exists a constant solution of (17), .

the compact convex setis positively invariant with respect to the system (17), and;

;

there existand a (real) eigenvectorofsuch thatfor all;

for all;

is the largest positively invariant set contained in.

Theorem 3

For system (11). When, there exists an alcohol free equilibriumis globally asymptotically stable in. When, there exists an alcoholism equilibriumis permanent in, that is to say, there exists anmsatisfies.

We will confirm that the system (11) satisfies all the hypotheses of Lemma 3.

Condition (i): Lemma 3 is satisfied if we suppose that .

Condition (ii): using the mean inequality and limit rule can validate the conclusion.

Condition (iii): notice that is an matrix with whenever , then from Lemma 2, there exists an eigenvector of and the associated eigenvalue is . Let , for , , therefore for all .

Condition (iv): we know and , so it is clearly satisfied.

Condition (v): let . If , then for . But since each term of the sum is nonnegative, then we get that . If this equation is established, then or . From system (11), if then for . That is . So, the only invariant set respect to (11) contained in H is , condition (v) is satisfied. This completes the proof.

Globally attractive of the alcoholism equilibrium

Theorem 4

When, the only alcoholism equilibriumin model (11) is globally attractive in.

We define the following functions, and for , where , are continuous and the right-hand derivative exists along solutions of (11). Let be a solution of (11), we may assume that , and . For a given and for sufficiently small from (11), if , we haveor for and , we have

According to the definition of M( y( t) ), we know

Then, if , for we haveor for and and since and , we conclude that . Therefore, if , then . Similarly, we can testify that if , . And if , then . If , then . Denote

Both Q( y ) and q( y ) are continuous and non-negative for . Notice that

Let and , then and . According to the LaSalle invariant set principle, any solution of (11) starting in will approach . But if , by Theorem 3 we know that . Then we conclude that any solution y(t) of (11), such that , satisfies , so is globally attractive in .

Control strategy

Timely stopping recuperator recurrence drinking alcohol and stopping the susceptible people to drink are two important and effective ways to prevent the spread of the alcoholism. Next, we give two kinds of control strategies for the treatment of recuperator for stopping relapsing and preventing the susceptible people to drink, respectively. Due to the fixed weight model and adaptive weight model have the same basic reproductive number, so we only study the control strategies of the fixed weight model.

Proportion treatment

Let be the treatment proportion for recuperator, , then system (11) becomes

By the next generation method, the basic reproductive number of (18) iswe can find that with the increase of the treatment proportion, is smaller.

Change the form of , we get that

is the basic reproductive number of (11). So, we can see that proportion treatment to recuperator is a very effective control strategy, and the bigger the proportion of recuperator to accept treatment, the alcoholism is more difficult to outbreak.

Proportion prevention

Let be the proportion of susceptible people who understand the harm of alcoholism and refuse to drink, , then system (11) becomes

By the next generation method, we get the basic reproductive number of (19) is

It is easy to know that , and the greater the proportion of susceptible people refuse to drink, the smaller the number of alcoholics.

Next, we are going to compare which strategy is more useful. Transform the form of and , we haveand

When , we know that . That is to say, when the proportion of recuperator to accept treatment is equal to the proportion of susceptible people to refuse drinking alcohol, the strategy of in proportion to prevent susceptible people drinking alcohol will be more effective.

Sensitivity analysis and numerical simulations

In this section, we perform some sensitivity analysis on the basic reproduction number in terms of the parameters. Our simulations take the scale-free networks with degree distribution is . Let , , and , , where , and are positive constants. Considering the influence of heavy alcoholics’ relapse and the weight between individuals. We focus on simulate the relapse parameter , the weight parameter and the nodes’ infectivity parameter .

From Fig. 2, it is clear that presents growth trend with the increase of the relapse parameter . It means that bigger alcoholism recurrence rate are easy to cause outbreaks of alcoholism. Figure 2a shows that the greater the weight parameter lead to greater . It means that the greater the link’s weight between two nodes, the easier the alcoholism broke out. Figure 2b shows that increases as the nodes’ infectivity parameter increases. That is to say, if a problem alcoholic has big“infectivity“, the alcoholism is more easy to broke out.

Fig. 2

The relationship between the basic reproduction number and the parameters on scale-free networks

Next, we perform some numerical simulations to illustrate our theoretical results of the models (11) and (12), so as to find better control strategies. The parameters are used as , , , , , , , It is clear that means no weight model, , mean fixed weight model and , mean adaptive weight model.

Figure 3a, b describe nodes without weight with and 0.09 respectively, Fig. 3a shows that when , the alcoholism dies out quickly. Figure 3b shows that the problem alcoholics’ population will maintain at a positive stationary level, which implies that the alcoholism will become endemic.

Fig. 3

The densities of alcoholics with different degrees and without weight when (a) and (b)

Figure 4a, b describe nodes with fixed weight with and 0.05 respectively, Fig. 4a shows that when , the alcoholism dies out quickly, which implies that the alcohol free equilibrium of (11) is stable. Figure 4b shows that when the problem alcoholics’ population will maintain at a positive stationary level of (11), which implies that the alcoholism will become endemic. Compared with Figs. 3a and 4b, for same parameters, basic reproduction number is 0.8981 and 1.6283, respectively. So we know that of model on networks with weights is larger than that on networks without weights.

Fig. 4

The densities of alcoholics with different degrees and fix weight when (a) and (b)

Figure 5a, b describe nodes with adaptive weight. We know that (12) and (11) have the same . So we choose the same parameters as Fig. 4. From Fig. 5a, b, we can see that since the adaptivity of weight, the alcoholics rapidly drop first and experience a valley then up to a small peak or dies out. This is caused by the behavior of people’s self-protection awareness. From this point, the adaptive model is more close to the actual situation.

Fig. 5

The densities of alcoholics with different degrees and adaptive weight when (a) and (b)

Figure 6a, b shows the densities of problem alcoholic individuals with different degrees and different the adaptive coefficient We know that the stronger the adaptive coefficient the greater the self-protection awareness of susceptible. It lead to alcoholics density maintain in a lower value.

Fig. 6

The densities of alcoholics with different adaptive coefficient in a and in b

Conclusion and discussions

In this paper, we proposed a modified SIRS alcoholism model with relapse on weighted networks to study the influences of individual’s contact patterns on drinking dynamics. We construct two complex network models with fixed weight and adaptive weight, respectively. We get the alcoholic propagation threshold which determines the propagation dynamics. We also obtain the existence of equilibria of (11). For the model with fixed weight, we prove that when , the model’s alcohol free equilibrium is global stable and alcoholism will disappear, otherwise, if , the alcoholism equilibrium is global attractivity and alcoholism will persistence. For the model with adaptive weight, we only make some numerical simulations. By comparing the alcoholism models with no weight, fixed weight and adaptive weight, we have:

Fixed weight model have larger propagation threshold than no weight model;

The adaptive weight cannot change the propagation threshold, but it can induce the alcoholism to decay quickly;

Strong adaptability can inhibit the alcoholism population reach a high level;

In order to eliminate alcoholism problem, first we should try to reduce the frequency of interaction between susceptible and problem alcoholic, that is to say decreasing the link’s weight between two nodes, this may be a effective measures. Second, through education or media to spread the dangers of alcohol abuse in order to decrease the relapse it is also a very effective measures.

The delay differential equations usually exhibit much more complicated dynamics than ordinary differential equations because the time delay can lead to instability, oscillation, or bifurcation phenomena (Bianca et al. 2015; Bianca and Guerrini 2014). There is a time delay during a susceptible individual becomes the problem alcoholic, so it is more realistic to consider the time delay in the modelling alcoholism process. We can modify (1) to the following model with delay

We leave this work for the future.