Mathematical Analysis of Four Fundamental Epidemiological Models for Monkeypox Disease Outbreaks: On the Pivotal Role of Human–Animal Order Parameters—In Memory of Hermann Haken

Abstract

Four fundamental models that describe the spread of Monkeypox disease are analyzed: the SIR-SIR, SEIR-SIR, SIR-SEIR, and SEIR-SEIR models. They form the basis of most Monkeypox diseases models that are currently discussed in the literature. It is shown that the way the model subpopulations are organized in disease outbreaks and evolve relative to each other is determined by the relevant unstable system eigenvectors, also called order parameters. For all models, analytical expressions of the order parameters are derived. Under appropriate conditions these order parameters describe the initial outbreak phases of exponential increase in good approximation. It is shown that all four models exhibit maximally two order parameters and maximally one human–animal order parameter. The human–animal order parameter firmly connects the outbreak dynamics in the animal system with the dynamics in the human system. For the special case of the SIR-SIR model, it is found that the two possible order parameters completely describe the dynamics of infected humans and animals during entire infection waves. Finally, a simulation of a Monkeypox infection wave illustrates that in line with the aforementioned analytical results the leading order parameter explains most of the variance in the infection dynamics.

1. Introduction

Monkeypox (Mpox) disease is an infectious disease that is endemic in several African countries [1]. During the last five decades, Mpox infection waves have repeatedly occurred in those countries [2,3]. In this context, understanding the initial outbreak phases of Mpox infection waves is of particular importance because these initial phases offer the opportunity for ad hoc interventions that may dramatically reduce the infection dynamics [4]. In particular, in line with previous works on COVID-19 waves [4], a vital step to understand the emergence of Mpox infection waves is to conduct model-based analyses that determine the initial organization of such waves. Insights obtained from such endeavors are not only relevant for the aforementioned African countries. Rather, due to international travel, Mpox infection waves can spread out to non-endemic countries as it was observed recently during the global 2022–2023 Mpox epidemic [5,6].

Mpox disease comes with symptoms like fever, rash, sore throat, and respiratory distress [7] and can lead to death [1,2,3]. The disease is caused by the Monkeypox virus that is transmitted from animals to humans [1,8]. In particular, rodents such as squirrels, rats, and mice can carry the virus but also certain monkeys (whence the name Monkeypox) [8]. However, once the virus has invaded a particular human population it can also spread out within the population by means of human to human transmissions [1]. Temporal limited outbreaks of Monkeypox disease have been observed since 1970 [1,2,3]. Such waves are typically observed in the human population, whereas the infection dynamics in the animal population is unobserved. During the period 2000–2009, when taking all observed disease outbreaks together, it is estimated that there had been about 10,000 cases of Monkeypox disease. For the period of 2010–2019, this estimate is higher at about 19,000 cases [3]. In particular, from September 2017 to April 2018 an infection wave spread out through several regions of Nigeria [2,8,9]. The wave in the human population was probably triggered by multiple, independent animal to human transmissions of the Monkeypox virus [8]. In addition, it was confirmed that human to human transmissions were involved in the disease outbreak [9]. Since January 2023, cases of Monkeypox disease started to rise dramatically in the Democratic Republic of the Congo (DRC). As of August 2024, the infection surge is still ongoing [10,11] and is about to spread out to several neighboring countries of the DRC [12]. Within the DRC the surge reached the 20,000 mark of suspected cases [13] that had never been reached before. In view of those developments, on 14 August 2024, the WHO announced its highest level of alert for Mpox and declared the Mpox outbreak as a health emergency of international concern [12]. Finally, as mentioned above, in several instances international travel brought the Mpox disease to non-endemic countries [2,3]. Most dramatically, during the period of 2022–2023 Mpox disease spread out globally (primarily in South, Middle, and North America and Europe [5,6]) and produced almost 100,000 infected cases worldwide [14]. While the 2022–2023 global epidemic eventually subsided due to intervention and prevention measures and in the absence of an animal reservoir, the situation in the endemic African countries and, in particular, the current DRC-centered outbreak are highly unpredictable.

For the time being, Mpox disease waves in endemic African countries triggered by waves in the respective animal reservoirs are likely to occur and they may or may not spill over to non-endemic countries [11].

Several efforts have been made to describe the emergence of Mpox disease in human–animal systems with the help of epidemiological models [15]. At the heart of these efforts are models that describe susceptible (S), infected (I), and recovered (R) individuals [4,16]. Such SIR models can describe both the human and animal populations of interest. Adding exposed (E) individuals that have been infected but are not yet infectious leads to SEIR models [4,16] that again may be used to describe both the human and animal populations of interest. Combining these two types of models, SIR and SEIR, we arrive at four fundamental models: SIR-SIR, SEIR-SIR, SIR-SEIR, and SEIR-SEIR models, where the first acronym describes the animal system, while the latter refers to the human system. An SIR-SIR Monkeypox model has been studied by Bhunu and Mushayabasa [17]. SIR-SIR Monkeypox models were also used by Emeka et al. [18] and Somma et al. [19] and have been slightly generalized to take vaccinated [18] and quarantined [19] individuals into account, respectively. Madubueze et al. [20] based their Monkeypox model on an SIR-SEIR system and fine-grained the human SEIR infection dynamics, among other things, by taking again the possibility of vaccination into account. Usman and Adamu [21], Peter et al. [22], Bankuru et al. [23], and Collins and Duffy [24] used Mpox SEIR-SEIR approaches featuring additional groups of quarantined [22] or vaccinated [21,23] individuals. A plenitude of highly detailed epidemiological models for the spread of Mpox disease in human–animal systems has been proposed in the literature. Such models include, for example, individuals with different degree of infectiousness [25,26], quarantined and isolated individuals [27], clinically ill or hospitalized individuals [25,26,28], and detected versus undetected cases [29]. All these models have in common that they are based on one of the aforementioned four fundamental models.

Despite these modeling efforts, what is missing is an analysis of the multi-compartmental components involved in Mpox outbreaks. Such multi-compartmental components describe compartments bound together to entities and are determined by stable and unstable eigenvectors of the respective human–animal systems. Multi-compartmental component analyses have been conducted for various epidemiological models describing COVID-19 outbreaks during the COVID-19 pandemic [4,30] and several virus dynamics models describing SARS-CoV-2 infections [31] and the human immune reaction [32]. In particular, the leading components that describe the initial organizations of COVID-19 outbreaks in populations have been determined in terms of so-called order parameters [33,34,35] for SIR systems, SEIR systems, and some higher-dimensional models [4]. As mentioned above, a comparable analysis is missing in the modeling literature on Mpox disease outbreaks. That is, the above reviewed studies on epidemiological Mpox disease models did not present any analysis of relevant Mpox disease order parameters. Therefore, the current study supplements the existing literature and adds a novel aspect to it. In the context of the COVID-19 pandemic, benefits of the order parameter perspective have been demonstrated. For example, it has been demonstrated that the impact of intervention measures can be conveniently analyzed when focusing on the leading building-blocks or order parameters of COVID-19 waves. Explicitly, using this approach, COVID-19 waves observed in the USA [36], Europe [37], China [38], Thailand [39], and Pakistan [40] have been analyzed. Interestingly, all systems investigated in those studies on COVID-19 waves exhibited only a single leading multi-compartmental component or order parameter. In contrast, as it will be shown below, Mpox infection dynamics in general is characterized by two order parameters. In anticipation of this novel aspect and in view of the absence of studies devoted to determine Mpox order parameters, the aim of the current study is to identify and compute the leading organizational elements or order parameters of the four fundamental models listed above and to interpret their qualitative and quantitative aspects. The aim is to demonstrate that they dominate and determine the initial phase dynamics of Mpox outbreaks. In doing so, it will be determined how compartments or subpopulations evolve relative to each other. In addition, the objective is to identify the remaining multi-compartmental components given in terms of (neutrally) stable eigenvectors and to explore their roles as well for the enfolding of Mpox infection waves.

The remainder of this study is structured as follows. The four fundamental models will be introduced in Section 2. The models will be analyzed in Section 3. In particular, key results regarding the multi-compartmental building-blocks in terms of order parameters and their amplitudes will be obtained in Section 3.1 and implications of those results will be discussed in Section 3.2. Section 3.3 will briefly exemplify that under certain conditions the dynamics along the aforementioned remaining eigenvectors makes essential contributions to the overall dynamics as well. In Section 3.4, some findings will be illustrated with the help of a simulated Mpox infection wave. Some conclusions will be drawn in Section 4. Certain limitations of the current study will be addressed in this section as well.

2. Methods

2.1. Four Fundamental Models

Let us define the four fundamental epidemiological models for Mpox infection dynamics.

2.1.1. Model A (SIR-SIR Model)

Let $S_a$, $I_a$, and $R_a$ denote the number of individuals in the subpopulations of susceptible, infected, and recovered animals, respectively. Likewise, let $S_h$, $I_h$, and $R_h$ denote the number of humans in the susceptible, infected, and recovered subpopulations, respectively. The evolution equations of the population variables read [17]

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{S}_a& =& -{\displaystyle \frac{{\beta }_a}{N_a}}{I}_a{S}_a,\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{I}_a& =& {\displaystyle \frac{{\beta }_a}{N_a}}{I}_a{S}_a-{\gamma }_a{I}_a,\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{R}_a& =& {\gamma }_a{I}_a,\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{S}_h& =& -\left({\displaystyle \frac{{\beta }_h}{N_h}}{I}_h+{\displaystyle \frac{{\beta }_{12}}{N_a}}{I}_a\right)S_h,\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{I}_h& =& \left({\displaystyle \frac{{\beta }_h}{N_h}}{I}_h+{\displaystyle \frac{{\beta }_{12}}{N_a}}{I}_a\right)S_h-{\gamma }_h{I}_h,\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{R}_h& =& {\gamma }_h{I}_h,\hfill \end{array}$$

(1)

where $N_a=S_a+I_a+R_a$ and $N_a=S_a+I_a+R_a$ denote the total animal and human populations, respectively. In Equation (1), ${\beta }_a$, ${\beta }_{12}$, and ${\beta }_h$ denote the effective contact rates (also called infection rates or transmission rates) of animal to animal, animal to human, and human to human transmissions, respectively. The parameters ${\gamma }_a$ and ${\gamma }_h$ denote recovery rates of infected animal and human individuals, respectively. The current study focuses on infection waves that take place on relatively short durations of a few months such that birth and death processes can be neglected. Moreover, deaths due to Mpox disease are neglected as well. The populations are assumed to be constant at least in good approximation such that variations in $N_a$ and $N_h$ can be neglected. Consequently, the evolution equations for $R_a$ and $R_h$ can be eliminated by putting $R_a=N_a-S_a-I_a$ and $R_h=N_h-S_h-I_h$. Moreover, in the current study, for the sake of simplicity, epidemiological models will be formulated in rescaled variables (i.e., with the help of fractions) as defined by [41]

$$s_j=S_j/N_j,i_j=I_j/N_j,r_j=R_j/N_j$$

(2)

for $j=a,h$. The SIR-SIR model (model A) defined by Equation (1) then becomes

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{s}_a& =& -{\beta }_a{i}_a{s}_a,\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{i}_a& =& {\beta }_a{i}_a{s}_a-{\gamma }_a{i}_a,\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{s}_h& =& -\left({\beta }_h{i}_h+{\beta }_{12}{i}_a\right)s_h,\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{i}_h& =& \left({\beta }_h{i}_h+{\beta }_{12}{i}_a\right)s_h-{\gamma }_h{i}_h\hfill \end{array}$$

(3)

with $r_a=1-s_a-i_a$ and $r_h=1-s_h-i_h$. The state vector of model A reads $\mathbf{x}=(s_a,i_a,s_h,i_h)$. As mentioned in the introduction, Bhunu and Mushayabasa [17], Emeka et al. [18], and Somma et al. [19] based their studies on SIR-SIR models as described by Equation (1).

2.1.2. Model B (SEIR-SIR Model)

