key: cord-0578088-s9h2sp2q
authors: Singh, Somya; Alajaji, Fady; Gharesifard, Bahman
title: A Finite Memory Interacting P'{o}lya Contagion Network and its Approximating Dynamical Systems
date: 2020-10-01
journal: nan
DOI: nan
sha: c70691c5c429538717317b692e050e818089c454
doc_id: 578088
cord_uid: s9h2sp2q

We introduce a new model for contagion spread using a network of interacting finite memory two-color P'{o}lya urns, which we refer to as the finite memory interacting P'{o}lya contagion network. The urns interact in the sense that the probability of drawing a red ball (which represents an infection state) for a given urn, not only depends on the ratio of red balls in that urn but also on the ratio of red balls in the other urns in the network, hence accounting for the effect of spatial contagion. The resulting network-wide contagion process is a discrete-time finite-memory ($M$th order) Markov process, whose transition probability matrix is determined. The stochastic properties of the network contagion Markov process are analytically examined, and for homogeneous system parameters, we characterize the limiting state of infection in each urn. For the non-homogeneous case, given the complexity of the stochastic process, and in the same spirit as the well-studied SIS models, we use a mean-field type approximation to obtain a discrete-time dynamical system for the finite memory interacting P'{o}lya contagion network. Interestingly, for $M=1$, we obtain a linear dynamical system which exactly represents the corresponding Markov process. For $M>1$, we use mean-field approximation to obtain a nonlinear dynamical system. Furthermore, noting that the latter dynamical system admits a linear variant (realized by retaining its leading linear terms), we study the asymptotic behavior of the linear systems for both memory modes and characterize their equilibrium. Finally, we present simulation studies to assess the quality of the approximation purveyed by the linear and non-linear dynamical systems.

Interacting urn networks are widely used in the field of applied mathematics, biology, and computer science to model the spread of disease [1] , [2] , consensus dynamics [3] , [4] , image segmentation [5] and the propagation of computer viruses in connected devices [6] and in social networks [7] . A general description for an interacting two-color urn network is as follows. We are given a network of N urns. At time t = 0, each urn is composed of some red and some black balls, where different urns can have different initial compositions, but no urn is empty.

At each time instant t, a ball is chosen for each urn with probability depending on the composition of the urn itself and of the other urns in the network, and then additional ("reinforcing") balls of the color just drawn are added to the urn. Letting U i,t denote the ratio of red balls in urn i at time t, the draw variable Z i,t , denoting the indicator function of a red ball chosen for urn i at time t, is governed by

where "w.p." stands for "with probability" and f : R N → (0, 1) is a real-valued function which accounts for the interaction in the network of urns. The stochastic process {Z i,t } ∞ t=1 is commonly known as a reinforcement process generated by an urn model.

Although, a variety of models are used in the literature depending on the application, Pólya urns are the most commonly used urn models (there are a few examples based on other interacting urn networks; e.g., see [8] , [9] for interacting Friedman urn networks). An interesting feature of Pólya urns is that their reinforcement scheme represent preferential attachment graphs in the sense that at each time step, we add balls of a particular color with a probability proportional to the number of balls of that color in the urn. This representation of Pólya urns as preferential attachment graphs makes interacting Pólya network a good choice to model the spread of infection in an interacting population. For more insight on similarities between Pólya urns and preferential attachment models see [6] , [10] , [11] .

In this paper, we study a network of N two-color Pólya urns. The objective is to model the spread of infection (commonly referred to as contagion) through this interacting Pólya urn network, where each urn is associated to a node (e.g., "individual") in a general network (e.g., "population") to delineate its "immunity" level. Each Pólya urn in the network contains red and black balls which represent units of "infection" and "healthiness" respectively.

The reinforcement of drawing a ball for each urn is mathematically formulated such that a weighted composition of other urns in the network (activated by an interaction row-stochastic matrix) affects the drawing process, hence capturing interaction between nodes.

