# Many-body methods in agent-based epidemic models

###### Abstract

The susceptible-infected-susceptible (SIS) agent-based model is usually employed in the investigation of epidemics. The model describes a Markov process for a single communicable disease among susceptible (S) and infected (I) agents. However, the disease spreading forecasting is often restricted to numerical simulations, while analytic formulations lack both general results and perturbative approaches since they are subjected to asymmetric time generators. Here, we discuss perturbation theory, approximations and application of many-body techniques in epidemic models in the framework for squared norm of probability vector , in which asymmetric time generators are replaced by their symmetric counterparts.

###### pacs:

02.50.-r, 03.65.Fd, 05.10.-a, 87.10.Ca^{†}

^{†}thanks: We have applied techniques from fundamental physics to epidemic models, inspired by the teaching of Roger Maynard, to whom we dedicate our findings.

Proper planning and management lie in the foundation of efficient health and sanitary policies Team (2016). The decision-making process usually relies on predictions from epidemic models and raw data to rule the best available policy to mitigate the disease spreading. Resource planning becomes even more relevant during the emergence of new communicable diseases, as improper actions may extend the duration or adversely impact health workers Heesterbeek et al. (2015). Despite the success attained by traditional epidemic models for large scale epidemics in the past, they have been unable to produce reliable predictions for small scale disease spreading in heterogeneous populations Heesterbeek et al. (2015); Bansal et al. (2007); Eames and Keeling (2002). This issue is further enhanced due to the stochastic nature of pathogen transmission mechanisms and patient care. As such, a considerable amount of epidemic models have been proposed to mimic the correct behavior for spreading of communicable infectious diseases.

The simplest susceptible-infected-susceptible (SIS) model considers the time evolution of a single communicable disease among susceptible (S) and infected (I) agents Pastor-Satorras et al. (2015). The infected agents may either transmit the disease to susceptible agents with constant probability , turning them into infected agents , or undergo the recovery process with probability rate , during a fixed time interval . Furthermore, two approaches are available to describe the disease spreading in the SIS model for a fixed population of size : compartmental and stochastic. In the compartmental approach relevant quantities are well-described by averages Keeling and Eames (2005), from which one derives non-linear differential equations. For instance, the number of infected agents , for fixed , in the compartmental SIS model is

(1) |

This is the expected behavior for large homogeneous populations, where fluctuations are negligible. Introduction of effective transmission () and recovery () probabilities rates contemplates effects due to heterogeneous population. This parameter estimations employ networks, with size , as substrate to display the varying degree of non-uniformity within a population. In this scheme, each vertex in the network contains a single agent, while the links among agents are given by the corresponding adjacency matrix . Thus, non-trivial topological aspects of are incorporated in the effective transmission and cure rates Bansal et al. (2007); Eames and Keeling (2002); Pastor-Satorras et al. (2015); Keeling and Eames (2005); Pastor-Satorras and Vespignani (2001).

The stochastic approach also requires networks to describe the population heterogeneity. However, contrary to the compartmental approach, the assumption about negligible fluctuations is removed, meaning averages alone are not sufficient to properly describe the disease spreading. In fact, fluctuations are intrinsic components in stochastic formulations and their relevance increases with vanishing van Kampen (1992), a much more realistic scenario in modeling emergence of small scale epidemics of communicable infectious diseases Heesterbeek et al. (2015). In this approach, transition probabilities among configurations of agents are expressed by the transition matrix , and take place during the time interval Reichl (1998); Alcaraz et al. (1994); Alcaraz and Rittenberg (2008). The disease transmission and agent recovery are modeled as probability vector of a Markov process with time interval taken to be small enough so that only one recovery or one transmission event takes place in the entire population. This is compatible with the Poissonian assumption Pastor-Satorras et al. (2015).

As a closing remark, one notices the transition matrix is analogous to the time evolution operator in quantum theories. As a result, the eigenvalues and eigenvectors of express the time evolution of agent-configurations. Despite the striking similarity between both operators, is often asymmetric, restraining its use to small values of for epidemic models in numerical simulations or introducing severe obstacles for analytic solutions. These issues hinder the systematic development of perturbation theories for agent-based epidemic models, often requiring fresh simulations to forecast the impact of small variations of the parameters of the model, contrary to the rationale behind perturbation theory.

