Exact solutions and superposition rules for Hamiltonian systems generalizing time-dependent SIS epidemic models with stochastic fluctuations

1.   Introduction

The mathematical theory of epidemics, although formally established as an independent discipline through the pioneering work of Brownlee, Hamer, Lotka and Ross at the beginning of the 20th century, can actually be traced back to the 18th century, when Bernoulli first proposed the use of differential equations to study the population dynamics of infectious diseases (see ^[1,2] and references therein). Among the various models proposed from the 1930s onward, because they combine both analytical and statistical methods, the so-called compartmental models play a relevant role within epidemiological dynamics, with each of the variables or compartments (hence suggesting the terminology) representing a specific stage of contagion. Basic important examples are given by the well-known SIR models with three compartments.^* The time evolution of these variables is described in terms of a non-autonomous system of ordinary differential equations ^[3]. There are many variants to this approach, as studied by e.g., Kermack and McKendrick ^[4] and Abbey ^[5], who introduced additional constraints, as in the SIRS model ^[6], in which immunity only lasts for a short period of time, or the MSIR model, which takes into account additional assumptions on immunity (see ^[1] and references therein for further details). The simplest assumptions are given by the SIS model, in which the individuals in the population do not acquire immunity after infection, implying that they remain susceptible to being infected again ^[7]. Albeit, the model appears to be quite simple, and the introduction of fluctuations leads to more realistic realizations, although these additional constraints are not entirely obvious. One recurrent way to take into account fluctuations is to consider stochastic variables and hence to describe the model in terms of stochastic differential equations ^[8,9,10]. For the particular case of SIS models, the use of the Hamiltonian formalism has led to an improved set of differential equations for the mean and variance of infected individuals. In ^[7], the spreading of the disease is treated through the use of Markov chains, in which at most one single recovery or transmission occurs in an infinitesimal interval, and by using the mean density of individuals and the variance as dynamical variables. The fluctuations that were included by following this approach correspond to a truncated stochastic expansion, where higher statistical moments are neglected. Other alternative approaches to epidemiological models based on the Hamiltonian machinery can be found in ^[11,12,13], where exact solutions of an SIR model have been obtained in ^[14,15]. An approach without a Hamiltonian perspective can be found in ^[16,17]. The Hamiltonian model in ^[7], in spite of the stochastic nature of the fluctuation, can still be treated as a deterministic system. Genuine stochastic (Hamiltonian) systems, i.e., stochastic differential equations defined on complete probability spaces, have been studied extensively ^[10,18,19] and successfully applied to SIS models (see ^[20,21] and references therein). In this context, the consideration of stochastic parameters/variables allows one to apply the formalism to other systems that can be modeled as a contagion-like process, such as information transmission models, networks or one-dimensional diffusion processes ^[21,22]. This fact legitimizes the idea that SIS models can be generalized by introducing time-dependent parameters.

^* Here, "S" stands for the individuals susceptible to the disease, "I" designates the infected individuals and "R" stands for the recovered ones.

In ^[13], the authors showed how the machinery of Lie-Hamilton (LH in short) systems ^{[23,24,25,26,27]} can be applied to the study of SIS models, providing a superposition principle for a time-dependent generalization of the model studied in ^[7]. From the perspective of LH systems, further generalizations of the underlying Hamiltonian system that preserve the property that one parameter is identified with a variable infection rate are conceivable, as they provide formal extensions of these epidemic models once the additional functions have been identified with relevant epidemiological parameters. In this sense, the time-dependent Hamiltonians considered in this work introduce some semi-stochatic effects by means of fluctuations. It is not our purpose to analyze the epidemiological validity of such extensions, but to present the formal analytical frame to which such models are subjected, and to study the existence of exact solutions and, alternatively, the existence of nonlinear superposition principles that can be computationally more efficient than a direct attempt to integrate the system.

We recall that Lie systems ^[27,28,29] are mainly characterized by the existence of (nonlinear) superposition rules, allowing a complete description of the solution of a system of differential equations in terms of a certain number of particular solutions and significant constants. The superposition principles are themselves deeply related to the existence of a finite-dimensional Lie algebra of vector fields associated with the system, which further leads to a realization of the system in terms of a sum of the Lie algebra generators with time-dependent coefficients ^[27]. As shown in ^[30], using the integrability properties of distributions generated by vector fields, the formalism can be extended to the case of stochastic differential equations, leading to a stochastic version of the Lie-Scheffers theorem and corresponding (local) superposition rules.

A particularly interesting case arises whenever the vector fields spanning the Lie algebra are Hamiltonian vector fields with respect to a symplectic structure, which is a feature that allows us to treat the system as a classical Hamiltonian system ^[31]. Systems of this type, called LH systems, have been extensively studied in low dimensions (see ^{[23,24,25,26,27]} and references therein) and constitute a powerful technique for the analysis of differential equations. One of the advantages of the compatibility with a symplectic structure resides in the fact that the so-called coalgebra formalism of (super-)integrable systems ^[32] can be applied to systematize the construction of time-independent constants of the motion ^[24], which finally provides a superposition principle through algebraic computations, thus avoiding the cumbersome integration of systems of ordinary or partial differential equations (ODEs and PDEs in short, respectively).

As mentioned, the main objectives of this work are to embed the SIS Hamiltonian models of ^[7,13] into larger Hamiltonian systems such that the underlying symplectic structure is preserved, and to determine to which extent exact solutions of these formal generalizations can be obtained, alternatively applying the techniques of LH systems to derive suitable superposition rules. The potential interpretation of the various new time-dependent coefficients introduced in this way is not considered, nor is the applicability to systems beyond epidemiological models, such as the spread of information over networks or numerical methods. In Section 2, we briefly review the SIS Hamiltonian model proposed by Nakamura and Martinez in ^[7], which possesses a constant infection rate. In Section 3, we reconsider the generalization proposed in ^[13] based on the introduction of a variable infection rate. The striking point is that the $t$-dependent Hamiltonian studied in ^[13] gives rise to a natural generalization by adding a second arbitrary $t$-dependent parameter in such a manner that the resulting Hamiltonian can be related to the so-called book LH algebra $\mathfrak{b}_2$ ^[26,33]. The classification of LH systems on the plane $\mathbb R^2$ ^[25,26] allows us to obtain a canonical transformation between the generic Cartesian coordinates $(x, y)$ and the mean density of infected individuals $\langle\rho\rangle$ and the variance $\sigma^2$ of this generalized Hamiltonian model. Using the geometric and algebraic properties of LH systems, we derive an exact explicit solution for this book Hamiltonian, which holds, in particular, for the model proposed in ^[13], where no exact solution was given.