Various other models have been proposed in the literature to portray contagion in networks using interacting Pólya processes. In [2] , the concept of "super urn" in networks is utilized to model spatial contagion, where drawing a ball for each urn consists of drawing from a "super urn" formed by the collection of the urn's own balls and the ones of its neighbours. In [1] , [12] , optimal curing and initialization policies were investigated for the network contagion model of [2] . In [13] , the authors introduce a symmetric type reinforcement scheme for interacting Pólya urns. Finally, an example of a more complicated interaction network is given in [14] , where a finite connected graph with each node equipped with a Pólya urn is considered and at any given time t, only one of the two interacting urns receive balls with probability proportional to a number of balls raised to some fixed power.

An important characteristic of most reinforcement processes generated via urn models is that they are non-Markovian in the sense that the composition of each urn at any given time affects its composition at every time instant thereafter. This property is not realistic when modelling the spread of infection as one should account for the possibility that infection is cured (or that the urn is removed, a possibility that we do not consider here). In this work, we consider an interacting Pólya urn network where each urn has a finite memory, denoted by M ≥ 1, in the sense that, at time instant t > M , the reinforcing balls added at time t − M are removed from the urn and hence have no effect on future draws. This notion of a finite memory Pólya urn was introduced in [15] (in the context of a single urn) to account for the diminishing effect of past reinforcements on the urn process, which is a realistic assumption when modelling contagion in a population. The resulting network draw variables

of the Pólya urn with memory M forms a Markov chain of order M . The memory parameter M gives a new degree of freedom to this interacting Pólya contagion network which makes it more suitable for the study of epidemics.

In particular, the draw process of the network with M = 1 represents the Markov process of the well-known susceptible-infected-susceptible (SIS) model [16] , [17] , [18] , [19] , [20] , [21] , [22] , as each urn i at time t > 1 exhibits two possible states via its draw variables, susceptible (Z i,t = 0) and infected (Z i,t = 1). Hence, our In [23] , a rescaled Pólya urn model with randomly fluctuating conditional draw probability is considered. Another

Markovian Pólya process is the Pólya-Lundberg process [24] , which was recently adapted in [25] to measure the dynamics of the SARS-CoV-2 pandemic, among many other models(e.g., see also [26] where the classical Pólya urn scheme is used).

The techniques used in the analysis of a finite memory Pólya process are quite different from the ones used for any general random reinforcement process. Standard techniques used for the latter case include the method of moments [27] , martingale methods [13] , [28] , stochastic approximations [29] , [30] , [31] and the embedding of reinforcement processes in continuous time branching processes [32] , [33] , [34] . A detailed discussion on these methods can be found in the survey [35] and in [36] . Unlike standard reinforcement processes, we are able to use Markovian properties in our analysis as our model of interacting urns with finite memory M yields an M th order Markov draw process. However, one drawback of working with a memory-M Markov chain over a network is that the size of its underlying transition probability matrix grows exponentially with both M and the network size. To account for this problem, after having introduced our interacting Pólya urn network and investigated its properties in detail, we formulate a dynamical system to tractably approximate its asymptotic behaviour. To obtain this dynamical system, we make the assumption that for any given time t > M , the joint probability distribution of draw processes at times t−1, · · · , t−M for any urn is equal to the product of marginals. This type of approximation is referred to as "mean-field approximation" and is commonly used in the literature on compartmental models, such as the SIS model. The key factor that distinguishes our treatment is the latitude provided by the consideration of memory M ≥ 1, in contrast to the SIS model which is based on a memory one (M = 1) Markov chain. In particular, as our simulations which are performed for both non-homogeneous and homogeneous networks verify that the nonlinear dynamical system that we obtain is a good approximation for the true (underlying) Markov process. We also characterize the equilibrium point of this dynamical system, when the nonlinear dynamical system is approximated by its linear part (the latter approximation is exact for the case with memory M = 1). More specifically, we show that when M = 1,the Markov process mimics precisely a linear dynamical system with a unique equilibrium (which can be exactly determined); while for M > 1, we note that the approximating linearized dynamical system has a unique equilibrium when its governing (block) matrix has a spectral radius less than unity.

