Epidemic models try to predict the spread of an infectious disease afflicting a specific population, see for example [14, 39, 42] . These models are rooted in the works of Bernoulli in the 18th century, when he proposed a mathematical model to defend the practice of inoculating against smallpox [27] . This was the start of germ theory.
At the beginning of the 20th century, the emergence of compartmental models was starting to develop. Compartmental models are deterministic models in which the population is divided into compartments, each representing a specific stage of the epidemic. For example, S represents the susceptible individuals to the disease, I designates the infected individuals, whilst R stands for the recovered ones. The evolution of these variables in time is represented by a system of ordinary differential equations whose independent variable, the time, is denoted by t. Some of these first models are the Kermack-McKendrick [28] and the Reed-Frost [1] epidemic models, both describing the dynamics of healthy and infected individuals among other possibilities. There are several types of compartmental models [25, 26, 40] , as it can be the SIS model, in which after the infection the individuals do not adquire immunity, the SIR model, in which after the infection the individuals adquire immunity, the SIRS model, for which immunity only lasts for a short period of time, the MSIR model, in which infants are born with immunity, etc. In this present work our focus is on the SIS model.
The SIS model. The susceptible-infectious-susceptible (SIS) epidemic model assumes a population of size N and one single disease disseminating. The infectious period extends throughout the whole course of the disease until the recovery of the patient with two possible states, either infected or susceptible. This implies that there is no immunization in this model. In this approach the only relevant variable is the instantaneous density of infected individuals ρ = ρ(τ ) depending on the time parameter τ , and taking values in [0, 1] . The density of infected individuals decreases with rate γρ, where γ is the recovery rate, and the rate of growth of new infections is proportional to αρ (1 − ρ) , where the intensity of contagion is given by the transmission rate α. These two processes are modelled through the compartmental equation
One can redefine the timescale as t := ατ and introduce the constant ρ 0 := 1 − γ/α, so we can rewrite (1.1) as
Clearly, the equilibrium density is reached if ρ = 0 or ρ = ρ 0 . Although compartmental equations have proven their efficiency for centuries, they are still based on strong hypotheses. For example, the SIS model works more efficiently under the random mixing and large population assumptions. The first assumption is asking homogeneous mixing of the population, that is, individuals contact with each other randomly and do not gather in smaller groups, as abstaining themselves from certain communities. This assumption is nevertheless rarely justified. The second assumption is the rectangular and stationary age distribution, which means that everyone in the population lives to an age L, and for each age up to L, which is the oldest age, there is the same number of people in each subpopulation. This assumption seems feasible in developed contries where there exists very low infant mortality, for example, and a long live expectancy. Nonetheless, it looks reasonable to implement probability at some point to permit random variation in one or more inputs over time. Some recent experiments provide evidence that temporal fluctuations can drastically alter the prevalence of pathogens and spatial heterogeneity also introduces an extra layer of complexity as it can delay the pathogen transmission [21, 49] .
It is needless to point out that fluctuations should be considered in order to capture the spread of infectious diseases more closely. Nonetheless, the introduction of these fluctuations is not trivial. One way to account for fluctuations is to consider stochastic variables. On the other hand, it seems that in the case of SIS models there exist improved differential equations for the mean and variance of infected individuals. Recently, in [43] , the model assumes the spreading of the disease as a Markov chain in discrete time in which at most one single recovery or transmission occurs in the duration of this infinitesimal interval. As a result [29] , the first two equations for instantaneous mean density of infected people ρ and the variance σ 2 = ρ 2 − ρ 2 are
where ∆ 3 = ρ 3 − ρ 3 . This system finds excellent agreement with empirical data [43] . Equations (1.2) and (1.3) are equivalent when σ becomes irrelevant compared to ρ . Therefore, a generalization of compartmental equations only requires mean and variance, neglecting higher statistical moments. The skewness coefficient vanishes as a direct consequence of this assumption, so that ∆ 3 := 3σ 2 ρ . For a big number of individuals (N 1), the resulting equations are
(1.4)
The system right above can be obtained from a stochastic expansion as it is given in [54] , as well.
Hamiltonian character of the model (1.4) . The investigation of the geometric and/or the algebraic foundations of a system permits to employ several powerful techniques of geometry and algebra while performing the qualitative analysis of the system [3, 5, 33] . This even results in an analytical/general solution of the system in our case. For example, we cite [32] and [44] for Lie symmetry approach to solve the classical SIS model. Particularly, the Hamiltonian analysis of a system plays an important role in the geometrical analysis of a given system. For instance, in compartmental theories we can cite an early work [45] that classifies the SIR model as bihamiltonian. We also refer to [23] for conformal Hamiltonian analysis of the Kermack-McKendrick Model. In a very recent study [7] , it is shown that all classical epidemic compartmental models admit Hamiltonian realizations.
In [43] , the SIS system (1.4) involving fluctuations has been recasted in Hamiltonian form in the following way: the dependent variables are the mean ρ and the variance σ 2 , and they both depend on time. Then, we define the dynamical variables q = ρ and p = 1/σ, so the system (1.4) turns out to be
We employ the abbreviation SISf for system (1.5) to differentiate it from the classical SIS model in (1.2) . The letter "f" accounts for "fluctuations".
We have computed the general solution to this system, finding a more general solution than the one provided by Nakamura and Martínez in [43] . Indeed, we have found obstructions in their model solution. We shall comment this in the last section gathering all our new results.
Our general solution for this system reads:
,
(1. 6) In order to develop a geometric theory for this system of differential equations, we need to choose certain particular solutions that we shall make use of. Here we present three different choices and their corresponding graphs according to the change of variables q = ρ and p = 1/σ. Let us turn now to interpret these equations geometrically on a symplectic manifold. The symplectic two-form ω = dq ∧ dq is a canonical skew-symmetric tensorial object in twodimensions. For a chosen (real-valued) Hamiltonian function h = h(q, p), the dynamics is governed by a Hamiltonian vector field X h defined through the Hamilton equation
where ι X h is the contraction operator in the tensor algebra, and dH is the (exterior) derivative of h. In terms of the coordinates (q, p), the Hamilton equations (1.7) become
It is possible to realize that the SISf system (1.5) is a Hamiltonian system since it fulfills the Hamilton equations (1.7). To see this, consider the Hamiltonian function
and substitute it into (1.8). A direct calculation will lead us to (1.5) . The skew-symmetry of the symplectic two-form implies that the Hamiltonian function is constant all along the motion. In classical mechanics, where the Hamiltonian is taken to be the total energy, this corresponds to the conservation of energy.
The goal of the present work. Our focus is the SISf model permitting fluctuations given in (1.5) . We shall perform both the Lie analysis and the Lie-Hamilton analysis of this system in order to explore the geometric/algebraic foundations of the model as well as to arrive at the analytical general solution. Further, we shall introduce a deformation (the so-called quantum deformation) of the model (1.5) to make it applicable to more complicated issues.
To get involved in the mathematical theories needs some special emphasis or/and some expertise. For that, we shall try to present both the theoretical and the practical results in the most accessible forms for a general audience as far as we can. Our plan to achieve this is to exhibit itinerary maps, graphs, and a comprehensive appendix.
The paper is organized into three main sections, and an appendix. Sections 2 and 3 concern the Lie and the Lie-Hamilton analysis of the SISf model (1.5), respectively. Section 4 proposes the so-called quantum deformation of the model and Appendix A is a brief review of some algebraic constructions, such as Lie algebras, Poisson algebras, Lie coalgebras and Poisson Hopf algebras, that we refer to in the main body of the paper.
A (nonautonomous) system of ordinary differential equations (ODEs) can be written as a timedependent vector field and vice versa. Such system is called Lie system if the associated timedependent vector field takes values in a finite-dimensional Lie algebra of vector fields [15, 16] . Equivalently, one can also define a Lie system as a system of ODEs admitting a(nonlinear) superposition principle, that is, a map allowing us to express the general solution of the system of ODEs in terms of a family of particular solutions and a set of constants related to initial conditions. Accordingly, in this section, we shall show that the SISf model (1.5) is a Lie system. For the completeness of this work, we start with a brief summary of some basic notions concerning the theory of Lie systems. We refer to Appendix A.1 for the definitions of Lie (sub)algebras, and to Appendix A.2 for the definitions of (co)tangent bundle and vector fields. For some further reading on Lie systems especially about its relevant role in Physics, Mathematics, Biology, Economics, and other fields of research, see [18, 35, 51] and extensive reference lists there in.
Let N be a manifold. A time-dependent vector field on N is a differentiable mapping
where pr 2 is the projection to the second factor in R × N . In this regard, we can consider a time-dependent vector as a set (family) of standard vector fields {X t } t∈R , depending on a single parameter. The converse of this assertion is also true, that is a set of vector fields depending smoothly on a single parameter can be written as a time-dependent vector field [18] . We plot the following diagram to summarize the discussions.
:
A time dependent vector field determines a non-autonomous ODE system
The suspension of a time-dependent vector field X on N is the vector field on R × N defined to be X + ∂/∂t, where t is the global coordinate for R. An integral curve of a time-dependent vector field X is an integral curve of the suspension [3] . Conversely, every system of first-order differential equations in normal form describes the integral curves of a unique time-dependent vector field.
Consider the case that N is a two dimensional manifold with local coordinates (q, p), and the extended space R × N with induced coordinates (t, q, p). A time dependent vector field is written as
where k and g are real valued functions on R × N . In this coordinate system, the dynamical equations governed by X determine a system of nonautonomous ordinary differential equations (2.5)q = k(t, q, p),ṗ = g(t, q, p).
Notice that, the SISf model (1.5) fits this picture with even the generalization ρ 0 = ρ(t). We shall turn to this discussion in the following subsections.
Vessiot-Guldberg algebra. Start with a time-dependent vector field X, and determine the set {X t } t∈R of vector fields. We denote the smallest Lie subalgebra containing the set {X t } t∈R by V X . The associated distribution D X for a time-dependent vector field X is the generalized distribution spanned by the vector fields in V X that is
The associated codistribution (D X ) • for a time-dependent vector field X is then defined to be the annihilator of D X that is
is a subbundle of the cotangent bundle T * N . Our interest relies in a distribution D X determined by a finite-dimensional V X and hence D X becomes integrable on the whole N . In this case, V X is called Vessiot-Guldberg Lie algebra. It is worth noting that even in this case, (D X ) • does not need to be a differentiable codistribution.
A function f on an open neighbourhood U is a local t-independent constant of motion for a system X if and only if df (x) in (D X x ) • | U for all x in N . This statement reads that (locally defined) time-independent constants of motion of time-dependent vector fields are determined by (locally defined) exact one-forms taking values in the associated codistribution. Therefore, (D X ) • is a crucial object in the calculation of such constants of motion for a system X.
Let X be a time-dependent vector field defined on a manifold N . A superposition rule depending on m particular solutions of X is a function
such that the general solution x(t) of X can be written as x(t) = φ(x (1) (t), . . . , x (m) (t), λ) for any generic family x (1) (t), . . . , x (m) (t) of particular solutions. Here, λ is a point of N related with the initial conditions. A system of equations admitting superposition rule is called Lie system.
A modern statement of this result is described in [16] . A first-order system (2.3) admits a superposition rule if and only if the associated time-dependent vector field X can be written as
for a certain family b 1 (t), . . . , b r (t) of time-dependent functions and a family X 1 , . . . , X r of vector fields on N spanning an r-dimensional real Lie subalgebra of vector fields. We refer to [35] for details, see [58, 50] for some further discussions and the first examples. The Lie-Scheffers Theorem yields that a system X admits a superposition rule if and only if V X is finite-dimensional, [18] .
An itinerary map to derive superposition rules. General solutions of Lie systems can also be investigated through superposition rules. There exist various procedures to derive them [4, 57] , but we hereafter use the method devised in [16] , which is based on the notion of diagonal prolongation [18] . Let us denote the m + 1-times Cartesian product of N by itself as N (m+1) . We denote by
the projection onto the first factor. Given a time-dependent vector field X on N , the diagonal prolongation X of X to the product space N (m+1) is a unique time-dependent vector field on N (m+1) projecting to X by the projection pr that is pr * X t = X t for all t. We also require X to be invariant under the permutations x (i) ↔ x (j) with i, j = 0, . . . , m. The procedure to determine superposition rules described in [16] goes as follows:
Step 1. Take a basis X 1 , . . . , X r of a Vessiot-Guldberg Lie algebra associated with the Lie system X.
Step 2. Choose the minimum integer m so that the diagonal prolongations X 1 , . . . , X r , of X 1 , . . . , X r to N m are linearly independent at a generic point.
Step 3. Obtain, for instance, by the method of characteristics, n functionally independent first-integrals F 1 , . . . , F n common to all the diagonal prolongations, X 1 , . . . , X r , to N (m+1) .
We require such functions to satisfy
Assume that these integrals take certain constant values, i.e., F i = k i where i = 1, . . . , n, and employ these equalities to express the variables (x 1 ) (0) , . . . , (x n ) (0) in terms of the variables of the other copies of N within N (m+1) and the constants k 1 , . . . , k n .
The obtained expressions constitute a superposition rule in terms of any generic family of m particular solutions and n constants. Let us apply these steps in the realm of the SISf model.
The model (1.5) can be generalized to a model represented by a time-dependent vector field
where the constitutive vector fields are computed to be
The generalization comes from the fact that ρ 0 (t) is no longer a constant, but it can evolve in time. Let us apply the steps introduced in the previous section one by one to arrive at the general solution.
Step 1. For the vector fields in (2.12), a direct calculation shows that the Lie bracket
is closed within the Lie algebra. This implies that the SISf model (1.5) is a Lie system. The Vessiot-Guldberg algebra spanned by X 1 , X 2 is an imprimitive Lie algebra of type I 14 according to the classification presented in [8] .
Step 2. If we copy the configuration space twice, we will have four degrees of freedom (q 1 , p 1 , q 2 , p 2 ) and we will archieve precisely two first-integrals in vinicity of the Fröbenius theorem. A first-integral for X t has to be a first-integral for X 1 and X 2 simultaneously. We define the diagonal prolongation X 1 of the vector field X 1 in the decomposition (2.13). Then we look for a first integral F 1 such that X 1 [F 1 ] vanishes identically. Notice that if F 1 is a first-integral of the vector field X 1 then it is a first integral of X 2 due to the commutation relation. For this reason, we start by integrating the prolonged vector field (2.14)
through the following characteristic system
Fix the dependent variable q 1 and obtain a new set of dependent variables, say (K 1 , K 2 , K 3 ), which are computed to be
Step 3. This induces the following basis in the tangent space
provided that q 1 is not zero. Introducing the coordinate changes exhibited in (2.16) into the diagonal projection X 2 of the vector field X 2 , we arrive at the following expression
To integrate the system once more, we use the method of characteristics again and obtain
We obtain two first integrals by integrating in pairs (K 1 , K 2 ) and (K 1 , K 3 ), where we have fixed K 1 . After some cumbersome calculations we obtain
By substituting back the coordinate transformation (2.16) into the solution (2.19) (please notice the difference between capitalized constants (K 1 , K 2 , K 3 ) and lower case constants (k 1 , k 2 ), we arrive at the following implicit equations
.
(2.20)
Let us notice that the equations (2.20) depend on a particular solution (q 2 , p 2 ) and two constants of integration (k 1 , k 2 ) which are related to initial conditions.
Let us show now the graphs and values of the initial conditions for which the solution reminds us of sigmoid behavior, which is the expected growth of ρ(t). As particular solution for (q 2 , p 2 ), we have made use of particular solution 2 given in Figure 2 through its corresponding values of q, p through the change of variables q = ρ and p = 1/σ. Notice that we have not included the error graph since it gives a constant zero graph because k 1 → 0
Step 3 revisited -linear approximation. Since the solution (2.20) is quite complicated, one may look for a solution of a linearized model. We first employ the following change of coordinates
In terms of these new variables, the system (2.18) reads
.
One can solve the system above by introducing a linear approximation
after which (2.22) reads
We can solve now v and w in terms of u and obtain
According to the algorithm for deriving superposition rule, we need to isolate the constants of integration k 1 and k 2 . Hence, the two first integrals read now
Now, if we substitute the coordinate changes in (2.21) and in (2.16) , we arrive at the following general solution
which can be written as
Notice that the solution depends on a particular solution (q 2 , p 2 ) and two constants of integration (k 1 , k 2 ), as in (2.20).
Let us show now the graphs and values of the initial conditions for which the solution reminds us of sigmoid behavior, which is the expected growth of ρ(t). As particular solution for (q 2 , p 2 ), we have made use of particular solution 2 given in Figure 2 through its corresponding values of q, p through the change of variables q =< ρ > and p = 1/σ. In this section, we shall show that the SISf model (1.5) is a Lie-Hamilton system [10, 17, 35] . Among the developed methods for Lie-Hamilton systems, we consider a very important recent method for the obtainance of solutions as superposition principles through the Poisson coalgebra method [8, 9] . The traditional method for the computation of superposition principles for Lie systems relies in the integration of systems of ordinary or partial differential equations [4, 57] , but in the case of Lie-Hamilton systems, the nonlinear superposition rule can be obtained straightforwardly through a Casimir of the Vessiot-Guldberg Lie algebra. We refer to Appendices A.3 and A.4 for informal summaries of the symplectic and the Poisson geometries, and to Appendix A.5 for Lie coalgebras. We refer Appendix A.6 for the definition of the symmetric algebra.
A Lie-Hamiltonian structure is a triple (N, Λ, h), where (N, Λ) stands for a Poisson manifold and h represents a t-parametrized family of functions h t : N → R such that the Lie algebra H := Lie({h t } t∈R ) generated by this family is finite-dimensional Poisson algebra. A timedependent vector field X is said to possess a Lie-Hamiltonian structure (N, Λ, h) if X t is the Hamiltonian vector field corresponding to h t , for each t ∈ R,
In this case, X is called a Lie-Hamilton system, and H is called a Lie-Hamilton algebra for X. Now we examine superposition rules for the case of Lie-Hamilton systems.
Itinerary map for superposition rules of the Lie-Hamilton Systems. Let X be a Lie-Hamilton system with a Poisson algebra H spanned by the linearly independent Hamiltonian functions {h 1 , . . . , h r }.
Step 1. Consider the injection g → H of a Lie algebra g into H that turns the basis v i of g into the Hamiltonian functions φ(v i ) := h i for i = 1, . . . , r. Referring to this inclusion, we can define Poisson algebra homomorphisms
where S(g) is the symmetric algebra.
Step 2. If C is a polynomial Casimir of the Poisson algebra S(g), say C = C(v 1 , . . . , v r ), then D(C) is a constant of motion for X, and so are the functions defined by . . . , v r ) ) , k = 2, . . . , m.
Let us notice, that each F (k) can naturally be considered as a function of C ∞ (M m ) for every m k.
Step 3. From the functions F (k) , we can obtain other constants of motion in the form
Indeed, since the prolongation X is invariant under the permutations
ij are also t-independent constants of motion for the diagonal prolongations X to M m . Let us now illustrate this procedure in the following example.
Over the plane, consider the following vector fields
Step 1. With respect to the canonical symplectic structure ω = dx ∧ dy, this corresponds Lie algebra, denoted by iso (2), determined by a basis
Step 2. The symmetric Poisson algebra S iso(2) has a non-trivial Casimir invariant given by
. Choosing the representation given in (3.6), we obtain a trivial constant of motion on (x, y) ≡ (x 1 , y 1 )
Introducing the coalgebra structure in S iso(2) through the coproduct (A. 34) , we obtain nontrivial first integrals.
(3.8)
Step 3. Furthermore, using the property of permutating subindices (3.4) , we find more first integrals
12 +F
13 +F
23 . We may consider as many first integrals as the number of degrees of freedom of the system in order to integrate it. The combination of these functions leads us to a superposition rule.
The SISf epidemic system in (1.5) admits a Hamiltonian formulation as in (1.9), where the corresponding Hamilton equations are (1.5) . In this section we will retrieve the Hamiltonian (1.9) using the theory of Lie systems, in particular, the theory of Lie-Hamilton systems. To retrieve the Hamiltonian, first we need to prove that equations in (1.5) form a Lie-Hamilton system.
We have already proven in (2.13) that (1.5) defines a Lie system. In order to see if it is a Lie-Hamilton system, we first need to check whether the vector fields in (2.13) are Hamiltonian vector fields. Consider now the canonical symplectic form ω = dq ∧ dp. It is easy to check that the vector fields X 1 and X 2 in (2.13) are Hamiltonian with respect to the Hamiltonian functions (3.10) R R, and it is isomorphic to the one defined by vector fields X 1 , X 2 . The Hamiltonian function for the total system is
and it is exactly the Hamiltonian function (1.9) proposed in [43] .
Lie-Hamilton systems can also be integrated in terms of a superposition rule, as it was explained in the preliminary section. In order to do that, we need to find a Casimir function for the Poisson algebra, but unfortunately, there exists no nontrivial Casimir in this particular case. It is interesting to see how a symmetry of the Lie algebra {X 1 , X 2 } commutes with the Lie bracket, i.e. the vector field
Notice too that Z is a conformal vector field, that is,
Since it is a Hamiltonian system, one would expect that a first integral for Z, let us say f , would Poisson commute with the Poisson algebra {h 1 , h 2 }, since Z = −Λ(df ). Nonetheless, this is not the case unless f = constant. This implies that the Casimir is a constant, hence trivial and the coalgebra method can not be directly applied. However, there is a way in which we can circumvent this problem by considering an inclusion of the algebra I r=1 14A as a Lie subalgebra of a Lie algebra to another class admitting a Lie-Hamiltonian algebra with a non-trivial Casimir. In this case, we will consider the algebra, denoted by I 8 iso(1, 1), due to the simple form of its Casimir. If we obtain the superposition rule for I 8 , we simultaneously obtain the superposition for I r=1 14A as a byproduct.
The Lie-Hamilton algebra iso(1, 1) has the commutation relations (3.14)
The Casimir associated to this Lie-Hamilton algebra is C = h 1 h 2 +h 3 h 0 . Let us apply the coalgebra method to this case. Mapping the representation without coproduct, the first iteration is trivial, i.e., F = 0. We could usethe second order coproduct and third order coproduct ∆ (2) and ∆ (3) , or the second order coproduct ∆ (2) together with the permuting subindices property. We need three constants of motion, this would be equivalent to integrating the diagonal prolongation X on (R 2 ) 3 . Using the coalgebra method and subindex permutation, one obtains
From them, we can choose two functionally independent constants of motion. Our choice is
23 = k 2 . The introduction of k 3 simplifies the final result, with expression
In the case that matters to us, I r=1 14A , the third constant k 3 is a function k 3 = k 3 (x 2 , y 2 , x 3 , y 3 ) and B 0. Notice though that this superposition rule is expressed in the basis (3.14), therefore, we need the change of coordinates between the present iso(1, 1) and our problem (3.10). See that the commutation relation {h 1 , h 3 } = −h 1 in (3.14) coincides with the commutation relation {h 1 , h 2 } = h 2 of our pandemic system (3.10). So, by comparison, we see there is a change of coordinates (3.18) x = −qp,
This way, introducing this change (3.18) in (3.16), the superposition principle for our Hamiltonian pandemic system reads
Here, (q 2 , p 2 ) and (q 3 , p 3 ) are two particular solutions and k 1 , k 2 , k 3 are constants of integration. Now, we show the graphics for < ρ >= q(t) and σ 2 = 1/p 2 using the two particular solutions in Figure 2 and Figure 3 provided in the introduction. Notice that we have renamed c = (k 2 − k 1 ± B) and k = (k 2 − k 1 ∓ B). 
A Lie-Hamilton system admits a Poisson-Hopf deformation, see [11] . The interest of Poisson-Hopf deformations of these models resides in the fact that many of the outcome deformed systems happen to be other systems that are already identified in the physics and mathematics literature. This is a highlight, since this basically implies that two different systems scattered in the mathematical physics literature can be merely related by a Poisson-Hopf transformation. This could be a start for the classification of the fuzzy number of integrable systems.
We have interpreted that the SISf model (1.5) is a Lie-Hamilton system. So, in this section we propose a (quantum) deformation of the SISf model [8, 12, 13, 19, 35] which does not rely on the complicated functional analysis considerations for unbounded quantum operators [30] . The introduction of the quantum approach will account for the existence of a possible interaction of the SISf model with a heat bath, which is effectively an infinite pool of thermal energy at a given constant temperature. This could correspond with centrally heated buildings. We refer Appendix A.6 for some details on Poisson-Hopf algebras.
To obtain a deformation of the Lie-Hamilton realization of the SISf model we make use of deformed Poisson-Hopf algebras. Following [11] , we summarize the procedure for the planar systems as follows:
Itinerary map to arrive at a deformation of a Lie-Hamilton system. Consider a Lie-Hamilton system given by a time-dependent vector field X t ∈ R 2 describing a Lie system.
There is a corresponding Lie algebra of Hamiltonian functions Lie({h t } t∈R ) satisfying (3.1) that defines the Lie system as a Lie-Hamilton system.
Step 1. Deform the space C ∞ ({h tz } * t∈R ) with a deformation parameter z ∈ R and accordingly the Poisson bracket (h z,1 , . . . , h z,l ) .
Here, F z,ij are certain smooth functions also depending smoothly on the deformation parameter z and such that
where ∇ is the gradient.
Step 2. Deform the vector fields X z,i as
Step 3. The deformed total dynamics is encoded in
Notice that when the parameter goes to the zero, we have
and the deformed dynamics reduces to the initial one.
For the SISf model (1.5), we start with the Vessiot-Guldberg algebra (2.13) labelled as I r=1 14A . To obtain a deformation of a Lie algebra I r=1 14A , we need to rely on a bigger Lie algebra, in this case, we make use of sl (2) . To this end, consider the vector fields X 1 and X 2 in (2.12), and let X 3 be a vector field given by (4.6)
where c ∈ R. Then, {X 1 , X 2 , X 3 } span a Vessiot-Lie algebra V isomorphic to sl(2) that satisfies the following commutation relations
This vector field X 3 admits a Hamiltonian function, say h 3 , with respect to the canonical symplectic form on R 2 , so that we have the family
Hence, {h 1 , h 2 , h 3 } span a Lie-Hamilton algebra H ω ; isomorphic to sl(2) where the commutation relations with respect to the Poisson bracket induced by the canonical symplectic form ω on R 2 are given by
Here is the deformation of the model by steps.
Step 1. Applying the non-standard deformation of sl(2) in [11] we arrive at the Hamiltonian functions
Here, shc(x) is the cardinal hyperbolic sinus function. Accordingly, the Poisson brackets are computed to be (4.11)
Step 2. The vector fields X z;1 and X z;2 associated to the Hamiltonian functions h z;1 and h z;2 exhibited in (4.10) are
(4.12)
We do not write explicitly the expression of the vector field X z;3 because it does not play a relevant role in our system. The deformed vector fields keep the commutation relations
Step 3. The total Hamiltonian function for the deformed model is
so that the deformed dynamics is computed to be
(4.15)
This system describes a family of z-parametric differential equations that generalizes the SISf model (1.5) , where the demographic interaction and both rates allow a more realistic representation of the epidemic evolution. According to the kind of deformation, this may be called a quantum family SISf model. Note that the SISf model can be recovered in the limit when z tends to zero.
For the present case, the constants of motion defined in (3.3) are computed to be (4.16)
after the quantization, the latter one becomes
where h (2) z;j := D (2) z (∆ z (v j )). This coproduct ∆ z can be described as a follows
More explictly, using the expressions given in (4.10), we have h (2) z;j =h z;j (q 1 , p 1 )e 2zh z;2 (q 2 ,p 2 ) + h z;j (q 2 , p 2 )e −2zh z;2 (q 1 ,p 1 ) , j = 1, 3
z;2 =h z;2 (q 1 , p 1 ) + h z;2 (q 2 , p 2 ).
(4.18)
So, to retrieve another constant of motion we can apply the trick of permuting indices. Then, here we have a second constant of motion, writing it implicitly,
z, (23) = shc 2zh
where the subindex (23) means that the variables (q 2 , p 2 ) are interchanged with (q 3 , p 3 ) when they appear in the deformed Hamiltonian functions h z;j and h (2) z;j(23) =h z;j (q 1 , p 1 )e 2zh z;2 (q 3 ,
z;2(23) =h z;2 (q 1 , p 1 ) + h z;2 (q 3 , p 3 ).
In (4.17), we have = k 12 , with k 23 , k 12 ∈ R, one is able to retrieve a superposition principle for q 1 = q 1 (q 2 , q 3 , p 2 , p 3 , k 12 , k 23 ) and p 1 = p 1 (q 2 , q 3 , p 2 , p 3 , k 12 , k 23 ). Notice that here (q 2 , p 2 ) and (q 3 , p 3 ) are two pairs of particular solutions and k 12 , k 23 are two constants over the plane to be related to initial conditions.
Here we present a summary of our new results, we conclude how they agree with the experimental data and how they clash with other models and results by other authors. Let us list our major results.
• We have achieved a more general solution (1.6) for the SISf system, a SIS model permitting fluctuations, presented in (1.5) than the one provided in [43] . • We have limited the range of C 1 and C 2 in the solution (1.6) for system (1.5), so that we obtain sigmoid-type or hyperbolic-type like behavior, which is the one that we could expect for the growth of infected individuals ρ(t). This range is 3C 2 C 1 10 14 C 2 . For simplification, we choose C 2 = 1. If one exceeds these limits, it is easy to spot singularities and unexpected behavior of the function. This is precisely what we have realized as we plotted the solutions given in [43] . We are still unsure whether their solution is an actual solution or if their choice of initial conditions are not the ideal to overcome singularities and find the real regime. • The error σ 2 decays with time, as the density of infected individuals stabilizes as constant. This is what one expects intuitively.
• We have been able to retrieve general solutions for equations (1.5) making use of geometric methods. The theory of Lie systems has permitted us to find a general solution in terms of a particular solution and a challenging choice of the constants of integration (k 1 , k 2 ), in both cases, the linear approximation and the unlinearized version. • We have been able to retrieve a general solution for (1.5) making use of the theory of Lie-Hamilton systems and it is expressed in terms of two particular solutions and a very challenging choice of two constants (c, k). One can check the choices of constants for every case in the legends of the graphics. • It is interesting to notice that using the geometric methods of the Lie-Hamilton theory, one retrieves the introduced hamiltonian (1.9) by [43] as performed in (3.11). • Furthermore, we have generalized this hamiltonian (1.9) by introducing a deformation parameter, know as the quantum deformation parameter using the Poisson-Hopf algebra method. The deformed hamiltonian (4.14) gives rise to a new integrable SIS model in which one includes a quantum deformation parameter (4.15).
Relation with Covid19 pandemic. One may wonder how the current pandemic of COVID19 could be related to a SISf-pandemic model. The SISf model is a very first approximation for a trivial infection process, in which there is only two possible states of the individuals in the population: they are either infected or susceptible to the infection. Hence, this model does not provide the possibility of adquiring immunity at any point. It seems that COVID19 provides some certain type of immunity, but only to a thirty percent of the infected individuals, hence, a SIR model that considers "R" for recuperated individuals (not susceptible anymore, i.e., immune) is not a proper model for the current situation. One should have a model contemplating immune and nonimmunized individuals. Unfortunately, we are still in search of a stochastic Hamiltonian model including potential immunity and nonimmunity.
Future work. As future work, we would like to extend our study to more complicated comparmental models, although at a first glance we have not been able to identify more Lie systems, at least in their current PDE form. We suspect that the Hamiltonian description of these compartmental models could nonetheless behave as a Lie system, as it has happened in our presented case. This shall be part of our future endeavors. Moreover, one could inspect in more meticulous detail how the solutions of the quantumly-deformed system (4.15) recover the nondeformed solutions when the introduced parameter tends to zero. We need to further study how this precisely models a heat bath and if this new integrable system could correspond to other models apart from infectious models. We would like to figure out whether it is possible to modelize subatomic dynamics with the resulting deformed hamiltonian (4.14).
Apart from that, there exists a stochastic theory of Lie systems developed in [31] that could be another starting point to deal with compartmental systems. In the present work we were lucky to find a theory with fluctuations that happened to match a stochastic expansion, but this is rather more of an exception than a rule. Indeed, it seems that the most feasible way to propose stochastic models is using the stochastic Lie theory instead of expecting a glimpse of luck with fluctuations.
As we have stated, finding particular solutions is by no means trivial. The analytic search is a very excruciating task. We think that in order to fit particular solutions in the superposition principle, one may need to compute these particular solutions numerically. Some specific numerical methods for particular solutions of Lie systems can be devised in [46] . Let A be a vector space over a field (in this work R). Consider the multiplication map m and a unit element map ι, given by
We say that the triple (A, m, ι) is an algebra. This algebra is said to be associative if the multiplication is associative, i.e.,
for all x, y, and z in A. Otherwise it is referred to as a nonassociative algebra. A mapping between two algebras is called homomorphism if it respects the multiplicatios and unit from one algebra to the other.
An algebra is called a Lie algebra, denoted by g, if the multiplication is skew-symmetric and satisfies the Jacobi identity [6, 24] . In this case, the multiplication operation is a Lie bracket [•, •]. The Jacobi identity is
for all x, x and x in g. If this identity is fulfilled, then the nonassociativity of the multiplication is followed. A mapping between two algebras is called a Lie algebra homomorphism if it preserves the Lie algebra structure.
Let B be an arbitrary subspace of a Lie algebra g. Let N be an n-dimensional manifold (a locally Euclidean space). A differentiable curve on N through the point x is a function γ : R → Q with, say, γ (0) = x. Let us define now an equivalence relation in the set of differentiable curves passing through x. Two curves γ andγ are equivalent if they take the same value at x and if the directional derivative of functions along them at x are the same, namely,
for all functions f : N → R. A tangent vector v(x) at x is an equivalence class of curves at x. The set of all equivalence classes of vectors, that is, the set of all tangent vectors at x is tangent space T x N at x in N . It is possible to show that T x N admits an n-dimensional vector space structure.
The union of all tangent spaces T x N in x ∈ N is
and it is a 2n-dimensional manifold that is called tangent bundle of N . From the point of view of applications, if a manifold is the configuration space of a physical system, then, the tangent bundle is the velocity phase space of the system. That is, it consists of all possible positions and all possible velocities. There is a projection, called tangent fibration , τ N mapping a tangent vector to its base point
Vector fields. A section of the tangent fibration is a vector field X on N mapping a vector to each point in N , that is a map Y : N → T N such that τ N • Y : id N where id N is the identity map on N . The set of vector fields X(N ) on a manifold Q admits a module structure over the ring of functions.
An integral curve of a vector field X with initial condition x, is a curve γ passing through γ(0) = x so that
From the theorem of existence and uniqueness of solutions for ODEs, we know that for every x there exists a unique integral curve γ satisfying (A.8). The flow of Y is a smooth one parameter group of diffeomorphisms defined by means of integral curves as follows
where γ(0) = x. This realization of vector fields enables us to define the directional derivative of a function f along a vector field Y as
where γ is an integral curve of Y . Notice that, if f is a constant of motion, then, Y (f ) = 0.
Lie Algebras. The space of all vector fields, denoted by X(N ), is a Lie algebra if it is endowed with a Lie bracket. For any two vector fields Y 1 and Y 2 on the base manifold N , the Lie bracket reads
for all smooth functions f on N . Here, the notation Y (f ) is simply the directional derivative of the function f in the direction of Y . A distribution generated by a finite set of vector fields, say Y 1 , ... Y r , is the collection of subspaces D x , generated by the vectors Y 1 (x), ... Y r (x), in all points. It is called an involutive distribution if the Lie algebra is closed for the set Y 1 , ... Y r .
Local realization in 2D. Assume now that N is a two dimensional manifold with local coordinates (q, p), then the tangent bundle T N admits the induced coordinates (q, p;q,ṗ). In this case, a vector field is written as
where k and g are real valued functions on N . In this coordinate system, the dynamical equations governed by X determine a system of nonautonomous ordinary differential equations (A.13)q = k(q, p),ṗ = g(q, p).
An integral curve to the vector field X in (A.12) is a solution of the system (A.13) given by q = q(t) and p = p(t).
A tangent space T q Q at point q admits a (finite dimensional) vector space structure. So, there exists the linear algebraic dual of T q Q, which we denote by T * q Q. We call T * x N a cotangent space. By collecting all of them, we define the cotangent bundle (A.14)
There is a projection π Q from T * Q to Q, mapping a covector α(q) to its base point q. Sections θ of the cotangent bundle are one-forms on Q satisfying π Q • θ = id Q .
2D Symplectic Manifold. Assume now that N is a two dimensional planar manifold with coordinates (q, p). For this case, we can understand N = T * Q for some one-dimensional Q.
There is a canonical symplectic two-form ω on N . A symplectic two-form is a closed, skewsymmetric and nondegenerate bilinear mapping taking two vector fields to the space of scalars, i.e., ω takes two vector fields Y 1 and Y 2 on N to a unique function. In this case, the symplectic two-form is computed to be (A. 15 ) ω = dq ∧ dp. for any z in g * , and any two function(als) h and g on g * , [37] . Here, the bracket on the right hand side is the Lie algebra bracket on g. Let us note that δh/δz and δg/δz denote the partial (for infinite dimensional cases Fréchet) derivatives of the function(al)s. Under the assumption of the reflexivity, they are elements of g.
Poisson bivectors. Being a derivation for each factor, a Poisson structure determines a unique bivector field Λ ∈ Γ(∧ 2 T N ), that we call a Poisson bivector, such that Referring to this mapping, we define a Hamiltonian vector field (A.18), as follows
where dh is the (exterior) derivative of h. The space of Hamiltonian vector fields induces an integrable generalized distribution F Λ on N associated to Λ. Pointwisely, each fiber F Λ x of this distribution is defined to be the image space of the mappingΛ, that is {X x | X ∈ ImΛ}. Here, the leaves are symplectic manifolds with respect to the restrictions of Λ [55] .
In order to arrive the formal definition of a coalgebra, one simply reverse the directions of the (multiplication and unit) arrows in (A.1). Accordingly, a vector space A is called a coalgebra if it admits a comultiplication and a counit given by
: A −→ R, (z) = 1 (A. 23) respectively. We ask that these operations must satisfy the relation
where Id is being the identity mapping. We denote a coalgebra by a triple (A, ∆, ). A mapping is called a coalgebra homomorphism if it respects the coalgebra structures. is both a Poisson and a Hopf algebra, then it is a Poisson-Hopf algebra and it is denoted by the tuple (A, m, ι, {•, •}, ∆, , γ) .
The space of smooth functions C ∞ (g * ) on the dual of a Lie algebra g is a Hopf algebra relative to its natural associative algebra with unit provided that for every z, z in g * , and f, g, h in C ∞ (g * ). Further, considering the Lie-Poisson bracket (A.19), C ∞ (g * ) turns out to be a Poisson-Hopf algebra.
Symmetric (co)algebra. The symmetric algebra S(g) of a (finite dimensional) Lie algebra g is the smallest commutative algebra containing g. To reach such algebra, we do the following. The second tensor power g ⊗ g of the Lie algebra is the space of real valued bilinear maps on the dual space. Iteratively, the kth tensor power g ⊗k is the space of real valued k-linear maps. Taking the direct sum of the tensor powers of all orders, we arrive at the tensor algebra Tg of g. Here, the multiplication is
We consider a basis {x 1 , . . . , x r } of the Lie algebra g. The space generated by the elements
is an ideal, denoted by R, of the tensor algebra Tg. The quotient space Tg/R is called a symmetric algebra and denoted by S(g). The elements of S(g) can be regarded as polynomial functions on g * , so, we can endow it with the Lie-Poisson bracket (A.19) that makes S(g) a Poisson algebra. One can show that S(g) can always be endowed with a coalgebra structure by introducing the comultiplication which, clearly, is also a Poisson algebra homomorphism.
Universal enveloping algebras. We consider once more the tensor algebra Tg of a Lie algebra g with a basis {x 1 , . . . , x r }. We now define the space generated by the elements
This space is an ideal, denoted by L, of the tensor algebra. The quotient space Tg/L is called universal enveloping algebra of g, and it is denoted by U (g). U (g) is the biggest associative algebra containing all possible representations of g. See that, if the Lie algebra bracket on g is trivial, then, the universal enveloping algebra U (g) is equal to the symmetric algebra S(g).