Using the subalgebra lattice of the classification of planar LH systems ^[25,26], a second extension of the book LH algebra $\mathfrak{b}_2$ to the oscillator algebra $\mathfrak{b}_2 \subset\mathfrak{h}_4$ is considered by introducing an additional $t$-dependent parameter. This case is considered in Section 4 in full detail. In particular, we present an exact solution for this new oscillator Hamiltonian, providing, alternatively, a superposition rule that, albeit not required because of the exact solution, will be convenient in the context of extensions of the system to higher-dimensional LH systems. In Section 5, we consider the maximal embedding chain (within the LH classification) $\mathfrak{b}_2 \subset\mathfrak{h}_4 \subset\mathfrak{h}_6$ that leads to the so-called two-photon Lie algebra $\mathfrak{h}_6$ ^[34,35]. Although, in this case, the first-order system of differential equations is still linear (in the Cartesian coordinates) with variable coefficients, a direct integration will be, generally, no more possible or computationally feasible due to its equivalence with a reduced system of Riccati type. This justifies the construction of a superposition principle, eventually allowing a description of the solution in terms of three particular solutions.

Finally, some open problems and future prospects are discussed in Section 6. For completeness of the exposition, in the Appendix, we briefly recall the fundamental properties of LH systems, the technical details of which can be found in ^{[23,24,25,26,27]}.

2.   SIS Hamiltonian model with stochastic fluctuations

As mentioned previously, one of the main premises of the SIS model is that all individuals are still susceptible to the infection after recovery, which implies that they do not acquire immunity. From this condition, we conclude that the model can be described by only using two compartmental variables. The first one, $I$, corresponds to the number of infected individuals, whereas the second compartment $S$ describes the number of individuals susceptible to the infection at a given time. The model assumes a large population of size $N$, a random mixture of the population, a regular and stationary age distribution and one single contaminating agent. The chances of infection for any individual extend through the whole course of the epidemic. In this context, one relevant variable is given by the density $\rho = \rho(\tau)$ of infected individuals, depending on the time parameter $\tau$ that takes values in the interval $[0, 1]$. The density of infected individuals decreases according to a rate $\gamma \rho$, with $\gamma$ denoting the recovery rate, while the growth rate of new infections is proportional to $\alpha \rho(1-\rho)$, with a transmission rate parameter $\alpha$. The equation that describes the variation of infected individuals, taking into account both contribution factors, is given by

One can redefine the timescale as $t : = \alpha \tau$ (where $\alpha\neq 0$) and introduce the constant $\rho_0: = 1-\gamma/\alpha$, so we can reformulate Eq (2.1) as

Clearly, the equilibrium density is reached if either $\rho = 0$ or $\rho = \rho_0$.

One major drawback of the model is that it assumes that each individual of the population lives to a certain maximal age $L$, and that, for each age $a < L$, the number of individuals of age $a$ is the same. While such a homogeneity assumption seems feasible in developed countries, where infant mortality is quite low, the hypothesis cannot be realistically assumed to hold for developing countries. To circumvent this difficulty, it is reasonable to introduce probability functions at some point, in order to allow a random variation over time in one or more of the inputs. Some experimental results show conclusively that temporal fluctuations can drastically modify the prevalence of pathogens and the spatial heterogeneity (see ^[36,37]).

In order to introduce fluctuations in the SIS model, we follow the general ansatz proposed by Nakamura and Martinez in ^[7], where the spreading of the disease is interpreted as a Markov chain in discrete time with at most one single recovery or transmission occurring in each infinitesimal interval. Therefore, the respectively equations for the instantaneous mean density of infected individuals $\langle\rho\rangle$ and the variance $\sigma^2 = \langle\rho^2\rangle-\langle\rho\rangle^2$ are ^[38]

where $\Delta_3 = \langle\rho^3\rangle-\langle\rho\rangle^3$ and the variance in the first equation slows down the growth rate of $\langle\rho\rangle$, recalling the Allee effect ^[39]. Equations (2.2) and (2.3) are equivalent whenever $\sigma$ becomes negligible as compared to $\langle\rho\rangle$. If, in addition, we ignore higher statistical moments, the resulting dynamical system describes a Gaussian variable evolving over time. This implies that $\Delta_3 \simeq3\sigma^2\langle\rho \rangle$; thus for a sufficiently large number of individuals $N$, the resulting equations read as follows

With the introduced fluctuations, Eq (2.4) corresponds to a stochastic expansion, according to ^[40]. In this situation, we assume that the density $\rho = \langle \rho\rangle +\eta$ is properly described by the instantaneous average, as well as some noise function $\eta$. For consistency, it is further assumed that $\langle \eta\rangle = 0$ and $\langle \eta^2\rangle = \sigma^2$. The latter system of Eq (2.4) allows for a Hamiltonian formulation, with the phase-space variables given by the mean density of infected individuals $\langle\rho\rangle$ and the variance $\sigma^2$. Considering

as dynamical variables ^[7], the system (2.4) adopts the following form:

It is straightforward to verify that these are the canonical equations associated with the Hamiltonian

As shown in ^[7], the system (2.4) can be solved exactly, providing the following solution:

where $\tilde c_1$ and $\tilde c_2$ are two constants depending on the initial conditions. In terms of the canonical variables $(q, p)$ in Eq (2.5), the solution of the system (2.6) adopts the following expression:

For a detailed analysis of all of the above results, we refer the reader to ^[7].

