Epidemic spreading is an important topic in complex-systems theory [1] and much research on its underlying dynamics has been conducted in recent years. Although a strong focus has been on the theoretical analysis of epidemic spreading [2, 3] , research has also included the control and prediction of disease outbreaks [4, 5] , the spread of rumors [6, 7] , and the propagation of computer viruses [8, 9] . As more and more infectious diseases such as Severe acute respiratory syndrome (SARS) [10] , Ebola virus [25] have brought disasters to humans, how to constrain the global pandemics has been one of most important and pressing challenges. In recent years, many immunization strategies have been proposed for containing and limiting epidemics. Traditional immunization strategies fall into two categories. The first category includes topology-based strategies, such as random immunization [11, 12] , targeted immunization [13, 14, 15] , acquaintance immunization [16] , and graph partitioning [17] . Recent successes have used a targeted destruction of the potential transmission network before an outbreak occurs. "Super-blockers" are identified and immunized to efficiently break network connectivity [18] . The second category includes those that focus on the dynamics of the diffusion of information about the disease, such as information-driven vaccination patterns [19, 20, 21] . Another research topic in epidemic spreading is developing optimal strategies of deploying limited resources such that the epidemic outbreak can be most efficiently suppressed [22, 23, 24] .
Most research on immunization strategies and optimal resource deployment assumes that available resources are fixed, static, and exist independent of the dynamic epidemic process, but in real-world scenarios the amount of such available resources as drugs, medical personnel, and financial support are strongly affected by the evolution of the disease. For example, a pandemic, e.g., the Ebola virus disease (EVD) [25] , can quickly become an enormous economic burden to a region [26] , and even after the disease has been brought under control the economic recovery of the region is slow [27] . Much recent research has examined how dynamic changes in resources affect the dynamics of epidemic spreading. Some research has focused on public resources [28, 29, 30] . For example, Ref. [28] describes how resource constraints caused by the outbreak of disease affect the dynamics of the epidemic. They assume that healthy individuals in the system provide the needed resources, and that the number of these healthy individuals decreases as the infection rate increases. Reference [30] finds that there is a critical amount of invested public resource needed to constrain the spread of a disease, and when that amount is larger than the critical value, the disease can be suppressed. If it is not, the fraction of infected individuals can quickly increase. Other researchers assume that real-world infected individuals cannot always receive public resources and must seek help from friends in their social circles [31] , and that understanding this phenomenon is important in controlling an epidemic. Reference [31] examines how social supports affect epidemic spreading in a double-layer multiplex network in which one layer is the pattern of resource allocation and the other is of epidemic spreading. They find a hybrid transition in the fraction of infected nodes that exhibits properties of both continuous and discontinuous phase transitions.
Although the above literature examines the dynamic evolution of resources and their influence on epidemic spreading, it overlooks the phenomenon of resource diffusion among individuals. Such resources as economic wealth constantly flow among individuals. An important topic for research involves the so-called "Matthew effect" [32] in which the flow of economic wealth tends to make the rich richer. This is relevant because infected individuals with wealth tend to receive better treatment and have a higher probability of recovering than those without.
To investigate the properties of resource diffusion and how it impacts disease spreading, we examine its multiplex structure [33, 34, 31] . We form a two-layer multiplex network of N nodes. Each node in one layer has a counterpart in the other layer. The structure of the two layers can differ. For example, a person may have one group of friends with whom they have regular face-to-face contact and another group of friends in the on-line world [35] .
Here we investigate how resource diffusion affects the dynamics of epidemic spreading in two-layer multiplex networks. We assume that resources diffuse among nodes in the social layer S, and that the disease spreads in the physical contact layer C. Because the diffusion of resources among nodes in layer S can be either preferential or random, we introduce a bias parameter α that controls the diffusion. When the nodes are healthy they can generate new resources. The recovery of infected nodes in layer C depends on the resources of their counterparts in layer S. Through simulations we find that the preferential diffusion of resources can change the phase transition type of the fraction of infected nodes at the steady state ρ(∞). When the degree of interlayer correlation r is below a critical value r c , and the initial fraction of infected nodes ρ(0) is large, i.e., ρ(0) = 0.99, the phase transition ρ(∞) changes from two continuous phase transitions to a single continuous transition as α increases. In addition, there are two hysteresis loops accompanying the two phase transitions when α is below a critical value α c , and one hysteresis loop when α > α c . When r > r c , the phase transition of ρ(∞) changes from multiple (when α is too large or too small) to discontinuous, and then to hybrid, with a initial continuous transition followed by a discontinuous transition. There is always a single hysteresis loop. Note that there is an optimal strategy of resource diffusion under which the disease can be most effectively suppressed, and the threshold reaches a maximum.
We model the coupling of the dynamics of disease spreading and resource diffusion in a double-layer multiplex network. Each individual has links with colleagues or coworkers in the physical contact layer and also with friends in the social relation layer. We construct the double-layer multiplex network model using the uncorrelated configuration model to independently generate layers S and C [36] . These two subnetworks have the same number of nodes N, and there is a one-to-one correspondence between nodes in the two layers. Each layer also has its own internal structure. In an uncorrelated double-layer network, the node degrees in the first layer are independent of the nodes degrees in the second. Thus a highdegree node in the first layer does not not necessarily have a corresponding high-degree node in the second. In contrast, in a correlated double-layer network the node degrees in one layer are somewhat dependent on the node degrees in the other layer. Quantitatively, we use the Spearman rank correlation coefficient r [37, 38] in which r ∈ [−1, 1] to characterize the degree correlation between the two layers. For example, when r > 0 the two layers are positively correlated. A larger r value indicates a higher probability that a high-degree node in the first layer matches a high-degree node in the second layer. In contrast, when r < 0 the two layers are negatively correlated. A smaller value of r indicates a higher probability that a high-degree node in the first layer matches a low-degree node in the second layer. The topological structure of the two layers are encoded in the two adjacency matrices A S = {a S ij } and A C = {a C ij }, respectively. If nodes i and j are connected by a link in layer S (C), a S ij = 1 (a C ij = 1), otherwise a S ij = 0 (a C ij = 0).
To examine how resource diffusion affects epidemic spreading we propose a resource-based susceptible-infected-susceptible (rSIS) model to describe the epidemic spreading in layer C. In the rSIS model, each node can be either susceptible or infected. The recovery process of the infected nodes depends on the resources of their counterparts in layer S. We denote ρ i (t) to be the the probability that node i is infected at time t, and ρ(t) the fraction of infected nodes at t, which is determined by averaging over the infection probability of all nodes
Here ρ(∞) is the fraction of infected nodes when t → ∞. We first randomly select a fraction of ρ(0) nodes to be seeds (infected nodes) and leave the remaining nodes in the susceptible state. At each time step the infected nodes transmit the disease to susceptible neighbors at an infection rate β. The recovery of infected nodes is dependent upon resources supplied by their counterparts in layer S.
Because resources can promote the recovery of infected nodes, we consider that when a node in layer S has greater resources the corresponding node in layer C will have a higher recovery rate. We denote µ i (t) the recovery rate of node i at time t, which is a monotonically increasing function of the resource quantity owned by the counterpart of i in layer S. Note that µ i (t) is a constant value for all nodes in the classical SIS model. Specifically, µ i (t) can be expressed
where µ 0 is the basic recovery rate, which we here fix at µ 0 = 0.1, and ω i (t) is the accumulated resources of the counterpart of node i in layer S at time t. The resource diffusion in layer S is dependent upon the state of nodes in layer C. At each time step, if node i in layer C remains in the S state, the corresponding node in layer S generates a new unit of resource. At the same time, depending on the sign of α, it preferentially transfers one unit of resource to one of its neighbors (the target neighbor). Note that the target neighbor is chosen independent of its state, but the target node does not transmit resources to neighbors if it is not in the S state.
We denote φ i→j the resource transfer probability from node i to j and assume that this transfer probability is related to the degree of j. Then φ i→j is
where δ ij = 1 if i = j, otherwise δ ij = 0. The parameter α allows us to tune the degree of preference. When α > 0, φ i→j is positively related to the degree of j and a high-degree neighbor has a high probability of being selected, but when α = 0, every neighbor of node i has the same probability of being selected. Note that when i = j node i retains the unit of resource during the current time step. The resources σ j (t) that node j acquires from healthy neighbors at time t, can be written
When node i in layer C is in the I state, the corresponding node in layer S does not generate a new resource unit nor does it transfer a resource unit to its neighbors. The accumulated resources of the counterpart of node i in layer S are consumed. For simplicity, we assume that infected nodes consume the all resources of their counterparts. Thus ω i (t) returns to 0 at the current time step. The susceptible nodes store the resources to distribute to neighbors or recover when they are infected in the following time.
We use synchronous updating [2] to simulate the coupled dynamic process of disease spreading and resource diffusion. At each time step with a probability β∆t a susceptible node is infected by one of its infected neighbors. Simultaneously, infected nodes recover with a probability µ i (t)∆t, where i = 1...N. We set a time step ∆t = 1 and run each simulation sufficiently long to ensure that the system enters a steady state in which either no nodes are infected or the number of infected nodes fluctuates within a small range.
Here we examine how preferential resource diffusion affects disease spreading in uncorrelated double-layer networks. We focus on networks with a heterogeneous degree distribution because many networked systems in both nature and technological applications are complex and have a heterogeneous degree distribution [39, 40] . We use an uncorrelated configurational model [41, 36] to build a double-layer network in which the degree distribution is P (k) ∼ k −γ S for layer S and is P (k) ∼ k −γ C for layer C, where γ S and γ C are the power exponents. We fix both values of the power exponential at γ S = γ C = 2.2, and both γ S and γ C are denoted to γ if there is no other special statement. To avoid degree correlations between two layers, each layer is made independent. Because the simulations are time consuming, we set the system size to N = N S = N C = 5000. For the maximum degree we use the structural cut-off k max ∼ √ N [42] and set the minimum degree at k min = 2 [43] . To determine the epidemic threshold, we use a susceptibility measurement [44, 45] where · · · is the ensemble averaging, and χ exhibits peaks at the transition points if they exist.
We first examine the fraction of infected nodes at the steady state ρ(∞) as a function of β with a small fraction of seeds ρ(0) = 0.01 and a large fraction of seeds ρ(0) = 0.99. ρ(0) = 0.01, and β per the persistence threshold when ρ(0) = 0.99 [48] . In addition, we denote β I per and β II per the first and the second invasion (persistence) thresholds.
We next examine the underlying mechanism of the hysteresis loop. Figures 1(c) , 1(f), and 1(i) show the ensemble average recovery rate at the steady state µ(∞) = 1/N µ i (∞), for α = −2.0, α = −1.0, and α = 1.0, respectively. We find that for these values of α, prior to the threshold the average recovery rate is µ(∞) = 1.0 and after the threshold it decreases continuously with β. When the spreading process begins with a low fraction of seeds, i.e., ρ(0) = 0.01, the recovery rate is higher than when there is a larger initial fraction of seeds, i.e., ρ(0) = 0.99 [see Figs. 1(c), 1(f), and 1(i)]. This is because when ρ(0) is small the fraction of susceptible nodes (1 − ρ(0)), is sufficiently high to generate a large number of resources. A lower recovery rate for ρ(0) = 0.99 delays the recovery of infected nodes and increases the infection rate λ = β/ µ(∞) [1] . Thus the disease breaks out at a lower threshold when ρ(0) = 0.99, and the value of ρ(∞) is larger than when ρ(0) = 0.01. Consequently there is a hysteresis loop. In addition, when α = −2.0 and α = −1.0 the two curves of µ(∞) for ρ(0) = 0.99 and ρ(0) = 0.01 overlap at some value of β that separates the parameter space of β into two regions. Thus there are two hysteresis loops in the separated regions.
To determine how preferential resource diffusion affects the dynamics of disease spreading, we examine ρ(∞) as a function of β and α ∈ [−2.0, 2.0]. Figures 2(a) and 2(b) show the phase diagrams with initial conditions ρ(0) = 0.01 and ρ(0) = 0.99, respectively, and Fig. 2(c) shows the difference between values of ρ(∞) in 2(a) and 2(b). Note that ρ(∞) increases continuously with β at each fixed α. In addition, when ρ(0) = 0.01 there is a single phase transition with one threshold β inv [circles in Fig. 2(a) ]. When ρ(0) = 0.99 there is a critical α c value below which there is a double phase transition with two transition points β I per and β II per [circles in Fig. 2(b) ]. Note that the thresholds in Figs. 2(a) and 2(b) are the peaks of susceptibility χ. We also find that when β is fixed, ρ(∞) first decreases and then increases with α when β is large, i.e., β > β inv (β > β I per ) if there are two thresholds), and we obtain the minimum value at the α opt where there is optimal resource diffusion that optimally suppresses disease spreading. Note also that the invasion threshold β inv and persistence threshold β I To further explore these results, we study the resource distribution (green circles) in layer S at the steady state when β = (β inv ) − for ρ(0) = 0.01, where β = (β inv ) − is the infection rate immediately below the threshold β inv [see Fig. 3 ]. When ρ(0) = 0.99 we see similar results. Here we denote ω(k S , ∞) the resource quantity of nodes with degree k S at the steady state, where k S is the degree of nodes in layer S, and ω(k S , ∞)/N the scaled value of ω(k S , ∞). Note that ω(k S , ∞) is shortened to ω(k S ). In addition, to determine how the resource distribution in layer S influences the recovery of nodes in layer C at each value of parameter α, we examine how resources are distributed in nodes whose counterparts in layer C have k C degrees, where k C is the degree of nodes in layer C. This allows us to observe the the change trend of recovery rate with α. Figure 3 (a) shows that when α = −2.0 resources move preferentially to low-degree nodes and ω(k S ) as expected decays rapidly with k S . In addition, most of the nodes in the two subnetworks with highly skewed degree distributions are low-degree and only a few are highdegree. Thus the counterparts of the high-degree nodes in layer C have a higher probability of being low-degree nodes in layer S because of the random correlation between the two layers. Thus most of the counterparts to the high-degree nodes in layer C have large values of ω(k S ) in layer S [yellow squares in Fig. 3(a) ], i.e., most high-degree nodes in layer C have a high recovery rate that delays outbreaks of the disease as β increases. When α = 1.0 resources move preferentially toward high-degree nodes in layer S and agglomerate on high-degree nodes at the steady state. When there is a random correlation between the two layers, most high-degree nodes correspond to low-degree nodes in layer C. Thus the resources of k-degree nodes in layer S increase with k S [see Fig. 3(c) ]. In contrast, the ω(k S ) decreases sharply with k C , which indicates that the recovery rate of the high-degree nodes in layer C rapidly declines when β increases and resources decrease. This in turn increases the effective infection rate λ = β/ µ(∞) in the system. Figure 2 shows that a severely skewed distribution of resources lowers the epidemic threshold and a large fraction of nodes when α is large, i.e., α = 1.0.
When α = −1.0, the diffusion of resources in layer S is less biased than when α = −2.0 or α = 1.0. We analyze Eq. (4) and find that although low-degree nodes still have a small advantage of acquiring resources, high-degree nodes can acquire approximately the same quantity of resource at each time step because they have more connections than low-degree nodes. Thus resources are distributed evenly for both high-degree and low-degree nodes [see Fig. 3(b) ]. When resource diffusion is optimal, all nodes in layer C have a rapid recovery rate [see Fig. 1 (f)] that reduces the infection probability between each pair of susceptible and infected nodes. Here the disease is suppressed to the greatest extent. Figure 2 shows that the highest epidemic threshold β inv (β per ) and lowest fraction of infected nodes ρ(∞) are obtained when resource diffusion is optimal, i.e., when α = −1.0.
There are extensive interlayer correlations in real-world multiplex systems [49, 50] . In social networks, for example, an individual with many daily face-to-face contacts with colleagues tends to also have many social network contacts [35] . In transportation networks, hub airports tend to correlate with hub rapid transit stations [51] . We here investigate how the degree correlations between the two layers impact the process of resource diffusion and the dynamics of disease spreading. To construct a double-layer correlated network with an adjustable degree of inter-layer correlation, we first generate two subnetworks of the same size N = 5000 and the same power exponent γ = 2.2 with a maximum positive or maximum negative correlation. We then rematch each pair of counterpart nodes with a probability q. Thus the interlayer correlation after rematching becomes [37, 38] r =| 1 − q | .
When the two layers are initially at maximum positive correlation r ≥ 0, otherwise r ≤ 0. Fig. 5(e) ]. In addition, when α = 1.0 the first threshold disappears when ρ(0) = 0.99. Later we will use a finite-size scaling analysis to demonstrate the discontinuous increase of ρ(∞) [52, 53, 54] . Note that, unlike when r = 0 or r = −0.8, there is single hysteresis loop for all values of α. We can obtain the same explanation for the hysteresis loops by analyzing the ensemble average recovery rate µ(∞) as a function of β, similar to when r = 0. To determine how preferential resource diffusion affects the dynamics of epidemic spreading when there is interlayer degree correlation, we use two-parameter (α, β) phase diagrams for r = −0.8 and r = 0.8 [see We obtain thresholds from susceptibility χ. To explain the optimization we examine the resource distribution of nodes in layer S and how resources are distributed on those nodes with counterparts in layer C that have k C degrees with initial ρ(0) = 0.01 when β = (β I inv ) − [β = (β inv ) − if it is a single phase transition]. Thus we obtain the scatter plots of ω(k S )/N versus k S (green circles) and ω(k S )/N versus k C (yellow squares). We obtain results similar to those when ρ(0) = 0.99. Figures 7(a) -7(c) show resource distributions for α = −1.5, 0, and 1.0, respectively, when r = −0.8. Note that when α = −1.5 the probability that resources move to low-degree nodes in layer S is high. Figure 7 (a) shows that ω(k S )/N decreases sharply when k S in layer S increases (green circles). In addition, when the correlation between the two layers is negative, highdegree nodes in layer C correlate with low-degree nodes. Because low-degree nodes are more numerous in a heterogeneous network, most low-degree nodes in layer C still have low-degree counterparts. Thus both high-degree and low-degree nodes in layer C can rapidly recover because there are adequate resources supplied by their counterparts in layer S [yellow squares in Fig. 7(a) ]. When this is the case, the disease is effectively constrained [see Fig. 6(a) ]. When α = 0 and α = 1.0, resources move preferentially to the few high-degree nodes in layer S and low-degree nodes receive little [green circles in Figs. 7(b) and 7(c) ]. When β increases, the recovery rate of high-degree nodes in layer C rapidly decreases because they cannot receive resources from their counterparts [yellow circles in Figs. 7(b) and 7(c)] and the disease is not constrained. Thus we see a small threshold and a large ρ(∞) when resources move preferentially to high-degree nodes in layer S.
When r = 0.8, to constrain disease spreading the recovery rate of both high and low degree nodes in layer C should maintain a high threshold. To achieve this, resources must diffuse to high-degree nodes in layer S, i.e., α ≃ −0.5, in a positive correlation between the two layers. Thus when α = α opt ≃ 0.5 there is a maximum threshold value and a minimum ρ(∞) value when β is fixed. Figure 7 (f) shows that when resources move only to high-degree nodes in layer S, i.e., when α = 1.0, there are no resources for the low-degree nodes in layer S. Figure 7 (d) shows that when resources move only to low-degree nodes, there are none for the high-degree nodes. In both of these extreme conditions, the node recovery rate in layer C declines rapidly as β increases, which causes an earlier outbreak of disease [see Figs. 6(d) and 6(e)].
We next use a finite-size scaling analysis to examine the discontinuous increase of ρ(∞) when α approaches α opt and the two network layers are positively correlated. We define ρ(N, ∞) the fraction of infected nodes at the steady state for a network with N nodes and ∆ρ(N, ∞) the maximum increase of ρ(N, ∞) during an infinitely small increase of β, which is expressed
where ∆β is an infinitesimal increment of β, set at ∆β = 0.001 in our simulations, and ρ(N, ∞, β) is the fraction of infected nodes at steady state when infection rate is β. When
there is a discontinuous increase in ρ(∞) [55, 51] . Note that we use α = 0 for the finitesize scaling analysis. Figure 8 (a) shows ∆ρ(N, ∞) as a function of N when α = 0.6 (orange triangles), α = 0.7 (green circles), and α = 0.9 (red squares). Note that when α = 0.6, ∆ρ(N, ∞) converges to 0 asymptotically. When α = 0.7 and α = 0.9, ∆ρ(N, ∞) asymptotically converges to a positive constant. Figure 8 (b) shows ρ(∞) as a function of β when ρ(0) = 0.01 for three typical values of interlayer correlation r = 0.6, r = 0.7, and r = 0.9 in a network of size N = 10000. Note that when r = 0.6, ρ(∞) increases continuously with β. When r = 0.7 and r = 0.9, ρ(∞) first increases slowly and continuously at β I inv , and then jumps discontinuously at β II inv , all of which are characteristics of a hybrid phase transition.
We next use extensive simulations to obtain the phase diagram of ρ(∞) in the twoparameter (r, β) plane with an initial condition ρ(0) = 0.01 when α = 0. When ρ(0) = 0.99 the results are similar. Figure 8 (c) shows that when the two layers are negatively correlated (r < 0), ρ(∞) increases continuously with β. When r > 0, there is a critical value point r c [point A in Fig. 8 (c) ]. When r ≥ r c there is a discontinuous change of ρ(∞) at the threshold. Note also that the epidemic threshold increases with r [see white line in Fig. 8 (c) ]. To explain the hybrid discontinuous phase transition, we plot the time evolution of total resources Ω(t) and infected fraction ρ(t) with the initial condition ρ(0) = 0.01 for r = 0.9 when α = 0 [see Fig. 9 (a)]. When ρ(0) = 0.99 the results are similar. Figure 9 (b) shows the corresponding resource distribution at the steady state. When β is immediately below (β II inv ) − , the scaled value of the total resources Ω(t)/N abruptly increases at the early stage of the diffusion process [green circles in Fig. 9(a) ] because almost all nodes in layer C are healthy and resources are constantly generated by the corresponding nodes in layer S. After a longer period of time t > 300 the system enters a steady state, and fluctuations stay within a small range (upper green circles). Here the resources of high-degree nodes are rapidly consumed, and the resource level for low-degree nodes remains high [see Fig. 9 (b)] indicating that the disease is localized around the high-degree nodes. We thus learn that before β II inv the system changes from a disease-free absorbing phase to a locally active phase (in which ρ(∞) reaches a finite small value) at β I inv [green line in Fig. 9(a) ]. For the sake of comparison, Fig. 9 (b) shows a plot of the resource distribution when β = 0.01.
When β = (β II inv ) + , the value of Ω(t)/N rapidly increases as the disease spreads from the local area of the seeds [red circles in Fig. 9(a) ]. As t increases ρ(t) slowly increases and Ω(t)/N reaches a peak value at a crossover time t * . After t * , Ω(t)/N drops rapidly, indicating that the newly-generated node resources in layer S are not sufficient to recover the infected nodes in layer C. The recovery rate of the infected nodes then declines as resources decrease, which induces an increase in the infection rate of the disease, especially in the hub nodes. Thus as the infection rate increases, the resources available in layer S further decrease and the node recovery rate in layer C decreases. Then a cascading effect appears that sharply increases ρ(t) from a small finite value to a value near 1.0 [red line in Fig. 9(a) ]. Figure 9 (a) shows ∆ρ(∞), which is the increase of ρ(∞) when β increases from (β II inv ) − to (β II inv ) + . This indicates a discontinuous jump in ρ(∞). Figure 9 (b) shows that all resources of all nodes in the network have been consumed, in contrast to when (β II inv ) − . Fig. 10(a) ]. When α < −1, β I inv decreases with r, but when α > −1, β I inv increases with r. When α < −1, resources move preferentially to low-degree nodes in layer S. To suppress the spreading, the nodes in layer S must supply enough resources to high-degree nodes in layer C. Thus negative interlayer correlation enhances the disease suppression. In contrast, when α > −1 high-degree nodes add resources in layer S. To constrain these high-degree nodes we must have high-degree counterparts in layer C. Thus we increase β I inv with r. Finally we explore the relationship among the optimal values of the bias parameter α opt at which the disease is maximally controlled. Figure 10(b) shows α opt as a function of r. Note that the value of α opt increases monotonically with r because, with the increase of the interlayer correlation, the probability that the large degree nodes in layer S have counterparts with large degrees also increases. To protect the large degree nodes in layer C, resources in layer S must diffuse preferentially to large degree nodes. Thus α opt increases with r.
We have explored how preferential resource diffusion affects the dynamics of disease spreading in correlated multiplex networks. We assume that resources diffuse in the social contact layer and that the disease is transmitted in the physical contact layer of the network. The two dynamical processes are coupled such that the generation and diffusion of resources in layer S are dependent on the state of nodes in layer C, and that the recovery of infected nodes in layer C are dependent on the resources of their counterparts in layer S. To model the disease spreading in layer C, we propose a resource-based susceptible-infected-susceptible (rSIS) model. Using extensive simulations we find that preferential resource diffusion can change the phase transition in ρ(∞), i.e, when the degree of interlayer correlation r is below a critical value, the transition ρ(∞) in ρ(0) = 0.99 changes from two continuous phase transitions to one single phase transition as the controlling parameter α increases. Note that when ρ(0) = 0.01 the transition of ρ(∞) is single and continuous throughout the parameter space of α. In addition, there are hysteresis loops in the continuous phase transitions. There are two hysteresis loops accompanied by two phase transitions and one single hysteresis loop accompanied by one single phase transition of ρ(∞). When r is above the critical value, the phase transition of ρ(∞) changes from multiple (α is too large or too small) to discontinuous, and then becomes hybrid and exhibits the properties of both continuous and discontinuous transitions (α is near the optimal value). Note that there is an optimal resource diffusion at each fixed value of r. When the diffusion of resources is optimal the threshold reaches a maximum and the disease can be maximally suppressed.
In recent years constraining disease epidemics in human populations has become a hot research topic and has attracted many workers across a variety of fields. Most research has focused on ways of optimally allocating limited public resources, but there has been little examination of how the resource diffusion among the individuals affects spreading dynamics. Our model fills this gap. There remain limits in our model. For example, because the model is complex we have not yet developed theoretical solutions, and thus theoretically obtaining an optimal solution α opt would be an interesting and important path for future research.