Analytical and computational study of an individual-based network model for the spread of heavy drinking

ABSTRACT Two simple models for the spread of heavy drinking among a network of individuals are re-introduced and analysed. We provide theorems on the spread of alcohol abuse for these models in cases involving simple connection schemes. Indicators for this spread that resemble the used in disease assessment are suggested and studied. We further provide computations with our models on general application networks and begin to study the reliability of the spread indicators.


Introduction
Over 22 million Americans suffer from substance abuse disorders, approximately 80% involve alcohol use disorder [22]. Alcohol consumption and substance abuse are often closely tied to social behaviours. Particularly among teens and young adults, drinking is typically due to social motivations and peer influences have been shown to be predictive of substance abuse [11][12][13]. The complexity and dynamic nature of social interactions highlight the potential usefulness of modelling and network-based approaches to understanding the interactions between behaviour patterns and alcohol use.
In the last decade there have been significant advances, mainly through the availability of vast amounts of data and fast computers, in our understanding of the structure of social networks (see [2,3,8,20,25]). The first author developed an integro-differential equation model for the spread of alcohol abuse [6] based on the network studies in [5,26]. In this report, we extend the research in these papers.
A wide range of mathematical modelling approaches have been used to examine social behaviour and alcohol use, including Markov models, differential equations, and individual-based networks [1,4,7,10,16,18,19,23,24]. Several of these efforts have examined the spread of alcohol abuse between social contacts including the SIR/SIS models presented in [4,17,21] as well as the data analysis in [23] and the agent-based model in [10].
In this paper we introduce two models based entirely on those in [6] and originally [5]. They each start from a network with N nodes which represent individuals. The first uses a cubic function to model the evolution of the drinking level index v m i of individual i on step m: N i is the list of nodes connected to i, its friends, and |N i | is the number of indices listed in N i . r i ∈ [0, 1] is the resilience and measures the ith individual's resistance, with respect to drinking, from peer pressure. We think of the v m i = 0 as being a healthy individual with respect to drinking; they could be light drinkers or sober. While v m i = 1 is an individual who is generally considered unhealthy as far as his/her drinking. The connection weights are 0 ≤ σ ij ≤ 1 for j ∈ N i indicating the strength of influence of person j on person i with so that n m i forms a σ -weighted average over the states j ∈ N i . Note that the weights σ ij are not necessarily symmetric. We also have a step or iteration with ending time M so m = 0, 1, 2, . . . , M. The parameter λ might be called a step size if we view (1) as a discretization of a differential equation (DE). We often think of λ ∼ 1/M but always insist on 0 < λ < 1.
The inspiration for the cubic model comes primarily from viewing it as an approximation to a DE where 0 and 1 are stable steady states and the r i is unstable. Thus, when individual i's friends are heavy drinkers (n m i > r i ) this individual tends toward unhealthy alcohol consumption (v m i → 1) and when the friends are healthier (n m i < r i ) individual i tends toward sobriety (v m i → 0). The second model we consider, referred to as the switch model, is v m+1 where H = H(s) is the Heaviside function, To see that (3) provides a reasonable model, note that on the step from m to m+1, we have v m+1 Thus, if n m i ≥ r i , then v m+1 i will be bigger than v m i and the opposite will occur if n m i < r i , This approach was introduced in [6] and, due to its similarity to certain neuronal network models, was especially easy to analyse. Initial configurations at m = 0 are provided for each of the schemes; We will prove in the next section that the v m i ∈ [0, 1]. The update formulas are set so the weighted influence of friends n m i when compared to the individual i 's resilience r i will give an appropriate adjustment to their drinking level v m i . That is, A common theme when examining transmission (whether behaviours, information or disease) is the evaluation of a bifurcation parameter, such as the basic reproduction number, which predicts whether the disease/information will spread or die out. Individual-based models often follow a natural threshold which suggests that it may be possible to find similarly predictable quantities in our case as well. A focus of this paper is on the accuracy of the following quantities: SF stands for scale free. More precisely, does χ 0 > 1 or χ SF 0 > 1 provide a good predictor on the spread of unhealthy drinking? Here, χ 0 is a ratio of average of the individual's initial drinking levels with the average of the individual's resiliences. The number χ SF 0 takes into account the connection scheme and thus each individual's influence.
The outline of this paper is as follows. In Section 2 we show for each model that the v m i ∈ [0, 1] for all i and m. We also prove theorems about the limiting behaviours in m for each model in the cases of all-to-all, star and leader-clique networks. In Section 3 we display our results on a variety of simulations where we look at generalizations of the simple networks studied in Section 2. We are focussed on the accuracy of χ 0 and χ SF 0 in this section. In Section 4 we experiment with our prediction scheme on two application networks. We will add some conclusions in Section 5.