In summary, our results provide a novel mathematical framework for the study of epidemics on networks in realistic scenarios where memory is a consideration.

This paper is organised as follows. In Section II, we describe our interacting Pólya contagion network with memory M , which generates a general M th order time-varying Markov network draw process. In Section III, we show that for the homogeneous case (i.e., when all urns have identical initial compositions and the same reinforcement parameters), the network Markov process is time-invariant, irreducible and aperiodic. We obtain the transition probability matrix of this Markov process, illustrate the calculation of its stationary distribution and establish its exact asymptotic marginal distributions. In Section IV, we derive the linear dynamical system for the general (non-homogeneous) network for M = 1 and develop the nonlinear dynamical system approximations for M > 1 using a mean-field approximation. We investigate the role of memory as well as the equilibrium of these dynamical systems (both for M = 1 and its linearized variant for M > 1). Simulation results are presented in Section V to assess the modeling quality of the linear and non-linear dynamical systems. Finally, conclusions and directions of future work are stated in Section VI.

We consider a network of N Pólya urns, where each urn can be associated to a node in an arbitrary network. At time t = 0, urn i contains R i red balls and B i black balls, i = 1, . . . , N . We let T i = R i + B i be the total number of balls in the ith urn at time t = 0, and assume that each urn contains nonzero red and black balls at time t = 0

i.e., R i > 0 and B i > 0. We also let U i,t denote the ratio of red balls in urn i at time t, with its initial value (at time t = 0) given by U i,0 = R i /T i .

We next define the reinforcement scheme, in the form of draw variables, Z i,t , associated with urn i at time t ≥ 1, for our proposed interacting Pólya contagion network:

1 if a red ball is drawn for urn i at time t 0 if a black ball is drawn for urn i at time t where the process of drawing a ball for urn i is governed by (1) and the function f is explicitly defined below.

The drawing mechanism (1) is applied simultaneously to all urns. If a red ball (respectively, a black ball) is drawn for urn i, we add ∆ r,i (t) red balls (respectively, ∆ b,i (t) balls) to urn i. This scheme, which we refer to as the urn scheme, is often captured by a matrix of the form:

We assume throughout that ∆ r,i (t) ≥ 0, ∆ b,i (t) ≥ 0, for all t ∈ Z ≥0 , and that there exist an urn i such that ∆ r,i (t) + ∆ b,i (t) = 0 at all times t.

To formulate the interaction part of the model, we start with defining an interaction matrix S as an N × N row-stochastic matrix with non-negative entries, i.e., each row in S sums to one. Entries of the interaction matrix S are denoted by s ij , where i, j ∈ {1, . . . , N }. The interaction matrix S can also be thought of as a weighted adjacency matrix of a directed graph with each node equipped with a Pólya urn.

Having defined the interaction matrix, we can now explicitly specify the function f used in the drawing mechanism (1). In particular, we set the probability of choosing a red ball from urn i at time t as follows:

Note that at time t > 1, the conditional probability of node i's draw variable Z i,t is a function of all past draw variables in the network, namely, {Z j,k } for j = 1, . . . , N and k = 1, . . . , t − 1. Furthermore, as all draws occur simultaneously, the draw variables Z i,t and Z i ,t are conditionally independent given all past draws in the network, for any i = i ; hence at any time t,

For ease of notation, we define the following (normalized) initial and reinforcement network parameters for i = 1, . . . , N and t ≥ 1:

We are now ready to describe the notion of finite memory for the above interacting Pólya contagion network. In particular, we consider the scenario where we keep the additional balls introduced in the urns after each draw only for a finite amount of time, M ≥ 1, which we call the memory of the network. In other words, the reinforcement process is altered such that, for all urns in the network, we remove the balls added at time t after the (t + M )th draw. This assumption, which was introduced in [15] in the context of a single Pólya urn, makes the model more realistic, because it accounts for the decrease in severity (or influence) of "infection" with time. Also, the balls in the urn at time t = 0 are never removed from the urn. In the epidemic setting, this initial composition of the urns can be referred to as the intrinsic or inherent immunity of the individuals against infection. We will next show that the 

