Suppressing epidemic spreading in multiplex networks with social-support

Although suppressing the spread of a disease is usually achieved by investing in public resources, in the real world only a small percentage of the population have access to government assistance when there is an outbreak, and most must rely on resources from family or friends. We study the dynamics of disease spreading in social-contact multiplex networks when the recovery of infected nodes depends on resources from healthy neighbors in the social layer. We investigate how degree heterogeneity affects the spreading dynamics. Using theoretical analysis and simulations we find that degree heterogeneity promotes disease spreading. The phase transition of the infected density is hybrid and increases smoothly from zero to a finite small value at the first invasion threshold and then suddenly jumps at the second invasion threshold. We also find a hysteresis loop in the transition of the infected density. We further investigate how an overlap in the edges between two layers affects the spreading dynamics. We find that when the amount of overlap is smaller than a critical value the phase transition is hybrid and there is a hysteresis loop, otherwise the phase transition is continuous and the hysteresis loop vanishes. In addition, the edge overlap allows an epidemic outbreak when the transmission rate is below the first invasion threshold, but suppresses any explosive transition when the transmission rate is above the first invasion threshold.


Introduction
An outbreak of such diseases as SARS [1] and H5N1 [2,3] puts at risk the lives of countless people. During the first nine months of the recent Ebola epidemic there were 4507 confirmed or probable cases of infection and 2296 deaths [4]. Increasing the investment of public resources to control a disease pandemic can be a serious economic burden, especially in developing countries [5,6]. Many researches have been done on how to optimize scarce public health care and immunization resources when attempting to control an epidemic [7,8,9,10], the goal being to minimize the number of infected individuals by determining that optimal allocation [11].
A complex network science approach is now being widely used to determine the impact of resource investment on spreading dynamics. Böttcher et al. [12] studied the impact of resource constraints on epidemic outbreaks and found that when the resources generated by the healthy population cannot cover the costs of healing the infected population the epidemics go out of control and discontinuous transitions [13,14,15,16,17] occur. Chen et al. [18] explored the critical influence of resource expenditure on constraining epidemic spreading in networks and found that public resources can affect the stability of the disease outbreak. At a certain disease transmission rate there is a critical resource level above which a discontinuous phase transition in the infected population occurs. Böttcher et al. [19] assumed that only the central nodes in a network can provide the necessary care resource, and they found that a discontinuous transition in infected nodes occurs when the central nodes are surrounded by infected nodes. All of these researches focus on how public resource investment affects the spread of disease.
In real-world scenarios only a small percentage of patients are assisted by public resources. The majority depend on help from family and friends who provide economic [20,21,22] and emotional support [23,24]. We thus study how social support from family and friends affects the dynamics of disease spreading. In a social network, a node has different connections in different settings. We can thus regard friendship ties (virtual contacts) and co-worker ties (physical contacts) as two different network layers. Although economic and medical resources and sources of information usually propagate through social relationships, diseases usually propagate through physical contacts. Thus we use a multiplex network of two-layers [25,26,27,28] to study how resource allocation in the social layer affects the spreading dynamics in the contact layer.
We use the susceptible-infected-susceptible (SIS) model in a multiplex network of twolayers to mimic the coupling dynamics between disease spreading and resource support. The disease propagates through the layer of physical contacts, but infected nodes seek help from their neighbors through the layer of social relations. Infected nodes receive resources from healthy neighbors and do not generate resources. We analyze the process using a dynamic message passing (DMP) approach [29,30,31,32]. We examine how degree heterogeneity affects the dynamical process and find that the infected density in the steady state (ρ) increases continuously at the first epidemic threshold and then jumps suddenly at the second threshold. Hysteresis loops exist in the phase transition of the infected density, and the size of the hysteresis region and the value of the invasion threshold decrease with the degree heterogeneity. Examining how edge overlap between the two layers affects the dynamics of spreading we find that the overlap has a critical value. When the overlap is below the critical value, the infected density first increases continuously and then discontinuously with disease transmission rate, and there are hysteresis loops. When the overlap is above the critical value, the phase transition of ρ is continuous and there is no hysteresis loop. We also find that when the transmission rate is below the first invasion threshold the disease outbreaks more easily for a large edge overlap, but when the transmission rate is above the first invasion threshold the edge overlap suppresses the disease spreading and the second invasion threshold increases as the overlap increases.