The SEIR-SIR model involves the class of exposed animals ($E_a$) that are infected but not yet infectious. The relative size of the exposed animal subpopulation is denoted by $e_a=E_a/N_a$. The evolution equations for the rescaled model variables read

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{s}_a& =& -{\beta }_a{i}_a{s}_a,\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{e}_a& =& {\beta }_a{i}_a{s}_a-{\alpha }_a{e}_a,\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{i}_a& =& {\alpha }_a{e}_a-{\gamma }_a{i}_a,\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{s}_h& =& -\left({\beta }_h{i}_h+{\beta }_{12}{i}_a\right)s_h,\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{i}_h& =& \left({\beta }_h{i}_h+{\beta }_{12}{i}_a\right)s_h-{\gamma }_h{i}_h\hfill \end{array}$$

(4)

and involve in addition to the SIR-SIR model parameters the parameter ${\alpha }_a$, which describes the transition rate of animals from being infected and non-infectious to being infected and infectious. For the fraction variables of recovered animal and human individuals, the following relations hold: $r_a=1-s_a-e_a-i_a$ and $r_h=1-s_h-i_h$. The state vector of the SEIR-SIR model defined by Equation (4) reads $\mathbf{x}=(s_a,e_a,i_a,s_h,i_h)$. Note that in what follows we will refer to the population variables $i_a$ and $i_h$ as infectious populations rather than infected populations in order to highlight the distinction between the variables $e_j$ and $i_j$ for $j=a,h$ that describe both infected individuals.

2.1.3. Model C (SIR-SEIR Model)

Like model B, model C involves a class of exposed individuals. Unlike model B, in model C this class ($E_h$) shows up in the human subsystem. The relative size is denoted by $e_h=E_h/N_h$. The evolution equations of the rescaled model variables read

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{s}_a& =& -{\beta }_a{i}_a{s}_a,\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{i}_a& =& {\beta }_a{i}_a{s}_a-{\gamma }_a{i}_a,\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{s}_h& =& -\left({\beta }_h{i}_h+{\beta }_{12}{i}_a\right)s_h,\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{e}_h& =& \left({\beta }_h{i}_h+{\beta }_{12}{i}_a\right)i_h-{\alpha }_h{e}_h,\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{s}_h& =& {\alpha }_h{e}_h-{\gamma }_h{i}_h\hfill \end{array}$$

(5)

and are complemented by the relations $r_a=1-s_a-i_a$ and $r_h=1-s_h-e_h-i_h$. In Equation (5) the parameter ${\alpha }_h$ describes the transition rate of infected but non-infectious humans (i.e., exposed humans) to become infectious. The model state vector reads $\mathbf{x}=(s_a,i_a,s_h,e_h,i_h)$. As mentioned in the introduction, Madubueze et al. [20] based their study on a Monkeypox model of the SIR-SEIR system type.

2.1.4. Model D (SEIR-SEIR Model)

The SEIR-SEIR model (model D) takes the possibility of exposed individuals into account both for the human and animal populations of interest. Accordingly, the model reads

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{s}_a& =& -{\beta }_a{i}_a{s}_a,\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{e}_a& =& {\beta }_a{i}_a{s}_a-{\alpha }_a{e}_a,\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{i}_a& =& {\alpha }_a{e}_a-{\gamma }_a{i}_a,\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{s}_h& =& -\left({\beta }_h{i}_h+{\beta }_{12}{i}_a\right)s_h,\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{e}_h& =& \left({\beta }_h{i}_h+{\beta }_{12}{i}_a\right)i_h-{\alpha }_h{e}_h,\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{s}_h& =& {\alpha }_h{e}_h-{\gamma }_h{i}_h\hfill \end{array}$$

(6)

with $r_a=1-s_a-e_a-i_a$ and $r_h=1-s_h-e_h-i_h$ and all parameters as defined for the previous models A, B, C. The state vector of the SEIR-SEIR model defined by Equation (6) reads $\mathbf{x}=(s_a,e_a,i_a,s_h,e_h,i_h)$. Various authors have used the SEIR-SEIR model defined by Equation (6) as a departure point to model the spread of Mpox disease [21,22,23,24] (see also the Introduction).

2.2. Admissible Model Parameters

In what follows it is assumed that all model parameters are positive like

$${\beta }_a,{\beta }_h,{\beta }_{12},{\gamma }_a,{\gamma }_h,{\alpha }_a,{\alpha }_h>0.$$

(7)

That is, degenerated special cases of models exhibiting vanishing model parameters are ignored.

3. Results

3.1. Potential Order Parameters

Order parameters are known to describe the emerging order of multi-component systems in various disciplines [33,35] and have recently been determined for various epidemiological models describing the COVID-19 pandemic [4]. In the context of the current study, it is useful to discuss order parameters in the context of potential order parameters—as they will be defined next.

Definition 1.

A potential order parameter is an eigenvector associated with an eigenvalue that for appropriate model parameters assumes positive values.

In other words, in general, dynamical systems exhibit two types of eigenvalues: eigenvalues that are negative or zero in any case and eigenvalues that for appropriate model parameters become positive (but may become negative or zero for other parameter values). Eigenvectors associated with the second type of eigenvalues will be called potential order parameters. The implication of the definition is that for those model parameters that make the eigenvalue in fact positive the potential order parameter is an order parameter or unstable eigenvector [33,35]. In contrast, for model parameters that make the eigenvalue negative the potential order parameter corresponds to a stable eigenvector [33,35].

Theorem 1.

All four models A–D defined by Equations (3)–(6) exhibit two potential order parameters for the admissible model parameters listed in Section 2.2.

Proof. 

It is sufficient to show that each model exhibits two eigenvalues that for appropriate, admissible model parameters can become positive.

The model A (SIR-SIR model) exhibits the disease-free fixed points $s_a^{*}\in [0,1]$, $i_a^{*}=0$, $s_h^{*}\in [0,1]$, $i_h^{*}=0$. Let us introduce the following two relative variables describing differences with respect to fixed point variables ${\delta }_j=s_j-s_j^{*}$ for $j=a,h$. Note that $i_a$ and $i_h$ already can be regarded as relative variables because the respective fixed point values read $i_a^{*}=0$ and $i_h^{*}=0$. The model A state vector of relative variables reads $\mathbf{u}=({\delta }_a,i_a,{\delta }_h,i_h)$. The eigenvalues of model A are obtained from a linear stability analysis for which it is assumed that all entries of the relative state vector are small quantities (i.e., the state of the epidemiological system is close to a fixed point of the system). Then, Equation (3) becomes a linear evolution equation of the form [4,33]