Model properties and threshold criteria
In this section we establish some basic properties of the methods as well as theorems on the thresholds for the spread of heavy drinking on simple networks.
The proof is by induction, we assume that v m i ∈ [0, 1] for all i but m is fixed. From this we have, n m i ∈ [0, 1]. Recalling that λ ≤ 1 and that all the r i ∈ [0, 1] give and so, v m+1 i ≤ 1 (as well as being greater than or equal to −1). Now, letting u m i = 1 − v m i we can show, using (2), and then use the same proof to conclude that |u m+1 i | ≤ 1 and thus |1 − v m+1 i | ≤ 1. This verifies that v m+1 i ≥ 0 and finishes the analysis for the cubic model case (1). For the switch model (3); we again use induction and iterate from the mth time level with the assumption that v m i ∈ [0, 1] for all i but fixed m. Again, we will show this in the m+1 case. Since v m i ≥ 0 and H(·) ≥ 0 it follows that, since all the terms on the right side of (3) are nonnegative, v m+1 This shows that 0 ≤ v m+1 i ≤ 1 for the switch model (3). Analysis of simple networks: We now introduce a set of simple networks ( Figure 1) that we see as possible prototypes for more complex ones. We study theoretically an all-toall network (Figure 1(A)) that we see as a place-holder for random networks. The star network (Figure 1(B)) provides some clues on the behaviour of scale-free systems and the leader-clique a further generalization of both the all-to-all and the star network. The idea of examining small networks is inspired by Golubitsky and Stewart [9].
In the star network ( Figure 1(B)) there is one leader or hub that is connected to each of the followers and thus has a significant influence. We see this as similar to a scale-free network where there tends to be a few nodes with high degrees of connectivity. We will see that it is important to take into account the numbers of connections in order to predict the spread of behaviour for both the star and eventually the scale-free networks.
In the leader clique network (Figure 1(C)), by varying the weight of the leader's connections to the follower clique members, we will be able to transition from an all-to-all situation to the star network.
All-to-all network: We now prove a theorem about the spread of heavy drinking in the simple case where the network is all-to-all as in Figure 1(A) with self-connections, uniform connection weights and uniform resiliences.
Then, for the models (1) and (3) with initial condition (4), Proof: For the cubic model, we have from (1) Note that since r i = r for all i, χ 0 =v 0 /r. When χ 0 > 1 we havev 0 − r > 0 and thus the sequence v m i will be increasing in m for i = 1, . . . , N. Since the sequence is bounded above by 1 it will converge to some limit between v 0 i and 1. But, since this sequence would increase beyond any limit in the open interval For the switch model, we sum equation (3) over all i, subtract r from each side and divide by N to find thatv Thus, ifv 0 − r > 0 thenv 1 − r > 0 and, by induction,v m − r > 0 for m = 1, 2, . . ..

Knowing this guarantees that
where we computed the geometric series sum on the last step. So then, since In the case when χ 0 < 1 ⇒v 0 < r we can argue thatv m+1 ≤v m and thus finishing the proof for the switch model case.
Remark: Theorem 2.2 shows χ 0 's relation to 1 is an exact predictor of the spread of unhealthy drinking in the all-to-all network with self-connections, uniform connection weights and constant resilience.
Star network: We now analyse the situation in Figure 1(B) which provides a highly reduced example of a network with a connection scheme that, like a scale free network, has a small number of nodes which have significantly higher degrees than the other.
The star network in the case of the cubic model would be v m+1 where the first equation is for a hub or super-connector with index H and the second is for each of the N F follower nodes numbered 1 through N F with constant resilience r F . The definition of the switch model on the star network is v m+1 where, as with the cubic model, the first equation is for a hub or super-connector with index H and the second set of equations is for each of the N F follower nodes numbered 1 through N F . We are now in a position to prove a theorem about this reduced model. This result will show that if the hub is a heavy drinker with low resilience and the followers are, on average, moderate drinkers, then as time goes on, all the members of the network will become alcohol abusers.

