Efficient method for comprehensive computation of agent-level epidemic dissemination in networks

Susceptible-infected (SI) and susceptible-infected-susceptible (SIS) are simple agent-based models often employed in epidemic studies. Both models describe the time evolution of infectious diseases in networks whose vertices are either susceptible (S) or infected (I) agents. Precise estimation for disease spreading is one of the major goals in epidemic studies but often restricted to heavy numerical simulations. Analytic methods using operatorial content are subject to the asymmetric eigenvalue problem, limiting the use of perturbative methods. Numerical methods are limited to small populations, since the vector space increases exponentially with population size N. Here, we propose the use of the squared norm of the probability vector to obtain an algebraic equation, which permits the evaluation of stationary states in Markov processes. The equation requires the eigenvalues of symmetrized time generators and takes full advantage of symmetries, reducing the time evolution to an O(N) sparse problem. The calculation of eigenvalues employs quantum many-body techniques, while the standard perturbation theory accounts for small modifications to the network topology.

One of the main goals in epidemics studies of communicable diseases is to correctly predict the time evolution of a given disease within a population [1].The forecasting procedure, which may take numerical or analytic formulations, often encounters obstacles due to heterogeneous populations and the disease spreading dynamics.For instance, ambiguous symptoms among distinct diseases may under or overestimate total reported infections, leading to incorrect estimates of transmission rates.Several epidemic models have been tailored to better grasp general behaviors in disease spreading [2,3].Among them, the simplest one is the susceptible-infected-susceptible (SIS) model.The SIS model is a Markov process and describes the time evolution of a single infectious disease in a population formed by susceptible (S) and infected agents (I).The infected agents carry the disease pathogens and may transmit them to susceptible agents with constant transmission rate β.The model also contemplates cure events for infected agents with constant cure rate γ and so does competition between cure and infection events.
There are two popular approaches often employed to mimic the disease spreading dynamic in populations with fixed size N: compartmental and stochastic ones [4].In the compartmental approach, relevant properties derived from either infected or susceptible agents are well-described by averages, a direct result from the random-mixing hypothesis [5].This enables one to derive non-linear differential equations to match the evolution of disease throughout the population.For instance, the number of infected agents in the compartmental SIS model, n(t), satisfies the following differential equation: with k = N − 1 and basic reproduction number [6] R 0 = β/γ .For homogeneous populations, this is the expected behavior.However, real agents differ from each other, leading to heterogeneous population, in disagreement with the random-mixing hypothesis [7].Stochastic approaches may also be further classified according to their descriptive variable.Similar to the compartmental model, the mesoscopic interpretation usually describes the time evolution of global variables [8,9], however, it allows fluctuations along time.Meanwhile, the microscopic approach describes the disease spreading of individual agents and their interactions, thus introducing fluctuations at the agent level over time.Both approaches mostly differ on how they treat fluctuations due to agent heterogeneity within a given population, in the epidemic processes.
Central to the microscopic stochastic approach is the underlying network used to reproduce the heterogeneity typically found within populations [10].In the network scheme, agents are represented by vertices and their connections are distributed according to the adjacency matrix A for the assigned network configuration (graph) [11].In this case, it is well-accepted that the mean number k of vertex connection in Eq. (1) describes the averaged process.Contrary to the random mixing hypothesis, non-trivial topological aspects of A may be incorporated in the effective transmission and cure rates, producing complex patterns in epidemics [3].The time evolution is dictated by the transition matrix T , whose matrix elements T µν are transition probabilities from network configuration ν to µ [12].In general, one often assumes Markovian behavior to describe disease transmission and cure events, in accordance with the Poissonian assumption [3].The difference between the compartmental and stochastic schemes leads to distinct evolution patterns for statistics as well.
For instance, Eq. (1) displays stable infected population for γ < β, power-law behavior for γ = β and exponential decay otherwise.While all three behaviors are also observed in the stochastic approaches, fluctuations become much more relevant when the number of infected agents, n(t) , is small compared to total population, N. Incidentally, this is the relevant regime to sanitary measures and health policies to contain real epidemics in early stages.
The Markovian approach produces accurate results if the infection transmission is known.
However, its usability is restricted to numerical simulations with small N since computational time is O(N 2 ).This weakness lies in the fact the T is generally non-hermitian [12].
Therefore, left and right eigenvectors are not related by transpositions, limiting the exact diagonalization only to small values of N or special transition matrices.One of the main goals in epidemic studies is the ability to correctly predict how small parameter or topological changes in the network affect the disease spreading.If such predictions are robust, preemptive actions to lessen the epidemic are also expected to achieve better results.This is exactly the subject of perturbation techniques, which make extensive use of scalar product between left and right eigenvectors.In epidemic models, however, one must deal with asymmetric transitions, prohibiting perturbative schemes based on normed scalar products.
Here, we have devised a method to avoid difficulties related to the non-hermiticity of T , where the Kronecker delta δ σ k ↑ = 1 if σ k =↑ and null otherwise.The set {C µ } spans a discrete Hilbert space À.For clarity's sake, we use the following notation: Latin indices run over vertices 1, . . ., N, while Greek indices enumerate 2 N configurations in À.
The next step required to assemble the transition matrix is the definition of operators and their actions on vectors in À.For instance, the operator σz k probes whether the agent located at vertex k is infected (↑) or not (↓): while the operator nk = (σ z k + 1)/2 extracts the number of infected agent at vertex k.Accordingly, the operator n = k nk extracts the total number of infected agents in the population.The k-th agent health status is switched by action of operators σ+ k and σ− k : Note that σ operators satisfy local fermionic anticommutation relations [13].However, their non-local algebraic commutation relations are bosonic: [σ β k , σγ k ′ ] = 0, for k = k ′ and β, γ = ±, z.
Let P µ (t) be the probability to find the system in the configuration |C µ , at time t.The collection of all P µ (t) forms the probability vector, |P (t) = µ P µ (t)|C µ , with µ P µ (t) = 1.For any Markov process, the transition matrix T describes allowed transitions among configurations such that |P (t + δt) = T |P (t) .Under Poissonian assumption [3], one only considers either a single cure or single infection event during a time interval δt.The Poissonian hypothesis tends to be more accurate for vanishing δt.
In the SIS model, any previously infected agent at vertex k is subjected to three distinct outcomes during the time interval δt: transmit the disease to one connected susceptible agent; cure itself; or remain unchanged.The operator σ− k nk produces the desired cure action, while A km σ+ m nk transmits the disease from the k-th agent to m-th agent, given the k-th agent is currently infected and the other is susceptible, as exemplified for the fully connected graph depicted in Fig. 2. If the cure and infection phases are independent from each other then T = Tcure Tinfec .Under this circumstances, the transition matrix is with Γ = γN/β.Once the explicit action of T is known, P µ (t) are readily evaluated.Fig. 3 exhibits numerical results for P µ (t) for µ = 0, 5, 2 N −1, parameter Γ/N = 0.0, 0.1, 0.3, 0.5, 1.2, P 1 (0) = 1 as initial condition and N = 12, in a fully connected network.For increasing Γ/N, the probability P 0 (t) to find the system without infected agents also increases, while the opposite holds true for P following: jk .The reason is the following: networks only assign distribution rules for connections, leaving the vertex distribution and, therefore, the Hilbert space unchanged.For each graph l = 1, . . ., M in the ensemble, one applies the associated transition matrix, T (l) , on the initial configuration |P (l) (0) , producing the probability vector |P (l) (δt) .In this way, one must also consider the ensemble averages.In particular, the average probability to find the system in configuration µ (t).Since the procedure is equivalent to the average of T over the graph ensemble -the network sample -one needs only to consider the network distribution of A. For clarity, we drop the bar symbol and always assume the average over graph ensemble.