Epidemic model with social-support
In a multiplex network of two-layers, each layer has N nodes and each node in the first layer has a counterpart in the second layer. Here the upper layer is the social relationship network (e.g., Facebook friends and family members) from which healthy nodes allocate resources to infected neighbors (see layer S in Fig. 1). The lower layer is the physical contact network through which the disease spreads (see layer C in Fig. 1). Variables A and B are the adjacency matrices of layer S and C with elements a ij and b ij . If nodes i and j are connected by one edge in layer S, a ij = 1, otherwise a ij = 0. The same is true in layer C. We denote by s υ the node state variable of node υ, and if it is in the susceptible state s υ = 0, otherwise s υ = 1. We assume that each healthy individual has a certain resource level r per unit time, which for simplicity we set at r = 1. Resources are distributed equally to infected neighbors. Figure 1 shows that node X distributes one resource unit to three infected neighbors in layer S, and that node Y distributes one resource unit to one infected neighbor. For the sake of analytical tractability, we assume that the total resource is not cumulative in the system, and if healthy nodes do not allocate their resources to neighbors they consume these resources themselves. In addition, infected nodes consume all of the received resources at the current time step, and each healthy individual generates a new one-unit resource at the next time step. Using this definition, the resources that node j gives to node i in layer S is Without resource support a node recovers spontaneously at a rate µ 0 [35], and for simplicity we assume µ 0 = 0. The recovery rate of i at time t is where µ i (t) ≡ µ(R i (t)), and R i (t) is the expected resources that node i receives from healthy neighbors. The µ r value is the coefficient that represents the efficiency of resource support from neighbors, µ r ∈ [0, 1], and k S i is the degree of i in layer S. The recovery rate of infected nodes is assumed to be positively related to the resource received from healthy neighbors in layer S. In real-world setting the cost of repairing a vital node in a complex system is much higher than the cost of repairing a common node. For example, because hub airports in airline networks play a vital role in connecting a large number of countries and regions, the repairing cost when they fail is much higher than that for lower-degree airports [33]. Similarly, the cost of repairing hub nodes in brain networks is much higher than the cost of repairing common nodes [34]. The same is true in epidemic spreading. Individuals exposed to viruses over a long period of time, e.g., medical staff members who are in constant contact with infected individuals, have large degrees in physical contact networks. Community leaders are also hub nodes in high-degree physical contact networks. In both cases the cost of curing these hub nodes being infected is much higher than other infected nodes in the contact networks. Thus we assume that the recovery rate of an infected node is negatively related to its degree.
We use the classical SIS model to investigate the spreading process in multiplex networks. Each individual can be either infected or susceptible. Susceptible individuals are healthy and are then infected by an infected neighbor at a rate β. Infected individuals recover at a rate µ i (t), which is assumed to be independent of the availability of social resources in previous researches [36,37].

Dynamic message-passing method
We use dynamic message-passing method to analyze the spreading dynamics. In this method a variable "message" passes through the directed edges of the network and does not backtrack to the source node. Our message is θ j→i , the probability that node j is infected by its neighbors other than i. In addition, ρ i (t) is the probability that node i is in the infected state at time t. The probability that an infected node i will connect to a healthy node j in layer S is a ij (1−θ j→i (t)), and the expected number of infected neighbors of node j is ℓ =i a jℓ θ ℓ→j (t) + 1, where the plus one takes into account that node i is infected. Thus the resource R i (t) that node i receives from healthy neighbors is Using this definition, the discrete-time version of evolution of ρ i (t) [38] is where ∆t is the time increment, which we set at ∆t = 1, and q i (t) is the probability that i is not infected by any neighbor in layer C, which is given by where N C i is the neighbor set of i in layer C. Note that to exclude any contribution of node i to the infection of j, we adopt θ j→i (t) instead of ρ j (t) in Eq. (5). Similarly, the discrete-time version of evolution of θ j→i (t) is Here (1 − φ j→i (t)) is the probability that j is infected by at least one neighbor other than i.
Here N C j \ i is the neighbor set of j excluding i, and the fraction of infected nodes at time t is where ρ i (∞) ≡ ρ i and µ i (t) ≡ µ i at the steady state t → ∞. Solving Eqs. (4) and (6) at the stationary state and we obtain the phase diagram of the model. We use iteration to numerically compute the evolution of the state of network nodes. Due to nonlinearities in Eqs. (3)-(7) they do not have a closed analytic form, and this disallows obtaining the epidemic threshold β c . If β > β c , ρ > 0, otherwise ρ = 0 in the steady state. When β → β c , ρ i → 0, θ j→i → 0, and the number of infected neighbors of each healthy node in layer S is approximately zero in the thermodynamic limit, prior to reaching the epidemic threshold (1 − θ j→i ) → 1. If we add these assumptions to Eq. (3) resource R i becomes R i → k S i , we will obtain the recovery rate µ i → µ r in the steady state [see Figs. 4(a) and 7(a)].
To solve Eq. (14) we define a 2E × 2E matrix J, where E is the number of edges and the elements of J are The system enters a global epidemic region in which the epidemic grows exponentially when the largest eigenvalue of J is greater than zero [31,37,32]. Thus we can obtain the epidemic threshold as where Λ J is the largest eigenvalue of J.