$${\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\mathbf{u}=L\mathbf{u}$$

(8)

with the linearization matrix L defined by

$$L=\left(\begin{array}{cccc}0& -{\beta }_a^{*}& 0& 0\\ 0& {\beta }_a^{*}-{\gamma }_a& 0& 0\\ 0& -{\beta }_{12}^{*}& 0& -{\beta }_h^{*}\\ 0& {\beta }_{12}^{*}& 0& {\beta }_h^{*}-{\gamma }_h\end{array}\right)$$

(9)

and ${\beta }_a^{*}={\beta }_a{s}_a^{*}$, ${\beta }_h^{*}={\beta }_h{s}_h^{*}$, and ${\beta }_{12}^{*}={\beta }_{12}{s}_h^{*}$. The matrix L exhibits the eigenvalues ${\lambda }_1=0$, ${\lambda }_2={\beta }_a^{*}-{\gamma }_a$, ${\lambda }_3=0$, and ${\lambda }_4={\beta }_h^{*}-{\gamma }_h$. The pair of eigenvalues ${\lambda }_1$ and ${\lambda }_2$ (and likewise the pair ${\lambda }_3$ and ${\lambda }_4$) are known as eigenvalues of SIR systems [4,42]. As will be shown below in the context of Theorem 2, the eigenvalues ${\lambda }_3$ and ${\lambda }_4$ can indeed be interpreted as eigenvalues of the SIR human subsystem. In contrast, while ${\lambda }_1$ and ${\lambda }_2$ formally correspond to SIR system eigenvalues such that it would be tempting to interpret them as SIR animal subsystem eigenvalues, as it will be shown below, at least for ${\lambda }_2$ such an interpretation is misleading. Rather, it is more appropriate to interpret ${\lambda }_2$ (in line with its derivation from matrix (9)) as an eigenvalue of the entire SIR-SIR system. ${\lambda }_2$ and ${\lambda }_4$ for appropriate, admissible model parameters (namely, ${\beta }_a^{*}>{\gamma }_a>0$ and ${\beta }_h^{*}>{\gamma }_h>0$, respectively) assume positive values. Therefore, the corresponding eigenvectors are the potential order parameters of the SIR-SIR system. They will be derived explicitly below.

For model B (SEIR-SIR model) the fixed points are given by $s_a^{*}\in [0,1]$, $e_a^{*}=0$, $i_a^{*}=0$, $s_h^{*}\in [0,1]$, $i_h^{*}=0$. The state vector of relative variables with respect to those fixed points reads $\mathbf{u}=({\delta }_a,e_a,i_a,{\delta }_h,i_h)$. For states sufficiently close to a given fixed point Equation (4) can be linearized and becomes Equation (8) with L defined by

$$L=\left(\begin{array}{ccccc}0& 0& -{\beta }_a^{*}& 0& 0\\ 0& -{\alpha }_a& {\beta }_a^{*}& 0& 0\\ 0& {\alpha }_a& -{\gamma }_a& 0& 0\\ 0& 0& -{\beta }_{12}^{*}& 0& -{\beta }_h^{*}\\ 0& 0& {\beta }_{12}^{*}& 0& {\beta }_h^{*}-{\gamma }_h\end{array}\right).$$

(10)

The matrix L exhibits the eigenvalues ${\lambda }_1=0$ and

$${\lambda }_{2,3}=-{\displaystyle \frac{{\gamma }_a+{\alpha }_a}{2}}\pm \sqrt{{\displaystyle \frac{{({\gamma }_a+{\alpha }_a)}^2}{4}}+{\alpha }_a({\beta }_a^{*}-{\gamma }_a)},$$

(11)

where the upper (lower) sign holds for ${\lambda }_2$ (${\lambda }_3$). These eigenvalues are known as eigenvalues of epidemiological SEIR models [4,42]. However, as will be shown in the context of Theorem 2 with respect to model B they actually describe eigenvalues of the entire SEIR-SIR model. The two remaining eigenvalues of L read ${\lambda }_4=0$ and ${\lambda }_5={\beta }_h^{*}-{\gamma }_h$ and denote SIR-human subsystem eigenvalues (as will be shown below). ${\lambda }_2$ assumes positive values for admissible parameters if ${\beta }_a^{*}>{\gamma }_a>0$ holds, whereas ${\lambda }_3$ is negative for all admissible model parameters [4]. ${\lambda }_5$ is positive for admissible parameters if ${\beta }_h^{*}>{\gamma }_h>0$ holds. In summary, ${\lambda }_2$ and ${\lambda }_5$ may assume positive values for admissible model parameters. By Definition 1, the corresponding eigenvectors describe the potential order parameters of the epidemiological SEIR-SIR system (4). They will be derived below.

In the case of model C (SIR-SEIR model) the fixed points are given by $s_a^{*}\in [0,1]$, $i_a^{*}=0$, $s_h^{*}\in [0,1]$, $e_h^{*}=0$, and $i_h^{*}=0$. Accordingly, the state vector of relative variables reads $\mathbf{u}=({\delta }_a,i_a,{\delta }_h,e_h,i_h)$ and satisfies Equation (8), which is again close to the aforementioned fixed points. For model C the linearization matrix L reads

$$L=\left(\begin{array}{ccccc}0& -{\beta }_a^{*}& 0& 0& 0\\ 0& {\beta }_a^{*}-{\gamma }_a& 0& 0\\ 0& -{\beta }_{12}^{*}& 0& 0& -{\beta }_h^{*}\\ 0& {\beta }_{12}^{*}& 0& -{\alpha }_h& {\beta }_h^{*}\\ 0& 0& 0& {\alpha }_h& -{\gamma }_h\end{array}\right).$$

(12)

A detailed calculation shows that the five eigenvalues of L are given in terms of three SEIR-human subsystem eigenvalues and two eigenvalues that formally look like SIR system eigenvalues. The SIR-system-like eigenvalues read ${\lambda }_1=0$ and ${\lambda }_2={\beta }_a^{*}-{\gamma }_a$. The SEIR-human subsystem eigenvalues read ${\lambda }_3=0$ and

$${\lambda }_{4,5}=-{\displaystyle \frac{{\gamma }_h+{\alpha }_h}{2}}\pm \sqrt{{\displaystyle \frac{{({\gamma }_h+{\alpha }_h)}^2}{4}}+{\alpha }_h({\beta }_h^{*}-{\gamma }_h)},$$

(13)

where the upper (lower) sign holds for ${\lambda }_4$ (${\lambda }_5$). ${\lambda }_2$ is positive if ${\beta }_a^{*}>{\gamma }_a>0$, ${\lambda }_4$ assumes positive values if ${\beta }_h^{*}>{\gamma }_h>0$, and ${\lambda }_5$ is negative for all admissible model parameters. In summary, only two eigenvalues may become positive: ${\lambda }_2$ and ${\lambda }_4$. The corresponding eigenvectors denote the potential order parameters of the SIR-SEIR system (5).

Finally, the fixed points of model D (SEIR-SEIR model) are described by $s_a^{*}\in [0,1]$, $e_a=0$, $i_a^{*}=0$, $s_h^{*}\in [0,1]$, $e_h=0$, and $i_h^{*}=0$. The model D state vector of relative variables reads $\mathbf{u}=({\delta }_a,e_a,i_a,{\delta }_h,e_h,i_h)$ and satisfies Equation (8) in the linear domain dynamics with L given by

$$L=\left(\begin{array}{cccccc}0& 0& -{\beta }_a^{*}& 0& 0& 0\\ 0& -{\alpha }_a& {\beta }_a^{*}& 0& 0& 0\\ 0& {\alpha }_a& -{\gamma }_a& 0& 0& 0\\ 0& 0& -{\beta }_{12}^{*}& 0& 0& -{\beta }_h^{*}\\ 0& 0& {\beta }_{12}^{*}& 0& -{\alpha }_h& {\beta }_h^{*}\\ 0& 0& 0& 0& {\alpha }_h& -{\gamma }_h\end{array}\right).$$

(14)

It can be shown that the six eigenvalues of L are given in terms of three SEIR eigenvalues for the human subsystem and three eigenvalues that again at least formally look like SEIR eigenvalues and only involve animal subsystem parameters. They are listed above already in the context of the SEIR-SIR and SIR-SEIR models. For the sake of clarity they are explicitly listed here. Those related to the animal subsystem parameters read ${\lambda }_1=0$ and

$${\lambda }_{2,3}=-{\displaystyle \frac{{\gamma }_a+{\alpha }_a}{2}}\pm \sqrt{{\displaystyle \frac{{({\gamma }_a+{\alpha }_a)}^2}{4}}+{\alpha }_a({\beta }_a^{*}-{\gamma }_a)},$$

(15)

where the upper (lower) sign holds for ${\lambda }_2$ (${\lambda }_3$). The three SEIR-human subsystem eigenvalues read ${\lambda }_4=0$ and

$${\lambda }_{5,6}=-{\displaystyle \frac{{\gamma }_h+{\alpha }_h}{2}}\pm \sqrt{{\displaystyle \frac{{({\gamma }_h+{\alpha }_h)}^2}{4}}+{\alpha }_h({\beta }_h^{*}-{\gamma }_h)},$$

(16)

where again the upper (lower) sign holds for ${\lambda }_5$ (${\lambda }_6$). In line with the previous discussion, it follows that ${\lambda }_3$ and ${\lambda }_6$ are always negative. In contrast, ${\lambda }_2$ and ${\lambda }_5$ are positive if ${\beta }_a^{*}-{\gamma }_a>0$ and ${\beta }_h^{*}-{\gamma }_h>0$ holds, respectively. In other words, for model D, ${\lambda }_2$ and ${\lambda }_5$ may assume positive values and the corresponding eigenvectors denote the potential order parameters of the SEIR-SEIR system (6).

In summary, all four fundamental models A, B, C, and D exhibit two potential order parameters. □

Corollary 1.

The four fundamental models A, B, C, and D defined by Equations (3)–(6) exhibit maximally two order parameters for the admissible model parameters listed in Section 2.2.

Proof. 

If the eigenvalues of both potential order parameters assume positive values, then the model under consideration exhibits two order parameters. This is the maximal number. □

Definition 2.

A human–animal order parameter is an order parameter (unstable eigenvector) that exhibits non-vanishing components both in the animal and human subsystems.

As we will see below there are eigenvectors that only exhibit components in either the animal subsystem or the human subsystem. Their coordinates (or amplitudes) describe dynamics in either of the two systems. In contrast, eigenvectors that exhibit components in both system connect the animal subsystem dynamics with the human subsystem dynamics. In particular, unstable eigenvectors (order parameters) with that property are of interest because they describe the emerging order of an Mpox outbreak in terms of multi-compartmental components that link both subsystems with each other. By Definition 2 these multi-compartmental components will be called human–animal order parameters.

Theorem 2.

All four fundamental models A, B, C, and D defined by Equations (3)–(6) exhibit maximally one human–animal order parameter for the admissible model parameters listed in Section 2.2.

Proof. 

In what follows, the eigenvectors of all four models will be derived and, in doing so, a constructive proof of Theorem 2 will be given. Eigenvectors will be denoted by ${\mathbf{v}}_j$ with $j=1,\dots ,4$ for model A, $j=1,\dots ,5$ for models B and C, and $j=1,\dots ,6$ for model D.

For model A, from matrix (9) it follows that eigenvectors related to ${\lambda }_1=0$ and ${\lambda }_3$ read ${\mathbf{v}}_1=(1,0,0,0)$ and ${\mathbf{v}}_3=(0,0,1,0)$. The eigenvector of ${\lambda }_4={\beta }_h^{*}-{\gamma }_h$ reads

$$\begin{array}{c}\hfill {\mathbf{v}}_4={\displaystyle \frac{1}{Z_4}}\left(\begin{array}{c}0\hfill \\ 0\hfill \\ -{\beta }_h^{*}\hfill \\ {\lambda }_4\hfill \end{array}\right)\end{array}$$

(17)

with $Z_4=\sqrt{{({\beta }_h^{*})}^2+{\lambda }_4^2}$. In order to derive ${\mathbf{v}}_2$ associated to ${\lambda }_2={\beta }_a^{*}-{\gamma }_a$ note that ${\mathbf{v}}_2$ satisfies

$$\begin{array}{c}\hfill \left(\begin{array}{cccc}-{\lambda }_2& -{\beta }_a^{*}& 0& 0\\ 0& 0& 0& 0\\ 0& -{\beta }_{12}^{*}& -{\lambda }_2& -{\beta }_h^{*}\\ 0& {\beta }_{12}^{*}& 0& {\lambda }_4-{\lambda }_2\end{array}\right){\mathbf{v}}_2=0\end{array}$$

(18)

Consequently, the first two components of ${\mathbf{v}}_2$ have the same structure as the two non-vanishing components of ${\mathbf{v}}_4$ such that ${\mathbf{v}}_2=(-{\beta }_a^{*},{\lambda }_2,a,b)/Z_2$ holds, where a and b (and $Z_2$) are still to be determined. Substituting this ansatz into Equation (18), and exploiting the third and fourth rows of the matrix equation, yields a and b (and $Z_2$). The result reads

$$\begin{array}{c}\hfill {\mathbf{v}}_2={\displaystyle \frac{1}{Z_2}}\left(\begin{array}{c}-{\beta }_a^{*}\hfill \\ {\lambda }_2\hfill \\ -{\beta }_{12}^{*}\left(1+{\displaystyle \frac{{\beta }_h^{*}}{{\lambda }_2-{\lambda }_4}}\right)\hfill \\ {\beta }_{12}^{*}{\displaystyle \frac{{\lambda }_2}{{\lambda }_2-{\lambda }_4}}\hfill \end{array}\right)\end{array}$$

(19)

with $Z_2=\sqrt{{({\beta }_a^{*})}^2+{\lambda }_2^2+{({\beta }_{12}^{*})}^2\xi }$ and $\xi ={[1+{\beta }_h^{*}/({\lambda }_2-{\lambda }_4)]}^2+{[{\lambda }_2/({\lambda }_2-{\lambda }_4)]}^2$. The eigenvectors ${\mathbf{v}}_2$ and ${\mathbf{v}}_4$ correspond to the potential order parameters of the SIR-SIR model (3). However, ${\mathbf{v}}_4$ exhibits non-vanishing components only in the human subsystem. In contrast, ${\mathbf{v}}_2$ exhibits non-vanishing components in both subsystems and, consequently, can describe an emerging order involving both subsystems. For ${\lambda }_2>0$, it follows from Definition 2 that ${\mathbf{v}}_2$ is the human–animal order parameter of the SIR-SIR system. Accordingly, the corresponding eigenvalue ${\lambda }_2$ should be interpreted as eigenvalue of the entire SIR-SIR system (as anticipated above). Note that the two non-vanishing components of ${\mathbf{v}}_4$ are known as components of the order parameter of SIR models [4]. Accordingly, for ${\lambda }_4>0$, on the one hand, the eigenvector ${\mathbf{v}}_4$ constitutes an order parameter of the epidemiological system (3). On the other hand, while it does not qualify as a human–animal order parameter, it may be regarded as the SIR order parameter of the human subsystem of Equation (3). Likewise, ${\lambda }_4$ may be regarded as an SIR human subsystem eigenvalue (as anticipated above)—in addition to its original role as eigenvalue of the model matrix (9) of model (3). In summary, the SEIR-SIR model (model A) maximally exhibits one human–animal order parameter.

For the SEIR-SIR model (model B) defined by Equation (4), the eigenvectors associated to ${\lambda }_1=0$ and ${\lambda }_4=0$ read ${\mathbf{v}}_1=(1,0,0,0,0)$ and ${\mathbf{v}}_4=(0,0,0,1,0)$. For ${\lambda }_5={\beta }_h^{*}-{\gamma }_h$ from the matrix (10) we obtain again an eigenvector with the components of an SIR order parameter:

$$\begin{array}{c}\hfill {\mathbf{v}}_5={\displaystyle \frac{1}{Z_5}}\left(\begin{array}{c}0\hfill \\ 0\hfill \\ 0\hfill \\ -{\beta }_h^{*}\hfill \\ {\lambda }_5\hfill \end{array}\right)\end{array}$$

(20)

with $Z_5=\sqrt{{({\beta }_h^{*})}^2+{\lambda }_5^2}$. The eigenvectors of ${\lambda }_{2,3}$ can be obtained using the ansatz ${\mathbf{v}}_{2,3}=(x,y,z,e,f)$. From Equation (10), it follows that the first three components $x,y,z$ can be determined independently from the remaining two. They address the animal SEIR subsystem. The solution reads

$$\begin{array}{c}\hfill \left(\begin{array}{c}x\hfill \\ y\hfill \\ z\hfill \end{array}\right)={\displaystyle \frac{1}{Z_j}}\left(\begin{array}{c}-{\beta }_a^{*}({\lambda }_j+{\alpha }_a)\hfill \\ {\beta }_a^{*}{\lambda }_j\hfill \\ {\lambda }_j({\lambda }_j+{\alpha }_a)\hfill \end{array}\right)\end{array}$$

(21)

for $j=2,3$, where $Z_j$ has still to be determined. Using the ansatz ${\mathbf{v}}_{2,3}=(x,y,z,e,f)$ in combination with the matrix (10) the remaining components $e,f$ can be determined. The result reads

$$\begin{array}{c}\hfill {\mathbf{v}}_j={\displaystyle \frac{1}{Z_j}}\left(\begin{array}{c}-{\beta }_a^{*}({\lambda }_j+{\alpha }_a)\hfill \\ {\beta }_a^{*}{\lambda }_j\hfill \\ {\lambda }_j({\lambda }_j+{\alpha }_a)\hfill \\ -{\beta }_{12}^{*}({\lambda }_j+{\alpha }_a)\left(1+{\displaystyle \frac{{\beta }_h^{*}}{{\lambda }_j-{\lambda }_5}}\right)\hfill \\ {\beta }_{12}^{*}({\lambda }_j+{\alpha }_a){\displaystyle \frac{{\lambda }_j}{{\lambda }_j-{\lambda }_5}}\hfill \end{array}\right)\end{array}$$

(22)

for $j=2,3$ and $Z_j$ defined such that $|{\mathbf{v}}_j|=1$. The eigenvectors ${\mathbf{v}}_2$ and ${\mathbf{v}}_5$ correspond to the potential order parameters of the SEIR-SIR model (4). However, ${\mathbf{v}}_5$ exhibits non-vanishing components only in the human subsystem. In contrast, ${\mathbf{v}}_2$ exhibits non-zero components in both subsystems. Consequently, for ${\lambda }_2>0$ it follows that ${\mathbf{v}}_2$ is the human–animal order parameter of the SEIR-SIR system. Furthermore, note that the subvector defined by Equation (21) is known as the order parameter of SEIR models [4]. This implies that the human–animal order parameter ${\mathbf{v}}_2$ exhibits in the animal subsystem the ordinary order parameter of epidemiological SEIR models. Finally, for ${\lambda }_5>0$ the eigenvector ${\mathbf{v}}_5$ is an order parameter of the SEIR-SIR model and may be regarded as the SIR order parameter of the human subsystem. ${\lambda }_5$ may be regarded as the corresponding eigenvalue of the SIR human subsystem (as anticipated above). In summary, the SEIR-SIR model (model B) maximally exhibits one human–animal order parameter.

The SIR-SEIR model (model C) exhibits the eigenvectors ${\mathbf{v}}_1=(1,0,0,0,0)$ and ${\mathbf{v}}_3=(0,0,1,0,0)$ associated to ${\lambda }_1=0$ and ${\lambda }_3=0$. From matrix (12), it follows that the eigenvectors ${\mathbf{v}}_{4,5}$ related to the eigenvalues ${\lambda }_{4,5}$ assume the form ${\mathbf{v}}_{4,5}=(0,0,x^{\prime },y^{\prime },z^{\prime })$. Since they address exclusively the human SEIR subsystem, the components should constitute an SEIR order parameter as shown in Equation (21). A detailed calculation shows that this is indeed the case:

$$\begin{array}{c}\hfill {\mathbf{v}}_j={\displaystyle \frac{1}{Z_j}}\left(\begin{array}{c}0\hfill \\ 0\hfill \\ -{\beta }_h^{*}({\lambda }_j+{\alpha }_h)\hfill \\ {\beta }_h^{*}{\lambda }_j\hfill \\ {\lambda }_j({\lambda }_j+{\alpha }_h)\hfill \end{array}\right)\end{array}$$

(23)

for $j=4,5$ with $Z_j$ such that $|{\mathbf{v}}_j|=1$. The derivation of ${\mathbf{v}}_2$ associated to ${\lambda }_2={\beta }_a^{*}-{\gamma }_a$ follows in part the derivation of ${\mathbf{v}}_2$ of the SIR-SIR model. As such, ${\mathbf{v}}_2$ satisfies

$$\begin{array}{c}\hfill \left(\begin{array}{ccccc}-{\lambda }_2& -{\beta }_a^{*}& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& -{\beta }_{12}^{*}& -{\lambda }_2& 0& -{\beta }_h^{*}\\ 0& {\beta }_{12}^{*}& 0& -({\alpha }_h+{\lambda }_2)& {\beta }_h^{*}\\ 0& 0& 0& {\alpha }_h& -({\gamma }_h+{\lambda }_2)\end{array}\right){\mathbf{v}}_2=0.\end{array}$$

(24)

Accordingly, the first two components of ${\mathbf{v}}_2$ constitute an SIR order parameter for the animal subsystem such that ${\mathbf{v}}_2=(-{\beta }_a^{*},{\lambda }_2,x^{″},y^{″},z^{″})/Z_2$ holds. A detailed calculation yields the remaining components $x^{″},y^{″},z^{″}$. The result reads

$$\begin{array}{c}\hfill {\mathbf{v}}_2={\displaystyle \frac{1}{Z_2}}\left(\begin{array}{c}-{\beta }_a^{*}\hfill \\ {\lambda }_2\hfill \\ -{\beta }_{12}^{*}{\displaystyle \frac{({\lambda }_2+{\alpha }_h)({\lambda }_2+{\gamma }_h)}{P_h({\lambda }_2)}}\hfill \\ {\beta }_{12}^{*}{\displaystyle \frac{{\lambda }_2({\lambda }_2+{\gamma }_h)}{P_h({\lambda }_2)}}\hfill \\ {\beta }_{12}^{*}{\displaystyle \frac{{\lambda }_2{\alpha }_h}{P_h({\lambda }_2)}}\hfill \end{array}\right)\end{array}$$

(25)

with the polynomial $P_h(\varphi )$ involving only human subsystem model parameters defined by

$$P_h(\varphi )=(\varphi +{\alpha }_h)(\varphi +{\gamma }_h)-{\alpha }_h{\beta }_h^{*}$$

(26)

and $Z_2$ chosen such that $|{\mathbf{v}}_2|=1$. Note that $P_h(\lambda )=0$ is the characteristic equation of ${\lambda }_{4,5}$ and yields the eigenvalues (13) of the SIR-SEIR model. However, in the context of ${\mathbf{v}}_2$, $P_h$ is applied to ${\lambda }_2$ rather than ${\lambda }_{4,5}$; see Equation (25). The eigenvectors ${\mathbf{v}}_2$ and ${\mathbf{v}}_4$ correspond to the potential order parameters of the SIR-SEIR model (5). ${\mathbf{v}}_4$ addresses only the human subsystem. In contrast, ${\mathbf{v}}_2$ addresses both subsystems. Consequently, for ${\lambda }_2>0$ we see that ${\mathbf{v}}_2$ is the human–animal order parameter of the SIR-SEIR system. For ${\lambda }_4>0$ the eigenvector ${\mathbf{v}}_4$ is an order parameter of the system and may be regarded as the SEIR order parameter of the human subsystem. However, it does not qualify as a human–animal order parameter. In summary, for the SIR-SEIR model (model C) there is maximally one human–animal order parameter.

The eigenvectors ${\mathbf{v}}_j$ of the SEIR-SEIR model (model D) defined by Equation (6) satisfy

$$\begin{array}{c}\hfill \left(\begin{array}{cccccc}-{\lambda }_j& 0& -{\beta }_a^{*}& 0& 0& 0\\ 0& -({\alpha }_a+{\lambda }_j)& {\beta }_a^{*}& 0& 0& 0\\ 0& {\alpha }_a& -({\gamma }_a+{\lambda }_j)& 0& 0& 0\\ 0& 0& -{\beta }_{12}^{*}& -{\lambda }_j& 0& -{\beta }_h^{*}\\ 0& 0& {\beta }_{12}^{*}& 0& -({\alpha }_h+{\lambda }_j)& {\beta }_h^{*}\\ 0& 0& 0& 0& {\alpha }_h& -({\gamma }_h+{\lambda }_j)\end{array}\right){\mathbf{v}}_j=0;\end{array}$$

(27)

see also Equation (14). The eigenvalues ${\lambda }_1=0$ and ${\lambda }_4=0$ are associated with the eigenvectors ${\mathbf{v}}_1=(1,0,0,0,0,0)$ and ${\mathbf{v}}_4=(0,0,0,1,0,0)$. The $3\times 3$ matrix in the right-bottom corner defines eigenvectors of the form ${\mathbf{v}}_j=(0,0,0,x^{‴},y^{‴},z^{‴})$. They have components like SEIR order parameter eigenvectors as in Equation (21) but describe the human subsystem. That is, for $j=5,6$ the eigenvectors read

$$\begin{array}{c}\hfill {\mathbf{v}}_j={\displaystyle \frac{1}{Z_j}}\left(\begin{array}{c}0\hfill \\ 0\hfill \\ 0\hfill \\ -{\beta }_h^{*}({\lambda }_j+{\alpha }_h)\hfill \\ {\beta }_h^{*}{\lambda }_j\hfill \\ {\lambda }_j({\lambda }_j+{\alpha }_h)\hfill \end{array}\right)\end{array}$$

(28)

with $Z_j$ such that $|{\mathbf{v}}_j|=1$. Not surprisingly, they resemble the eigenvectors of the SEIR system of the SIR-SEIR model (compare Equations (23) and (28)). In view of the left-upper $3\times 3$ matrix in Equation (27), the remaining eigenvectors ${\mathbf{v}}_{2,3}$ assume the form

$$\begin{array}{c}\hfill {\mathbf{v}}_j={\displaystyle \frac{1}{Z_j}}\left(\begin{array}{c}-{\beta }_a^{*}({\lambda }_j+{\alpha }_a)\hfill \\ {\beta }_a^{*}{\lambda }_j\hfill \\ {\lambda }_j({\lambda }_j+{\alpha }_a)\hfill \\ a^{″}\hfill \\ b^{″}\hfill \\ c^{″}\hfill \end{array}\right),\end{array}$$

(29)

where $a^{″},b^{″},c^{″}$ still need to be determined. Substituting Equation (29) into Equation (27) allows one to determine those components, which leads to

$$\begin{array}{c}\hfill {\mathbf{v}}_j={\displaystyle \frac{1}{Z_j}}\left(\begin{array}{c}-{\beta }_a^{*}({\lambda }_j+{\alpha }_a)\hfill \\ {\beta }_a^{*}{\lambda }_j\hfill \\ {\lambda }_j({\lambda }_j+{\alpha }_a)\hfill \\ -{\beta }_{12}^{*}({\lambda }_j+{\alpha }_a){\displaystyle \frac{({\lambda }_j+{\alpha }_h)({\lambda }_j+{\gamma }_h)}{P_h({\lambda }_j)}}\hfill \\ {\beta }_{12}^{*}({\lambda }_j+{\alpha }_a){\displaystyle \frac{{\lambda }_j({\lambda }_j+{\gamma }_h)}{P_h({\lambda }_j)}}\hfill \\ {\beta }_{12}^{*}({\lambda }_j+{\alpha }_a){\displaystyle \frac{{\lambda }_j{\alpha }_h}{P_h({\lambda }_j)}}\hfill \end{array}\right),\end{array}$$

(30)

with $P_h$ defined by Equation (26) and $Z_j$ such that $|{\mathbf{v}}_j|=1$. The last three components in Equation (30) related to the human subsystem resemble the human subsystem components of ${\mathbf{v}}_2$ of the SIR-SEIR model. As indicated in Equation (30) they are identical except for the pre-factor ${\lambda }_j+{\alpha }_a$ (compare Equations (25) and (30)). The eigenvectors ${\mathbf{v}}_2$ and ${\mathbf{v}}_5$ correspond to the potential order parameters of the SEIR-SEIR model (6). ${\mathbf{v}}_5$ addresses only the human subsystem, while ${\mathbf{v}}_2$ addresses both subsystems. Consequently, ${\mathbf{v}}_2$ for ${\lambda }_2>0$ is the human–animal order parameter of the SEIR-SEIR model (model D). For ${\lambda }_5>0$ the eigenvector ${\mathbf{v}}_5$ is also an order parameter of the SEIR-SEIR model. It may be regarded as the SEIR order parameter of the human subsystem (but not as a human–animal order parameter). In summary, for the SEIR-SEIR model (model D) there is maximally one human–animal order parameter.

In conclusion, all four models exhibit maximally one human–animal order parameter. □

3.2. Implications: Amplitude Dynamics and the Role of Human–Animal Order Parameters

Let us define the amplitudes $A_1,\dots ,A_m$ of the model eigenvectors ${\mathbf{v}}_j$ implicitly by the expansion [4]

$$\begin{array}{c}\hfill \mathbf{u}=\sum _{j=1}^m{\mathbf{v}}_j{A}_j\end{array}$$

(31)

with $m=4$ for model A, $m=5$ for models B and C, and $m=6$ for model D, where it is assumed that the eigenvectors are linearly independent from each other and form a complete vector basis (i.e., degenerated, special cases may be discussed separately). The expansion (31) holds for arbitrary $\mathbf{u}$ (i.e., it holds beyond the linear initial phase dynamics that will be defined below) [4]. In line with Equation (31), the state vector of the model under consideration can be expressed like

$$\begin{array}{c}\hfill \mathbf{x}={\mathbf{x}}^{*}+\sum _{j=1}^m{\mathbf{v}}_j{A}_j,\end{array}$$

(32)

where ${\mathbf{x}}^{*}$ denotes the fixed point of interest (see also Section 3.1). Explicitly, the amplitudes $A_j$ can be computed from either $\mathbf{u}$ or $\mathbf{x}$ like

$$\begin{array}{c}\hfill A_j={\mathbf{w}}_j·\mathbf{u}={\mathbf{w}}_j·(\mathbf{x}-{\mathbf{x}}^{*}),\end{array}$$

(33)

where the dot denotes the scalar product and ${\mathbf{w}}_j$ the bi-orthogonal vector associated to ${\mathbf{v}}_j$ [4]. In general, the bi-orthogonal vectors ${\mathbf{w}}_j$ can be determined numerically with the help of the analytical expressions for ${\mathbf{v}}_j$ [4]. As mentioned in the introduction, initial phases are crucial phases of infection waves. The expansion defined by Equation (32) can be discussed for such initial phases in the context of the four fundamental models A, B, C, and D. To this end, it is helpful to make the following definition.

Definition 3.

The linearized initial phase dynamics is the dynamics of the state vector $\mathbf{x}$ as defined by the linearized evolution Equation (8) with $\mathbf{x}={\mathbf{x}}^{*}+\mathbf{u}$ and an initial state $\mathbf{x}(t=0)$ in an ϵ-environment of the fixed point ${\mathbf{x}}^{*}$ such that $|\mathbf{u}(t=0)|<ϵ$.

The idea here is that infection waves typically start with a small number of infected individuals such that the human–animal system initially is close to a fixed point. Mathematically, this property of the initial state to be in a close vicinity of a fixed point can be expressed by requiring that the initial state is in an $ϵ$-environment of a fixed point and by choosing a small value for $ϵ$. If $ϵ$ is sufficiently small, the linearized model defined by (8) describes a good approximation of its original nonlinear model (either A, B, C, or D). It is beyond the current study to define precisely what is meant by a good approximation. It is sufficient to note that on the one hand the accuracy of the linear approximation solution as measured by reasonably defined quantities typically improves when $ϵ$ is made smaller and smaller. On the other hand, for Mpox infection waves, at a certain point in time, the nonlinear aspects of the models A, B, C, and D will become relevant. At that point in time, the linear approximate model (8) will fail to give an accurate description of the infection dynamics. In summary, Equation (8) is tailored to describe the initial phase dynamics, as it is also pointed out in Definition 3. Solutions of the initial phase dynamics as defined in Definition 3 and by Equation (8) are given in terms of the superposition

$$\begin{array}{c}\hfill \mathbf{x}={\mathbf{x}}^{*}+\sum _{j=1}^m{\mathbf{v}}_j{A}_j(0)exp\{{\lambda }_{j}t\},\end{array}$$

(34)

where $A_j(0)$ denotes the initial amplitudes at the initial time point $t=0$. The initial amplitudes can be computed from the initial state $\mathbf{x}(0)$ with the help of Equation (33) [4]. For example, for the SIR-SIR model (3) the superposition solution (34) reads

$$\begin{array}{c}\hfill \mathbf{x}={\mathbf{x}}^{*}+{\mathbf{v}}_2{A}_2(0)exp\{{\lambda }_{2}t\}+{\mathbf{v}}_4{A}_4(0)exp\{{\lambda }_{4}t\}+{\mathbf{h}}_0,{\mathbf{h}}_0={\mathbf{v}}_1{A}_1(0)+{\mathbf{v}}_3{A}_3(0).\end{array}$$

(35)

Theorem 3.

For the admissible model parameters listed in Section 2.2 the initial phase dynamics of any Mpox outbreak dynamics as described by the four fundamental models A, B, C, and D defined by Equations (3)–(6) and the linearized evolution Equation (8) corresponds to one of three qualitatively different scenarios.

In other words, despite the differences across the four models, there exist only three qualitatively different scenarios how Mpox disease outbreaks (as described by those models) initially evolve.

Proof. 

Let ${\lambda }^{(h,X)}$ and ${\mathbf{v}}^{(h,X)}$ denote the eigenvalue and eigenvector of the potential order parameter of the human subsystem of the model X with $X=A,B,C,D$. Likewise, let ${\lambda }^{(12,X)}$ and ${\mathbf{v}}^{(12,X)}$ denote the eigenvalue and eigenvector of the potential human–animal order parameter of the model X. Let us next turn to Equation (34). All four models exhibit two neutrally stable eigenvectors that are associated with zero eigenvalues and point into the directions of $s_a$ and $s_h$, respectively. Since we consider an Mpox disease outbreak, the disease-free fixed point under consideration must be unstable [4], which implies that at least one of the two eigenvalues ${\lambda }^{(h,X)}$ and ${\lambda }^{(12,X)}$ must be positive. Furthermore, there are maximally two positive eigenvalues. In total, these considerations lead to the following three Mpox disease outbreak scenarios:

$$\begin{array}{ccc}\hfill Scenario(i)& :& {\lambda }^{(h,X)}>0,{\lambda }^{(12,X)}<0,\hfill \\ \hfill Scenario(ii)& :& {\lambda }^{(h,X)}<0,{\lambda }^{(12,X)}>0,\hfill \\ \hfill Scenario(iii)& :& {\lambda }^{(h,X)}>0,{\lambda }^{(12,X)}>0.\hfill \end{array}$$

(36)

In scenario (i), the outbreak dynamics is characterized by two neutrally stable directions (see above) and one unstable direction given by the SIR ($X=A,C$) or SEIR ($X=B,D$) human subsystem order parameter ${\mathbf{v}}^{(h,X)}$. The remaining directions are given in terms of stable eigenvectors. In scenario (ii), the outbreaks dynamics is characterized again by two neutrally stable directions. There is one unstable direction given by the human–animal order parameter ${\mathbf{v}}^{(12,X)}$. The remaining directions are given in terms of stable eigenvectors. Although scenario (i) and (ii) have in common that they both feature only a single unstable direction, scenarios (i) and (ii) differ from each other qualitatively because they involve different types of order parameters. Further details about this difference will be discussed below. Finally, scenario (iii) describes an outbreak dynamics that is characterized by two neutrally stable directions and two unstable directions given in terms of the two maximally possible order parameters discussed in Section 3.1. The remaining directions for the models B, C, and D are given by stable eigenvectors. Scenario (iii) differs from scenarios (i) and (ii) qualitatively by the number of unstable directions. □

In what follows, the scenarios will be discussed in more detail. In scenario (i) we have ${\lambda }^{(h,X)}>0$ and ${\lambda }^{(12,X)}<0$. The state $\mathbf{x}$ of the system under consideration evolves away from the fixed point along the direction ${\mathbf{v}}^{(h,X)}$, which has only non-vanishing components in the space of the human system (e.g., for model A: ${\mathbf{v}}^{(h,X)}={\mathbf{v}}_4$ with $v_{4,1}=v_{4,2}=0$, $v_{4,3}=-{\beta }_h^{*}/Z_4$ and $v_{4,4}={\lambda }_4/Z_4$). That is, $A^{(h,X)}(t)=A^{(h,X)}(0)exp\{{\lambda }^{(h,X)}t\}$ holds and describes an exponential increase in the amplitude related to ${\mathbf{v}}^{(h,X)}$. In contrast, let $\xi (X)$ denote the index of the human subsystem order parameter ${\mathbf{v}}^{(h,X)}$ of scenario (i) (e.g., $\xi (A)=4$ for model A). Then, components of the initial state $\mathbf{x}(0)$ in all other directions as measured by $A_j(0)={\mathbf{w}}_j(\mathbf{x}(0)-{\mathbf{x}}^{*})$ for $j\ne \xi$ and $j=1,\dots ,m$ either decay in magnitude or remain constant. More precisely, let $s(i)$ denote the indices of negative eigenvalues associated with stable eigenvectors with $i=1$ for model A, $i=1,2$ for models B and C, and $i=1,2,3$ for model D. Then $A_{s(i)}(t)=A_{s(i)}(0)exp\{{\lambda }_{s(i))}t\}$ describes a dynamics towards the fixed point. Let us split the overall dynamics into three components: an outwards dynamics (${\mathbf{x}}_{out}$) describing the dynamics away from the fixed point along unstable eigenvectors (i.e., order parameters), an inwards dynamics (${\mathbf{x}}_{in}$) describing the dynamics towards the fixed point along stable directions, and ${\mathbf{h}}_0$ describing the constant part of the dynamics related to the two neutrally stable directions. Accordingly, Equation (34) for all three scenarios (i), (ii), and (iii) becomes