3.   Generalization of the SIS Hamiltonian model with a variable infection rate

Quite recently, a generalization of the Hamiltonian (2.7) was proposed in ^[13] by considering a $t$-dependent infection rate through the used of a smooth function $\rho_0(t)$, which amounts to introducing a $t$-dependent basic reproduction number $R_0(t)$ that is actually observed in more accurate epidemic models ^[41,42,43]. The remarkable point, as was explicitly shown in ^[13], is that the resulting $t$-dependent Hamiltonian inherits the structure of an LH system ^{[23,24,25,26,27]}. In what follows, we summarize the main results.

If we consider that $\rho_0 = \rho_0(t)$ and insert it into Eq (2.6), we obtain the following system:

The dynamics of the resulting system is easily described in terms of a $t$-dependent vector field (see (A.7) in the Appendix):

and such that the commutation rule

holds. This implies that the generalized SIS model (3.1) determines a Lie system ^[27,28,44], the Vessiot-Guldberg algebra of which is isomorphic to the so-called book algebra, here denoted by $\mathfrak{b}_2$. In this context, ${\bf X}_A$ can be regarded as a dilation generator, while ${\bf X}_B$ corresponds to a translation generator.

In addition, as $(q, p)$ are canonical variables, we can consider the usual symplectic form in the phase-space:

and compute the corresponding Hamiltonian functions $h_A, \, h_B$ associated with the vector fields ${\bf X}_A, \, {\bf X}_B$ through the contraction or inner product of $\omega$, that is,

The functions $h_A, \, h_B$ are given by

and satisfy the Poisson bracket

Thus ${\bf X}_A, \, {\bf X}_B$ are Hamiltonian vector fields, and the $t$-dependent Hamiltonian of the generalized SIS model (3.1) reads as follows (see (A.11)):

reproducing the Hamiltonian function (2.7) proposed in ^[7] for a constant $\rho_0$. This shows that both cases can be studied by means of a unified geometrical approach. In consequence, the generalized SIS model (3.1) is not merely a Lie system, but also an LH one ^{[23,24,25,26,27]}, as was already proved in ^[13].

3.1.   Generalized SIS Hamiltonian from the book algebra

By taking into account the LH theory, the SIS Hamiltonian (3.8) can be further generalized through the introduction of a second arbitrary $t$-dependent parameter $b(t)$ in the form

leading to the following system of differential equations:

Clearly, the associated $t$-dependent vector field is just

with ${\bf X}_A$ and ${\bf X}_B$ given in Eq (3.2).

Summing up, for any choice of the parameters $\rho_0(t)$ and $b(t)$, the differential equations described by Eq (3.10) always determine an LH system with an associated $\mathfrak{b}_2$-LH algebra. For this reason, we can say that Eq (3.10) is a generalized (time-dependent) Hamiltonian from the book LH algebra that includes the SIS epidemic model (3.1) studied in ^[13] for the values $b(t)\equiv 1$, as well as the model (2.6) introduced in ^[7] for the values $\rho_0(t)\equiv \rho_0$ and $b(t)\equiv 1$.

We recall that the $\mathfrak{b}_2$-LH algebra appears within the classification of planar LH systems presented in ^[25,26] as the class I$_{14A}^{r = 1}\simeq \mathbb{R} \ltimes \mathbb{R}\simeq \mathfrak{b}_2$. Let us consider the Hamiltonian vector fields spanning the $\mathfrak{b}_2$-LH algebra in terms of the usual Cartesian coordinates $(x, y)\in \mathbb R^2$, as given in ^[33]:

These verify the Poisson bracket (3.7) with respect to the standard symplectic form

The vector fields associated with (3.12) are obtained by using Eq (3.5) and given by

satisfying the commutator relation given by Eq (3.3). From the classification of LH systems, we deduce the existence of a canonical transformation between the variables $(q, p)$ in Eqs (3.2) and (3.6) and $(x, y)$ in Eqs (3.12) and (3.14), respectively:

preserving the symplectic form

It is worthy to be mentioned that the canonical transformation (3.15) is quite different from the change of variables considered in ^[13] (cf. Eq (34)), as the latter does not preserve the symplectic form. Observe that, in the absence of the canonical character of such a transformation, the variables $q$ and $p$ used in ^[13] do not correspond to canonical dynamical variables; hence their direct interpretation as the mean density of infected individuals $\langle\rho\rangle$ and the variance $\sigma^2$ given by Eq (2.5) is no more entirely obvious.

3.2.   Exact solution for the book SIS Hamiltonian

The remarkable fact is that the system (3.10) becomes separable in the coordinates $(x, y)$, as can be routinely verified:

The equations are uncoupled and can be solved in a straightforward way by quadrature. Considering again the canonical transformation (3.15), the corresponding solution of the initial system (3.10) is obtained.

The explicit solution of the linear system (3.17) is given by

where $c_1$ and $c_2$ are the two constants of integration, determined by the initial conditions, while $a$ is a real number that ensures the existence of the integrals over the interval $[a, t]$. By introducing the transformation (3.15), we arrive at the general solution of Eq (3.10) for arbitrary $t$-dependent parameters $\rho_0(t)$ and $b(t)$:

Recall that the resulting mean density of infected individuals and the variance of the model are just $\langle\rho(t)\rangle = q(t)$ and $\sigma^2(t) = 1/p^2(t)$, respectively (see Eq (2.5)). Observe also that, as a byproduct of the general result above, the exact solution of the system (3.1) considered in ^[13] is easily obtained by setting $b(t)\equiv 1$ and keeping a variable $\rho_0(t)$, a fact that complements the results of that work.

Finally, we fix $b(t)\equiv 1$ and consider a constant $\rho_0(t)\equiv \rho_0$. Then, we can choose, for instance, $a = 0$ in the integral for ${\Theta}(t)$ of Eq (3.18), that is, ${\Theta}(t) = \rho_0 t$. Hence, the solution given by Eq (3.19) reduces to

recovering, as expected, the solution given by Eq (2.9) of the system (2.6) studied in ^[7], provided that the constants of integration are redefined as