.

(4) 

Similarly, the total number of balls in urn i at time t is

The result then follows by dividing (5) by (6), and by normalizing both numerator and denominator by T i .

We now establish the Markov property for the network draw process {Z t } ∞ t=1 . Proposition 1: For an IPCN(M, N ) system, the stochastic process given by {Z t } ∞ t=1 is a time-varying Markov chain of order M .

Proof 2: Let a t = (a 1,t , · · · , a N,t ) ∈ {0, 1} N . Using (2) and by virtue of the conditional independence of the draw variables Z i,t and Z i ,t given all past draws in the network for all i = i , we have for t ≥ M that

As a result, we have that

Hence the process {Z t } ∞ t=1 is a time-varying M th order Markov chain.

In many settings of contagion propagation, the "individuals" in the network, being the urns in our setting , are "identical" in the sense of having similar initial parameters. Motivated by this, in this section we develop further the stochastic properties of the IPCN(M, N ) system for the case where the underlying parameters given in (3) are uniform across all urns, by setting ∆ r, 

By a proof similar to the one given in Proposition 1, we next note that for a homogeneous IPCN(M, N ), the draw process {Z t } ∞ t=1 is a time-invariant Markov chain of order M . Indeed, using (2) and (8), or by directly simplifying (7), we have

where the conditional probabilities do not depend on time; hence for the homogeneous

We next examine more closely the properties of this time-invariant Markov chain. Starting with the case of M = 1, i.e., for the homogeneous IPCN(1, N ) system, the transition probabilities of {Z t } ∞ t=1 are given by (9) .

Also, the probability of going from state a = (a 1,t , · · · , a N,t ) to state b = (b 1,t+1 , · · · , b N,t+1 ) can expressed as

with d ∈ {1, · · · , N }. We denote the transition probability matrix of this Markov process with memory M = 1 and

, whose entries are given above. Note that for a memory M = 1 Markov process, it is possible to go from any state to any state in one time step with a positive transition probability. Hence, the Markov chain is irreducible and aperiodic.

For 

in one time step is non zero if and only if a ij = b i(j−1) for i ∈ {1, · · · , N } and j ∈ {2, · · · , M }. If the transition probability is nonzero, it is given by

where d ∈ {1, · · · , N }. Similar to the case of memory one, we denote the transition probability matrix for memory is irreducible and aperiodic. For memory M > 1, to prove irreducibility of the Markov chain, we show that given any two states, it is possible to go from one state to another in finitely many time steps with a positive probability. Let us fix two arbitrary states, a = ((a 11 , a 21 , · · · , a N 1 ), · · · , (a 1M , a 2M , · · · , a N M )) and b = ((b 11 , b 21 , · · · , b N 1 ), · · · , (b 1M , b 2M , · · · , b N M )). We next construct an M -step path (which occurs with a positive probability) between states a and b.

• Suppose the Markov chain is in state W t = a at time t. At time t + 1, we go from state a to state,

Since a ij = a (0) i(j−1) for i ∈ {1, 2, · · · , N } and j ∈ {2, 3, · · · , M }, the transition probability of going from state a to a (0) is nonzero and can be obtained using (11).

• At time t + 2 we go from state a (0) to state a (1)

Following this pattern of adding one N -tuple from state b at each time step, we will reach state b in M time steps.

In summary, choosing any initial state, we can reach any other state of the Markov chain in at most M steps.

Hence, the Markov chain is irreducible. Also, note that the period of the state with all zeros is one. Since all the states of an irreducible Markov chain have the same period, we obtain that this Markov chain is aperiodic.

We now have an explicit formula for the entries of the transition probability matrix Q (M,N ) . Since the timeinvariant Markov chain {W t } ∞ t=1 is irreducible and aperiodic, it has a unique stationary distribution and it is ergodic.