II. SQUARED NORM
Up to O(δt 2 ), P µ (t) obeys the following system of differential equations, where Ĥ ≡ (½− T )/δt is the time generator.For time independent Ĥ, is the solution of Eq. ( 6).The operator Ĥ governs the dynamics with eigenvalues {λ µ } and the respective left {χ µ } and right {φ µ } eigenvectors.The eigenvalues satisfy λ µ ≥ 0, vanishing for stationary states [12].Statistics for observable Ô(t) are calculated according to . Among the relevant observables in disease spreading models, the mean number of infected agents, n(t) , and variance, σ 2 (t), exemplified in Fig. 4, are often relevant variables.Formally, they admit eigendecomposition: n(t) = µν γ µν e −λν t and σ 2 (t) = µν ξ µν e −λν t − n(t) 2 , with   derivative of |P (t)| 2 is obtained from the hermitian operator Ĥ = (1/2)( Ĥ + ĤT ), 1 2 Unlike Ĥ, the operator Ĥ has eigenvalues {Λ µ } but the left eigenvectors are computed from right eigenvectors {ψ µ } by simple Hermitian conjugation.The trade-off is that Λ µ may assume negative values, as shown in Fig. 6, and the coefficients ψ µ |P (t) = g µ (t) are complex numbers.As such, the coefficients g µ are not probabilities.Despite this shortcoming, the coefficients g µ are used to evaluate configuration probabilities: An important expression is derived from Eq. ( 7), subjected to the constraint µν C µ |ψ ν g ν (t) = 1.Now, Eq. ( 9) takes a simpler form if is constant, which is the expected outcome whenever the system reaches at least one stationary state.In such case, Eq. ( 9) reads where the collection of coefficients gµ,l ≡ lim t→∞ g µ (t) describes the l-th stationary state.
Table I displays Eq. ( 10) non-trivial solution for SI model with N = 3 and β/N = 0.1.
This simple example is chosen since the solution can be evaluated by brute force and tested against the correct answer.Of course, the trivial solution gµ,0 = δ µ,0 and Λ 0 = 0 also satisfy Eq. (10).
In addition to stationarity, |P (t)| 2 may also assume maximal or minimal values at time instants t c , leading again to Eq. ( 10), the difference being only the evaluation of coefficients g µ (t) at t = t c .Numerical examples are shown in Fig. 7.The time instant t c is important for dynamics as the extremal condition |P (t c )| 2 informs us when the disease spreading rate changes its growth pattern.Accordingly, t c may also be used to estimate the maxima for narrow peaked statistics.For instance, the nonexistence of cure creates a rapid transient   phase in SI model, with all agents infected as stationary state.During the transient, the variance n 2 (t) − n(t) 2 is well described by a narrow function with peak near t c .The estimation improves as N increases.Therefore, by solving the constrained algebraic Eq. ( 10), either directly or via functional minimization, one also evaluates crucial statistics.
We note Eqs. ( 9) and ( 10) introduce a novel way to tackle stochastic problems: asymmetric operators are replaced by symmetric operators and the eigenspectra are used to evaluate stationary states in Eq. (10).Furthermore, the stationary states obtained in this way carry the network topological information as the adjacency matrix determines the eigenvalue distribution.In the large N ≫ 1 regime, the eigenspectrum becomes dense and it is convenient to analyze Eq. ( 9) using the continuous variable Λ.For completeness sake, we briefly discuss this regime in Appendix A. Alternatively, one may consider the extremal |P (t c )| 2 and obtain t c .In turn, t c may be employed to estimate the time at which statistics develop maxima, as long as they are narrow peaked functions.Since the method is valid for any Markov process, it can be employed for more realistic epidemic models.