Here, we first briefly review results Nakamura et al. (2016) derived from the squared norm of the probability vector . The formalism proposed therein allow us to further explore the operatorial content of the corresponding Markov process. More specifically, we demonstrate the time evolution is also achieved in a framework that only requires eigenvalues derived from Hermitian operators. This leads to a constrained multivariate equation, which can be solved by standard optimization techniques. We explore the fact the time evolution heavily relies on the eigendecomposition and formulate perturbative corrections to epidemic models. Additionally, we also discuss methods usually employed in quantum many-body problems and Statistical Physics, such as the the Bethe-Peierls approximation Bethe (1935) and the Holstein-Primakoff transformations Holstein and Primakoff (1940); Emary and Brandes (2003a, b).

## I SIS model

One of the key aspect of agent-based models is the use of networks to describe heterogeneity among the distinct agents Eames and Keeling (2002); Keeling and Eames (2005); Bansal et al. (2007); Pastor-Satorras et al. (2015), as Fig. 1 depicts. By taking into account the individuality of each agent, the web of connections between them creates disease spreading patterns, due to the development of characteristic pathogen mobility within the underlying population. Epidemic time duration is shorter for populations consisting of loose connected agents, whereas the potential of disease dissemination is expected to be stronger for hub-like agents.

Networks are characterized by their topological quantities such as degree, connectivity, centrality, etc., and we label the collection of these descriptive quantities as . As long as they share the same set of topological characteristic values , several distinct objects may in fact represent the same network . Graphs are natural candidates to represent networks Albert and Barabási (2002) since they are mathematical constructions formed by interconnected vertices (). For each graph, the adjacency matrix () describes all present connections among vertices; the matrix elements are , if vertices and are connected or vanish otherwise. In this context, an ensemble of graphs is a convenient way to represent a network, i.e., one graph is a single realization of the network , whose set of topological quantities are .

The correspondence between networks and agent-based models requires that each vertex contains exactly one agent, which state belongs a discrete domain. More precisely, holds the discrete state -th agent and reproduce the connection among agents. Once agents and their interconnections are properly written, we address the -agent configuration health state. Let the current health status of the -th agent be , which for may acquire two values, namely, (susceptible) or (infected). The vector

(2) |

with , describes the health status configuration of agent. Since there are two health states available per agent, the total number of distinct configurations is . We enumerate each configuration using binary arithmetic:

(3) |

For instance, the configuration containing only healthy individuals is , whereas the configuration where only the agent is infected is . Henceforth, we set the following notation: Latin integer indices run over agents , while Greek integer indices enumerate configurations .

The set formed by configurations spans the finite vector space . Within , one may define the relevant operators for epidemic spreading. The operator shows whether the -th agent is infected () or not (), namely,

(4) |

The number of infected agents at vertex is obtained via the operator

(5) |

while the total number of infected agents in the population is

(6) |

Agent health status are flipped by the action of and :

(7) | ||||

(8) |

null otherwise. They are combined to form another . The localized operators and satisfy well-known algebraic relations. For each , form local algebra, with the following structure constants: and . However, they also exhibit local fermionic anticommutation rules, , and non-local bosonic relations, , for and . The dual fermionic-bosonic behavior is a familiar occurrence in quantum spinchains Alcaraz et al. (1988); Nakamura et al. (2016). Usually, it is advisable to select either the fermionic behavior via the Jordan-Wigner transformation Batista and Ortiz (2001) or, alternatively, the bosonic behavior via the Holstein-Primakoff transformation Holstein and Primakoff (1940). We postpone the behavior-selection as our intention in this section concerns general aspects.

For any Markov process, the probability vector represents the system and is written as

(9) |

where is the probability to find agents in the configuration , at time , subjected to probability conservation constraint,

(10) |

Another integral part of the Markov process is the transition matrix , which weight transitions among configurations in a fixed time interval , creating the temporal succession:

(11) |

The details concerning disease transmission or recovery in epidemic model are included in by considering operators that act over the -agent configurations . In the SIS model, an infected agent at vertex is subjected to three distinct outcomes after the action of : transmit the disease to a susceptible connected agent; recover to the susceptible state; or remain unchanged.

