{'fraction_non_alphanumeric': 0.031932213424110116, 'fraction_numerical': 0.010323392288906979, 'mean_word_length': 4.561111111111111, 'pattern_counts': {'":': 1, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
How network properties and epidemic parameters influence stochastic SIR dynamics on scale-free random networks

With the premise that social interactions are described by power-law distributions, we study a SIR stochastic dynamic on a static scale-free random network generated via configuration model. We verify our model with respect to deterministic considerations and provide a theoretical result on the probability of the extinction of the disease. Based on this calibration, we explore the variability in disease spread by stochastic simulations. In particular, we demonstrate how important epidemic indices change as a function of the contagiousness of the disease and the connectivity of the network. Our results quantify the role of starting node degree in determining these indices, commonly used to describe epidemic spread.


Introduction
Real systems of interactions between humans are commonly studied as interaction networks. In this setting, an edge in a network represents a possible interaction between the connected nodes. The mechanism of disease transmission upon contact provides a strong case for the use of these networks in modelling epidemic. Spreading processes of infectious diseases take place over these networks, thus their structures can be essential for understanding disease transmission. In particular, knowledge on the interplay between network structure and transmission parameters can improve predictions and prevention and control strategies.
A key factor in epidemic spread is the population network structure where a disease may spread. However, a complete picture, in principle, may require the knowledge of every individual in a population and its relationships, for example, as in contact tracing Ferretti et al. (2020). Networks used in different fields have some common geometric characteristics regarding the distribution of the nodes and edges, for example, a few nodes that may act as hubs and the vast majority of nodes with a few neighbours Pastor-Satorras et al. (2015). A classification in this context is given by network generation models, used as generators of synthetic networks, with controlled topological properties. Several types of networks have been proposed over time Keeling and Eames (2005); de Arruda et al. (2018). A simple model was built by Erdös and Rényi Erdös and Rényi (1959) where nodes are connected according to a uniform probability without any preference. However, some structural properties observed in real-world world networks cannot be reproduced by this model as empirically verified by Newman Newman and Park (2003). It has been observed in real networks that the degree distribution of individuals is far from homogeneous and only a few individuals have several connections and the majority have a few connections Albert et al. ( , 2000. Barabási and Albert proposed a model to generate scale-free networks  with a connection mechanism that mimics the natural formation of social contacts. However, their algorithm cannot control the value of the exponent of the power-law distribution Keeling and Eames (2005). The configuration model can overcome this issue, generating a network with a given degree sequence. This model is used, for example, in social dynamics as it captures connectivity features of this class of networks Boccaletti et al. (2006). The study of how epidemic evolve on networks has been tackled with various approaches, theoretical and computational. Theoretical results on this topic can be found, for example, in the monographs Estrada (2011); Kiss et al. (2017). In López-García (2016), the author derives results on a small heterogeneous population before validating them via numerical simulations. In Strona et al. (2018), the authors explore the interplay between network properties and disease characteristics on the spread of the latter. Many other articles are devoted to the study of the impact of the network structure on the evolution of an epidemic Ball et al. (2010); Draief et al. (2006); Ganesh et al. (2005); Zhang et al. (2013). Another interesting feature, studied in Rocha et al. (2011), is the implementation of dynamic contacts, meaning a network in which nodes may delete and create edges in time.
To this end, in this paper we investigate how the interplay between the connectivity of different configuration model networks and the contagiousness of the disease affect the magnitude of the epidemic. Following a summary of the relevant mathematical background, we derive a lower bound and a closed formula for the probability of extinction of the disease in a network, as a function of all the parameters involved and the degree of the initially infected node. We introduce a stochastic model on scale-free random network, which we use to run simulations on different instantiations of the model. The simulations, carried out via a specialized modification of the Gillespie algorithm, validate our theoretical results regarding the probability of extinction of the disease.
We provide a thorough analysis via simulations that highlights the influence of network connectivity and infectiousness on three key epidemic indices: the maximum number of simultaneously infected individuals, the total number of eventually infected individuals, and the duration of the epidemic. In particular, we focus on the role of the position of the nodes from which the disease starts spreading as well as the number of initially infected nodes, exploring simulations for different instantiations of the model parameters. We compare the evolution of the epidemic, measured by the aforementioned indices, in four different cases for the initial infected nodes, categorized by their position in the respective network: hub (degree in the tail of the distribution), mean degree, peripheral (low degree), and randomly chosen. The stochastic simulations with our model, in agreement with the deterministic lower bound for the probability of extinction, provide a quantification of the role of model parameters in the spread of the epidemic on scale-free networks. Overall, our results demonstrate how the combination of deterministic considerations and stochastic simulations can be used to provide insights on epidemic spread.

Power-law networks
In this Section, we describe scale-free networks, which can be used to mathematically evaluate a given model or describe a real network structure, mimicking social contacts between individuals. We then describe the classic algorithm used to generate networks with a certain degree distribution, that is, configuration model (CM).

Scale-free networks
Many real-world graphs follow power law degree distributions, that is, degree distributions with the probability of a node to have k direct neighbours given by p k ∼ k −α . These kind of networks are also called scale-free random graphs. The preferential attachment model of  produces a degree distribution with α = 3. However, there are many examples in which α ∈ (2, 3] (see e.g. (Durrett, 2006, Sec. 1.4)). To explain the scale-free property, Barabási defined the power-law distribution both with a discrete and a continuum formalism; we refer to Barabási (2015) for a more in-depth introduction to this topic, and recall here the definitions and results we need for our analysis.
In a scale-free network, the probability of having i neighbours is where C is the normalization constant such that where k min is the smallest degree we allow a node to have.
The main difference between a random and a scale-free network comes in the tail of the degree distribution, representing the high-k region of p k : highdegree nodes, called hubs, are naturally present in scale-free networks, contrary to random networks. Since all real networks are finite, we may expect that nodes assume a maximum degree, k max , called the natural cut-off of the degree distribution p k (see Boguñá et al. (2004)). This quantity represents the expected size of the largest hub in a network. To calculate k max , we assume that in a network G(N, V ) of N nodes we expect at most one node in the (k max , ∞) regime. We proceed as in Barabási (2015), and set ourselves in the continuum formalism. In this case, the probability to observe a node whose degree exceeds Since the degree distribution has the form of the natural cut-off follows Hence, the highest attainable degree is directly proportional to a power of N between 1 2 (corresponding to α = 3) and 1 (α = 2). Once we have obtained this values for k max , we can compute again the normalization constant C such that kmax kmin p k = 1.

Configuration model
In the case of large networks, the adjacency matrix is typically not available. However, in the case of real networks of which we know the degree sequence, we can generate a graph with precisely the same degrees. Given the number of nodes N and the sequence of degrees {k i } 1≤i≤N of length N (we omit the subscript from now on, for ease of notation), the aim is to construct an undirected graph with N vertices in which the i-th vertex has precisely degree k i . We denote such graphs with G(N, {k i }). Given a degree distribution obtained from observing a stochastic network, the algorithm used to fit this distribution is called configuration model (CM). To construct the network, we start by assigning to each vertex i in the set of nodes a random degree k i , drawn from the chosen probability distribution p(k). Clearly, k max ≤ N − 1 since no vertex can have a degree larger than N − 1. The degrees of the vertices are represented as half-links or stubs, thus we impose the constraint that the sum N i=1 k i must be even, that is First, two stubs are connected to form an edge. After that, another pair of stubs are chosen from the remaining 2m − 2 stubs and connected, respecting the preassigned degrees. The network is completed by repeating this procedure until the stubs run out. The result of this construction is a random network whose degrees are distributed according to p k Newman (2010). If L denotes the numbers of degrees assumed in the network and N 1 , N 2 , . . . , N L the number of nodes of each degree, the average degree in the network is given by Kiss et al. (2017) wherek i are the different degrees of the nodes in the graph, listed from the sequence {k i }. This formula is equivalent to k = k kp k , for this specific realization. Note that this linkage procedure does not exclude self-loops or multiple edges, but their expected number is bounded (see e.g. (Hofstad, 2016, Prop. 7.1)). When the size of the graph N → +∞ with a fixed degree distribution, self-loops and multiple edges become less and less apparent in the global dynamics (see e.g (Durrett, 2006, Th. 3.1.2)).
Let us remark the fact that a giant component, a connected component of the network that contains a significant proportion of the entire nodes, exists if and only if k 2 − 2 k > 0. (2) In power-law distributions, the second moment k 2 diverges for any chosen exponent α ∈ (2, 3) whereas the first moment k is finite. This means that (2) is always satisfied for any configuration model with a power-law tail to its degree distribution as long as 2 < α < 3. Hence, there will always be a giant component. However, the giant component is not necessarily connected to the whole network, and the remaining nodes may form smaller components Newman (2010).

Lower bound for the probability of extinction
We work on an SIR epidemiological dynamics model on a network built with the CM model algorithm, described in Section 2.2. In this network, we introduce a number of infected individuals in an otherwise fully susceptible population; the exact number of infected individuals is specified each time it changes. Later, in Section 4, we study how the introduction of a different number of initially infected individuals affects the evolution of the epidemic. In this section, we consider an epidemic which starts with only one infected individual, and we derive a deterministic formula for the probability of extinction, that is, the probability that one infected individual in the network does not cause a major outbreak. We use this probability as a benchmark to compare with our simulations. As usual, we assume that an individual remains infected for an exponentially distributed time of parameter γ. During its infectious period, an individual infects each of its neighbours (independently of the others) according to a Poisson process of parameter β. Note that modellers often assume that β does not depend on the number of contacts; however, we notice that in most epidemic models (and data that are collected) the infection probability decreases with the number of contacts. The basic reproduction number R 0 , in this setting, is given in Iannelli and Pugliese (2014), in which the expression is computed as where C v is the coefficient of variation, defined as C v = k 2 k . Using this definition, the expression of R 0 can be written also as We focus on this slightly unrealistic, but analytically tractable assumption. We approximate the initial phase of an epidemic by a branching process where all contacted individuals are susceptible. This is a standard approximation in epidemic models; thus, cliques and triangles are neglected, and the population is assumed to be large enough. First, we compute the probability that an infected with i neighbours infects k of them. We start by considering one of them, conditioning on the length of the infectious period: where R := β/γ. We can not now use the binomial distribution, because infections of different contacts are not independent, but correlated by the length of the infectious period. This is due to the fact that if the infectious period is short, it is likely that no neighbour will be infected, while if it is long, most of them will. Indeed, if Q is the number of infected neighbours, The other expressions are more complicated: We can clarify the expression in (4) as Then, by induction, we obtain Although the infections of different contacts are not independent, the expected number of contacts infected, E(Q) can be easily obtained as iR R + 1 because the expected value of a finite sum is the sum of the expected values. We need these computations in order to estimate the probability of extinction of the branching process, to which we refer as a minor epidemic. Indeed, we have a multitype branching process (types are the number of neighbours) and its probability of extinction is given by the smallest positive solution of s = f (s) where s = (s 1 , . . . , s N ) (if N is the number of types) and To be precise, once we found the smallest solution s * , s * i = f i (s * ) represents the extinction probability starting with one individual of type i. For the derivation of these formulae, we refer to Harris (1963); Mode (1971). Note that f i (1, . . . , 1) = 1; hence (1, . . . , 1) is always a solution of f (s) = s. If we find a smaller solution, then the probability of extinction is smaller than 1; otherwise, the probability of extinction is 1. We then need to compute p i (k). An individual with i neighbours will have m 1 of type 1, . . . , m N of type N (where We could then compute the probability of infecting k 1 out of m 1 , k 2 out of m 2 , . . . , k n out of m n and sum over all possible combinations. We do not need to go through all that because the probability of infecting one neighbour is independent of its properties (i.e. of its degree). Hence, if k = k 1 + . . . + k N , We compute the probability of infecting k neighbours, and then that of the k infected neighbours, k 1 were of type 1,. . . k N of type N . We can go a further step in the computation of f i (s); indeed One could go beyond this, but the formulae would become increasingly cumbersome. In order to compute the relevant solution numerically, one can start with a vector s 0 ≤ s * (for instance s 0 = 0) and then compute s n = f (s n−1 ). Iterating this, we converge to the required fixed point.
Remark 1. However, one may notice one property; obviously s 1+Ri . For instance, if R = 0.002, as in the first line of the example, and i ≤ 50, we get s * i > 1 1.1 ≈ 0.91, thus we immediately see that the extinction probability is very high. Hence, a lower bound for the probability of extinction is given by We recall that R = β/γ and i is the degree of the initially infected node.
In order to see if the probability of extinction is 1 or lower, we can resort to R 0 , the spectral radius of the N × N matrix M , whose elements m ji are the expected number of infected of type j generated by an infected of type i. This is very easy to compute, conditioning on the number of neighbours: M is a matrix of rank 1; hence its spectral radius is easy to compute. From Since exists C > 0 such that v j = Cjp j , thus This is exactly the expression in (3).

The stochastic model
We build a network following the CM algorithm 2 3 , choosing k min = 1, and fixing N = 10000; k max then varies with α as described by We focus on the interval 2 ≤ R 0 ≤ 3; its value varies with the parameters α ∈ [2, 3) and β as shown in Figure 1, following (3). The simulation algorithm implements a specialised version of the Gillespie algorithm Gillespie (1977) that makes use of the network structure for enhancing the simulation efficiency. This refinement provides an equivalent algorithm for interactions on a network as in our model, while drastically reducing the simulation times in comparison to an implementation that accommodates the standard Gillespie algorithm. The main difference in our implementation is in the computation of the second order reaction propensities: at any time, an infection can occur as an instantiation of a second order reaction. The propensity of an infection is proportional with the number of possible infection events. However, on a network, the number of possible infection events is given by the the number of edges that connect the susceptible and infected nodes in the network. By exploiting this fact, our algorithm computes the propensity of the second order reactions as the product of the number of these edges and the reaction rate constant β. The propensities of the first order reactions are computed as usual as the product of the number of infected nodes and the rate constant γ.
For each pair or parameters listed above, we simulated the epidemic 100 times. In each simulation, we recorded the maximum number of simultaneously infected nodes and the number of eventually infected, that is, the total number of recovered at the end of the epidemic with I = 0 in the SIR model. For each pair of parameters, we computed the average of these values, depicted below.

Deterministic vs. stochastic "hub case"
In this section, we verify our stochastic model with respect to the deterministic results in Section 3 for the case in which the infection starts from a hub of the network. The first infected node is chosen as the node with the maximum degree in the network whereby the degree is inversely proportional to α: this follows from expression (1).
We first evaluated the probability of extinction in 100 simulations for each pair of parameters α and β and then we considered different threshold values to have a minor epidemic where the number of eventually infected remains small; the threshold values considered in this study are 1, 2, 3 and 4. We then computed the difference between this data and the deterministic probability of extinction. Finally, we computed the mean and the standard deviations in order to choose the "best threshold value" for the comparison.  Table 1, we observe the best agreement when the threshold is set to 3. After choosing the threshold to define the minor epidemic, we compared the deterministic and simulation results to verify our model. This comparison is depicted in Figure 2. The simulation results are in good agreement with the deterministic results, however they do not coincide exactly. This can be explained by the fact that theoretically we consider a network with an infinite number of nodes, but practically we use a finite network with a large number of nodes. In Figure 2 we compare, for the probability of extinction of the epidemic, the deterministic value (blue) obtained in Section 3 with the stochastic realizations (red) with simulations. For each (α, β) couple, we performed 100 simulations. We observe that for the smallest value of β, the probability of extinction seems to be almost independent of the connectivity of the network. Moreover, when β = 0.002, even in the cases in which the epidemic does not immediately die out, the number of eventually infected individuals in each of the 100 iterations do not exceed 15, which can be seen in detail in the supplementary material. However, when the transmission rate increases, the probability of extinction increases for greater values of α. This behaviour can be noted also in Fig. 8 ("Hub node").

Different numbers of initially infected nodes
First, we considered as separate cases, epidemic starting with three different values of initially infected nodes, that are I(0) = 1, I(0) = 5 and I(0) = 10, with four possible initial positions in the network: random, hub, intermediate, and peripheral; we compare the results in the heatmaps in Figures 3, 4, 5 and 6.  If the position of the initially infected node changes, the dynamics of the spread of the disease mutates accordingly. Indeed, when the first infected nodes are on the periphery, the epidemic is diffusing slowly compared to the case in which the epidemic starts in a node with more contacts. The random case is qualitatively intermediate between the peripheral case and the mean-degree case: this follows from the distribution of the degree of nodes in the network, since we have few nodes in the tail of the distribution (hub). Ten Infected  The resulting numbers in the heatmaps do not seem to change significantly between the I(0) = 5 and I(0) = 10 case, suggesting a saturation effect already for a small (compared to the total population) number of initially infected nodes. Moreover, in the random, mean-degree and peripheral cases the number of infected (both simultaneously and eventually) is increasing as the initial number of infected nodes grows.  The hub case represents a special case, as shown in Figure 6. There, the number of the infected nodes is greater in the case of only one initially infected node than the other two cases. This behaviour can be explained in two distinct ways: first of all, the position of the first infected nodes are assigned in a "hierarchical manner", that is, since the hubs live in the tail of the power-law distribution, the second infected node inserted in the network has a degree lower than the first one. Moreover, if two neighbouring nodes are surrounded by several infected at the same time, they will be infected and the will recover approximately simultaneously, creating a certain saturation in the number of available susceptible nodes.
The hub case is a "worst case scenario": the epidemic spreads even if it starts from only one infected node, if it has enough strength in terms of its parameters to do so. Moreover, it is very unlikely that more than one infected node enters a susceptible population simultaneously. Hence, in the next section, we analyze the hub scenario with I(0) = 1 in greater detail.  As expected, the epidemic spreads more with an increase in β, and less with an increase in α, recalling that a greater α indicates an overall less connected network. Moreover, the more central the initial infection is, the more the disease is able to spread through the nodes of the network.
We notice a noisiness of the final time, compared to the other quantities, which is known in literature and to be expected, since the final part of the epidemic is a subcritical branching process, which is characterized by a great variability in its duration.
We notice that the average of maximum simultaneously and eventually infected nodes is increasing in β, and decreasing in α; this trend becomes more evident the higher the number of initially infected nodes. This is due to the fact that the higher number of initially infected nodes, the smaller the probability of all of them being in a small hubs, disconnected from the giant cluster.  (3)), thus causing a greater spread of the epidemic for fixed values of α: this is evident by looking at columns of the following heatmaps. The opposite is true for α: greater values of α mean a less connected network, which impairs the spread of the disease as it can be noticed by looking at rows of the heatmaps; it is less evident, from equation (3), at least, that greater values of α imply smaller values of R 0 , but our simulations confirm the intuition and what was known in the literature for this model Barabási (2015).
Clearly, the more central and connected the initially infected node is, the greater the magnitude of the infection, as the heatmaps labelled "hub", "mean degree" and "peripheral" illustrate; the "random" heatmaps give, qualitatively, an average of the other three.
The only quantity which remains noisy is the final time, visualized in 9: the tendencies are the same as in the two previous figure, i.e. values decreasing with α and increasing with β, but the resulting heatmap is quite noisy.

Discussion
We have shown how an epidemic effects a synthetic network, following a cutoff power law distribution for the degree of the nodes through the use of the configuration model. We provided a deterministic result on the probability of extinction of the infectious disease, which is based on the initial condition of the network and on the initially infected node. Our deterministic results and stochastic simulation results are in good agreement; the gap in the values stems from the fact that we have considered an infinite network in the deterministic framework, and a finite network, yet with a large number of nodes, in the stochastic one. We have analysed various possible initial conditions for the networks we simulated via heatmaps, and compared them with the available deterministic insight. In particular, we explored how the position of the initially infected individual influences the whole epidemic, measured through three indices: eventually infected individuals, maximum simultaneously infected individuals and overall duration of the epidemic. Our analysis and numerical exploration confirm that the infectiousness of the disease is directly proportional to the spread of the epidemic; moreover, we conclude that the same disease (i.e., characterized by the same parameters β and γ) infects more individuals in networks generated with a lower exponent power-law, as one would expect intuitively.
We have simulated the epidemic as a SIR stochastic dynamic using a spe-cialised version of the Gillespie algorithm. Our algorithm is versatile; many different additional features can be implemented, for example, dynamics edges, contact tracing, quarantine. We leave these as inspiration for future work.