Numerical and simulation results
To examine how resource support affects epidemic dynamics, we perform numerical computations and stochastic simulations in the networks. Because many real-world complex networks have a highly skewed degree distribution, e.g., Facebook [40] and the World Wide Web [41], we focus on networks with a heterogenous degree distribution. We assume that the two layers of the network have the same degree sequences (k S i = k C i ). Thus for simplicity we denote k i to be the degree of node i in both layers S and C.
To build our multiplex network we use an uncorrelated configuration model (UCM) [42] with a given degree distribution P (k) ∼ k −γ in which γ is the degree exponent. Here a smaller γ implies a more heterogeneous degree distribution. The maximum degree is determined by the structural cut-off k max ∼ √ N [43] and we set the minimum degree at k min = 3. In addition we disallow multiple and self-connections and set the network size as N = 10000. When studying resource support from neighbors, we eliminate any possibility of spontaneous recovery, i.e., µ 0 = 0, and assume that node recovery is solely dependent on the amount of resources received. Here we set the efficiency parameter at µ r = 0.6 and the µ r value does not affect the result [44,37].
To determine the epidemic threshold, we use a susceptibility measure [45,46] where . . . is the ensemble averaging, and χ exhibits peaks at the transition points. We now examine how degree heterogeneity and edge overlap between the two layers of the network affect its dynamic features.