4.   Generalized Hamiltonian from the oscillator algebra

The classification of LH systems ^[25,26] shows that the two-dimensional book LH algebra appears as a Lie subalgebra of other classes. This allows us to further extend the results of the previous section by considering higher-dimensional LH algebras that entail the introduction of additional $t$-dependent parameters $b_i(t)$, and such that the book SIS system (3.10) is recovered for the vanishing values of these parameters. Among the various LH algebras into which the book algebra can be embedded, the natural candidate is the class I$_{8}$. As a Lie system, this only requires the addition to $X_A$ and $X_B$ in Eq (3.14) of a new vector field. Keeping the same notations of ^[25,26], I$_{8}$ can be expressed in terms of the Cartesian coordinates $(x, y)$ such that the the three vector fields are given by

which obeys the commutation rules

and thus closing on a Lie algebra isomorphic to the (1+1)-dimensional Poincaré algebra $\mathfrak{iso}(1, 1)$. The Lie system ${\bf X}$ determined by the $t$-dependent vector field

leads to the following (linear) system of differential equations:

The generators ${\bf X}_i$ (4.1) then become Hamiltonian vector fields $h_i$ with respect to the symplectic form given by Eq (3.13), namely,

Observe that, in order to ensure that these functions close with respect to the Lie-Poisson bracket, it is necessary to add the central generator $h_0$:

It follows that the resulting LH algebra is isomorphic to the centrally extended Poincaré algebra $\overline{\mathfrak{iso}}(1, 1)$, which is also isomorphic to the oscillator algebra $\mathfrak{h}_4$. Within this identification, $h_3$ can be interpreted as the number operator, while $h_1$, $h_2$ can be seen as ladder operators. Alternatively, the differential equation (4.4) can also be obtained through the Hamilton equations with canonical variables $(x, y)$ from the Hamiltonian given by

The canonical transformation (3.15) allows us to express the vector fields described by Eq (4.1) and the Hamiltonian vector fields described by (4.5) in terms of the canonical variables $(q, p)$ in Eq (2.5):

This procedure leads to a genuinely coupled nonlinear first-order system of differential equations that generalizes the book SIS system (3.10) with an additional parameter $b_1(t)$:

This system can be regarded as a generalized time-dependent Hamiltonian from the oscillator LH algebra that extends the SIS epidemic model (3.10) for the embedding chain $\mathfrak b_2\subset \mathfrak h_4$, or, alternatively, as a formal Hamiltonian system into which the SIS epidemic model has been embedded.

4.1.   Exact solution for the oscillator time-dependent Hamiltonian

In spite of the apparently complicated structure of the system of oscillatory type (4.9), an exact solution can still be found by taking into account the previous cases and the fact that, by means of the canonical transformation (3.15), the system decouples in the new variables, while still preserving the symplectic form $\omega$ and, hence, the underlying geometric structure. The system (4.4), being linear, can be easily integrated and gives

where, again, $c_1$ and $c_2$ are the two constants of integration determined by the initial conditions. By applying the canonical transformation (3.15), we obtain the exact solution of Eq (4.9), which adopts the form

It can be trivially verified that this solution reduces to the solution given by Eq (3.19) of the book SIS system for $b_1(t)\equiv 0$ and $b_2(t)\equiv b(t)$.

It should be mentioned that the integrability by quadratures of the systems (3.10) and (4.9) can actually be inferred from the structural properties of the book and oscillator algebras as Vessiot-Guldberg algebras of a Lie system. As shown in ^[45] (see also ^[46,47,48] and references therein), if the Vessiot-Guldberg algebra is solvable, then the integrability by quadratures is guaranteed.

4.2.   Superposition rule for the oscillator time-dependent Hamiltonian

As both the generalized book and oscillator systems determined by the differential equations comprising Eqs (3.10) and (4.9) are classical Lie systems, they always admit a (nonlinear) superposition rule ^[44] (see (A.2)). As these equations also determine LH systems, an explicit superposition rule can be found in terms of $t$-independent constants of the motion constructed by means of the coalgebra formalism presented in ^[24] (see also ^[26,33]), and as based on the Casimir operators (see the Appendix). This formalism solely fails for the book algebra $\mathfrak{b}_2$, as this Lie algebra does not admit nonconstant Casimir invariants. However, using appropriate embeddings of the book algebra, such as $\mathfrak{b}_2\subset \mathfrak{h}_4$, the exact solutions given by Eqs (3.19) and (4.11) for both the book and oscillator SIS systems can be obtained. According to the general features of LH systems, for the oscillator Hamiltonian (4.7), we can find a superposition rule, which, particularized to the book SIS Hamiltonian (3.9), would provide a superposition rule for the latter. Albeit, deriving a superposition principle for the system (4.9) is redundant; as the system can be solved explicitly, it is advantageous to compute it for comparison purposes. As will be shown in the next section, the oscillator time-dependent $\mathfrak{h}_4$-Hamiltonian (4.7) appears itself as a special case of another and last extension, that is, the two-photon $\mathfrak{h}_6$-system for which, although an explicit solution of the associated system is formally conceivable, its computational implementation is probably too cumbersome due to a nontrivial coupling.

We now proceed to compute $t$-independent constants of the motion for the oscillator SIS system (4.9) and deduce a superposition rule. Let us consider the oscillator LH algebra expressed in a basis $\{ v_1, v_2, v_3, v_0 \}$ that formally satisfies the same Poisson brackets described by Eq (4.6):

There exists a non-trivial quadratic Casimir element given by

In order to derive a superposition rule for this case, we require $m = 2$ particular solutions and $n = 2$ significant constants (see (A.2)), and, hence, the indices $k = 2, 3$ in the application of the coalgebra formalism ^[26] (see (A.13)). From (4.5), we construct the Hamiltonian vector fields $h^{({k})}_i$ in Cartesian coordinates $(x, y)\in \mathbb{R}^2$:

Each set satisfies the Poisson brackets described by Eq (4.6) with respect to the symplectic form