We next illustrate this Markov process via a simple example. It is easy to see that π = [π 00 , π 01 , π 10 , π 11 ] satisfies the equation πQ (1,2) = π. We thus have • lim t→∞ P (Z 1,t = 1) = π 10 + π 11 = ρ, lim t→∞ P (Z 2,t = 1) = π 01 + π 11 = ρ,

• lim t→∞ P (Z 1,t = 0) = π 00 + π 01 = σ, lim t→∞ P (Z 2,t = 0) = π 00 + π 10 = σ.

Also using (8) with M = 1, we have for i = 1, 2 that

Hence, irrespective of the used interaction matrix, the asymptotic marginal (one-fold) distributions and urn compositions for the IPCN(1, 2) system are the same as for the single (memory one) Pólya urn studied in [15] ;

this result is proved in general in Theorem 1 below. We, however, next observe that the asymptotic 2-fold draw distributions for the IPCN(1, 2) urns do not match their counterparts for the single Pólya urn process of [15] .

Indeed, the 2-fold (joint) distribution vector of the single-urn (stationary) Pólya Markov chain in [15] is given bỹ

Furthermore, for the homogeneous IPCN(1, 2) system, the joint probability P (Z 1,t = a 1 , Z 1,t+1 = b 1 ) for urn 1 of the homogeneous IPCN(1, 2) system is given by

Thus noting the conditional independence of Z 1,t+1 and Z 2,t+1 given (Z 1,t , Z 2,t ), and using the IPCN(1, 2) matrix Q (1, 2) along with the fact that lim t→∞ P (Z 1,t = a 1 , Z 2,t = a 2 ) = π a1,a2 , we obtain

which shows explicitly by how much the asymptotic 2-fold draw distribution for urn 1 deviates fromπ (2) . Note that by setting s 11 = 1, the error term π 01 (1 − s 11 )δ/(1 + δ) reduces to zero, making the two distributions match, as expected (since when s 11 = 1, urn 1 only interacts with itself).

• Note that it is much harder to derive in closed-form the stationary distribution for the homogeneous IPCN(M, N ) system with M > 1 and N > 2 but we have the following asymptotic marginal probabilities for a homogeneous IPCN(M, N ) system. 

for all urns i in the network.

Using (11), we obtain Setting γ := γ − ρ1 N in the above equation, we have that

Since the eigenvalues of S have absolute values less than or equal to one (as S is a row-stochastic matrix), we obtain that

As seen in the earlier section, in general it is not easy to obtain the stationary distribution for the IPCN(M, N )

Markov chain characterized in (7), which has 2 M N states. Due to this exponential increase in the size of the transition probability matrix with the number of urns N and memory M , it is difficult to analytically solve for the stationary distribution in terms of the system parameters. In this section, we present the main core of our paper, i.e., a class of dynamical systems whose trajectory approximates the infection probability at time t for any urn i in the network. Notably, for the case of M = 1, we observe that the process naturally leads to an exact linear dynamical system without using any approximation. For M > 1, we use a mean-field approximation to obtain an approximating dynamical system. This is in the same theme as the classical SIS model, and its variations, where a dynamical system is often used instead of the original Markov chain. A few examples in which an approximate dynamical system is constructed for a Markov chain are given in [16] , [17] , [18] , [22] . Throughout the section, we consider δ r,i (t) = δ r,i and δ b,i (t) = δ b,i for all time instances t, i.e., we remove the time dependence from the reinforcement parameters of the Markov process.

Our main objective is to analyze the behaviour of the draw variables {Z i,t }, when memory is one. In particular, we obtain a dynamical system for the evolution in time of P (Z i,t = 1), which we outline next.

For ease of notation, given an IPCN(M, N ) and an urn i, we denote the infection probability at time t by

Recall from Lemma 1 and (7) with M = 1 that the conditional infection probability of urn i at time t, given all the draw variables at time t − 1, is given by

Now taking expectation with respect to (Z 1,t−1 , · · · , Z N,t−1 ) on both sides of (13), we get