Effects of degree heterogeneity
To investigate how degree heterogeneity affects spreading dynamics, we disallow any edge overlap between the two layers, i.e., nodes are randomly connected by edges in layer S and layer C, and the amount of edge overlap m e is approximately 0 in the thermodynamic limit.
To examine ρ as a function of β, we randomly select one percent of the nodes to be seeds (ρ(0) = 0.01). Figure 2(a) shows the epidemic spreading for γ = 2.4 and γ = 3.2. Note the hybrid phase transition in ρ that exhibits properties of both continuous and discontinuous phase transitions. As β increases ρ grows continuously at β I inv .
Then an infinitely small increase in β induces an sudden jump of ρ at β II inv , where β I inv and β II inv are the first and second invasion thresholds. The ρ transition type indicates that there are three possible system states, (i) completely healthy, (ii) partially infected, and (iii) completely infected. This differs significantly from the classical SIS model. In addition, we find hysteresis loops in the phase transition of ρ when γ = 2.4 and γ = 3.2 [see Fig. 2(a)]. When the seed density is initially low, e.g., ρ(0) = 0.01, the disease breaks out at the invasion threshold β I inv , but when it is initially high, e.g., ρ(0) = 0.9, the disease breaks out at the persistence threshold β per . The arrows in Fig. 2(a) indicate the direction of the hysteresis loops. We determine critical points β I inv and β II inv and persistence threshold β per using the susceptibility χ shown in Fig. 2(b). The theoretical results obtained from the numerical iterations agree with the simulation results [see the lines in Fig. 2(a)].
We next determine how degree heterogeneity (i.e., parameter γ) influences the spreading dynamics. Figure 3(a) shows the two-parameter (β, γ) phase diagram. The parameter space is partitioned into three regions according to ρ value. When β < β I inv , the system falls into the no-epidemic regime, i.e.,the green and part of the purple area below β I inv . When β I inv ≤ β < β II inv , it falls into the low-epidemic regime (bounded by two critical lines) in which ρ increases slowly with β. Finally, above β II inv , ρ suddenly jumps to the high epidemic regime (red) in which approximately all nodes are infected. The regime between invasion threshold β II inv and persistence threshold β per is the hysteresis region (purple). The values of β I inv and β II inv both increase with γ. Although we can obtain the theoretical value of β I inv from Eq. (16), we cannot obtain the theoretical value of β II inv and β per by linearizing the equations around ρ i → 0 and θ j→i → 0 and thus we must apply numerical methods using Eqs. (4) and (6). We first define a judgment value ǫ that is linear with system size N. Without loss of generality we set ǫ = 0.3. We then define the jump size ∆ρ to be where ∆β is an infinitesimal increment in β, which we set at ∆β = 0.001, and ρ(β) is the infected density in the steady state when the transmission rate is β. We obtain the threshold The high epidemic region with large value of ρ (denoted by red color) in steady state, the no epidemic region with zero value of ρ (denoted by green color), and the low epidemic region with small value of ρ (part of the purple region bounded by the two critical lines). The hysteresis region (denoted by purple color) is bounded within β II inv and β per (denoted by red squares). The two invasion thresholds β I inv (denoted by lower blue circles), β II inv (denoted by upper blue circles) and persistence thresholds β per are determined by the susceptibility measure χ. Theoretical results obtained from the DMP method are denoted by dotted lines in the figure. (b) The thresholds interval β II inv − β I inv is plotted as a function of system size N for three different values of γ: γ = 2.0 (red triangles), γ = 2.2 (blue circles), and γ = 2.8 (dark grey squares). Error bars are smaller than the symbols used for the data points.
when ∆ρ ≥ ǫ is at a certain β value in the thermodynamic limit [47,16]. Using the numerical method, we obtain the second invasion threshold β II inv and the persistence threshold β per . Figure 3 shows that the theoretical values marked by dotted lines agree with the simulation results. The change in the system state among the three regions indicates that the phase transitions of ρ are hybrid. Figure 3(a) shows that the low epidemic and hysteresis regions expand as γ increases.
To demonstrate that there are two invasion thresholds in networks with heterogeneous degree distribution, we use a finite-size scaling analysis [48]. Figure 3 as a function of N for γ = 2.0, γ = 2.2, and γ = 2.8, where • is the norm operator. Figure 3(b) shows the values of β II inv − β I inv converging asymptotically to positive constant values in the thermodynamic , which implies the two invasion thresholds do not merge when γ ≤ 2.2 and the two are always present in networks with a heterogeneous degree distribution.
To analyze the sudden jump of ρ and the hysteresis loops, we examine the transmission process analytically using mean-field approximation in random regular networks (RRNs), which corresponds to the limit γ → ∞. Through a bifurcation analysis we account for the existence of the sudden jump of ρ and the hysteresis loops (see appendix information). Note that the first threshold β I inv disappears in the RRNs and the transition of ρ is discontinuous when it is not hybrid [see Fig. 9(a)].
To explain the hybrid transition when γ is finite, i.e., when γ ≤ 3.2, we investigate the number of susceptible neighbors around each infected node in layer S and their recovery rates as a function of β. In the steady state the number of each infected node's susceptible neighbors in layer S is n s i and their fraction n s i /k i . Here the recovery rate is µ i . To evaluate the collective state, we examine the average quantity n s /k of n s i /k i and the average quantity µ of the recovery rate. Figure 4(a) shows plots of n s /k and µ as functions of β for γ = 2.8. We find that both n s /k and µ are constant when β < β I inv , which implies zero values for ρ. They then slowly decrease until they reach the β II inv , at which point an infinitesimal increase in β causes a jump in n s /k and µ . Figures 4(b)-4(d) show the time dependence near β I inv and β II inv . Figure 4(b) shows the time evolution of the infected density ρ(t) around β I inv ≃ 0.023 for γ = 2.8. The difference in ρ(∞) for β just below and above threshold β I inv is ∆ρ [see Eq. (18)]. Note that ρ(∞) increases slowly at β I inv , i.e., a small increment ∆ρ ≃ 0.022. We next examine the time evolutions of the average resources of the infected nodes R(t) and the hub nodes R h (t) . Note that without loss of generality we can assign hub node status to nodes with a degree larger than k = 30. Note also that when β is just below β I inv both R(t) and R h (t) increase until t = t * , which implies that all infected nodes have acquired sufficient resources to recover and ρ(t) drops to zero. In contrast, when β is just above β I inv and the transmission rate is low, the promotion effect of the hub nodes allows the disease to spread on a finite scale and infected nodes have sufficient resources for recovery. Thus infection and recovery processes are balanced, and the values of R(t) and R h (t) fluctuate around a finite value when t → t ∞ [see Fig. 4(b)]. As β smoothly increases   Figures 4(c) and 4(d) show a critical time t * ≃ 220 at which β is approximately β II inv ≃ 0.033. At the early stage of the propagation process, i.e., when t < t * , the disease spreads through the local seed nodes. Because most of neighbors of the infected nodes in layer S remain healthy, they have a sufficient resource level to recover. Here the infection and recovery processes are balanced. As the ρ(t) value increases slowly the available resources levels R h (t) and R(t) for β ≃ β II inv slowly decrease [see Figs. 4(c) and 4(d)]. When β < β II inv the infection and recovery processes remain balanced when t → ∞, thus the density of infection fluctuates around a small finite value when t → t ∞ (ρ(∞) ≃ 0.18) [see Fig. 4(c)]. Because infection and recovery processes are balanced in β I inv ≤ β < β II inv , the value of ρ increases slowly. Note that as hub nodes disappear in the RRNs the disease is suppressed until β reaches a threshold at which point it jumps discontinuously, the balance disappears (see Appendix), and only one threshold remains. When β > β II inv the transmission rate is relatively large and the balance between infection and recovery is broken. Infecting the healthy nodes in layer C decreases the resources available to the nodes in layer S and delays the recovery of infected nodes. This recovery delay increases the effective transmission probability in layer C, more healthy nodes are infected, and both the available resources and the recovery rate decrease. This causes a cascading infection in system nodes that is accelerated when hub nodes are surrounded by infected nodes, and this can cause total system failure. Figure 4(d) shows an abrupt drop of R h (t) and R(t) at t * when β > β II inv . Figure 4(c) shows a rapid increase in the density of infection from a small value ρ(t * ) ≃ 0.18 to a high value ρ(∞) ≃ 1.0. Note that in the steady state the large difference ∆ρ between β < β II inv and β > β II inv causes explosive transitions. This explains the hybrid transition in networks with a heterogeneous degree distribution. Figure 4(d) shows the evolution of the resource level in the hub nodes. This explains the decrease in the two invasion thresholds and the gap that appears between the two thresholds with the increase of degree heterogeneity. A more heterogeneous network has more hub nodes and is more sensitive to increases in β. Thus increasing the degree heterogeneity reduces the gap between the two thresholds [see Fig. 3].
These numerical and simulation results differ greatly from the classical SIS model. In the multiplex networks with a heterogeneous degree distribution, degree heterogeneity enhances disease spreading and the phase transition is hybrid. Besides, there are hysteresis loops in the phase transition of ρ, and the interval between the two invasion thresholds and the hysteresis region decreases as degree heterogeneity increases. When γ → ∞ the network is approximately a RRN, β I inv disappears as hub nodes disappear, and the transition is discontinuous.