From this result, using the Casimir invariant given by Eq (4.13), we obtain two $t$-independent constants of the motion $F^{(k)}$ (A.14) for the diagonal prolongation $\widetilde{ {\bf X}}^3$ to $(\mathbb{R}^2)^3$ ^[26,33]:

$F^{(2)}$ and $F^{(3)}$ are functionally independent and in involution in $C^\infty\bigr((\mathbb{R}^2)^3 \bigl)$ with respect to the symplectic form (4.15), that is,

We observe, in particular, that the function $F = C\left(h^{(1)}_0, h^{(1)}_1, h^{(1)}_2, h^{(1)}_3\right)$ vanishes identically (see (A.14)). Two additional constants of the motion can be deduced by taking into account the permutation symmetry of the indices in $F^{(2)}$ (see (A.16)):

so that $F^{(3)} = F^{(2)}+ F_{13}^{(2)}+ F_{23}^{(2)}$.

The formalism thus provides us with four constants of the motion described by Eqs (4.16) and (4.18) for the oscillator Hamiltonian (4.7), which can be identified by using some of the significant constants $k_i$ within the superposition rule (A.2) of interest. We set

There are several alternatives for the choice of the two necessary constants $k_i$ in order to obtain a superposition rule ^[33]. For computational reasons, we make the choice $F^{(2)} = k_1$ and $F^{(3)} = k$ in order to express $(x_1, y_1)$ in terms of $(x_2, y_2, x_3, y_3)$, along with the two constants $k_1, k$. After some routine algebraic manipulations, we arrive at the following nonlinear superposition rule:

where

Finally, we apply the canonical transformation (3.15) to the above results while keeping the canonical symplectic form given by (4.15):

As a result, we get $t$-independent constants of the motion described by Eq (4.16) for the oscillator SIS Hamiltonian (4.7), expressed in terms of the Hamiltonian vector fields described by Eq (4.8) with the canonical SIS variables given by (2.5):

The corresponding superposition rule for Eq (4.9) turns out to be

with $B$ given in Eq (4.21) and

Summarizing, we have obtained the general solution $(q_1, p_1)$ of the oscillator system of differential equations given by Eq (4.9) that generalize the Hamiltonian (3.8) in terms of two particular solutions $(q_2, p_2)$ and $(q_3, p_3)$ and two significant constants $k_1$ and $k$ in an explicit form. We remark that the superposition rule given by Eq (4.24) holds for any $t$-dependent parameters $b_1(t)$, $b_2(t)$ and $\rho_0(t)$ within the oscillator SIS Hamiltonian (4.7); thus this result can also be regarded as a superposition rule for the book SIS system (3.10) as the particular case with $b_1(t) = 0$.

5.   Generalized Hamiltonian from the two-photon algebra

Following the classification of LH systems ^[25,26] on the plane $\mathbb R^2$, the embedding $\mathfrak{b}_2\subset \mathfrak{h}_4$ can still be extended up to a maximal chain $\mathfrak{b}_2\subset \mathfrak{h}_4\subset \mathfrak{h}_6$, with $\mathfrak{h}_6$ (class P$_5$ in ^[25]) being the so-called two-photon algebra ^[34,35], which corresponds to the highest-dimensional non-solvable LH algebra on the plane. This embedding of the oscillator algebra hence leads to a maximal extension of the book SIS system through this chain. Specifically, we have to add two vector fields ${\bf X}_4$ and ${\bf X}_5$ to the oscillator basis described by Eq (4.1), which, in the coordinates $(x, y)$, read as follows:

These generators satisfy the following commutation rules:

which determine a 5-dimensional Lie algebra isomorphic to $\mathfrak{sl}(2) \ltimes \mathbb{R}^2$, where $\mathfrak{sl}(2) = {\rm{span}}\{ {\bf X}_3, {\bf X}_4, {\bf X}_5\}$ and $\mathbb{R}^2 = {\rm{span}}\{ {\bf X}_1, {\bf X}_2 \}$. The corresponding Lie system ${\bf X}$ is now determined by the $t$-dependent vector field

giving rise to the following system of differential equations:

Therefore, these equations generalize the oscillator Lie system (4.4). Albeit, the system is still linear for $b_4(t) b_5(t)\ne 0$, so an explicit integration is quite cumbersome, and, from the computational point of view, rather ineffective. Explicitly, from the first equation in Eq (5.4), we find that

which leads to the second-order inhomogeneous linear equation

where

In order to solve the equation, a particular solution of the inhomogeneous equation and the general solution of the homogeneous part must be computed. The homogeneous part itself can be reduced by the change of variables

to a Riccati equation

which, in general, cannot be solved explicitly (i.e., the solution may not be expressible in terms of the classical functions). Alternatively, Eq (5.6) can be reduced to a free second-order equation $z^{\prime\prime}(\tau) = 0$ by using the Lie symmetry method ^[49]. However, the general solution is still far from being easily obtained due to the complicated form of the linearizing point transformation. These computational difficulties justify the use of a superposition principle for the system (5.4), with the additional advantage that the results can be easily translated to the coordinates $(q, p)$ via the canonical transformation (3.15). We further observe that the system in these coordinates is far from being linear.

Starting from the realization of Eq (5.1), it can be easily verified that the generators ${\bf X}_i$ are Hamiltonian vector fields $h_i$ with respect to the symplectic form described by Eq (3.13), with the following functions:

Again, a central generator $h_0$ must be added to obtain a Lie-Poisson algebra, where the Poisson brackets are given by

The resulting 6-dimensional LH algebra is a central extension of $\mathfrak{sl}(2) \ltimes \mathbb{R}^2$, isomorphic to the two-photon Lie algebra $\mathfrak{h}_6$ ^[34,35] and further equivalent to the (1+1)-dimensional centrally extended Schrödinger Lie algebra ^[50]. Clearly, $\mathfrak{h}_6$ contains the book $\mathfrak{b}_2 = {\rm{span}}\{ h_2, h_3\}$ and the oscillator $\mathfrak{h}_4 = {\rm{span}}\{h_1, h_2, h_3, h_0\}$ LH algebras, as well as the simple Lie algebra $\mathfrak{sl}(2) = {\rm{span}}\{ h_3, h_4, h_5\}$. Hence, we have the inclusions