$$\begin{array}{c}\hfill \mathbf{x}={\mathbf{x}}^{*}+{\mathbf{x}}_{out}+{\mathbf{x}}_{in}+{\mathbf{h}}_0.\end{array}$$

(37)

Specifically, for scenario (i), we obtain

$$\begin{array}{ccc}\hfill {\mathbf{x}}_{out}& =& {\mathbf{v}}^{(h,X)}{A}^{(h,X)}(0)exp\{{\lambda }^{(h,X)}t\},{\lambda }^{(h,X)}>0,\hfill \\ \hfill {\mathbf{x}}_{in}& =& \sum _{i=1}^{m-3}{\mathbf{v}}_{s(i)}{A}_{s(i)}(0)exp\{{\lambda }_{s(i)}t\},{\lambda }_{s(i)}<0.\hfill \end{array}$$

(38)

The component ${\mathbf{x}}_{out}$ is most relevant for the disease outbreak. As mentioned above, ${\mathbf{v}}^{(h,X)}$, addresses only the human subsystem. Consequently, scenario (i) describes an outbreak due to human to human transmissions (while the infection dynamics in the animal subsystem subsides). In line with earlier studies on COVID-19 outbreaks [4], ${\mathbf{v}}^{(h,X)}$ describes the organization of this type of Mpox outbreak. For example, for the SIR-SIR model (3) (model A), from the order parameter ${\mathbf{v}}_4$ defined by Equation (17) it follows that during the initial phase of an Mpox outbreak as a result of the outwards dynamics ${\mathbf{x}}_{out}$ changes $\Delta s_h$ and $\Delta i_h$ in the relative sizes of the susceptible and infected populations satisfy