Theorem 2.3:
If drinking behaviour is governed by either the cubic or switch models above on the star network, v 0 Proof: First, we show this for the cubic method. Averaging the second equation (7) Since We now show this for the switch network; averaging over the follower nodes we havē v m+1 Using an induction argument beginning with our assumptions at As in the previous proofs we have that {v m F } is an increasing sequence bounded by 1. Again, we can argue that its limit must be 1 by showing that any other limit betweenv 0 F and 1 would lead to a contradiction.
In the same way, we can prove v m+1 and then, since {v m H } is increasing we must have v m H → 1 as m → ∞.
Remark: Of course, the network structure is crucial here. Suppose r F >v 0 Notations for averages for the entire star network are as follows: So, we have These results show that χ 0 is not a good predictor of the emergence of heavy drinking in this (star) network.
At the same time, for this star network, we have If we continue to suppose r H <v 0 F < r F , v 0 H > r F and N F large (So all v's → 1 and χ 0 < 1.), then So χ SF 0 > 1 is consistent with the spread of heavy drinking in this situation.
Leader and clique network: We now consider a leader and clique model with the connection scheme as in Figure 1(C). Here, the network equations for the cubic model would look like the following: The influence function n m LC involves a weighted sum with parameter β L > 0; For the switch model we have v m+1 Note that the connection weights are In this case we adjust the influence of the leader through the parameter β L . If β L ≥ 1 is large, the leader has a strong influence and the network becomes somewhat similar to the star network except the leader is influenced only by himself. While if β L = 1 then the leader's influence is the same as the other members of the clique and n m LC is just the average over all the nodes including the leader If, further, r L = r C then this is just the all-to-all network studied earlier.
Note that here our quantities χ 0 and χ SF 0 may not good indicators of behaviour since we have manipulated the weights with this β L parameter.
An interesting situation is where the leader is a heavy drinker with low resilience and either the clique has moderate drinkers or the leader has a significant influence (β L large). The theorem below shows that all members will become alcohol abusers in these situations.

Theorem 2.4:
then for either the cubic or switch models we have v m L → 1 and v m i → 1 for i = 1, . . . , N C as m → ∞.
Proof: As in the other proofs, the key is to show that all members have an increasing drinking index. From the monotonicity and the fact that the limits must end up at a fixed point (in this case value 1) we can complete the proof. Ifv 0 C > r C then, with our assumption v 0 L > r C , So, we have n 0 LC > r C > r L . Thus {v m L } and each of the {v m i } will be increasing sequences and they all must tend to 1.
If, instead, the clique members are healthy at the start sov 0 C < r C ; we need to consider n 0 LC more carefully. Since r C −v 0 C is positive, we use the right value in the 'Max' in Equation (11). So, we have LC > r C > r L . This statement shows the sequences will again be increasing and finishes the proof. To investigate the leader clique network further, we carried out a simple computation with the cubic algorithm on the network in Figure 1(C) (with N C = 9 nodes instead of the 6 nodes displayed in Figure 1(C)) and the weighting scheme (10). For β L values ranging from 5 to 50 we ran the cubic algorithm for M = 150 steps with r C = 0.65, r L = 0.25, λ = 0.20, v 0 L = 0.75 and the v 0 i were chosen randomly on [0, 1/2] and their mean wasv 0 C = 0.323. We examine the longer time behaviour of the average drinking level of the whole network (including the leader)v M . Figure 2(A) shows the final outcomes,v M as the leader's influence is increased in the case when the leader is already a heavy drinker (0.75 = v 0 L > r L = 0.25) while all members of the clique are healthy (v 0 i ≤ 1/2 < r C = 0.65). Of course, heavy drinking ensues as the leader's influence is increased.
If we knew the connections in our network (the knowledge we assume for χ 0 and χ SF 0 ) as well as the influence which includes the weights of the connections, we could compute the following parameter:  whose values, as they depend on β L , are plotted in Figure 2(B). Observe that the point (in β L ) where χ LC 0 = 1 (where β L ∼ = 18) is an accurate predictor of the emergence of heavy drinking (Figure 2(B) -members v i -index getting close to 1). 5

Simulations
In this section we extend the χ 0 analysis developed in Theorems 2.2-2.3 to general graphs using computer experimentation. We have developed MATLAB codes that create either random or scale-free networks with N nodes and then implement either the cubic or switch network models for M steps. For these simulations a symmetric adjacency matrix W was generated where each entry w ij is 0 or 1 depending on whether there is a connection between nodes (individuals) i and j. Thus, N i is just the nonzeroes in row i, |N i | is the number of nonzeroes in that row and the σ ij = 1/|N i | if w ij = 1. Figure 3 has our results for many runs on a fixed random network with the cubic model in order to assess how well the χ 0 -criterion derived in Theorem 2.2 predicts actual behaviour. The initial conditions and resiliences are randomly chosen for each run and uniformly distributed on [0, 1]. We find that χ 0 =v 0 /r is well correlated with the final averages for the indexv at end step M in the case of this random network (correlation r Corr ∼ = 0.93). In Figure 4 we display results from computer runs on a scale-free network with the cubic model. Here, we find that χ 0 =v 0 /r is not a good predictor of spread (correlation r Corr ∼ = 0.77). However, if we use χ SF 0 then we typically find a better correlation (r Corr = 0.93).
In Figure 5 we continue the investigations whose results we showed in Figures 3 and 4. In Figure 5(A), we show how the correlation between χ 0 and spread is enhanced as more and more connections are added to the random network.
In Figures 5(B) and (C) we look at histograms of χ 0 and χ SF 0 for scale-free networks that enhance the results observed in Figure 4. In each of these cases we did 15 experiments where each consisted of 25 runs on a particular fixed scale-free network. The results reinforced the predictability of χ SF 0 for such scale-free networks. Figure 4(B) has the correlation of χ 0 andv M as 0.737, while Figure 4(C) has the correlation of χ SF 0 andv M as 0.908.
In Figure 6 we display results from computer runs on a scale-free network with the switch model. Here, we find that χ 0 is again not a good predictor of spread (correlation r Corr ∼ = 0.47). However, if we use χ SF 0 then we typically find a better correlation r Corr = 0.78. Unfortunately, the correlation between χ SF 0 is not as solid with the switch model as it was with the cubic.

Applications -karate club and multiple clique networks
In this section we investigate the accuracy of χ 0 and χ SF 0 for our cubic model on the spread of heavy drinking for two example networks.
In the karate club [27] there are two leaders who are not on friendly terms. In Figure 7 we display the network and its degree distribution. In Figure 8 we provide our results from applying our cubic model on this application network. The results in Figure 8 where the correlation to χ SF 0 is much better than to χ 0 was typical of our computations. The karate network has two hubs which suggest it is closer to a scale-free network than a random network. The computations show χ SF 0 is a better predictor of spread than χ 0 which is consistent with our new ideas. The correlation between predicted χ SF 0 values and outcomes was 0.89 while it was only 0.61 for χ 0 .
We also considered a multiple clique network as displayed in Figure 9(A). In this case there are several individuals who know members outside their own clique as well as others who only have friends (connections) within their clique. In our computations we found that both the χ 0 and χ SF 0 quantities are good predictors for the multiple clique network (see Figures 9 and 10). Here, the correlation between final outcomes for drinkingv and χ 0 was 0.86 and was 0.81 for χ SF 0 .

Discussion
In this paper we have examined two related simple difference models (cubic and switch) for the spread of alcohol abuse on networks. Each node or individual i in the network is identified by their drinking index v m i ∈ [0, 1] at each time step m. Inspired by the idea of reproduction number or R 0 in disease spread, we investigated three parameters χ 0 , χ SF 0 and χ LC 0 that depended on initial members drinking levels, resiliences and, for the latter two, connections and influence. We found that χ 0 was a good predictor of behaviour for the all-to-all and random networks. The χ SF 0 and χ LC 0 were better predictors for the behaviour of networks where the influence of individual nodes can be vastly different either through connections or weights.
In the future we intend to investigate treatment programmes such as alcoholics anonymous with our models. We find the question of finding optimum interventions for given M−Clique Network:  networks and initial connection schemes very appealing and would look to [14,15] for inspiration.