To this end, defining the vector P (t) as P (t) = [P 1 (t), P 2 (t), · · · , P N (t)] T , we obtain the following dynamical system for the IPCN(1, N ) network. IPCN(1, N ) system, the infection vector satisfies the equation

where J N,1 ∈ R N ×N , C N,1 ∈ R N ×1 are matrices with respective entries:

Proof 5: Follows from (14) .

We next examine the equilibrium of this linear dynamical system. Proof 6: It is enough to show that the spectral radius of the matrix J N,1 is less than one; since the spectral radius is less than, or equal to, the row sum norm of the matrix, it is enough to show that the row sum norm of J N,1 is strictly less than 1, see [37] . Note that

since 0 ≤ ρ j ≤ 1 and δ r,j , δ b,j ≥ 0 for all j ∈ {1, 2..., N }. Hence, the sum of absolute values of entries in ith row of the matrix J N,1 satisfies

which yields the result.

As an illustration, we find the equilibrium of the linear dynamical system (15) for a much simpler IPCN(1, N ) system.

Given an IPCN(1, N ) system with S = I,

Proof 7: For an IPCN(1, N ) system with S = I, we have that

In this case, since S = I and hence the draw variables of urns are independent of each other. The asymptotic value of P i (t) for i ∈ {1, · · · , N } is given by the equilibrium point of the linear dynamical system

which is given by P * ∈ R N whose ith component is given by (17) .

Another way to find this equilibrium point is to write the transition probability matrix for a single urn using (7) and solving for stationary distribution to obtain lim t→∞ P i (t). The transition probability matrix for a single nonhomogeneous urn i is given by

On solving for the stationary distribution, [π 0 , π 1 ]Q (1,1) = [π 0 , π 1 ], we obtain that π 1 indeed equals the right-handside (R.H.S.) of (17).

We also illustrate Theorem 3 by examining the special homogeneous case. This aligns with the result in Theorem 1. IPCN(1, N ) system, the equilibrium of (15) is given by P * = ρ1 N , where 1 N is vector of ones of size N .

Proof 8: By Theorem 3, the equilibrium P * is given by

Note that the row sums in (I − J N,1 ) are given by 1 1 + δ , i.e.,

Therefore,

We now construct a class of dynamical systems which approximates the IPCN(M, N ) Markov chain in (7) .

Unlike the memory M = 1 case, we need to resort to approximations here to obtain dynamical systems. We use the following mean-field approximation here:

• We assume that for every time instant t > M , for each urn i, Z i,t−1 ,· · · ,Z i,t−M are approximately independent of each other; i.e., at any given time instant t > M , we assume that

for all j ∈ {1, 2, · · · , N }.

For the IPCN(M, N ) system, we have from (7) that

Now, taking expectation with respect to on both sides of (19) and using the linearity property of expectation, we obtain

where B M := {(a 1 , a 2 , · · · a M ) | a k ∈ {0, 1} for k ∈ {1, 2, · · · , M }}.

Now we use the mean-field approximation (18) in (20) to obtain the following class of approximating nonlinear dynamical systems

).

For simplicity of notation, we write (22) in the following way:

where

We next give a useful rearranged form of (23). In particular, even though (23) appears to be complicated, after some simplifications, the coefficients of the nonlinear terms follow a binomial pattern. To give an idea of this binomial pattern, we will first present a few examples and then give a proof formula for the rearranged form of (23).

Example 4.1: We note that for the IPCN(2, 2) system, the approximating dynamical system is given by

Next, by expansion, we observe that for IPCN(3, 2) system, the approximating dynamical system is given by

where one can already observe the binomial pattern that we hinted at. •

We will now obtain a rearrangement of (23) for a general IPCN(M, N ) system.

Theorem 4: For the IPCN(M, N ) system, the approximating dynamical system (23) can be written as

where H n,M := {(d 1 , d 2 , · · · , d n ) d i ∈ {1, · · · , M }, d i = d j ∀i, j ∈ {1, · · · , n}}.