$$\begin{array}{c}\hfill {\mathbf{x}}_{out}\Rightarrow {\displaystyle \frac{\Delta i_h}{\Delta s_h}}={\displaystyle \frac{v_{4,4}}{v_{4,3}}}=-{\displaystyle \frac{{\lambda }_4}{{\beta }_h^{*}}}.\end{array}$$

(39)

Accordingly, a decay of susceptibles by $\Delta s_h<0$ comes with an increase in infectious individuals by $\Delta i_h=-{\lambda }_4\Delta s_h/{\beta }_h^{*}>0$ and vice versa an increase in infectious individuals by $\Delta i_h>0$ is associated with a decrease in susceptibles of $\Delta s_h=-{\beta }_h^{*}\Delta i_h/{\lambda }_4<0$. Similar considerations can be made for models B, C, and D based on the human subsystem order parameters defined by Equations (20), (23), and (28), respectively. For example, for the SIR-SEIR model (5) (model C), due to the outwards dynamics ${\mathbf{x}}_{out}$ the populations $e_h$ and $i_h$ change relative to each other like

$$\begin{array}{c}\hfill {\mathbf{x}}_{out}\Rightarrow {\displaystyle \frac{\Delta i_h}{\Delta e_h}}={\displaystyle \frac{v_{4,5}}{v_{4,4}}}={\displaystyle \frac{{\lambda }_4+{\alpha }_h}{{\beta }_h^{*}}}.\end{array}$$