The same system (5.4) of differential equations can alternatively be deduced from the Hamilton equations with canonical variables $(x, y)$ by using the Hamiltonian associated with the vector field ${\bf X}$ (5.3):

Similar to the previous section, the canonical transformation (3.15) leads to the expression of the vector fields described by Eq (5.1) and the Hamiltonian vector fields described by Eq (5.10) in terms of the canonical variables $(q, p)$ given by Eq (2.5), namely,

The corresponding first-order system that generalizes the oscillator SIS system (4.9) adds two new parameters $b_4(t)$ and $b_5(t)$:

This system is a formal generalization of the SIS epidemic time-dependent Hamiltonians associated with both the book and oscillator SIS systems (3.10) and (4.9). Whether this system truly corresponds to a realistic epidemic model (or other processes that can be modeled, like contagions) depends on whether the added parameters can be properly identified by using some fluctuation parameters. Our purpose is not to analyze these possibilities, but rather to provide a procedure that allows one to obtain solutions of such systems by using a superposition principle determined by the LH structure. It is worthy to be observed that, as a byproduct of the embeddings given by Eq (5.12), this method also provides a $t$-dependent system for the $\mathfrak{sl}(2)$-LH algebra whenever $b_1(t) = b_2(t) = 0$, although this case cannot be considered as a realistic generalization of the SIS epidemic Hamiltonian of Eq (2.6) due to the vanishing of $b_2(t)$.

In contrast to the previous book and oscillator Hamiltonian systems, formerly written in Cartesian coordinates $(x, y)$ in Eqs (3.17) and (4.4), the two-photon system (5.4) is nontrivially coupled, which prevents us from finding a solution by quadratures for the case that $b_4(t)b_5(t)\neq 0$ (see Eq (5.6)). However, as already mentioned, we can deduce a superposition rule and then apply the canonical transformation (3.15) to obtain a superposition rule for the system with the variables $(q, p)$ given by Eq (2.5). To this extent, we follow the same steps as used in Section 4.2 for the oscillator time-dependent Hamiltonian (the more technical details are described in the Appendix).

We start by considering the $\mathfrak{h}_6$-LH algebra in a basis $\{ v_1, \dots, v_5, v_0 \}$ that satisfies the Poisson brackets described by Eq (5.11). Besides the trivial central generator, the two-photon Lie-Poisson algebra possesses a third-order Casimir invariant given by (see ^[26,35])

This implies that a superposition rule will be given in terms of $m = 3$ particular solutions and $n = 2$ significant constants; hence, the relevant indices for the application of our formalism (see (A.13)) are $k = 2, 3, 4$. The Hamiltonian vector fields $h^{({k})}_i$ are easily obtained from Eq (5.10), being explicitly given by

These functions satisfy the Poisson brackets given by Eq (5.11) with respect to the symplectic form

Now, introducing Eq (5.17) into the Casimir (5.16), we find $t$-independent constants of the motion described by (A.14) for the diagonal prolongation $\widetilde{ {\bf X}}^4$ to $(\mathbb{R}^2)^4$. These are $F = F^{(2)} = 0$ and

where we have made use of the permutation symmetry of the variables for $F^{(3)}$.

In this case, a convenient choice for the $n = 2$ significant constants intervening in the superposition rule is given by the functions $F^{(3)} = k_1^2$ and $F_{34}^{(3)} = k_2^2$, where $k_1$ and $k_2$ are constants. In terms of these, after some algebraic manipulation, we can write $(x_1, y_1)$ as the general solution of the system (5.4) in terms of three particular solutions $(x_2, x_2)$, $(x_3, x_3)$ and $(x_4, x_4)$, together with the two constants $k_1$ and $k_2$, in the following form:

where, in shorthand notation,

At this stage, we apply the canonical transformation (3.15) to the above results, always preserving the symplectic form given by Eq (5.18):

Then, we directly deduce from Eq (5.19) the $t$-independent constants of the motion in the canonical variables described by Eq (2.5). In particular, the constant $k_1$ such that $F^{(3)} = k_1^2$ now becomes

The resulting superposition rule for the two-photon system (5.15), although cumbersome, determines the general solution $(q_1, p_1)$ in terms of rational functions of the three particular solutions $(q_2, p_2)$, $(q_3, p_3)$ and $(q_4, p_4)$ and the significant constants $k_1$ and $k_2$, being explicitly given by

where

The relevant observation at this point is that this superposition principle remains unaltered for any choice of the $t$-dependent parameters $b_i(t)$ and $\rho_0(t)$ appearing in the Hamiltonian (5.13); thus, it holds for all of the Hamiltonian systems considered previously, and, in particular, for the SIS models (2.6) and (3.1). This means specifically that, if any of the generalized systems turns out to be experimentally identified with a realistic model, the main properties of the original systems remain unaltered.

6.   Conclusions

In the context of SIS epidemic models, the Hamiltonian formulation has shown to be an effective technique to unveil the main features concerning the dynamics of such systems, even in those cases where a closed expression of the solutions is not possible. One of such models ^[7] has recently been reconsidered from this point of view, which has further led to a generalization involving a time-dependent infection rate ^[13]. Both systems correspond to a single class of LH systems, allowing a unified approach. Taking into account that the LH algebra that underlies these models can be embedded into higher-dimensional LH algebras on the real plane, in this work, we have systematically analyzed the formal extensions of the SIS models to first-order systems of differential equations depending on additional time-dependent parameters. Moreover, it has been shown that any Hamiltonian system (whether identifiable with an SIS model or not) based on the book algebra $\mathfrak{b}_2$ admits an exact solution in closed form, which, particularized to the above-mentioned models, reproduces the exact solution given in ^[7] and completes the analysis in ^[13], providing the general solution of the latter. The general $\mathfrak{b}_2$-LH system already formally generalizes these systems, as it introduces a second parameter, the interpretation of which has, however, not been considered. In addition, using the embedding chain $\mathfrak{b}_2\subset \mathfrak{h}_4\subset \mathfrak{h}_6$ of LH algebras ^[25], further formal extensions of the models have been proposed, up to a maximal case depending on five arbitrary time-dependent parameters. While the systems based on the oscillator LH algebra $\mathfrak{h}_4$ have been shown to admit an exact solution (this case additionally completing the results in ^[13], where a superposition rule was given but with non-canonical variables), the computational obstructions to effectively integrate the corresponding LH system for $\mathfrak{h}_6$ suggests the derivation of a superposition principle, the existence of which is guaranteed by the properties of LH systems. This case also represents the maximal possible extension of the time-dependent Hamiltonian models based on the assumptions of ^[7,13]. This means that, once the additional parameters have been identified with suitable characteristics of an epidemic model, the formalism developed in this work immediately provides either the general solution of the system, or a superposition rule that depends at most on three particular solutions.

