Finite symmetries in agent-based epidemic models

We present an algorithm which explores permutation symmetries to describe the time evolution of agent-based epidemic models. The main idea to improve computation times relies on restricting the stochastic process to one sector of the vector space, labeled by a single permutation eigenvalue. In this scheme, the transition matrix reduces to block diagonal form, enhancing computational performance.

In recent years, the emergence of Zika and Ebola viruses have attracted much attention from scientific community after reports of their aggressive effects, microcephaly in newborns [1] and high mortality rate [2][3][4], respectively. Despite their intrinsic transmission differences, both viruses spread in a population starting from a single infected individual based on her geographic localization and relationship network. Contact tracing and proper clinical care planning are key parts of the World Health Organization (WHO) strategic plan [5] to mitigate on-going transmissions and incidence cases, requiring the correct spatiotemporal dissemination of the disease. This assertion has renewed the interest in agent-based epidemic models (ABEM).
ABEM are mathematical models that describe the evolution of infectious diseases among a finite number N of agents, along time. For that purpose, agents are labeled using integer numbers k = 0, 1, . . . , N − 1 whereas contacts between agents are mapped via an adjacency matrix A. The matrix elements are A ij = 1 if the j-th agent connects to the i-th agent and vanishes otherwise. Accordingly, the set formed by agents and their interconnection is expressed as a graph as depicted in Fig. 1. In this way, heterogeneity arises naturally since the individuality of agents is taken into account, distinguishing ABEM from compartmental epidemic models [6]. The simplest ABEM, the susceptible-infected-susceptible model (SIS), considers only two health states for agents, infected |1 or susceptible |0 , and the occurrence of the following events during a time interval δt [7]. An infected agent may undergo a cure event and The time interval δt is often chosen so that sequential cure-cure or transmission-cure events are unlikely within the available time window. This is the so-called Poissonian hypothesis [8].
Following Ref. [9], any configuration of N agents is obtained by direct composition of individual agent states. Let µ be an integer that labels the µ-th configuration so that |µ ≡ |n N −1 · · · n 1 n 0 , with n k = 0, 1 and µ = n N −1 2 N −1 + · · · + n 0 2 0 . A simple example for N = 4 is |8 ≡ |1000 , which represents the configuration where only the third agent is currently infected.
From this scheme, it is already clear that there exists 2 N configurations in total since there are two available states for each agent. In what follows, we employ the notation: Latin indices enumerate agents 0, 1, . . . , N − 1 while Greek indices enumerate configuration states 0, 1, . . . , 2 N − 1.
Despite the existence of this exact solution, the applicability of Eq. (3) at this stage is limited to small N ∼ O(20). The reason is the exponential growth of the underlying vector space as 2 N . Here we show algorithms to generate the operatorsT andĤ using finite symmetries or, equivalently, permutation symmetries via Cayley's theorem [13]. These algorithms are usually applied to condensate matter physics [14,15] but due to recent developments in the disease spreading dynamics [9], they may also be employed in epidemiology studies. For pedagogical reasons, we first show how to build the complete 2 N vector space and the corresponding transition matrix. Next, we explore cyclic permutations to construct the cyclic vector space, in whichT is broken down into N smaller blocks. Lastly, we consider the most symmetric cases, which reduce the problem to O(N). These instances correspond to the mean field or averaged networks. The iteration of sparseT over |π(t) produces the desired disease evolution among agents. Relevant steps are shown in Algorithm 1. Numerical codes are shown in pseudocode and Python. 1