Proof 9: We show that (25) is obtained by a rearrangement of (23). The R.H.S. of (23) is given by:

The constant term can be extracted from (26) by setting a j = (0, 0, · · · , 0) in B M for 1 ≤ j ≤ N and is given by

Note that the order of (27) is M . In order to get the nth degree term (where 1 ≤ n ≤ M in (27)), we need to

and the rest M − n chosen terms have to be 1. We then look at the coefficient of the chosen nth order term. Note that the coefficients of the chosen P j (t − k)'s are either 1 or −1, depending on the tuple a j . Given a tuple a j , there are exactly v M , P j (t − k)'s with coefficients 1 and the rest n − v M of them have coefficient −1.

Summing over all the possible coefficients of the nth degree term of (26) we get

Finally, we can obtain the nth degree terms (1 ≤ n ≤ M ) for the other N − 1 urns in exactly the same way as above.

The analysis of the nonlinear dynamical systems given in (25) is clearly more intricate than the one in the case with memory one, where the evaluations were given by a linear dynamical system, namely (15) . This being said,

given that the presence of nonlinearity is due to the product of probabilities, we can use a further approximation by considering the leading linear terms.

Corollary 3: The linear part of the dynamical system (25) is given by

Proof 10: Setting n = 1 in (25) we obtain

Equation (28) gives an approximate linear dynamical system for the IPCN(M, N ) system. For M ≥ 1 network of N urns, we define

Using (25) and dropping the nonlinear terms, we can write, We next present a few simulations to assess how close these class dynamical systems are to our Markov process.

We provide a set of simulations * to illustrate our results. For this purpose, we have considered four different setups which are aimed at demonstrating the impact of memory, as well as initial urn compositions and reinforcement parameters. In particular, for the first two networks with N = 10 (i.e., fig. 1 and fig. 2 ), we use δ r values that are significantly larger than the δ b values in fig. 1 and δ b values significantly larger than δ r values in fig. 2 . In fig. 3 , we consider larger size non-homogeneous networks with N = 100. We simulate the IPCN(M, N ) system for M = 1, 2, 3 and their corresponding approximating (nonlinear) dynamical systems given by (25) . We also simulate the linear approximation (28) of the nonlinear dynamical system for each M = 2, 3. Recall that for M = 1, the linear dynamical system in (15) exactly characterizes the underlying Markov draw process. Finally, in fig. 4 , we simulate a homogeneous IPCN(M, N ) system. Throughout, for the given IPCN(M, N ) system, we plot the average empirical sum at time t, which is given by

For each plot, the average empirical sum is computed 100 times and the mean value is plotted against time. For the dynamical systems, we plot the average infection rate at time t, which is given by 1

We first note from the simulations that for the network with M = 1, the linear system in (15) matches the empirical sum of the draw process, as expected since in this case the linear system is exact.

We next observe that the nonlinear dynamical system (25) is always a good approximation for the IPCN(M, N )

system. Note that in (25) , the order of the approximating nonlinear dynamical system is equal to the memory of the IPCN(M, N ) system and therefore, when we drop nonlinear terms from (25) to obtain the linear approximation (28) as we expect, the approximation gets worse. For M > 1, we can see this worsening of linear approximation in fig. 1 and fig. 3 . However, in some exceptional cases, the linear approximation performs well. An example of this behavior is presented in fig. 2 , where the linear approximations perform as well as the nonlinear ones. An important aspect of these simulations is that the reinforcement parameters play a major role in determining the asymptotic value of the probability of infection. For example in fig. 1 , since the δ r parameters are significantly larger than the δ b parameters (i.e., infection is much more likely than recovery), the asymptotic value of the plots is higher (i.e., the urns tend towards having a larger composition of red balls). Similarly in fig. 2 , since the δ b values are significantly larger than the δ r values, the asymptotic value of the plots are lower (i.e., the urns tend towards having a larger proportion of black balls). Furthermore, the better performance of the linear system observed in fig. 2 relative to fig. 1 and fig. 3 is attributed to the fact that the constant term in the linear approximation (given by (28)) increases when δ r is increased and decreases when δ b is increased. Depending on how large δ r is, the probability of infection as approximated by (28) can exceed 1 and hence the linear approximation does not perform well for these cases. Whereas, no matter how large δ b gets, the probability of infection never gets smaller than 0 and hence the linear approximation performs comparatively better in this case. Lastly, we observe from the simulations for the homogeneous IPCN(M, N ) system in fig. 4 that the empirical sum as well as the linear and nonlinear dynamical approximations converge to ρ irrespective of the memory of the system. This phenomenon is indeed shown in Theorem 1 for any homogeneous IPCN(M, N ) system.