(40)

Any increase in exposed humans by $\Delta E_h=100$ individuals implies an increase in infectious individuals by $\Delta I_h=({\lambda }_4+{\alpha }_h)/{\beta }_h^{*}·100$ individuals (where we have used that $\Delta I_h/\Delta E_h=\Delta i_h/\Delta e_h$ holds). In summary, during the initial phase of a scenario (i) outbreak the component ${\mathbf{x}}_{out}$ that drives the outbreak establishes rigid relationships between the dynamics of the subpopulations $s_h$ and $i_h$ (all models) and $e_h$ (models C,D). These relationships, in turn, are determined by the model-specific human subsystem order parameter ${\mathbf{v}}^{(h,X)}$.

The Mpox outbreak scenario (ii) is characterized by ${\lambda }^{(h,X)}<0$ and ${\lambda }^{(12,X)}>0$. The system state $\mathbf{x}$ evolves away from the fixed point along the direction of the human–animal order parameter ${\mathbf{v}}^{(12,X)}$. Equation (37) holds with

$$\begin{array}{ccc}\hfill {\mathbf{x}}_{out}& =& {\mathbf{v}}^{(12,X)}{A}^{(12,X)}(0)exp\{{\lambda }^{(12,X)}t\},{\lambda }^{(12,X)}>0,\hfill \\ \hfill {\mathbf{x}}_{in}& =& \sum _{i=1}^{m-3}{\mathbf{v}}_{s(i)}{A}_{s(i)}(0)exp\{{\lambda }_{s(i)}t\},{\lambda }_{s(i)}<0.\hfill \end{array}$$

(41)

By Definition 2, the human–animal order parameter ${\mathbf{v}}^{(12,X)}$ exhibits components both in the animal and human subsystems (e.g., see Equation (19) for model A). Consequently, scenario (ii) describes Mpox outbreaks that involve infection outbreaks in the animal subsystems that drive Mpox outbreaks in the corresponding human subsystems. While in the previously discussed scenario, scenario (i), the infection dynamics in an animal subsystem immediately subsides and the initial phase of a wave is caused by human to human virus transmissions, in the scenario (ii) the infection dynamics in a human subsystem would subside immediately if the system would be decoupled from its animal reservoir. More precisely, if we would put ${\beta }_{12}=0$, then due to the fact that in scenario (ii) we have ${\lambda }^{(h,X)}<0$ the disease-free fixed points of the human subsystem under consideration are neutrally stable. In other words, due to the coupling with ${\beta }_{12}>0$ the infection wave of the animal subsystem under consideration drives an infection wave in the corresponding human subsystem. Moreover, ${\mathbf{v}}^{(12,X)}$ describes the organization of scenario (ii) outbreaks caused by the outwards dynamics ${\mathbf{x}}_{out}$. For example, for model A from Equation (19) it follows that during the initial phase of such outbreaks changes $\Delta s_a$ and $\Delta i_a$ in the relative sizes of the susceptible and infectious animal populations satisfy

$$\begin{array}{c}\hfill {\mathbf{x}}_{out}\Rightarrow {\displaystyle \frac{\Delta i_a}{\Delta s_a}}={\displaystyle \frac{v_{2,2}}{v_{2,1}}}=-{\displaystyle \frac{{\lambda }_2}{{\beta }_a^{*}}}.\end{array}$$

(42)

Accordingly, a decrease in susceptibles animals as measured by $\Delta s_a<0$ comes with an increase in infectious animals as measured by $\Delta i_a$ with $\Delta i_a=-{\lambda }_2\Delta s_a/{\beta }_a^{*}>0$ and vice versa an increase in infectious animals $\Delta i_a>0$ is associated with an decrease in susceptible animals like $\Delta s_a=-{\beta }_a^{*}\Delta i_a/{\lambda }_2<0$. Importantly, the human–animal order parameter ${\mathbf{v}}_2$ of model A also describes the coupling between the animal and human subsystems. For example, changes $\Delta i_h$ and $\Delta i_s$ due to the outwards dynamics component ${\mathbf{x}}_{out}$ are given by

$$\begin{array}{c}\hfill {\mathbf{x}}_{out}\Rightarrow {\displaystyle \frac{\Delta i_h}{\Delta i_a}}={\displaystyle \frac{v_{2,4}}{v_{2,2}}}={\displaystyle \frac{{\beta }_{12}^{*}}{{\lambda }_2-{\lambda }_4}}>0.\end{array}$$

(43)

That is, the outwards dynamics exhibits the property that an increase in infectious animals, say, by 1% is associated with an increase in the population of infectious human individuals by $1\%·{\beta }_{12}^{*}/({\lambda }_2-{\lambda }_4)$. This example and Equation (43) illustrate that there is a rigid coupling between the animal and the human subsystems, which in the case of SIR-SIR systems (model A) is described in detail by the human–animal order parameter ${\mathbf{v}}_2$. With the help of the previously derived human–animal order parameters defined by Equations (22), (25), and (30) for models B, C, and D, respectively, similar explicit conclusions can be drawn about the emerging order involved in Mpox waves as described by those models.

The third scenario, scenario (iii), is characterized by two positive eigenvalues, ${\lambda }^{(h,X)}>0$ and ${\lambda }^{(12,X)}>0$, and describes Mpox infection waves established by an interplay (or co-existence) of two order parameters: ${\mathbf{v}}^{(h,X)}$ and ${\mathbf{v}}^{(h,12)}$. Accordingly, the outwards dynamics away from the fixed point does not take place along a single direction. Rather, it takes place in a plane spanned by the vectors ${\mathbf{v}}^{(h,X)}$ and ${\mathbf{v}}^{(h,12)}$. More precisely, Equation (37) holds with

$$\begin{array}{ccc}\hfill {\mathbf{x}}_{out}& =& {\mathbf{v}}^{(h,X)}{A}^{(h,X)}(0)exp\{{\lambda }^{(h,X)}t\}\hfill \\ & & +{\mathbf{v}}^{(12,X)}{A}^{(12,X)}(0)exp\{{\lambda }^{(12,X)}t\},{\lambda }^{(h,X)}>0,{\lambda }^{(12,X)}>0,\hfill \\ \hfill {\mathbf{x}}_{in}& =& \sum _{i=1}^{m-4}{\mathbf{v}}_{s(i)}{A}_{s(i)}(0)exp\{{\lambda }_{s(i)}t\},{\lambda }_{s(i)}<0,\hfill \end{array}$$

(44)

where $A^{(h,X)}$ and $A^{(12,X)}$ denote the amplitudes of the eigenvectors ${\mathbf{v}}^{(h,X)}$ and ${\mathbf{v}}^{(h,12)}$. According to the initial phase dynamics described by the linearized equation (8), the amplitudes $A^{(h,X)}$ and $A^{(12,X)}$ measuring distances along the order parameter directions increase exponentially in magnitude over time. In doing so, the state of the system evolves further and further away from the disease-free fixed point in the 2D plane spanned by the two order parameters. The precise trajectory depends on the model parameters ${\lambda }^{(h,X)}>0$, ${\lambda }^{(12,X)}>0$ and initial conditions $A^{(h,X)}(0)$, $A^{(12,X)}(0)$. As indicated in Equation (44), an inward dynamics does not exist for model A, while for models B and C we have $i=1$ and for model D we have $i=1,2$.

For all three scenarios and all four models, during an initial interval $[0,T]$ Equation (37) may be used to compute approximative solutions to the exact solutions of the nonlinear models. That is, let $\mathbf{x}$ denote the solution of one of the models A, B, C, or D. Then, Equation (37) may be used to describe an approximative relationship like

$$\begin{array}{c}\hfill \mathbf{x}\approx {\mathbf{x}}^{*}+{\mathbf{x}}_{out}+{\mathbf{x}}_{in}+{\mathbf{h}}_0.\end{array}$$

(45)

The power of this approximation comes for situations in which the inwards dynamics is negligible. Such situations may arise when the eigenvalues of the inwards dynamics are relatively large in the amount such that ${\mathbf{x}}_{in}$ decays rapidly towards zero or when the initial amplitudes $A_s(i)(0)$ are relatively small compared to the order parameter amplitudes $A^{(h,X)}(0)$ and $A^{(12,X)}(0)$. If ${\mathbf{x}}_{in}$ can be neglected, then Equation (45) simplifies to yield

$$\begin{array}{c}\hfill \mathbf{x}\approx {\mathbf{x}}^{*}+{\mathbf{x}}_{out}+{\mathbf{h}}_0.\end{array}$$

(46)

For the three scenarios, Equation (46) reads explicitly

$$\begin{array}{ccc}\hfill (i)& :& \mathbf{x}\approx {\mathbf{x}}^{*}+{\mathbf{v}}^{(h,X)}{A}^{(h,X)}(0)exp\{{\lambda }^{(h,X)}t\}+{\mathbf{h}}_0,\hfill \\ \hfill (ii)& :& \mathbf{x}\approx {\mathbf{x}}^{*}+{\mathbf{v}}^{(12,X)}{A}^{(12,X)}(0)exp\{{\lambda }^{(12,X)}t\}+{\mathbf{h}}_0,\hfill \\ \hfill (iii)& :& \mathbf{x}\approx {\mathbf{x}}^{*}+{\mathbf{v}}^{(h,X)}{A}^{(h,X)}(0)exp\{{\lambda }^{(h,X)}t\}+{\mathbf{v}}^{(12,X)}{A}^{(12,X)}(0)exp\{{\lambda }^{(12,X)}t\}\hfill \\ & & +{\mathbf{h}}_0.\hfill \end{array}$$

(47)

The discussion so far focused on the initial phase dynamics. This discussion was centered around the order parameters and their amplitudes. As mentioned above in the context of the Definition 3, when looking at the later stages of an infection wave, then, in general, nonlinear aspects of epidemiological models become relevant. At those later stages, all amplitudes may make considerable contributions to the infection dynamics. This issue will be illustrated in Section 3.3.

In closing this section, let us point out that for the SIR-SIR model there exists a peculiarity that does not exist for the higher-dimensional models B, C, and D. The two potential order parameters ${\mathbf{v}}_2$ and ${\mathbf{v}}_4$ of the SIR-SIR model and their amplitudes completely describe the populations of infectious individuals $i_a$ and $i_h$. That is, let $P_i$ the projection of the state $\mathbf{x}$ in the subspace of $i_a$ and $i_h$, then

$$\begin{array}{c}\hfill \left(\begin{array}{c}{i}_a(t)\hfill \\ i_h(t)\hfill \end{array}\right)=P_i\left({\mathbf{x}}^{*}+{\mathbf{v}}_2{A}_2(t)+{\mathbf{v}}_4{A}_4(t)\right)\end{array}$$

(48)

holds. Equation (48) holds for any time point t and is not an approximation. Equation (48) follows from the fact that ${\mathbf{v}}_1$ and ${\mathbf{v}}_3$ do not have any components in the subspace of $i_a$ and $i_h$.

3.3. Role of the Neutrally Stable Eigenvectors and Their Amplitudes

Let ${\mathbf{v}}^{(a,0,X)}$ and ${\mathbf{v}}^{(h,0,X)}$ denote the eigenvectors of model X that exhibit only an $s_a$ or $s_h$ component, respectively. For example, for model A we have ${\mathbf{v}}^{(a,0,X)}={\mathbf{v}}_1=(1,0,0,0)$ and ${\mathbf{v}}^{(h,0,X)}={\mathbf{v}}_3=(0,0,1,0)$ (see Section 3.1). Let ${\lambda }^{(a,0,X)}$ and ${\lambda }^{(h,0,X)}$ denote the zero eigenvalues associated with those eigenvectors. Likewise, let $A^{(a,0,X)}$ and $A^{(h,0,X)}$ denote the amplitudes of ${\mathbf{v}}^{(a,0,X)}$ and ${\mathbf{v}}^{(h,0,X)}$.

Theorem 4.