The recovery event for the -th agent is executed by the operator . Brief inspection shows the action is quite simple: if the -th agent is currently infected, the health status flips to susceptible. On contrary, if the -th agent is susceptible, it returns the null vector. Although the recovery event does not involve the underlying network, the disease transmission event does. As a result, the corresponding operator transmits the disease from the -th agent to -th agent. Similarly to the recovery process, the operator only checks whether the -th agent is infected. The difference appears due to the adjacency matrix and the infection operator . Note that if the -th agent is already infected, the operatorial action returns the null vector as well. Finally, the event to remain unchanged requires diagonal operators and accounts for all the other possible non-diagonal events. The operator which provides the number of available outcomes of disease transmission for -th agent is ; whereas the operator that accounts for the chance to not recover is simply .

Under the Poissonian assumption, one only considers a single recovery or a single infection event, per time interval . Under this circumstances, the transition matrix reads

(12) |

with . The diagonal components are the probabilities for the configuration to remain unchanged after one time step. Disease spreading explicitly carries the network topology due to the contribution of . More importantly, the transition matrix contains non-Hermitian operators and, hence, its left and right eigenvectors, and , respectively, are not related by Hermitian transposition.

We emphasize that the construction of Eq. (12) considers only a single graph. In general, agent-based models are built under the hypothesis of . The reasoning behind this choice lies in the network averaging process. If the graph is large enough , one expects to recover the network topological quantities within a single realization. This statement is the equivalent to the ergodic hypothesis, where the ensemble average over graphs is replaced by the average within a single graph (self-averaging/annealing case). This expectation is reasonable but cannot hold for moderate . In what follows, we explicitly consider the network ensemble containing graphs (quenched). For that purpose, let be the -th graph in the network ensemble. Furthermore, for each graph there is a corresponding adjacency matrix . For fixed initial condition, Eq. (11) tells us the time evolution is a linear transformation so the average over the network ensemble is estimated by , where is the transition matrix using the graph (). For the SIS model, the averaging process is tracked down to

(13) |

with corresponding to the adjacency matrix of and the bar symbol represents the ensemble average. Notice that is not restricted to the integers 0 or 1 any longer. In practice, Eq. (13) claims the estimator is a real matrix and we can safely drop to bar symbol. However, sudden changes in the network topological properties must be investigated using Eq. (13), from which one derives perturbative operators.

Since in Eq. (12) is time independent, the Taylor expansion of Eq. (11) leads to the following system of differential equations:

(14) |

where are the matrix elements of the time generator

(15) |

The operator governs the dynamics of disease spreading and whose normal modes are labeled after the eigenvalues . The formal solution to Eq. (14) is

(16) |

Clearly, the eigenvalues must satisfy , for any , vanishing only for stationary states Reichl (1998). In addition, the corresponding spectral density function depends on the couplings present in Eq. (12), namely, the disease transmission and recovery probabilities as well as the network average adjacency matrix .

## Ii Theoretical approach

One of the main goals in epidemic models is the development of methods to predict the way epidemics change when parameters are subjected to small variations. If such predictions are robust, preemptive actions to lessen the epidemic are also expected to achieve better results. In physical theories, small changes in couplings or substrate are often investigated under perturbative schemes based on simpler models, which often have known orthonormal modes. In epidemic models, however, one must work with asymmetric operators and their left and right eigenvectors in Eq. (16). Despite the complexities related to the operatorial content of , the stochastic nature of the problem ensures the conservation of total probability for any . Another relevant descriptive variable derived from is the squared norm,

(17) |

As noted in Ref. Nakamura et al. (2016), probability conservation does not enforce conservation of along time, occurring only after the system reaches complete stationarity. Therefore, may be used to assess general properties of the stochastic model during both transient and stationary phases.

Since the squared norm can only assume values in the interval , one may consider a single differential equation to investigate the time evolution of . Taking the time derivative of Eq. (17) and using Eq. (14) results in

(18) |

Unlike , the symmetrized time generator

(19) |

is Hermitian with orthonormal basis and corresponding eigenvalues , for . The eigenvalues differ from their counterparts , since are not positive semi-definite, while the complex coefficients

(20) |

are not probabilities, even though they are used to evaluate the configurational probabilities

(21) |

Figure 2 illustrates the time evolution of for arbitrary Markov process.

The spectral decomposition of in Eq. (18) produces:

(22) |

subjected to the constraint . Notice that Eq. (22) is valid for any time instant . As such, one may also integrates the expression in Eq. (22) taking into account the probability constraint with help of one Lagrange multiplier ,

(23) | ||||

(24) |