III. PERTURBATION THEORY
The eigenvalues Λ µ are crucial to Eq. ( 10) whereas the eigenvectors |ψ µ are required to ensure the probability conservation constraint.In this section, we consider small perturbations to the network link distribution and their corresponding effects on disease spreading in the SIS model.

IV. REGULAR NETWORK
Fully connected networks are the simplest instances of a larger set known as regular networks.Other relevant element in the same set is obtained when the connection patterns among vertices are periodic.Lattices are their spatial representation and are widely employed to describe translation invariant systems.Their natural eigenset contains long and short range modes, allowing analytical tools to inspect long and short range disease spreading behavior, their characteristic frequencies and long range correlations.Since perturbative effects are our main concern here, we only consider a network with single period, or equivalently, a one-dimensional lattice of size N with periodic boundary condition, as Fig. 9a) illustrates.The adjacency matrix is A P and the matrix elements are (A The perturbative scheme to the network topology adds connections with probability δp ≪ 1 among vertices not previously connected, as shown in Fig. 9b).The perturbation creates shortcuts throughout the network, favoring rapid disease dispersion, in an attempt to mimic the relevant aspects found in small-world networks [16].For a single graph realization, translation invariance breaks and the important expression N −1 k Ôk = Ô1 is no longer valid for a general observable Ôk .However, for a large ensemble of graphs, the average transition matrix recovers translation invariance.The reasoning behind this claim lies in the fact all vertices would have 2 + δp(N − 2) neighbors, on average.Let p k,k ′ = δp be the probability to create a single link between V k and V k ′ , including nearest-neighbor vertices.Clearly, the idea is to emphasize the emergence of translation invariance and to interpret the perturbation operator as the meanfield disease spreading operator V in Eq. ( 14).Under this assumption, the contributions to the adjacency matrix due to perturbations are δp(1 − δ kk ′ ).One must be careful to subtract contributions from links already accounted by A P , resulting in the symmetric time generator ĤP i.e., the perturbation operator is proportional to δp ′ .The solution for δp ′ = 0 is obtained using techniques from strongly correlated systems and spinchains, in momentum space [13,17].