For the admissible model parameters listed in Section 2.2, any Mpox infection wave as described by the four fundamental models A, B, C, or D defined by Equations (3)–(6) exhibits the following final stage properties. All amplitudes converge to zero expect for the amplitudes $A^{(a,0,X)}$ and $A^{(h,0,X)}$ associated with the zero eigenvalues ${\lambda }^{(a,0,X)}={\lambda }^{(h,0,X)}=0$. The final fixed point values of $A^{(a,0,X)}$ and $A^{(h,0,X)}$ for $t\to \infty$ correspond to the decay of the respective susceptible subpopulations $s_a$ and $s_h$ over the entire course of the infection wave like $A^{(a,0,X)}(t\to \infty )=s_a^{*}(t\to \infty )-s_a^{*}(0)<0$ and $A^{(h,0,X)}(t\to \infty )=s_h^{*}(t\to \infty )-s_h^{*}(0)<0$.

In Theorem 4, $s_a^{*}(0)$ and $s_h^{*}(0)$ denote the fixed point values $s_a^{*}$ and $s_h^{*}$, respectively, considered in Section 3.1. The specifier $(0)$ has been added to distinguish more clearly between the final and initial fixed point values of a wave. Theorem 4 illustrates that while the order parameters and their amplitudes are the key building-blocks that describe the initial phase of an Mpox infection wave, the remaining amplitudes can play crucial roles at later stages during the time course of an infection wave. Theorem 4 highlights the role of the neutrally stable eigenvectors and their amplitudes.

Proof. 

As worked out in Section 3.1 in the context of Theorem 2, the neutrally stable eigenvectors are the only eigenvectors that do not exhibit components of infected individuals whether exposed or infectious, human or animal (i.e., they do not exhibit components such as $e_a$, $i_a$, $e_h$, $i_h$). They only exhibit components in the direction of susceptible populations. This implies that all other eigenvectors feature at least one component that describes an infected subpopulation. Furthermore, as discussed in Section 3.1 in the context of Theorem 1, the models A, B, C, and D only exhibit disease-free fixed points. Consequently, any wave dynamics eventually approaches a disease-free fixed point that features zero infected individuals [4,41]. This implies that in the long term

$$\underset{t\to \infty }{lim}{A}_j(t)=0\forall j:j=1,\dots ,m,j\ne n(a,0,X),j\ne n(h,0,X),$$

(49)

where $n(a,0,X)$ and $n(h,0,X)$ correspond to the indices j of the eigenvalues ${\lambda }^{(a,0,X)}$ and ${\lambda }^{(h,0,X)}$, respectively. That is, all amplitudes vanish expect for those related to the neutrally stable eigenvectors. From Equations (32) and (49), it then follows that the state vector in the long term satisfies

$${\mathbf{x}}^{*}(\infty )={\mathbf{x}}^{*}(0)+{\mathbf{v}}^{(a,0,X)}{A}^{(a,0,X)}(\infty )+{\mathbf{v}}^{(h,0,X)}{A}^{(h,0,X)}(\infty ),$$

(50)

where ${\mathbf{x}}^{*}(0)$ denotes the initial fixed point vector under consideration and ${\mathbf{x}}^{*}(\infty )$ denotes the new fixed point vector that the wave assumes when it has completely subsided (i.e., for $t\to \infty$). From Equation (50) and the definition of ${\mathbf{v}}^{(a,0,X)}$ and ${\mathbf{v}}^{(h,0,X)}$ it follows that

$$\begin{array}{c}\hfill \Delta s_{a,\infty }=s_a^{*}(t\to \infty )-s_a^{*}(0)=A^{(a,0,X)}(t\to \infty ),\\ \hfill \Delta s_{h,\infty }=s_h^{*}(t\to \infty )-s_h^{*}(0)=A^{(h,0,X)}(t\to \infty ),\end{array}$$

(51)

which is the statement made in Theorem 4. □

Equation (51) in combination with the results presented in Section 3.2 points out that changes in susceptibles $\Delta s$ can be determined by different types of eigenvectors and their amplitudes. During initial stages, changes $\Delta s$ in relation to changes of other populations are determined by order parameters (see, e.g., Equations (39) and 42)). During these initial stages, the exponential decay of susceptible populations again is determined by order parameters (see Equation (47)). In contrast, the final stage values of susceptibles are determined by the amplitudes of the neutrally stable eigenvectors. This also implies that during the course of an infection wave at some point in time the amplitudes $A^{(a,0,X)}$ and $A^{(h,0,X)}$ of the neutrally stable eigenvectors make essential contributions to the infection wave dynamics. This issue will be illustrated in Section 3.4 below.

3.4. Simulation

In this section some aspects of the aforementioned results will be illustrated by means of a simulation. For the sake of brevity, only a simulation for the simplest model, the SIR-SIR model defined by Equation (3), will be presented. The following model parameters were used: ${\beta }_a=40/$y, ${\gamma }_a=12/$y, ${\beta }_h=32.85/$y, and ${\gamma }_h=28.08/$y [23], where “y” stands for one year. The goal was to simulate a wave with a peak at about 2 months after wave onset. Such a 2-months-peak has been observed during the 2017 Monkeypox outbreak in Nigeria [8,9]. To this end, ${\beta }_{12}$ was assumed to be ${\beta }_{12}=10/$y, which produced the intended 2-months peak (see below). For the selected model parameters the two non-zero eigenvalues of the SIR-SIR model were found to be ${\lambda }_2=24.0/$y and ${\lambda }_4=4.44/$y. That is, the model described a scenario (iii) outbreak involving two order parameters: the human–animal order parameter ${\mathbf{v}}^{(12,A)}={\mathbf{v}}_2$ associated with ${\lambda }_2^{(12,A)}=24.0/$y and the human subsystem order parameter ${\mathbf{v}}^{(h,A)}={\mathbf{v}}_4$ associated with ${\lambda }_4^{(h,A)}=4.4/$y.

Equation (3) was solved numerically (using a Euler forward method with a time step $\tau$ of $\tau =0.01$days $=2.74·{10}^{-5}$years) to obtain trajectories for the state variables $s_a$, $i_a$, $s_h$, $s_i$. From the trajectories of the state variables thus obtained the trajectories of the amplitude variables $A_1,\dots A_4$ were computed. To this end, the explicit expression for ${\mathbf{v}}_2$ and ${\mathbf{v}}_4$ (see Equations (19) and (17), respectively) were used and ${\mathbf{w}}_j$ for $j=1,\dots ,4$ were computed numerically [4]. The amplitudes $A_j$ were then obtained from Equation (33) for $j=1,\dots ,4$.

A noted in Section 3.2 and Section 3.3, while the order parameter amplitudes $A_2^{(12,A)}$ and $A_4^{(h,A)}$ initially increase exponentially, the increase eventually is stopped and they decay to zero when the wave eventually subsides. This is consistent with the fact that the evolution equations of amplitudes of epidemiological models in general are nonlinear [4]. In this context, it is important to point out that even when order parameter amplitudes (such as $A_2^{(12,A)}$ and $A_4^{(h,A)}$) stop to increase exponentially and nonlinear effects become relevant, then still for some period order parameter amplitudes continue to make the main contributions to the infection dynamics. The reason for this is that initially the remaining amplitudes either decayed in magnitude or remained constant. In order to illustrate the role of the two order parameters ${\mathbf{v}}_2^{(12,A)}$ and ${\mathbf{v}}_4^{(h,A)}$ and their amplitudes of the SIR-SIR model for the entire duration of the simulated scenario (iii) infection wave, an approximation ${\mathbf{x}}_{app}$ of $\mathbf{x}$ was used that was based on the two order parameters like

$$\begin{array}{c}\hfill {\mathbf{x}}_{app}={\mathbf{x}}^{*}+{\mathbf{v}}_2^{(12,A)}{A}_2^{(12,A)}(t)+{\mathbf{v}}_4^{(h,A)}{A}_4^{(h,A)}(t)+{\mathbf{h}}_0\end{array}$$

(52)

with ${\mathbf{h}}_0$ as defined in Equation (35). Note that this approximation goes beyond the initial phase approximations given by Equation (47) for scenario (iii) outbreaks of models A, B, C, and D. In Equation (52), the amplitudes do not necessarily increase in an exponential manner.

Finally, in order to quantify the contributions that the order parameter amplitudes $A_2^{(12,A)}$ and $A_4^{(h,A)}$ as well as the neutrally stable amplitudes $A_1^{(a,0,A)}$ and $A_3^{(h,0,A)}$ make towards the infection dynamics of the simulated wave, an amplitude space perspective was taken with the amplitude space defined by the four-dimensional space spanned by the amplitude variables $A_1,\dots ,A_4$ [4]. For each amplitude at each time point t the variance explained by that amplitude at that time point t was determined. More precisely, explained variance scores were computed like ${\mathrm{Score}}_j(t)=\mathrm{var}(A_j)(t)/{\sum }_{k=1}^4\mathrm{var}(A_k)(t)$, where $\mathrm{var}(A_j)(t)$ denotes the variance of the amplitude trajectory up to time point t (i.e., $\mathrm{var}(A_j)(t)={(T-1)}^{-1}{\sum }_{k=1}^{k^{*}}{[A_j(t_k)-M_{j,t}]}^2$ with $M_{j,t}$ being the mean $M_{j,t}=T^{-1}{\sum }_{k=1}^{k^{*}}{A}_j(t_k)$ with $t=t_{k^{*}}$ and $T=k^{*}$). Note that this is a time series framework where mean values and variances are computed from samples that consists of data taken from trajectories at discrete time points $t_k$.

Figure 1 presents some of the simulation results. Panels (a) and (b) show the four state variables $s_a$, $i_a$, $s_h$, and $i_h$ from top to bottom as solid black lines. Panel (a) shows the first 60 days. This period describes the simulated initial outbreak and the increase in the size of human infectious population towards its peak value. In contrast, panel (b) shows the total simulation period of 180 days and includes later stages of the simulated infection wave that describe the subsiding of the infection dynamics. As expected from the model equations of $s_a$ and $s_h$, the populations $s_a$ and $s_h$ decayed monotonically. In contrast, $i_a$ and $i_h$ formed infection waves. For the selected parameters, both populations reached peak values at approximately the same time.

Panel (c) shows the amplitudes $A_1,\dots ,A_4$ as functions of time during the entire 180 days simulation period with $A_2^{(12,A)}$ and $A_4^{(h,A)}$ given as solid black and gray lines, respectively, and $A_1^{(a,0,A)}$ and $A_3^{(h,0,A)}$ given as dotted black and gray lines, respectively. As can be seen in panel (c), during the first 60 days the human-environment order parameter amplitude $A_2^{(12,A)}$ (solid black) increased monotonically and played the dominant role among all four amplitudes. The human subsystem order parameter amplitude $A_4^{(h,A)}$ (solid gray) also varied over time during that 60 days interval but its variations were relatively small as compared to $A_2^{(12,A)}$. This is not surprising, because for the selected model parameters ${\lambda }_2$ was about 5.5 times larger than ${\lambda }_4$. The remaining two (neutrally stable) amplitudes $A_1^{(a,0,A)}$ and $A_3^{(h,0,A)}$ (dotted lines) stayed almost constant during the initial 60 days period. After 60 days, the amplitudes $A_1^{(a,0,A)}$ and $A_3^{(h,0,A)}$ started to make essential contributions to the infection dynamics. The values of $A_1^{(a,0,A)}$ and $A_3^{(h,0,A)}$ at the simulation stop of $t=180$ describe in good approximation the drop in susceptibles $s_a$ and $s_h$ as discussed in the context of Theorem 4 (compare panels (b) and (c)).