where and are the fixed initial and final time instants, respectively. Equations (23) and (24) share striking similarity with the classical action of Mechanics Goldstein (1950). Since we are only interested in stationarity and extrema states so we can neglect Eqs. (23-24) for now. In fact, the condition of vanishing derivative in Eq. (22) results in the following algebraic equation:

(25) |

where either are the coefficients corresponding to stationary states or local extrema.

Formally, Eq. (25) may be solved using optimization algorithms for constrained quadratic problems. Of course, optimization problems are still formidable numerical problems and also involves the numerical approach taken for each optimization algorithm as well. Failure to converge to correct solutions or only walk in a particular neighborhood in the solution space are among common sources of problems. Furthermore, the derivation of Eqs. (22) and (25) assumes the eigenvalues and the corresponding eigenvectors are known, which again may be a complex diagonalization problem depending on the algebraic form of .

However, the aforementioned hardships are the crucial aspects one must consider to decide whether to use Eqs. (22) and (25) or Eq. (14). The answer is very simple: is only useful if symmetries are present in but not in Nakamura et al. (2016). Additional symmetries simplify the diagonalization problem and also introduce explicit bounds in root-finding procedures. It is easy to find examples where such symmetry increase occurs. For instance, consider so that the corresponding Hermitian generator is , whose eigenvectors are grouped according to the quantum number . Therefore, if additional symmetries are available for , traditional optimization algorithms become valuable resources to solve Eq. (25) and, thus, the stationary states of Markov processes and their corresponding occurrence probabilities.

As a practical example to verify the results in Eq. (25), consider the SIS model with agents and and , in the fully connected network . This set of parameters and network topology reproduce the SI model. From Eqs. (12) and (15), the Hermitian time generator for the SIS model is

(26) | ||||

(27) | ||||

(28) |

In this case, there are four eigenvalues of relevant to the description of stationary states, namely, , , and . The trivial stationary state (none infected) is obtained setting . The stationary state where all agents are infected is obtained using , and .

## Iii Perturbative methods and approximations

The ability to predict causal effects in the disease spreading is surely desirable for any epidemic model, as it allows decision-makers to select adequate strategies to mitigate new incidence cases and funding priorities. Among them, effects caused by small perturbations in the underlying contact network are specially important for agent-based models. As Ref. Pastor-Satorras et al. (2015) states, heterogeneous population hinders analytical insights about perturbative effects. Nonetheless, Eq. (25) provides an alternative way to introduce perturbative methods and approximations to epidemic models, since it relies on the Hermitian generator . This is relevant because the standard time independent perturbation theory may be used to evaluate corrections to the quantities relevant to Eq. (25), namely, the eigenvalues and the coefficients .

For the sake of clarity, we consider a finite number of topological values () to characterize the network. In this context, a simple perturbative scheme is attained by considering the change , with for any and fixed . The quantities are expected to be complicated functions so that the variations are not completely independent. Nonetheless, they are still calculated from estimators derived from the average adjacency matrix . Therefore, one expects a corresponding perturbative matrix to be added to the average adjacency matrix:

(29) |

The details of the matrix representation are specific for each perturbation set adopted, but the relevant information lies in the coupling , as it provides a natural perturbative parameter. Now, it is a simple task to identify the perturbation in the time generator,

(30) | ||||

(31) |

A few selected cases merit further attention. First, the special case , with , recovers the effective coupling formulation in random networks. Another interesting case occurs if and are periodic regular networks with distinct periods, , respectively. Depending on the initial condition and the ratio , the perturbation may either connect all available states, or lock the time evolution in a periodic cycle.

Hermiticity is sufficient to warrant Rayleigh-Schrödinger perturbation theory and produces first order corrections to eigenvalues and eigenvectors, respectively,

(32) | ||||

(33) |

The prime indicates the sum excludes degenerate states with eigenvalue . Substituting these results into Eq. (25) and discarding second order corrections, one arrives at

(34) |

While the perturbative corrections are given by Eqs. (32) and (33), the relation in Eq. (34) shows they might not be independent. Similarly, the perturbative corrections for configurational probabilities read

(35) |

Alternatives to perturbation theory are readily available as well, since the only requirement for Eq. (25) are eigenvalues and eigenvectors of . This means analytical and numerical techniques, usually available only for quantum many-body theories, are now available to epidemic models such as the Bethe-Peierls meanfield approximation (BPA) Bethe (1935) and the bosonification Holstein and Primakoff (1940).