V. BETHE-PEIERLS APPROXIMATION
In general, perturbations to topology are not required to affect all vertices in the same manner.For instance, consider a network whose links are distributed according to a parametric probability density function p(ω).If the network undergoes a parameter change ω → ω + δω, one may expect A ij → A ij + δω G ij .The matrix G carries all modifications experienced by the network under the change.Accordingly, the symmetric time generator is Ĥ = (β/N)( Ĥ0 + Ĥ1 + δω V G ).The perturbation V G is Hermiticity is sufficient to warrant RayleighSchrödinger perturbation theory.Furthermore, Eq. ( 10) requires first order perturbative corrections Λ (1) µ and g (1) µ must satisfy µ 2Re(g * µ g (1) Here, as usual, Λ (1) µ = ψ µ | V |ψ µ and g (1) ).In addition to perturbative methods, analytical and numerical techniques from manybody theories are now available to epidemic models.This is also true for approximations, such as Bethe-Peierls [18].In this approximation, nk is replaced by global average n.Application to the SIS model in an arbitrary network produces the effective time generator 2 Ĥ′ β/N = ΓN + n j κ j + j Ω j (cos θ j σz j − sin θ j σx j ), where κ j = k A kj is the degree of j-th vertex, Ω j = 2(Γ 2 + n2 κ 2 j ), cos θ j = (Γ − nκ j )/Ω j and sin θ j = (Γ + nκ j )/Ω j .The effective generator in Eq. ( 19) is diagonalized by rotations around the y-axis.

VI. CONCLUSION
In compartmental approaches to epidemics, the role of fluctuactions is underestimated when the population of infected agents is scarce.Disease spreading models using agent based models are limited to small population sizes due to asymmetric time generators and their large O(2 2N ) dimensions.Our findings show |P (t)| 2 is sufficient to avoid the mathematical hardships that accompany asymmetric operators.The squared norm provides a novel way to obtain stationary states and extremal configurations in general Markov processes, including epidemic models.Once stationary states are secured, the standard Rayleigh-Schrödinger perturbative technique becomes available to epidemics, making use of symmetrized operators and their eigenvalues and eigenvectors.The method paves the way for evaluation of corrections to configuration probabilities caused by perturbations in complex topologies, where analytical results are scarce.