We formulated an interacting Pólya contagion network with finite Markovian memory. We showed that for the homogeneous case, i.e., when all urns have identical initial conditions and reinforcement parameters, the underlying Markov process is irreducible and aperiodic and hence has a unique stationary distribution. We also derived the exact asymptotic marginal infection distribution. For the non-homogeneous interacting Pólya contagion network, we constructed dynamical systems to evaluate the network's infection propagation. We showed that when memory M = 1, the probability of infection can be exactly represented by a linear dynamical system which has a unique equilibrium point to which the solution asymptotically converges. For memory M > 1, we used meanfield approximations to construct approximating dynamical systems which are nonlinear in general; we obtained a linearization of this dynamical system and characterize its equilibrium. We provided simulations comparing the corresponding linear and approximating nonlinear dynamical systems with the original stochastic process. Notably, we demonstrated that the approximating nonlinear dynamical system performs well for all tested values of memory and network size. Future work include analyzing the stability properties of the nonlinear model, studying the scaling of the approximations with the size of the network, and designing curing strategies for the proposed model with systematic comparisons with the SIS model. 

Curing epidemics on networks using a Polya contagion model

A Polya contagion model for networks

On consensus in a correlated model of network formation based on a Pólya urn process

Opinion dynamics under social pressure

Image segmentation and labeling using the Polya urn model

On the spread of viruses on the internet

A dynamic model of social network formation

Synchronization and fluctuation theorems for interacting Friedman urns

Interacting urns on a finite directed graph

On a preferential attachment and generalized Pólya's urn model

Asymptotic behavior and distributional limits of preferential attachment graphs

Initialization and curing policies for Pólya contagion networks

Synchronazation via interacting reinforcement

A generalized Pólya's urn with graph based interactions

A communication channel modeled on contagion

Improved bounds on the epidemic threshold of exact SIS models on complex networks

Epidemic spreading in real networks: An eigenvalue viewpoint

Virus spread in networks

Epidemic processes over time-varying networks

On the dynamics of deterministic epidemic propagation over networks

Analysis and control of epidemics: A survey of spreading processes on complex networks

Stochastic Epidemic Models and their Statistical Analysis

The rescaled Pólya urn: Local reinforcement and Chi-squared goodness of fit test

A non-homogeneous Markov early epidemic growth dynamics model. application to the SARS-CoV-2 pandemic

Asymptotic incidence rate estimation of SARS-COVID-19 via a Pólya process scheme: A comparative analysis in Italy and European countries

Bernard Friedman's urn

Martingale functional central limit theorems for a generalized Pólya urn

Stochastic Approximation: A Dynamical Systems Viewpoint, ser. Texts and readings in Mathematics

Stochastic approximation with random step sizes and urn models with random replacement matrices

Randomized urn models revisited using stochastic approximation

Branching Processes, ser

Embedding of urn schemes into continuous time Markov branching processes and related limit theorems

Functional limit theorems for multitype branching processes and generalized Pólya urns

A survey of random processes with reinforcement

Synchronization of reinforced stochastic processes with a network based interaction

Nonlinear Dynamics and Chaos: With Applications to Physics

Statistical mechanics of complex networks

We sincerely thank Yanglei Song for providing the proof of Theorem 1 and improving the results in Section IV.

Here, the diagonal blocks of matrix, J N,M (i, i) is given by