In the BPA, the operator is replaced by global average , producing the effective time generator

(36) |

where is the degree of -th vertex, , and . The effective generator in Eq. (36) is diagonalized by rotations around the -axis:

(37) | ||||

(38) |

Due to the main requirement , the BPA rules out its applicability during transient. For stationary states, however, the BPA provides a convenient coarse particle picture accompanied by all eigenvalues and corresponding eigenvectors.

Symmetries are crucial ingredients to reduce the complexity associated with the spectral decomposition in quantum many-body problems. In what follows, we investigate finite networks with permutation invariance to better understand the role played by finite symmetries. This is possible due to Cayley’s theorem Tinkham (2003). In the SIS model, the fully connected network exhibits the desired symmetry. This simple case is used as training grounds for non-trivial networks. The first step is to unravel the role played by quantum angular operators

(39) | ||||

(40) |

The algebraic relations are and so that form a compact Lie algebra with Casimir operator .

The SIS symmetric time generator is obtained from Eqs. (27) and (28) and expressed using Eqs. (39) and (40):

(41) |

It should be noted the appearance of global angular operators is a direct consequence of the network choice adopted here, as it captures important global properties, including rotations. For non-trivial network topologies, one must consider localized angular momentum operators as usual in many-body problems. A brief inspection of Eq. (41) shows

(42) |

meaning the eigenvalues are conserved quantities and the number of infected agents may only assume the following values for fixed value . Under these circumstances, one introduces the Holstein-Primakoff transformations, which exchange the set of angular operators for the harmonic oscillator destruction and creation operators, and , respectively, with . The transformations for the sector are

(43) | ||||

(44) |

Usually, the Holstein-Primakoff transformations are most useful when the average occupation number satisfies for . Under this condition one may expand the square root and keep the linear order in . Consequently, the analysis of epidemics suggests the employment of coherent states,

(45) |

since they satisfy the eigenvalue equation with . Here, we only consider . Another remarkable property of coherent states is that several observables are derived from the Poisson distribution. Under this scheme,

(46) | ||||

(47) | ||||

(48) | ||||

(49) |

and is the incomplete Gamma function for integer . As a closing remark, we emphasize the rationale behind this approach: since coherent states are eigenvectors of , they remain unchanged under successive actions of , making them suitable candidates to characterize epidemics as .

## Iv Conclusion

Fluctuations are integral part in stochastic processes. In compartmental approaches to epidemics, their role are underestimated when the population of infected agents is scarce and heterogeneous. However, disease spreading models describing Markov processes are limited to small population size due to asymmetric time generators, which hinder the development of novel analytical insights. This issue is addressed by employing , which provides a single multivariate equation, requiring only eigenvalues and eigenvectors of symmetric time generators. One way to exploit this fact is to divert efforts in solving the symmetric spectral decomposition. Our finding shows the development of a perturbation theory in epidemic models, with emphasis in the aspects produced by topological perturbations, as show in Eq. (31). In addition, the Bethe-Peierls approximation provides the complete eigenspectrum and eigenvectors. Within this approximation, one focus in particle-like normal modes. Finally, the Holstein-Primakoff transformations exploits the network rotation symmetry to uncover a framework with quantum oscillators. From this result, one concludes coherent states are suitable candidates to study large scale epidemics. Since coherent states can also be used to study losses, we expect them to provide further hints about the epidemic decay times.

###### Acknowledgements.

We are grateful for T. J. Arruda and F. Meloni comments during manuscript preparation. A.S.M acknowledges grants CNPq 485155/2013 and CNPq 307948/2014-5.## References