I. TRANSITION MATRIX
The transition matrixT for SIS model [9] considering N two-state agents iŝ where Γ = γ/β, δ kl is the Kronecker delta, is the localized number of infected operator (n k = 0, 1) and are Pauli raising and lowering localized operators, respectively. Local algebraic relationships are [n k ,σ ± k′ ] = ±δ k′k and [σ + k ,σ − k′ ] = δ k′k (2n k − ½). Inspection of Eq. (4) readily showsT is not Hermitian. This means left-and right-eigenvectors are not related by Hermitian conjugation. In this scenario, the correct time evolution of π µ (t) using Eq. (3) requires the complete eigendecomposition, i.e. 2 N eigenvalues accompanied by 2 N right-eigenvectors and 2 N left-eigenvectors. This is often the main criticism against ABEM [8].
However, the scenario described above is not entirely correct. The rationale assumed all eigenstates are equally relevant, which is incorrect whenever A exhibits invariance upon the action of groups (sets of transformations). Symmetries allow for the matrix representation ofT in block diagonal form, as shown in Fig. 3. Eigenvectors related to each block share the same eigenvalue (degeneracy), as usual in quantum mechanics [16]. Therefore, the trick lies in selecting the appropriate base in respect to a given symmetry, redirecting computational efforts towards smaller blocks, which is always more efficient than working directly with the full matrix.
In the SIS model, cure events result from actions of one-body operators,σ − kn k ≡σ − k , on configuration vectors. Infection events are two-body operators: one infected agent may transmit the communicable disease to a susceptible agent after interaction between them, in the time interval δt. Interestingly, the resulting interaction also depends on symmetries available to the adjacency matrix A. The symmetries available to A may be further explored to assemble the initial vector space, reducingT to its block diagonal form.
Group operations over A are always finite transformations. One may explore the isomorphism between finite groups and the permutation group via Cayley theorem [13] to build permutation invariant subspaces. To that end, one must select the finite group and the corresponding symmetry. For graphs, the circular representation provides a convenient context to explore the existing symmetries as where N µ is the normalization. For clarity sake the integer p = 0, 1, . . . , N − 1 is used to label the eigenvalue sector. The representative vector |µ p describes the linear combination of N-agent configurations related to |µ by cyclic permutations. For instance, 3 corresponds to the representative vector for µ = 3, with N = 3 in the p = 0 sector. By construction, the vectors |µ p satisfy the eigenvalue equation The eigenvalues e −2ıπp/N are derived fromP N = ½. Since cyclic permutations never change link distributions, only node labels, cyclic permutation eigenvectors are suitable candidates to reduceT to block diagonal form whenever [P ,T ] = 0.

II. CYCLIC VECTOR SPACE
The complete picture of infection dynamics generated by SIS model requires the utiliza-  In what follows, we address four relevant points regarding the permutation eigenvectors |µ p , namely, the criteria used to label eigenvectors; normalization; number of infected agents; and the permutation operation.
The order convention is necessary to calculate the relative phase between configurations related by permutations, in non-trivial linear combinations. Consider Since |µ andP |µ are related by a single cyclic permutation, |φ = e −2ıπp/N |µ p . Note that the linear combinationP |µ p + |µ p = (1 + e −2ıπp/N )|µ p vanishes for p = N/2. Despite the simplicity of the previous example, it already illustrates the relevance of phase difference among cyclic vectors.
Normalization. According to Eq. (8), the squared norm of representative vectors is The evaluation of the scalar product µ|µ p follows directly from Eq. (8). One notices the configuration |µ may appear only once for several linear combinations |µ p , so that For instance, this is the case of 1|1 p . However, a given configuration |µ may contribute more than once if there exist an integer 1 ≤ r ≤ N such thatP r |µ = |µ , i.e., after r cyclic permutations the configuration repeats itself. SinceP N = ½, it follows N/r is the number of times the configuration |µ appears in |µ p . Each contribution adds . This result is conveniently summarized using the repetition number where the primed sum indicates N/r in the upper limit is an integer number. Therefore, µ|µ p = R µ,p /N µ and one obtains from Eq. (15).
We now show two examples to consolidate the discussion around R µ,p and N µ , for N = 4 and two infected agents. The configuration state |3 = |0011 requires N cyclic permutations to repeat itself, so that R 3,p = 1 for any p and the corresponding normalization for |3 p is  Table I for further reference.
Number of infected agents. The number of infected agents using representative vectors is   Fig. 6b).
Next, we focus attention only to sector p = 0, as it holds both all-infected and all-cured representative vectors. This route allows for the exact evaluation of π 2 N −1 , the probability the disease has reached every agent in the system; or the evaluation of π 0 , the disease eradication probability. Roughly speaking, the p = 0 sector also holds the largest dimensionality.
Consider each integer µ in [0, 2 N ) as a potential candidate to assemble the symmetric vector spaces for fixed p. By performing N − 1 cyclic permutations over |µ , one determines the representative state |µ p in Eq. (8) as well as the number of repetitions R µ,p , hence the norm N µ . Algorithm 3 calculates the representative vector |µ p associated with configuration |µ .
Due to the order convention adopted here, the string representation must be converted to the integer representation at the if -clause test. The representative configurations are then stored either in a list or dictionary. As additional benefit, since vector spaces are independent on the problem at hand, the set of representatives may also be stored in a database for further use in different problems, as long as they are subjected to the same symmetry.