Effects of edge overlap
In social networks two individuals can be friends in the social relation layer and coworkers in the physical contact layer. In transportation networks two cities can be connected by both an expressway and a railway. Thus edge overlap is essential in the science of complex networks, especially when studying percolation in multiplex networks [49]. Here we examine how the amount of edge overlap m e between the two layers affects the spreading dynamics. To eliminate the effect of structure, we fix the values γ = 2.2 and < k >= 9. We then use UCM to build a multiplex network with two identical layers m e = 1.0. To generate a variety of m e values, with a probability q = 1 − m e we rewire pairs of links in layer S. Figure 5(a) shows a plot of ρ as a function of β with two typical values m e = 1.0 and m e = 0.2. Note that when the edges between the two layers overlap completely (m e = 1.0) the infected density ρ smoothly increases from 0 to 1 and there is no hysteresis loop. When the rate of edge overlap between two layers is lowered, i.e., when m e = 0.2, a hybrid phase transition appears. The infected density ρ smoothly increases at β = β I inv and then the system acquires a low epidemic region (β I inv < β < β II inv ) in which ρ slowly increases. Subsequently at β = β II inv an infinitesimally small increase in β causes an abrupt jump in ρ and the disease suddenly spreads throughout the entire system. Hysteresis loops appear in the transition process and the arrows indicate their direction. Figure 5(b) shows that the invasion thresholds (i.e., β I inv ) and β II inv and the persistence threshold β per are determined by the susceptibility χ. Note that the hysteresis loop disappears when m e = 1.0, and it no longer satisfies the definition of ǫ > 0.3 at β L . Thus β L is an inflection point at which the increase in ρ accelerates. The theoretical results from the DMP method agree with the simulation results.
To determine how the amount of edge overlap between the two layers affects the spreading dynamics, we perform simulations for values of m e from 0 to 1 and obtain the space in the plane (β, m e ) shown in Fig. 6. The parameter space is separated into phase regions I and II by a critical value of edge overlap m c e ≃ 0.8. When m e < m c e the system falls into phase I in which the phase transition of ρ is hybrid and the space is again separated into three regions by two invasion thresholds β I inv (lower blue circles) and β II inv (upper blue circles). When β < β I inv the system has a no-epidemic region (green) in which all nodes are healthy and in a steady state. When β I inv ≤ β < β II inv the system has a low-epidemic region (orange) in which the infected density ρ increases continuously from 0 to a finite value until it reaches the second invasion threshold β II inv . Figure 3(b) shows a small low-epidemic region when m e ≤ 0.2 such that when γ = 2.2 and m e = 0.0 the value of β II inv − β I inv converges to a non-zero constant value when N → ∞. When β ≥ β II inv the system jumps abruptly to a high epidemic region (red) in which the disease spreads throughout the entire system. The hysteresis loops (purple) appear in phase I. In contrast, when m e ≥ m c e the system falls into phase II in which the phase transition of ρ is continuous. The value of ρ smoothly increases from 0 to 1 and the hysteresis loops disappear. Figure 6 shows that when β < β I inv the value of β I inv decreases as the amount of edge overlap increases. Here edge overlap promotes disease spreading. When β ≥ β I inv the value of β II inv increases as the amount of edge overlap increases. Here edge overlap suppresses disease spreading. We obtain the theoretical value of β I inv using Eq. (16) and β II inv and β per using the method in Section IV.A. Figure 6 shows that the theoretical values marked by the dotted lines agree with simulation results.
To clarify these results, Figs. 7(a) and (b) show a plot of the average recovery rate µ and the number of susceptible neighbors around each infected individual n s /k as functions of β. Note that when the two layers overlap completely (m e = 1.0), n s /k and µ decrease at the first threshold β I inv to a certain value and then decrease continuously to zero, indicating that the infected density in the steady state increases continuously up to 1 as β increases. Time evolution of ρ(t) for β close to β II inv , ∆ρ is the jump of ρ in the steady state for β is close below β II inv and close above β II inv . (c). Time evolution of average resource of all infected nodes R(t) , hub nodes R hub (t) , and infected density ρ(t). t * is the moment when all the neighbors of the infected nodes are in healthy state (there is no definition of R hub (t) and R(t) ). (d). The average resource of R(t) and R hub (t) as functions of t.
In contrast, when m e = 0.5 there are two abrupt jumps of n s /k and µ at β I inv and β II inv , respectively. Here µ jumps sharply to zero at β II inv [see Fig. 7 Figure 7(c) shows the average resource of all infected nodes R(t) and the average resource of hub nodes R h (t) as a function of t for m e = 1.0. Note that when β is immediately below β I inv both R h (t) and R(t) increase continuously until t = t * . When t > t * there is no definition of resource because all infected nodes recover [see ρ(t) for β < β I inv at t = t * ]. When β ≥ β I inv the infection and recovery rates are balanced, and both R h (t) and R(t) fluctuate around a finite value when t → t ∞ . Thus all the infected nodes recover with a certain probability and ρ(t) also fluctuates around a finite value when t → t ∞ [see ρ(t) for β > β I inv ]. With an increase in β, resource availability decreases continuously as the number infected nodes increases until the disease spreads throughout the system and no available resources remain [see Figs. 7(a) and 7(c)]. This accounts for the continuous increase in ρ when m e = 1.0. Figure 7(d) shows the time evolution of the infected density and available resource level for m e = 0.5. Note that when β is immediately below β II inv at the early stage the disease propagates within the local range of seed nodes, and there are sufficient healthy neighbors in layer S to temporarily suppress the spread. This causes a brief increase in available resources at the beginning of the propagation process and a slight decline in the density of infection. Subsequently the disease rapidly spreads along the edges in layer C. When edges in layer S link out (m e = 0.5), with a high probability that infected nodes in layer S infect their neighbors, R h (t) and R(t) rapidly decline, and ρ(t) rapidly increases. Eventually infection and recovery become balanced, and ρ(t), R h (t) , and R(t) converge to finite values. When β > β II inv there is also a temporary increase in both the available resources and the density of infection. However, when the propagation begins, unlike when m e = 1.0 [see Figs. 4(c) and 4(d)] there is no balanced period at the beginning of the process. The infection of the S-state nodes reduces the resource available to a large number of infected nodes in layer S and delays their recovery. This recovery delay further increases the transmission probability in layer C. Thus R h (t) and R(t) ) decline sharply to zero, the density of infection rapidly increases to 1, and cascading infection occurs.
An increase in the overlap between two layers indicates an increase in the local social circle of an individual. When an individual's colleagues (those frequently in contact, defined as the contact layer) and friends (the social relations, defined as the social layer) are the same group of people, the links in these two layers largely overlap. When β < β I inv , seed nodes initially transmit the disease only to immediate neighbors with whom they are in frequent contact. This high-value local effect causes infected nodes to have a higher probability of linking with other infected nodes in layer S and lowers the level of resources available from neighbors. Thus the overlap between two layers increases network fragility against the invasion of the disease, and increases the probability of an epidemic breakout, and thus lowers the epidemic threshold β I inv . In contrast, a lower value of overlap rate between the two layers indicates a more global social circle, neighbors of nodes in the social layer differ from neighbors in the contact layer. The infected nodes in the contact layer can acquire resources from healthy neighbors in the social layer. Thus the network is more robust against the invasion of the disease, and there is a relatively high epidemic threshold β I inv . This is the reason β I inv decreases as m e increases, as shown in Fig. 6. When β I inv ≤ β < β II inv , hub nodes promote disease transmission, the disease breaks out in a finite range, a sufficient number of healthy neighbors are present in layer S to help infected nodes to recover, and infection and recovery remain balanced.