Figure 2 .
Figure2.Markov process during the transmission phase.The probability P ↑↓↑↓ (t + δt) to find the system in the configuration |↑↓↑↓ , at time t + δt, depends on probability P ↑↓↓↓ (t) that the system was previously in the configuration |↑↓↓↓ and then transitioned with conditional probability p(↑↓↑↓ | ↑↓↓↓) = β/N to state |↑↓↑↓ .Analogous rationale applies to the configuration |↓↓↑↓ .The other possibility is that the system was already in the state |↑↓↑↓ at time t and remains unchanged during the time interval δt.As such, the probability to remain unchanged equals to one minus the probability to change to any other state.In this example, the graph is fully connected and there are 4 such transitions.

Figure 3 .
Figure 3. Configuration probabilities P µ (t) for µ = 0, 5 and 2 N −1 with N = 12 in a fully connected network.In a), probability P 0 (t) to observe all-cured configuration at time t for various couplings Γ/N .In b), P 5 (t) refers to the probability of transient configuration |↑ 1 ↑ 2 ↓ 3 • • • ↓ N , while c) exhibits the probability with all-infected agents, P 2 N −1 (t), in log scale.The coupling Γ/N = 0 represents the SI model, whose stationary state is described by all-infected configuration.For Γ/N = 0.1, the stationary state is a linear combination of distinct C µ , including all-cured C 0 and all-infected C 2 N −1 configurations.For intermediate couplings Γ/N = 0.3 and Γ/N = 0.5, the stationary state is C 0 with large transient ∆τ ∼ o(N 4 ).
Although left and right eigenvectors are expected to decompose the identity, their actual computation is rather cumbersome, doubling the computational effort and are specially prone to convergence errors.They also lack a clear analytical interpretation.Here we consider the squared norm, |P (t)| 2 = µ P 2 µ , which remains invariant under unitary transformations.First, total probability conservation µ P µ (t) = 1 does not warrant |P (t)| 2 conservation over time.Examples are found in Markov processes that evolve to unique stationary states since more configurations are available during the transient.The epidemic models considered in this study fall into this category, as shown in Fig.5.Second,

Figure 4 .
Figure 4. Standard deviation σ(t) and mean n(t) for SIS model with N = 12 in the fully connected network.The statistics σ(t) and n(t) are shown in a) and b), respectively.Intermediate cure/infection rates Γ/N = 0.3 and Γ/N = 0.5 eradicate the disease after very large time intervals: σ(t) exhibits initial rapid growth, develops a maximum at t ′ c ≡ t ′ c (Γ/N ) and then decays as powerlaw.

Figure 5 .
Figure 5. |P (t)| 2 in SIS model with N = 12 and P 1 (0) = 1.For each coupling parameter Γ/N , |P (t)| 2 always develops a global minimum followed by constant value at stationary state, being unity only for single state configurations.Logarithmic scale is employed to emphasize extremal points at t c ≡ t c (Γ/N ).

Figure 9 .
Figure 9. One-dimensional periodic lattice with N = 8 vertices.a) Vertex V k connects with vertex V k+1 and V k−1 with periodic boundary conditions, V N +1 = V 1 and V 0 = V N .Perturbative link addition with probability δp = δp ′ /(1 + δp ′ ) increases the mean degree d(k) by p(N − 2), allowing long range transmissions.In b), the graph shows the regular connections for V 3 and an additional connection to V 8 .
′ if A kk ′ = 1, and 0 otherwise.Within our framework, each vertex V k holds a single agent, whose current status is identified by σ k .The variable σ k may acquire two values, namely, σ k =↓ (susceptible) or σ k =↑ (infected), fulfilling the two-state requirement.