As mentioned, we have focused primarily on the analytical and geometrical properties of LH systems that contain the model (3.1), without pursuing further the viability of such models at the epidemiological or biological level. As the method developed allows one to derive solutions regardless of the particular significance of the parameters, for any appropriate epidemiological extension, the main properties are directly obtained by insertion of the values of the parameters in the corresponding solutions. In this context, the relevant question that emerges toward realistic applications is whether these additional arbitrary time-dependent functions can effectively be identified by using some compartmental parameters subjected to fluctuations, in the line of the first generalization proposed in ^[13]. In this sense, each possible choice of the parameters must be studied separately, with a proper contextualization of the parameters, based on experimental data. The epidemiological applicability of the generalized LH systems remains unknown.

It may be asked whether the LH formalism can be applied to other epidemic models, specifically in order to obtain exact solutions. The models considered in ^[51,52] lead, for certain identifications of the parameters, to differential equations very similar to either the matrix or projective Riccati equations. Although these are known to correspond to Lie systems (with semisimple Vessiot-Guldberg algebra) and hence admit a superposition rule (see ^[27,53] and references therein), in general, they do not admit an LH structure; thus their integrability by quadratures must be analyzed by other means. Actually, among the 28 isomorphism classes of Lie algebras of vector fields over $\mathbb{R}^2$, only 12 of them correspond to algebras of Hamiltonian vector fields ^[25], and further applications to models in biology or information theory should be searched among these LH algebras. Concerning the integrability by quadratures of LH systems, it can be inferred from the general theory of Lie systems for the case of solvable Vessiot-Guldberg algebras ^[45,48]. For a non-solvable Vessiot-Guldberg algebra this property does not hold in general, although for certain Lie systems additional integrability conditions, appropriately combined with a reduction procedure ^[54], may lead to an integration by quadratures. It is currently an open problem whether, for the case of non-solvable Vessiot-Guldberg algebras, the compatibility with a symplectic structure, beyond providing a systematic procedure to construct superposition rules, also determines additional integrability conditions that allow explicit integration of a system.

Although the extension of SIS models described by Eq (2.6) as a particular case has been exhausted in terms of the LH formalism, there is an additional possibility of further extending the systems by applying the so-called Poisson-Hopf algebra deformation of LH system formalism recently introduced in ^[33,55,56], which is based on quantum groups ^[32,57]. The main idea behind this procedure is to combine LH systems with quantum deformations, leading to genuinely nonlinear systems that are no more described by a finite-dimensional LH algebra, but in terms of involutive functional distributions. As shown in ^[55,56], for such systems, a kind of deformed superposition rule and deformed systems of differential equations can also be found, with the additional feature that the quantum deformation parameter can be identified with a perturbation parameter. In this frame, it is conceivable to analyze the generic deformed systems resulting from the embedding chains given by Eq (5.12), i.e., $\mathfrak{b}_2\subset \mathfrak{h}_4\subset \mathfrak{h}_6$ and $\mathfrak{sl}(2) \subset \mathfrak{h}_6$ (a special case with a quantum deformation of $\mathfrak{sl}(2)$ has already been considered in ^[13]), and to study to which extent the corresponding equations can be explicitly solved. In this respect, we recall that all possible quantum deformations for the oscillator Lie algebra $\mathfrak{h}_4$ can be found in ^[58], while those corresponding to the two-photon $\mathfrak{h}_6$ are presented in ^[50] in a Schrödinger Lie algebra basis. Work along the various research lines commented above is currently in progress.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

R. C. Stursberg and F. J. Herranz have been partially supported by the Agencia Estatal de Investigación (Spain) under grant PID2019-106802GB-I00/AEI/10.13039/501100011033. F.J.H. also acknowledges support from the Regional Government of Castilla y León (Junta de Castilla y León, Spain) and the Spanish Ministry of Science and Innovation MICIN and the European Union NextGenerationEU (PRTR C17.I1). The authors acknowledge the referees for their valuable comments and suggestions that have greatly improved the presentation, as well as for pointing out additional pertinent references.

Conflict of interest

The authors declare no conflicts of interest.

Appendix

Basic properties of Lie-Hamilton systems

We briefly recall the main definitions and properties of LH systems. For additional details, the reader is referred to ^{[23,24,25,26,27]} and the references therein.

Given a system of first-order ordinary differential equations on a manifold $M$ of dimension $n$ with coordinates ${\bf x} = \{x_1, \dots, x_n \}$,

we say that it admits a fundamental system of solutions if the general solution can be written in terms of $m\leq n$ functionally independent particular solutions $\{ {\bf x}_{1}, \dots, {\bf x}_{m} \}$, where ${\bf x}_{s} = \bigl((x_{1})_{s}, \dots, (x_{n})_{s}\bigr)$, and there are $n$ constants $\left\{ k_{1}, \dots, k_{n}\right\}$:

The relation (A.2) is called a superposition rule for the system (A.1), while the set $\{ {\bf x}_{1}, \dots, {\bf x}_{m} \}$ provides a fundamental system of solutions. It is straightforward to verify that this system can be described in equivalent form by means of the vector field