III. MATRIX ELEMENTS
The next step is the evaluation ofT in the p = 0 sector. Infection and cure dynamics are the main actors in this context as they inform the way representative vectors |µ 0 interact with each other,T |µ 0 = ′ {ν} T νµ |ν 0 . The prime indicates the sum runs over all eigenvectors in the p = 0 sector, while cyclic permutation invariance implieŝ Eq. (19) tell us the action ofT on the linear combination |µ 0 is calculated from the simpler operationT |µ . The resulting vectors are then permuted, producing the corresponding matrix elements. For instance, considerT |7 0 for N = 3: The relevant data structure forT are the off-diagonal transitions, which are further subdivided into two categories: one or two-body contributions. This is illustrated in Fig. 7 for the SIS model . The finite set of transition rules are passed as a lookup table or, if available, a dictionary. Data is organized as follows: each entry represents one or two-body configuration whose value corresponds to one tuple. Each tuple holds two immutable values: the configuration to which the entry transitions to and the assigned coupling strength. Two-body operators differ from their one-body counterparts due to the fact they require two agent loops and information from the adjacency matrix A, as seen in the Algorithm 5.

IV. CASIMIR VECTOR SPACE
The recent advances in the disease spreading dynamics in realistic populations are intimately linked to network theory [8]. In this context, a network corresponds to an ensemble of graphs sharing common characteristics, whose adjacency matrix occurs according to the network probability distribution function (NPDF). In this sense, a graph is one sample or realization of the network. Statistical properties of networks are derived for each graph and then taking the ensemble average and deviation. In practice, when graphs in the ensemble are large enough (N ≫ 1) and representatives, statistics may also be evaluated for each graph and extrapolated as those of the network.
Two cases hold particular importance for applications of network theory in epidemic models: the mean field and random networks. In the first case, all agents are connected to each other, meaning one infected agent may potentially infect anyone. Hence, the disease tends to spread faster than in constrained networks. Furthermore, all graphs in the mean field ensemble share the same adjacency A MF . In the other case, the connection between agent i and j occurs with fixed probability ρ. However, graphs in the random network ensemble differ from each other. Here, we only consider ensemble averages as a way to extract statistical properties, which is equivalent to set A random Here, our main concern is to employ the cyclic permutation eigenvectors |µ p to generate the eigenvectors of the complete permutation group, |s, m; p . The eigenvectors |s, m; p reduceT in mean field or random networks to block diagonal form with dimension O(N).
Generalization for q number of states is readily available by changing to the integer representation, with a k = 0, 1, . . . , q − 1, concomitant with additional off-diagonal transitions. For instance, the SIRS ABEM generalizes the SIS model as it introduces the removed (R) health state for agents. This additional state often means the agent has recovered from the illness and developed immunity, has been vaccinated or has passed away. In any case, once removed, the agent takes no part in the dynamics of disease transmission, hindering infection events [8]. As such, cure with immunization or death events produce the transition I → R, with probability γ while vaccination S → R occurs with probability ξ. To reflect the updated one-body off-diagonal transition, the following transition rules dictionary is used: Without loss of generality, the initial condition can always be written as where the primed sum runs only over the indices µ, which also labels the representative vectors. The cyclic permutationP k generates the remaining configurations related to |µ whereas the coefficients π µk are the corresponding initial probabilities. From the eigenvalue equation Eq. (9), one calculates the scalar product whereπ µp = N −1/2 k π µk e 2ıπpk/N is the discrete Fourier transform of π µk . Using the previous example, with one infected among N = 3 agents, with π µk = δ µ1 /3 so that R 1p = 1,π 1p = δ p0 / √ 3, and the previous result is recovered.
We now address the case where the evaluation of the desired statistics requires several permutation sectors. In the worst case scenario, every permutation sector contributes equally to the computation. Therefore one must diagonalize each block in order to obtain the relevant eigenvalues and eigenvectors. As a crude approximation, one may consider the N for each block [17], whereas the complexity range for full diagonalization is Thus diagonalization of N blocks reduces the total complexity from N −1 up to N −2 . More importantly, because blocks are disjointed, they can be diagonalized in different processors.
As the closing remark, the algorithms presented here are most suitable for networks with invariance by cyclic permutations. However, they are also convenient whenever the algebraic commutator can be approximated by [T ,P ] =Ô, where the operatorÔ is symmetric under cyclic permutations, [Ô,P ] = 0. In particular,Ô = q 0 ½+q 1P y + β=z,± q βŜ β , with constant q j (j = 0, 1, z, ±) and y ∈ Ê, creates interesting disease spreading dynamics such as localized disease source for q β = qδ β,0 .  The symbol % stands for modulo integer operation while double forward slashes stands for integer division. for k = 0 to k < N do ⊲ Loop over agents for j = 0 to j < N do 3: for i = 0 to i < N do 4: q ← (L j L i )