The relative importance of the four amplitudes during different stages of the simulated infection wave can also be explained with the help of their explained variance scores shown in panel (d). Panel (d) shows the explained variance scores for $A_2^{(12,A)}$ and $A_4^{(h,A)}$ given as solid black and gray lines, respectively, and $A_1^{(a,0,A)}$ and $A_3^{(h,0,A)}$ given as dotted black and gray lines, respectively. That is, in panel (d) the same color coding is used as in panel (c). Accordingly, during the first 60 days, the human–animal order parameter amplitude $A_2^{(12,A)}$ (solid black) explained most of the variance. During the first 20 days, the human subsystem order parameter amplitude $A_4^{(h,A)}$ (solid gray) also played a role such that during that period the two order parameter amplitudes taken together explain almost 100 percent of the variance in the infection dynamics. After 60 days, the explained variance score of $A_2^{(12,A)}$ decayed sharply, indicating that $A_2^{(12,A)}$ stopped playing the dominant role. For that later stages of the wave the neutrally stable amplitudes $A_1^{(a,0,A)}$ and $A_3^{(h,0,A)}$ explained most of the infection dynamics variance.

The gray dashed lines plotted in panels (a) and (b) show the state variable approximation ${\mathbf{x}}_{app}$ defined by Equation (52). As can be seen, the solutions for the infectious populations $i_a$ and $i_h$ as described by ${\mathbf{x}}_{app}$ were found to be identical to the exact solutions $i_a$ and $i_h$ (i.e., the gray dashed lines run exactly on top of the solid black lines). This illustrates the peculiarity of the SIR-SIR model expressed by Equation (48), namely, that the two potential order parameters (which are both actual order parameters for the simulated scenario (iii) outbreak) describe exactly the dynamics of $i_a$ and $i_h$. That is, as far as state variables $i_a$ and $i_h$ are concerned, ${\mathbf{x}}_{app}$ is not an approximation but an exact description of the infection dynamics. Moreover, during the initial period of 60 days ${\mathbf{x}}_{app}$ is a fair approximation of the state dynamics of $s_a$ and $s_h$ (see panel (a)). However, the simulation revealed that after that period ${\mathbf{x}}_{app}$ became a poor approximation of the dynamics of $s_a$ and $s_h$. This is consistent with the results presented in panels (c) and (d) and also illustrates what has been discussed previously in the context of Theorem 4. Panels (c) and (d) demonstrate that for the simulated Mpox infection wave the neutrally stable amplitudes $A_1^{(a,0,A)}$ and $A_3^{(h,0,A)}$ indeed became important after the initial phase passed by. Since ${\mathbf{x}}_{app}$ neglects these amplitude contributions and $A_1^{(a,0,A)}$ and $A_3^{(h,0,A)}$ are associated with eigenvectors that point into the direction of the susceptible populations $s_a$ and $s_h$, it is not surprising that ${\mathbf{x}}_{app}$ did not adequately capture the dynamics of the susceptibles of the simulated wave during the entire course of the wave.

In summary, it was found that the two order parameters and their amplitudes provided an exact description of the infectious human and animal populations (which is a peculiarity of the SIR-SIR model that does not hold for the remaining models B, C, and D). Furthermore, for the selected model parameters the order parameter approximation ${\mathbf{x}}_{app}$ also provided a good fit of the susceptible population dynamics during the initial increasing phase of the wave. For the selected model parameters, the human–animal order parameter amplitude $A_2^{(12,A)}$ made the main contribution, while the human subsystem order parameter amplitude $A_4^{(h,A)}$ made only a secondary contribution. At later time points the neutrally stable amplitudes $A_1^{(a,0,A)}$ and $A_3^{(h,0,A)}$ became important.

Finally, as mentioned above, for the simulated wave the human–animal order parameter amplitude dominated over the human subsystem order parameter amplitude. Therefore, changes in the population sizes of $i_a$ and $i_h$ should be determined approximately by ${\mathbf{v}}_2^{(12,A)}$ as described by Equation (43). Graphically speaking, the phase curve $i_h(i_a)$ in the 2D subspace spanned by $i_a$ and $i_h$ should follow the projection of the order parameter ${\mathbf{v}}_2^{(12,A)}$ into that subspace. Figure 2 shows the phase curve $i_h(i_a)$ of the simulated infection wave for the first 60 days (panel (a)) and for the entire simulation period (panel (b)) as solid black lines. The phase curves shown in Figure 2 were drawn from the solutions $i_a$ and $i_h$ presented in panels (a) and (b) of Figure 1. The phase curves were also drawn from the solution given by ${\mathbf{x}}_{app}$ (see the dashed gray lines). As expected (see Equation (48) again), the phase curves computed from ${\mathbf{x}}_{app}$ were identical to the phase curves obtained directly from the SIR-SIR model solutions. Importantly, panels (a) and (b) present the projection of the order parameter ${\mathbf{v}}_2^{(12,A)}$ into the 2D plane of $i_a$ and $i_h$ as red dotted lines. As can be seen in panel (a), the phase curve initially followed closely the order parameter ${\mathbf{v}}_2^{(12,A)}$. However, when the infection wave was about to reach the infection peak (i.e., towards the end of the 60-day period) the phase curve started to deviate from ${\mathbf{v}}_2^{(12,A)}$. The initial part during which the phase curve followed ${\mathbf{v}}_2^{(12,A)}$ is consistent with the scenario (ii) approximation shown in Equation (47). The deviation from ${\mathbf{v}}_2^{(12,A)}$ at the end of the 60-day period demonstrates the role of the second order parameter ${\mathbf{v}}_4^{(h,A)}$. In this context, note again that when taking this second order parameter ${\mathbf{v}}_4^{(h,A)}$ into account, we obtain the gray dashed line that runs on top of the exact solution. That is, the difference between a phase space dynamics along the red straight line and the actual dynamics indicated by the dashed gray line was entirely due to the contribution of the secondary order parameter: the human subsystem order parameter ${\mathbf{v}}_4^{(h,A)}$. Panel (b) demonstrates that for the subsiding part of the infection wave from 60 days to 180 days in crude approximation the phase curve $i_h(i_a)$ followed the direction defined by ${\mathbf{v}}_2^{(12,A)}$ (dotted red line). However, $i_h(i_a)$ deviated clearly from ${\mathbf{v}}_2^{(12,A)}$. As argued above, this deviation was due to the the term ${\mathbf{v}}_4^{(h,A)}{A}_4^{(h,A)}(t)$. In the interval from 60 to 180 days, the amplitude $A_4^{(h,A)}(t)$ formed a very shallow U-shaped curve (see the solid gray line in panel (c) of Figure 1). In line with this U-shape curve, the subsiding branch of the phase curve $i_h(i_a)$ in panel (b) of Figure 2 running from the top-right corner towards the fixed point $i_a^{*}=i_h^{*}=0$ first deviated slightly, subsequently reached a maximal deviation from ${\mathbf{v}}_4^{(12,A)}$, and finally approached again the direction specified by ${\mathbf{v}}_4^{(12,A)}$. In this context, note that ${\mathbf{v}}_4^{(12,A)}$ is attached in the 2D space at $\mathbf{x}(0)$ with $i_a(0)\ne 0$, $i_h(0)\ne 0$. Consequently, close to simulation stop at $t=180$ days when the phase curve $i_h(i_a)$ was about to approach the fixed point $i_a^{*}=i_h^{*}=0$ the phase curve crossed the ${\mathbf{v}}_4^{(12,A)}$-line.

4. Conclusions

Mpox outbreaks in endemic countries are typically initiated by animal to human transmissions of the Monkeypox virus [1,8]. However, once the Monkeypox virus has arrived in a particular human population, human to human transmissions may make a crucial contribution to the infection dynamics as well [1]. The current study is the first study that takes an order parameter approach that has been developed in the wake of the COVID-19 pandemic [4] to the field of Mpox infection outbreaks in order to address these two mechanisms. In four fundamental epidemiological models for Mpox outbreaks, two types of multi-compartmental components or building-blocks have been identified that shape and determine Mpox outbreaks. They are given in terms of human subsystem order parameters, on the one hand, and human–animal order parameters, on the other hand (see Section 3.1). The latter multi-compartmental components link the dynamics in animal subsystems rigidly with the dynamics in human subsystems during initial phases of Mpox outbreaks (see Section 3.2). The mathematical analysis revealed that each of the four fundamental models exhibits maximally only one human–animal order parameter and maximally only one human subsystem order parameter (see Theorems 1 and 2).

Previous work by Ma [42] pointed out the role of eigenvalues describing the temporal aspects of the initial phases of infectious disease outbreaks. Among other things, in the study by Ma [42] it has been assumed that the observable at hand, $X(t)$, increases exponentially according to a simple exponential function like $X(t)=X(0)exp(\lambda t)$. This situation holds for scenarios (i) and (ii) assuming that the approximations shown in Equation (47) hold. In contrast, as such, the single-parameter exponential increase model does not hold for scenario (iii) even under the simplified conditions described in Equation (47). Having said that, the simulation study presented in Section 3.4 showed that a scenario (iii) outbreak may effectively look like a scenario (i) or (ii) outbreak when one of the two relevant positive eigenvalues is relatively large with respect to the other. That is, the current study supplements the study by Ma [42] by identifying conditions under which single-parametric laws like $X(t)=X(0)exp(\lambda t)$ hold.

The current study focused on four fundamental Mpox models that effectively exhibit between four (SIR-SIR model) and six (SEIR-SEIR model) state variables. As mentioned in the introduction, more sophisticated models based on those four fundamental models have been proposed that feature relatively high-dimensional state spaces. In this context, Al-Shomrani et al. [25] proposed a 12-dimensional model based on the SEIR-SEIR model (6). Importantly, at the bifurcation point when the disease-free fixed point becomes unstable the authors were able to derive an analytical expression for the leading eigenvector (which in this case is the eigenvector associated with a zero eigenvalue). While the study by Al-Shomrani et al. was not tailored to discuss order parameters, it nevertheless shows that discussions as presented in the current study and calculations as carried out in Section 3 are not limited to the fundamental models defined by Equations (3)–(6): high-dimensional models may be analyzed in a similar way as the fundamental models have been analyzed in the current study.

Various studies as mentioned in the introduction have studied the impact of prevention and intervention measures on the spread of Mpox disease (see, e.g., [19,20,21,22,23]). In the context of early COVID-19 outbreaks in China [38], USA [40], and Thailand [39] during the year 2020, model-based analyses have found some evidence that intervention measures resulted in changes in the order parameters that shaped the respective outbreaks. Order parameters either changed their orientations [38,39] and/or qualitatively switched from unstable eigenvectors to stable eigenvectors [38,39,40] resulting in a subsiding of the COVID-19 surges of interest. The current study provides a basis to discuss similar impacts of prevention and intervention measures on the order parameters of Mpox infection models such as the fundamental models defined by Equations (3)–(6). While a detailed discussion is beyond the scope of this study, it should only be pointed out that in the context of Mpox waves intervention measures may result in a switch between different types of outbreak scenarios. More precisely, it is plausible to assume that in an endemic country witnessing a scenario (iii) outbreak, interventions may considerably reduce human to human transmission such that the eigenvalue ${\lambda }^{(h,X)}$ turns from a positive to a negative value. Since ${\lambda }^{(12,X)}$ would be not affected by such intervention measures, the outbreak may continue in the human subsystem (as long as ${\beta }_{12}>0$) in terms of a scenario (ii) outbreak. In doing so, a switch from a scenario (iii) to a scenario (ii) would occur. The remaining scenario (ii) outbreak could only be entirely stopped in the human subsystem by de-coupling the human population under consideration completely from the infectious animal reservoir (i.e., by administering intervention measures that lead to ${\beta }_{12}=0$).

In the four fundamental models that were examined in the current study, demographic terms were neglected. The reason for this was that the primary aim of the current study was to study the emergence of relative short-lived Mpox infection waves during which variations of population sizes due to birth and deaths can be assumed to be negligibly small. In future work, our analysis may be generalized to take the vital dynamics of populations into account. However, when demographic terms and/or variations in the total sizes of populations are taken into account, then the resulting epidemiological models typically increase in complexity, that is, they become higher-dimensional as compared to models that perform without those features. The aim of the current study was to discuss the main idea, namely, the existence of certain multi-compartmental components or building-blocks that shape and determine Mpox infection outbreaks with the help of analytical expressions. Such analytical expressions can be derived conveniently for relatively low-dimensional models as those discussed in the current study. In contrast, taking demographic terms and/or variations in population sizes into account may come at the cost that considerations have to be based on numerical approaches, that is, at the cost of losing concreteness.