Conclusions
We have investigated how the level of social support affects spreading dynamics using the susceptible-infected-susceptible model in social-contact coupled networks. Links in the social layer represent relationships between friends or families through which healthy nodes allocate recovery resources to infected neighbors. Links in the contact layer represent daily physical contacts through which the disease can spread. Infected nodes do not have resources, and their recovery depends on obtaining resources in layer S from healthy neighbors. We assume the recovery rate of an infected node to be a function of the resources received from healthy neighbors. We use the DMP method to analyze the spreading dynamics. We first examine how degree heterogeneity impacts disease spreading. We find that degree heterogeneity enhances disease spreading, and due to the existence of hub nodes there is a balanced interval β I inv < β < β II inv in which the infection and recovery processes remain balanced. The value of ρ increases continuously from 0 to a finite value at the first invasion threshold β I inv , increases slowly in β I inv < β < β II inv , then suddenly jumps at β II inv . Thus the transition of ρ is hybrid. In addition, increasing the degree exponent γ in the network increases the gap between the two thresholds and the hysteresis region. To analyze the sudden jump of ρ and the hysteresis loops, we examine the spreading process analytically using mean-field approximation in RRNs. Through a bifurcation analysis we account for the existence of the sudden jump of ρ and the hysteresis loops. In addition, in the RRNs the balanced interval disappears when there is a lack of hub nodes. The first invasion threshold β I inv thus disappears. We next fix the degree heterogeneity and investigate the effect of edge overlap between the two layers. We find that there is a critical value m c e . When m e < m c e there is a second invasion threshold β II inv that increases with m e . The value of ρ smoothly increases at β I inv and then suddenly jumps at β II inv , revealing the transition of ρ to be hybrid with the presence of hysteresis loops in this region [see Fig. 6]. In contrast, when m e > m c e the phase transition of ρ is continuous and the hysteresis loops disappear. In addition, when β < β I inv seed nodes can only transmit the disease locally at the early stage. Here an increase in global connectivity with a lower rate of overlap in the social layer (layer S) increases the probability of linking to healthy neighbors and increases the probability that infected nodes will recover. Thus the first invasion threshold β I inv decreases as the overlap rate m e increases. When β > β I inv , increasing the transmission rate increases the fraction of infected nodes, and an increase in global connectivity in layer S increases the probability of linking to infected neighbors and lowers the recovery rate. Thus the second invasion threshold β II inv increases with m e when m e < m c e .
Although researchers in different scientific fields have focused on ways of constraining disease epidemics in human populations, most scientific literature has been devoted to questions concerning the optimum allocation of public resources or the impact of government investment on spreading dynamics. There has been little examination of how social supports affect spreading dynamics, and our novel model fills this gap. In future research on the impact of social supports on spreading dynamics we will focus on such elements as how degree correlation and the clustering coefficient affect epidemic spreading. Other related topics will include the effect of preference-driven resource allocation on spreading dynamics and the interplay between disease dynamics and resource dynamics. The steady state of the spreading process corresponds to conditions dρ(t)/dt = 0 and dθ(t)/dt = 0. We denote θ(∞) as θ and obtain We also define g(θ) as the function of θ in the steady state, which is Here g(θ) is tangent to the horizontal axis at θ c (∞), which is the critical value in the limit t → ∞. The critical condition is dg(θ) dθ | θc = 0.
Solving Eq. (25) we also obtain the critical transmission rate. From Eq. (23) we see that θ = 1 and θ = 0 are two trivial solutions. Figure 8 shows that the number of solutions for Eq. (24) is dependent on β and there is a critical value of β at which three roots of Eq. (24) emerge, implying that a cusp bifurcation occurs. A bifurcation analysis [50] of Eq. (24) indicates that the physically meaningful stable solution of θ will suddenly increase, and there is an alternate outcome-explosive growth in ρ. Whether the unstable state stabilizes to an outbreak state (θ > 0,ρ > 0) or an extinct state (θ = 0,ρ = 0) depends on the initial infection density ρ(0), thus a hysteresis loop emerges. To distinguish the two thresholds of the hysteresis loop, we denote β per as the persistence threshold corresponding to the nontrivial solution θ c > 0 of Eq. (24) at which the disease initially has a large ρ(0) value. Here β inv is the invasion threshold corresponding to the nontrivial solution θ c = 0 of Eq. (24) at which the disease initially has a small ρ(0) value. The interval [β per , β inv ) is the hysteresis region. Figure 8 shows an example illustrating the relationship between ρ and β when k = 10. Note that g(θ) is tangent to the horizontal axis at θ = 1.0 when β per ≃ 0.0066 and at θ = 0.0 when β inv ≃ 0.066 respectively. When 0.0066 ≤ β < 0.066 three roots of Eq. (24) emerge, indicating that a saddle-node bifurcation occurs and the physically meaningful stable solution of θ increases suddenly to 1. If the disease initially has a relatively small infection density, e.g., ρ(0) = 0.01, the system converges to the stable state ρ = 0, which corresponds to θ = 0. On the other hand if the disease initially has a relatively large infection density, e.g., ρ(0) = 0.9, the system converges to the stable state ρ = 1, which corresponds to θ = 1. When β ≥ 0.0066 and β < 0.0066, ρ(0) has no effect on the stable state of the system. Thus β = 0.0066 is the persistence threshold β per and β = 0.066 is the invasion threshold β inv . Figure 9(a) shows the numerical and simulation results in RRNs with a degree k = 10. In RRNs the first invasion threshold β I inv disappears and the transition of ρ is discontinuous, i.e., not hybrid, due to the lack of the hub nodes. In addition the hysteresis loops exist in the transition of ρ. The orange dashed line and the blue line correspond to the theoretical results for ρ(0) = 0.01 and ρ(0) = 0.9, respectively, obtained from Eqs. (21) and (22). Figure 9(b) shows the susceptibility measurement χ vs β for ρ(0) = 0.01 and ρ(0) = 0.9. From these results we find that the theoretical results obtained from the mean-field approximation agree with the simulation results in RRNs.