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] , image segmentation [4] and the propagation of computer viruses in connected devices [5] and in social networks [6] . 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 [7] , [8] 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 [5] , [9] , [10] .
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] , [11] , optimal curing and initialization policies were investigated for the network contagion model of [2] . In [12] , 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 [13] , 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 [14] (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 {Z t } ∞ t=1 of the Pólya urn with memory M forms a Markov chain of order M , see [14] . There are other examples in the literature where a Markovian version of the Pólya process is studied. In [15] , 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 [16] , which was recently adapted in [17] to measure the dynamics, among many other models, the SARS-CoV-2 pandemic.
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 [18] , martingale methods [12] , [19] , stochastic approximations [20] , [21] , [22] and the embedding of reinforcement processes in continuous time branching processes [23] , [24] , [25] . A detailed discussion on these methods can be found in the survey [26] and in [27] . 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 well-known susceptible-infectious-susceptible (SIS) model [28] , [29] , [30] , [31] , [32] , [33] , [34] . 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 (of small-to-medium size) verify, the nonlinear dynamical system that we obtain approximates the true (underlying) Markov process with improved accuracy as we increase the memory order M . 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 (linear) dynamical system always admits a unique equilibrium (which can be exactly determined); while for M > 1, we note that the 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 and derive its stationary distribution and asymptotic marginal distributions for some specific settings. In Section IV, we develop the dynamical system approximation of our model (which is linear when M = 1 and non-linear when M > 1) using a mean-field approximation and investigate the role of memory, as well as the equilibrium of this dynamical system (both for M = 1 and its linearized variant for M > 1). Simulation results are presented in Section V to justify the choice of our approximate 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 no urn is empty at time t = 0, i.e., T i > 0 for all i. 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) , where 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. The ratio U i,t of red balls in urn i at time t can be expressed in terms of the draw variables as follows:
where
To formulate the interaction part of the model, we start with defining an interaction matrix S to be 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,
Remark 2.1: An example of an interacting Pólya urn network is given in [12] , where the authors studied a symmetric interaction matrix where all its entries are positive (i.e., all urns are connected to each other) with identical diagonal entries and identical off-diagonal entries; hence a complete graph is associated to this interaction matrix. In this paper, however, we deal with more general interactions, which are possibly sparse; moreover, as we will describe shortly, we are interested in finite memory models. To this end, it is worth stating a few remarks about the so-called super-urn scheme proposed in [2] which gives a time-varying interaction on a network of urns. In [2] , a super urn is associated to urn i and is formed by combining the balls of urn i with the balls of all its neighbouring urns in the network. In this set-up, the drawing mechanism for each urn i consists of randomly selecting a ball from the super urn associated to urn i. The probability of drawing a red ball for urn i at time t is thus given by
where N i is defined as the set of all the neighbours of urn i and the urn i itself [2] , and the second equality follows from the definitions of X j,t and U j,t . Comparing with (2), one can observe that it is not possible to define a constant interaction matrix that allows to capture this model as a special case of the one proposed here. Nevertheless, if we allow the interaction matrix to change in time, by comparing N j=1 s ij U j,t to the above expression, we conclude that by setting
we obtain a time-dependent interaction matrix for the super-urn model introduced in [2], making the system of [2] a special case of our interaction model (under a time-dependent interaction matrix). This being said, we only consider time-independent interaction matrices in this paper for the sake of analytical tractability. •
We next note that we can write U i,t recursively as follows:
For ease of notation, we define the following (normalized) initial and reinforcement network parameters:
We are ready now to describe a finite memory adaptation of the above interacting Pólya contagion network, which is the focus of this work. 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 urn. 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 [14] 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. We will next show that the sequence of N -tuple draw variables {(Z 1,t , . . . , Z N,t )} ∞ t=1 for our finite memory interacting Pólya network forms an M th order Markov chain. For brevity, from now on, we denote our finite memory interacting Pólya contagion network by IPCN(M, N ), where M is the memory of the network and N is the number of urns in the network; unless we state otherwise, we assume that the underlying parameters are given by (4) . Before showing the Markov property induced by the network-wide draw variable Z t := (Z 1,t , . . . , Z N,t ),
we present a useful characterization of the ratio of red balls at time t in terms of the IPCN(M, N ) parameters in (4).
For an IPCN(M, N ) system, we have that for t ≥ M + 1,
.
Proof: Recall that U i,t is the ratio of the red balls in urn i at time t. However, the expression of U i,t for a finite memory urn, will have additional terms due to the removal of balls along with addition of the balls at every time t ≥ M + 1. More precisely, given a finite memory M for the network, at every time instant t ≥ M + 1, we remove the balls added at time t − M ; hence for t ≥ M + 1, the number of red balls in urn i at time t is
Similarly, the total number of balls in urn i at time t is
The result then follows by dividing (6) by (7), 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 2.3: 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: 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
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 states. 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 (4) are uniform across all urns, by setting ∆ r,i (t) = ∆ b,i (t) = ∆ > 0, R i = R and T i = T , for all i and t. We refer to this setting as the "homogeneous case" throughout. Using (3), we have that
Also by Lemma 2.2, for a homogeneous IPCN(M, N ), we have that
for t ≥ M + 1, where ρ = R T = 1 − σ and δ = ∆ T . By a proof similar to the one given in Proposition 2.3, 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 (10), or by directly simplifying (8), 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 (11) . Also, the probability of going from state a = (a 1,t , · · · , a N,
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 the case of M > 1, since {Z t } ∞ t=1 is an M th order Markov process with 2 N states, the process {W t } ∞ t=1 , defined by W t := (Z t , Z t+1 , · · · , Z t+M ), becomes an (equivalent) Markov chain of order one with an expanded alphabet of 2 N M states. For the Markov chain {W t }, the transition probability of going from state a = ((a 11 , a 21 , · · · , a N 1 ), · · · , (a 1M , a 2M , · · · , a N M ))
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 M and N urns by the matrix Q (M,N ) = q (M,N ) ab whose entries are given above. In the next lemma, we extend the above irreducibility and aperiodicity properties for the Markov process {W t } ∞ t=1 with M ≥ 1. 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 =
. 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, W t+1 = a (0) = Z t+1 = (a 12 , a 22 , · · · , a N 2 ), · · · , Z t+M = (a 1M , a 2M , · · · , a N M ), Z t+M +1 = (b 11 , b 21 , · · · , b N 1 ) .
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 (13).
• 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 consider a few examples of this Markov process. 
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 (10) 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 [14] .
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 [14] . Indeed, the 2-fold (joint) distribution vector of the single-urn (stationary) Pólya Markov chain in [14] is given bỹ
Furthermore, for the homogeneous IPCN(1, 2) system, the joint probability P (Z 1,t = a 1 , Z 2,t = 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. To see this, note that for any urn i in the given interacting Pólya urn network,
where j is some urn in the network. Taking t → ∞ in the above equation, we obtain lim t→∞ P (Z i,t = 1) = ρ 1 + δ (1 + δ) = ρ.
•
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 (13), 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. 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 [28] , [29] , [30] , [34] . 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 P i (t) := P (Z i,t = 1).
Recall from Lemma 2.2 that the conditional infection probability of urn i at time t, given all the draw variables at time t − 1, is given by
Thus
× P (Z 1,t−1 = a 11 , · · · , Z N,t−1 = a N 1 )
where B N,1 := {(a 11 , a 21 , · · · , a N 1 ) | a i,1 ∈ {0, 1} for i ∈ {1, 2, · · · , N }}.
Note that in the R.H.S. of (15), we can obtain P (Z i,t |Z 1,t−1 = a 11 , · · · , Z N,t−1 = a N 1 ) using the transition probability matrix given in (12) , but the joint probability P (Z 1,t−1 = a 11 , · · · , Z N,t−1 = a N 1 ) is hard to determine explicitly. We thus consider a version of the so-called mean-field approximation, which has been commonly used in the literature on epidemics; for instance, see [28] , [29] . We now describe this approximation in more detail. For a given time instant t, we let
In other words, we assume that at any given time t, the draw variables {Z i,t−1 } are approximately independent of each other at time t − 1 and hence the joint probability P (Z 1,t−1 = a 11 , · · · , Z N,t−1 = a N 1 ) is approximately equal to the product of marginals. Using this approximation and (14), we have that
Note that (18) already gives us a dynamical system for the evaluations of P i (t). We find it useful to first demonstrate how we can arrive at a linear discrete-time dynamical system representing the evolution of P i in time for the case where N = 2. We formulate this derivation in the next example.
Example 4.1: Consider an IPCN(1, 2) system, and suppose that we wish to study an approximation to the evolution of P (Z i,t = 1) in time. By expanding (18), we have that
After rearrangement and cancellation of some terms, we obtain that
Similarly
Using (19) and (20), we obtain
1 (0) + s 11 β 
1 (0)
• Our main idea now is to generalize the procedure that we introduced in the above example, by establishing a linear dynamical system for the IPCN(1, N ) system. To this end, let us define the vector P (t) as P (t) = [P 1 (t), P 2 (t), · · · , P N (t)] T . For the IPCN(1, N ) system, the infection vector satisfies the equation
where J N,1 ∈ R N ×N , C N,1 ∈ R 1×N are matrices with respective entries:
The proof of this result relies on the following lemma.
For real numbers x 1 , x 2 , · · · , x N ,
where the set B N,1 is defined in (16) .
Proof: We prove (22) using induction on N . First, note that the base case N = 2 trivially holds, since
Suppose now that (22) holds for N . We proceed to prove that it also holds for N + 1. To this end, note that
Before giving the proof of Theorem 4.1, for ease of displaying the equations, we introduce the following notation,
We are now in a position to prove Theorem 4.1.
Proof of Theorem 4.1: We show that for all i ∈ {1, 2, ..., N },
We use induction on N to prove Theorem 24. The base case for N = 2 is already proved. Assume that (24) holds for N . Then, for a network of N + 1 urns, we can use (14) to write:
After simplifications, we obtain that
where we have used Lemma 4.2 to conclude the last equality. The result now follows by induction.
We next look at the equilibrium of this linear dynamical system. Proof: 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 [35] . 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
which yields the result.
As an illustration, we find the equilibrium of the linear dynamical system (21) for a much simpler IPCN(1, N ) system.
Given an IPCN(1, N ) system with S = I,
Proof: For an IPCN(1, N ) system with S = I, we have that
In this case, since S = I, the draw variables of urns are independent of each other and hence the mean-field approximation given in (17) holds true with an equality sign. Hence, the 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 (26) .
Another way to find this equilibrium point is to write the transition probability matrix for a single urn using (8) 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 R.H.S. of (26).
We also illustrate Theorem 4.3 by examining the special homogeneous case. For a homogeneous IPCN(1, N ) system, the equilibrium of (21) is given by P * = ρ1 N , where 1 N is vector of ones of size N .
Proof: By Theorem 4.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 (13) .
Unlike the memory M = 1 case, the dynamical systems that we obtain here are inherently nonlinear. Also, we will use mean-field approximations in two different ways:
• We assume that for every time instant t > M , the draw variables of urns {Z i,t } are independent of each other for times t − 1, t − 2, · · · , t − M , where M is the memory of the network.
• We assume that for every time instant t, for each urn i, Z i,t−1 , · · · , Z i,t−M do not depend on each other.
Clearly, the latter approximation was not present in the case where M = 1. For the IPCN(M, N ) system, we have
Now letting, a iM = (a i1 , a i2 , · · · , a iM ) and Z i,t−1 = (Z i,t−1 , Z i,t−2 , · · · , Z i,t−M ), we have
We now use the mean-field approximation
i.e., the past network draw variables (from time t − 1 to t − M ) are independent from each other (spatially and in time). Then, using approximation (29) in (28) and using (27) , we obtain
where Θ t (a j ,k , P j ) is defined in (23) and B N,M := {((a 11 , a 12 , · · · , a 1M ), (a 21 , a 22 , · · · , a 2M ), · · · , (a N 1 , a N 2 , · · · , a N M )) | a ik ∈ {0, 1} for i ∈ {1, 2, · · · , N }, k ∈ {1, 2, · · · , M }},
Note that (30) is a nonlinear discrete-time dynamical system when we consider the case of memory M ≥ 2. The main aim of this section is to give a useful rearranged form of (30). In particular, even though (30) 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 (30).
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 (30) for a general IPCN(M, N ) system. Before moving on to the main theorem, we need the following technical result.
Lemma 4.6: For positive real numbers x jk , j ∈ {1, · · · , N } and k ∈ {1, · · · , M }, the following holds true:
where the set B N,M is defined in (31) .
Proof: We re-index the numbers in (32) such that:
Now, we can use Lemma 4.2 to obtain that
Hence the proof follows.
We now move on to our main theorem of this subsection, in which we write the nonlinear dynamical system for general memory M with each order term separated.
Theorem 4.7: For the IPCN(M, N ) system, the approximating dynamical system (30) 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: We show that (33) is obtained by a rearrangement of (30) . First, note that we can write the R.H.S of (33) as:
Now, using Lemma 4.6 we have
Hence the R.H.S of (33) is given by
The constant term can be extracted from (34) by setting a j = (0, 0, · · · , 0) in B 1,M for 1 ≤ j ≤ N and is given by
Now, fixing j ∈ {1, 2, · · · , N }, we expand the term
Note that the order of (35) is M . In order to get the nth degree term (where 1 ≤ n ≤ M in (35)), 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 j,M , P j (t − k)'s with coefficients 1 and the rest n − v j,M of them have coefficient −1.
Summing over all the possible coefficients of the nth degree term of (34) we get
Finally, we can obtain the nth degree terms (1 ≤ n ≤ M ) for the other N − 1 urns in the exactly the same way as above.
The analysis of the nonlinear dynamical systems given in (33) is clearly more intricate than the one in the case with memory one, where the evaluations were given by a linear dynamical system, namely (21) . This being said,
given that the presence of non linearity is due to the product of probabilities, we can use a further approximation by considering the leading linear terms.
Corollary 4.8: The linear part of the dynamical system (33) is given by
Proof: Setting n = 1 in (33) we obtain
Equation (36) gives an approximate linear dynamical system for the IPCN(M, N ) system. For M ≥ 1 network of N urns, we define P (t) := [P 1 (t), · · · , P 1 (t − M ), , P 2 (t), · · · , P 2 (t − M ) · · · , P N (t), · · · , P N (t − M )] T .
Using (33) and dropping the nonlinear terms, we can write,
where, J N,M is a block matrix with N 2 blocks of size M × M . 
Here, the diagonal blocks of matrix, J N,M (i, i) is given by 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 three different setups. The first two examples are aimed to demonstrate the impact of memory, as well as urn properties. In particular, for two networks with N = 100 urns (highly non-homogeneous) and N = 10 urns (non-homogeneous, but with a smaller range of parameters), depicted in Figure 1 and Figure 2 , respectively, we simulate the IPCN(M, N ) system for M = 1, 2, 3 and their corresponding approximating (nonlinear) dynamical systems given by (33) . We also simulate the linear approximation (36) of the nonlinear dynamical system for each M = 2, 3. Finally, in Figure 3 ,
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 We observe from the simulations that the (nonlinear) dynamical system (33) is a good approximation for the IPCN(M, N ) system. This can be seen in particular in Figure 1 , where the parameters are highly non-homogeneous across the network. We also observe that as memory increases, the quality of this approximation increases. This quality also depends on the values of ρ's and δ's. The dynamical system approximation is better in Figure 2 as compared to Figure 1 . The main reason for the approximation to be better is that the values of δ r 's and δ b 's are close together in the scenario considered in Figure 2 as compared to Figure 1 . We observe again in Figure 2 that the nonlinear dynamical system is a better approximation as we increase the memory. Note that in (33) , 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 (33) to obtain the linear approximation (36), as we expect, the approximation gets worse. The worsening of linear approximation with increasing memory is clearly visible in Figure 1 . Lastly, we observe from the simulations for the homogeneous IPCN(M, N ) system in Figure 3 that the empirical sum as well as linear and nonlinear dynamical approximations converge to ρ irrespective of the memory of the system.
This phenomenon is indeed shown in Theorem 4.5 for any homogeneous IPCN(1, N ) system.
We have formulated an interacting Pólya contagion network with finite memory. We have shown that for the homogeneous case, i.e., all urns have identical initial conditions and reinforcement parameters, the underlying Markov process is irreducible and aperiodic and hence has a unique stationary distribution. For our interacting Pólya contagion network, we have constructed an approximating dynamical system to evaluate the probability of infection in each urn using a mean-field approximation. We show that for memory M = 1, the approximating linear dynamical system has a unique equilibrium point to which the solution asymptotically converges. For memory M > 1, the approximating dynamical system is nonlinear in general; we obtain a linearization of this dynamical system and characterize its equilibrium. We have provided simulations comparing the corresponding linear and nonlinear approximating dynamical systems with the original stochastic process to show the quality of our approximation.
Notably, we have demonstrated that the approximating dynamical system performs well, particularly, when the urn characteristic are not far from each other. Also, as the memory of the system increases, the quality of the approximation improves. Future work include analyzing the stability properties of the nonlinear model and designing curing strategies for the proposed model with systematic comparisons with the SIS model.