- Team (2016) WHO Ebola Response Team, ‘‘After ebola in west africa — unpredictable risks, preventable epidemics,” New England Journal of Medicine 375, 587–596 (2016), http://dx.doi.org/10.1056/NEJMsr1513109 .
- Heesterbeek et al. (2015) H. Heesterbeek, R. M. Anderson, V. Andreasen, S. Bansal, D. De Angelis, C. Dye, K. T. D. Eames, W. J. Edmunds, S. D. W. Frost, S. Funk, T. D. Hollingsworth, T. House, V. Isham, P. Klepac, J. Lessler, J. O. Lloyd-Smith, C. J. E. Metcalf, D. Mollison, L. Pellis, J. R. C. Pulliam, M. G. Roberts, and C. Viboud, “Modeling infectious disease dynamics in the complex landscape of global health,” Science 347 (2015).
- Bansal et al. (2007) S. Bansal, B. T. Grenfell, and L. A. Meyers, “When individual behaviour matters: homogeneous and network models in epidemiology,” J. R. Soc. Interface 4, 879–891 (2007).
- Eames and Keeling (2002) K. T. D. Eames and M. J. Keeling, “Modeling dynamic and network heterogeneities in the spread of sexually transmitted diseases,” Proc. Natl. Acad. Sci. USA 99, 13330–13335 (2002).
- Pastor-Satorras et al. (2015) R. Pastor-Satorras, C. Castellano, P. Van Mieghem, and A. Vespignani, “Epidemic processes in complex networks,” Rev. Mod. Phys. 87, 925–979 (2015).
- Keeling and Eames (2005) M.J. Keeling and K.T.D Eames, “Networks and epidemic models,” J. R. Soc. Interface 2, 295–307 (2005).
- Pastor-Satorras and Vespignani (2001) Romualdo Pastor-Satorras and Alessandro Vespignani, “Epidemic spreading in scale-free networks,” Phys. Rev. Lett. 86, 3200–3203 (2001).
- van Kampen (1992) N. G. van Kampen, Stochastic Processes in Physics and Chemistry, Vol. 1 (Elsevier Science, Amsterdam, 1992).
- Reichl (1998) L.E. Reichl, “A modern course in statistical physics,” (1998).
- Alcaraz et al. (1994) F.C. Alcaraz, M. Droz, M. Henkel, and V. Rittenberg, “Reaction-diffusion processes, critical dynamics, and quantum chains,” Annals of Physics 230, 250 – 302 (1994).
- Alcaraz and Rittenberg (2008) F. C. Alcaraz and V. Rittenberg, “Directed abelian algebras and their application to stochastic models,” Phys. Rev. E 78, 041126 (2008).
- Nakamura et al. (2016) G. M. Nakamura, A. C. P. Monteiro, G. C. Cardoso, and A. S. Martinez, “Fast computation method for comprehensive agent-level epidemic dissemination in networks,” ArXiv e-prints (2016), arXiv:1606.07825 [physics.soc-ph] .
- Bethe (1935) H. A. Bethe, “Statistical theory of superlattices,” Proc. R. Soc. A 150, 552–575 (1935).
- Holstein and Primakoff (1940) T. Holstein and H. Primakoff, “Field dependence of the intrinsic domain magnetization of a ferromagnet,” Phys. Rev. 58, 1098–1113 (1940).
- Emary and Brandes (2003a) Clive Emary and Tobias Brandes, “Chaos and the quantum phase transition in the dicke model,” Phys. Rev. E 67, 066203 (2003a).
- Emary and Brandes (2003b) C. Emary and T. Brandes, “Quantum chaos triggered by precursors of a quantum phase transition: The Dicke model,” Phys. Rev. Lett. 90, 044101 (2003b).
- Albert and Barabási (2002) R. Albert and A.-L. Barabási, “Statistical mechanics of complex networks,” Rev. Mod. Phys. 74, 47–97 (2002).
- Alcaraz et al. (1988) Francisco C Alcaraz, Michael N Barber, and Murray T Batchelor, “Conformal invariance, the xxz chain and the operator content of two-dimensional critical systems,” Annals of Physics 182, 280 – 343 (1988).
- Nakamura et al. (2016) G. M. Nakamura, M. Mulato, and A. S. Martinez, “Spin gap in coupled magnetic layers,” Physica A 451, 313 – 319 (2016).
- Batista and Ortiz (2001) C. D. Batista and G. Ortiz, “Generalized jordan-wigner transformations,” Phys. Rev. Lett. 86, 1082–1085 (2001).
- Goldstein (1950) H. Goldstein, Classical mechanics (Addison-Wesley, 1950).
- Nakamura et al. (2016) G. M. Nakamura, A. C. P. Monteiro, G. C. Cardoso, and A. S. Martinez, “Finite symmetries in epidemic models,” ArXiv e-prints (2016), arXiv:1609.06980 [q-bio.PE] .
- Tinkham (2003) Michael Tinkham, Group theory and quantum mechanics (Courier Corporation, 2003).