called the $t$-dependent vector field associated with the system (A.1). The existence of fundamental systems of solutions was first analyzed by Lie ^[44], who established a characterization in terms of finite-dimensional Lie algebras. The Lie-Scheffers theorem asserts that the system (A.1) admits a fundamental system of solutions if and only if it can be represented in terms of $\ell$ $t$-dependent parameters $b_i(t)$ in the form

and such that the vector fields

span an $\ell$-dimensional Lie algebra $\mathfrak{g}$, where the numerical constraint $n\, m\geq \ell = \dim\mathfrak{g}$ is satisfied. Thus, the generators given by Eq (A.5) obey the generic commutation relations given by

Under these conditions, the $t$-dependent vector field described by Eq (A.3) is reformulated as

The Lie algebra $\mathfrak{g}$ generated by the vector fields ${\bf X}_i ({\bf x})$ is usually called a Vessiot-Guldberg Lie algebra of the system (A.1), while either the system itself or ${\bf X} (t, {\bf x})$ is referred to as a Lie system ^[27,28].^†

^† It should be observed that Vessiot-Guldberg algebras are not uniquely determined, and hence do not constitute an invariant of the system. It has thus sense to speak about minimality, in the sense that $\mathfrak{g}$ is minimal if no proper subalgebra of $\mathfrak{g}$ is a Vessiot-Guldberg-Lie algebra of the system (A.1).

A Lie system ${\bf X} (t, {\bf x})$ is called an LH system if it admits a Vessiot-Guldberg Lie algebra $\mathfrak{g}$ of Hamiltonian vector fields with respect to a Poisson structure ^[27], with the compatibility condition of the generators ${\bf X}_i ({\bf x})$ and the symplectic form $\omega$ being given by the Lie derivative:

The corresponding Hamiltonian functions $h_i ({\bf x})$ associated with ${\bf X}_i({\bf x})$ are obtained through the inner product (see ^[31]):

The Hamiltonian functions $h_i({\bf x})$ span (eventually adjoining a constant function $h_0$) an $\ell$-dimensional Lie algebra ${\mathcal H}_\omega$ called an LH algebra, with Poisson brackets given by

where $C_{ab}^{c}$ are the same structure constants written in Eq (A.6). Therefore, $\mathcal{H}_{\omega}$ is spanned by a set of functions $\left\{h_1, \ldots, h_\ell\right\}\subset C^\infty({{\cal M}})$, where ${{\cal M}}$ is a suitable submanifold of $M$ that ensures that each $h_i$ is well defined. The $t$-dependent Hamiltonian function $h (t, {\bf x})$ associated with the $t$-dependent vector field described by Eq (A.7) leading to the same Lie system (A.1) is given by

As LH systems correspond to a particular class of Lie systems, a superposition rule (A.2) is guaranteed. In general, standard methods ^[28,29] for deriving superposition rules for Lie systems require the integration of PDEs or ODEs, which is often a cumbersome task. Nevertheless, if the Lie system is also an LH one, there is an algebraic approach ^[24] (see also ^[26,33]) that allows a systematic obtainment of $t$-independent constants of the motion, from which a superposition rule can be deduced. Such a formalism is based on the so-called coalgebra symmetry approach for superintegrable systems (see ^[32] and references theirein for technical details). The only formal requirement is that the LH algebra possesses a non-trivial Casimir invariant.

Explicitly, let us consider the LH algebra ${\cal H}_\omega$ of an LH system ${\bf X}$ with Hamiltonian (A.11), expressed as a Lie-Poisson algebra with generators $\{ v_1, \dots, v_\ell \}$ that formally satisfy the same Poisson brackets given by Eq (A.10) and assume that it admits the non-constant Casimir function

The underlying coalgebra structure of the LH algebra leads to the definition of the following Hamiltonian vector fields:

where ${\bf x}_{s} = \bigl((x_{1})_{s}, \dots, (x_{n})_{s}\bigr)$ denotes the coordinates in the $s^{{\rm th}}$-copy submanifold ${{\cal M}}\subset M$ within the product manifold $\mathcal{M}^k$. It can be proved (see ^[24,26,33]) that the functions defined through the Casimir (A.12) and the Hamiltonian vector fields given by Eq (A.13) as

are $t$-independent constants of the motion for the diagonal prolongation $\widetilde { {\bf X}}^{m+1}$ of the LH system ${\bf X}$ to the product manifold ${{\cal M}}^{m+1}$, with $\widetilde { {\bf X}}^{m+1}$ given by

Observe that $F^{(k)}$ denote functions of $C^\infty({{\cal M}}^{m+1})$, and that they can be considered as $t$-independent constants of the motion for the LH system (A.11). Furthermore, if every $F^{(k)}$ is non-constant, they constitute a set of $m$ functionally independent functions in $C^\infty({{\cal M}}^{m+1})$ that are in involution (i.e., they Poisson-commute). Moreover, the functions $F^{(k)}$ can be used to generate other constants of the motion in the following form:

where $S_{ij}$ denotes the permutation of the variables ${\bf x}_i\leftrightarrow {\bf x}_j$ in ${{\cal M}}^{m+1}$.

We finally recall that, in order to obtain a superposition rule (A.2) that depends on $m$ particular solutions, one should find a set $\{I_1, \ldots, I_n\}$ of $t$-independent constants of the motion on $\mathcal{M}^{m+1}$ for $\widetilde{ {\bf X}}^{m+1}$ such that ^[27]

This allows us to express the coordinates ${\bf x}_{m+1} = \left\{ (x_1)_{m+1}, \dots, (x_n)_{m+1}\right\}$ in terms of the remaining coordinates in $\mathcal{M}^{m+1}$ and the constants $k_1, \ldots, k_n$ defined by the conditions $I_1 = k_1, \ldots, I_n = k_n$. The constants of the motion described by Eqs (A.14) and (A.16) are generally sufficient to determine the set $\{I_1, \ldots, I_n\}$, and hence to deduce a superposition rule for the LH system in a direct algebraic way.