On the Integrability of the SIR Epidemic Model with Vital Dynamics

In this paper, we study the SIR epidemic model with vital dynamics _ S = −βSI + μðN − SÞ, _ I = βSI − ðγ + μÞI, _ R = γI − μR, from the point of view of integrability. In the case of the death/birth rate μ = 0, the SIR model is integrable, and we provide its general solutions by implicit functions, two Lax formulations and infinitely many Hamilton-Poisson realizations. In the case of μ ≠ 0, we prove that the SIR model has no polynomial or proper rational first integrals by studying the invariant algebraic surfaces. Moreover, although the SIR model with μ ≠ 0 is not integrable and we cannot get its exact solution, based on the existence of an invariant algebraic surface, we give the global dynamics of the SIR model with μ ≠ 0.


Introduction and Statement of the Main Results
Over the past one hundred years, the mathematical modelling of epidemics has been the object of a large number of studies. The Susceptible-Infected-Recovered (SIR) model is one of the most interesting and best understood nonlinear epidemic models [1,2]. The first SIR model was proposed by Kermack and McKendrick in 1927 [2]. After that, many different epidemic models including time delay, age structure, space factor, white noise, multigroup, and seasonality have been proposed and studied (see [1] and the references therein). Observing the spread of marketing message is analogous to an epidemic, Rodrigues and Fonseca [3] used a SIR model to study the effects of a viral marketing strategy. We consider the SIR epidemic model with vital dynamics which is given by where S is the number of healthy individuals who are susceptible to the disease, I is the number of infected individuals who can transmit the disease to the healthy ones, and R is the number of individuals who have been infected and then recovered from the disease and the parameters β, γ, μ, and N denote the average number of contacts per infective per day, the recovery rate, the death rate, and the initial total fixed number of the individuals, respectively. It is assumed that the birth rate is equal to the death rate in this model (1). For μ = 0, it is called the SIR model without vital dynamics and was proposed by Kermack and McKendrick [2]. The aim of this paper is to study system (1) from the integrability point of view. The integrability of differential equations has been an old and important problem and has attracted much attention. Many scholars have developed a lot of ways to deal with the integrability for both partial differential equations (see [4][5][6][7][8][9][10][11][12][13] for instance) and ordinary differential equations (see [14][15][16] for instance). In particular, the existence of first integrals plays a crucial role in the integrability of ordinary differential equations [17][18][19][20]. If the considered system of ordinary differential equations admits a straight line solution, by the differential Galois method, one can usually prove that this system has no rational first integrals for almost all the parameters [21][22][23]. However, this method cannot tell whether this system is integrable for the remaining parameters. The SIR epidemic model (1) is exactly the case. Another tool that deals with the integrability of 3D differential systems is the Darboux theory of integrability [15,16], which is a useful tool to find first integrals for polynomial ordinary differential equations, and has been successfully applied in many nonlinear models [24][25][26][27]. This theory can also help us make a more precise analysis of the global dynamics of a system topologically (see [28,29] and the references therein). In the framework of the Darboux theory of integrability, we will give a complete classification of the irreducible Darboux polynomials, of the polynomial first integrals, of the proper rational first integrals, and of the algebraic integrability for the SIR model.
In the case of the constant population, i.e., μ = 0, the SIR model (1) admits a first integral and can be reduced into a planar system with respect to the variables ðS, IÞ. The integrability of the reduced planar system has been investigated extensively by different methods, namely, the Adomian decomposition method [30], the homotopy analysis method [31], and the variational iteration method [32]. Recently, Bohner et al. investigated a rational SIR model with the constant population and the time-dependent coefficients and present an alternative solution method to Gleissner's approach [33]. Based on a quantum mechanical method, Williams et al. provided an exact analytical solution of the stochastic SIR model with the constant population [34]. However, in the case of the varying population, i.e., μ ≠ 0, the integrability of the full SIR model is an area where little research has been done. Harko et al. reduced the SIR model (1) into the Abel equation and provided its form series solution by using a perturbation approach [35]. To our knowledge, the first integrals of the SIR model (1) with μ ≠ 0 have not been investigated previously. The aim of this work is to cover this gap.
Recall that a real polynomial FðS, I, RÞ ∈ ℝ½S, I, R is called a Darboux polynomial of system (1) if it satisfies for some polynomial KðS, I, RÞ ∈ ℝ½S, I, R, called a cofactor of F. Clearly, if F is a Darboux polynomial, then F = 0 is invariant with respect to the flow of system (1). Hence, we call F = 0 an invariant algebraic surface of system (1). It is well known that a polynomial function f = f n 1 1 , ⋯, f n m m is a Darboux polynomial iff each irreducible factor f i , i = 1, ⋯, m is also a Darboux polynomial. Hence, for simplicity, we only focus on the irreducible Darboux polynomials of system (1).
Our main result on the integrability of system (1) is as follows, which characterizes all irreducible Darboux polynomials of system (1). In addition, Darboux polynomials can help us construct the algebraic (polynomial or rational) first integrals. By Theorem 1, we can easily obtain the following result.

Corollary 2.
The following statements hold for system (1).
(1) It has a polynomial first integral if and only if μ = 0, and in this case, the polynomial first integral is F = S + I + R (2) It has no any proper rational first integral

(3) It is not algebraically integrable
It is not surprising that the SIR model without vital dynamics, i.e., μ = 0, has a first integral F = S + I + R which is the total number of individuals in the given population. Furthermore, system (1) with μ = 0 has another first integral GðS, I, RÞ = γ ln S + βR, which can help us compute the number Sð+∞Þ of individuals that will never contract the infection. Let us mention that the first integral G = γ ln S + βR can be built by the classical Jacobi last multiplier method, observing system (1) admits a Jacobi last multiplier M = 1/S I. In a word, the system with μ = 0 is a completely integrable system with two functionally independent first integrals F, G. Based on this fact, we have the following remarks on its integrability.
(i) The orbits of system (1) with μ = 0 are contained in the curves and its general solutions are given by the following implicit functions: where c 1 , c 2 , c 3 are constants (ii) System (1) with μ = 0 has two Lax formulations _ L i = ½L i , N i , i = 1, 2, where the matrices L i and N i are given by 2 Advances in Mathematical Physics (iii) System (1) with μ = 0 has infinitely many Hamilton-Poisson realizations parameterized by the group S Lð2, ℝÞ; that is, ðℝ 3 , f·, · g ab , H cd Þ is a Hamilton-Poisson realization where the Poisson bracket f·, · g reads for any functions f , g ∈ C ∞ ððℝ + Þ 3 , ℝÞ, v = 1/ðSIÞ, the Casimir function the Hamiltonian function and the coefficients a, b, c, d ∈ ℝ such that ad − bc = 1 As mentioned above, the SIR model with μ = 0 is integrable with two first integrals, which implies it is orbitally equivalent to a linear differential system. Meanwhile, the SIR model with μ ≠ 0 has no polynomial/rational first integrals. However, based on the existence of an invariant algebraic surface, we can characterize the global phase portraits of the SIR model with μ ≠ 0, which helps us understand the final evolutions of this model and the spread of the disease. Theorem 3. The following statements hold for system (1) with μ ≠ 0.
(a) All orbits with initial points not on the invariant algebraic surface F 1 = 0 and not at infinity are heteroclinic ones, which all positively approach the surface F 1 = 0 and negatively go to infinity (b) The dynamics of system (1) at the infinity S 2 is topologically equivalent to the one described in Figure 1 (c) System (1) on the invariant algebraic surface F 1 = 0 has four topologically different phase portraits, which are described in Figure 2 Theorem 3 provides the final evolutions of this model. As we have shown, there are only two possible final evolutions for this model. In the first type, both infective and recovered individuals tend to zero; that is, the disease fails to spread, and in the second type, the infective individuals cannot tend to zero; that is, the disease is endemic.
The paper is organized as follows: in Section 2, we prove Theorem 1 and Corollary 2. The proof of Theorem 3 will be given in Section 3. In the last section we draw our conclusions, including some discussions on the biological meaning of our results.

Proof of Theorem 1 and Corollary 2
2.1. Proof of Theorem 1. When μ = 0, it is easy to see that S + I + R − N is a Darboux polynomial of system (1) with zero cofactor. In what follows, we deal with the case μ ≠ 0. For simplicity, we introduce the change of variables and a time rescaling t =t/ðμ + γÞ. Then, system (1) becomes where the prime ′ denotes a derivative with respect to the new timet and the coefficients a = μNβ/ðμ + γÞ 2 > 0 and b = μ/ðμ + γÞ ∈ ð0, 1Þ. Clearly, to complete the proof of Theorem 1, we need only to prove the next result. Suppose f ðx, y, zÞ is a Darboux polynomial of system (10) with the cofactor kðx, y, zÞ, that is,

Advances in Mathematical Physics
Comparing the degree of both sides of (11), we get that the cofactor kðx, y, zÞ is a polynomial with degree less than two. Without loss of generality, we set kðx, y, zÞ = k 0 + k 1 x + k 2 y + k 3 z. Then, we claim k 3 = 0. In fact, we write f as powers in the variable z, i.e., f ðx, y, zÞ = ∑ n j=0 f j ðx, yÞz j , where each f j is a polynomial in the variables x and y and f n ≠ 0. Computing the coefficient in (11) of z n+1 , we have 0 = k 3 f n ðx, yÞ which implies k 3 = 0.
We make the change of the variables Let n be the highest weight degree in the weight homogeneous components of f in x, y, z with the weight exponent ð1, 1, 0Þ. Set where F i is a weight homogeneous polynomial of weight degree ðn − iÞ in ðX, Y, ZÞ. Since f is a Darboux polynomial of system (10) with the cofactor k, it is not difficult to check that F is a Darboux polynomial of system (12) with the cofactor K, that is, where we still use x, y, z instead of X, Y, Z.
Equating the terms with α 0 in (14) yields where L stands for a linear partial differential operator To solve (15), we introduce the change of variables and transform (15) into where F 0 is F 0 written in terms of u, v, w. Solving (18) yields F 0 = G 0 ðv, wÞðv − uÞ k 1 u −k 2 with G 0 being an arbitrary smooth function in v, w. Clearly, in order for F 0 ðx, y, zÞ = F 0 ðu, v, wÞ = G 0 ðx + y, x exp ðzÞÞy k 1 x −k 2 to be a weight homogeneous polynomial of degree n in x, y, z, we must have k 1 = m, k 2 = −s and for some nonzero number c 0 ∈ ℝ/f0g and nonnegative integers s, m, l ∈ ℕ ∪ f0g. Similarly, equating the terms with α in (14) leads to Substituting F 0 = a 0 x s y m ðx + yÞ l into (20) yields Proceeding as above, by (17) we rewrite (21) as x y z q p Figure 1: Phase portrait at infinity on the Poincaré ball of system (1). Note that system (1) has two closed curves of equilibria fx = 0, y 2 + z 2 = 1g ∪ fy = 0, x 2 + z 2 = 1g and two nodes q, p. 4 Advances in Mathematical Physics Integrating the above equation with respect to u, we obtain with G 1 being an arbitrary smooth function. In order for F 1 to be a weight-homogeneous polynomial of weight degree n − 1, we must have G 1 = c 1 ∈ ℝ which is a constant, k 0 = −bs − bl − m, and Equating the terms with α j in (14) for j = 2, 3, ⋯, n, we have In particular, for j = 2 we substitute F 0 , F 1 into (25) and obtain or equivalently

Advances in Mathematical Physics
Solving the above linear differential equation, we get In order for F 2 to be a weight-homogeneous polynomial of weight degree n − 1, we have bc 1 + al + as = 0, Proceeding as above, from (25) with j = 3 and (29), we get Working in a similar way to solve F 2 , we can prove that where dilogð·Þ is the dilogarithm function defined by and the coefficients M i , i = 1, ⋯, 8 are given by In order for F 3 to be a weight-homogeneous polynomial of weight degree n − 3, we have M 5 = M 6 = M 7 = 0 and the function G 3 is a constant, which implies s = 0, with c 3 being a constant. Using the same argument similar to that above, we solve F 4 and get Advances in Mathematical Physics with c 4 ∈ ℝ. Furthermore, by recursive calculations, we can prove that Therefore, the Darboux polynomial of system (10) and the cofactor k = Kj α=1 = mx − ðb + mÞ, which completes the proof.  (2) of Corollary 2. Assume system (1) has a proper rational first integral f = P/Q with P, Q being relative prime; then, P, Q are two different Darboux polynomials with the same nonzero cofactor; see [16] for instance, which contradicts Theorem 1. (1) and (2) and a well-known result that a polynomial vector field in ℝ n is algebraically integrable if and only if it has n − 1 functionally independent rational first integrals.

Proof of Statement (b) of Theorem 3.
In order to understand the global dynamics of system (1) at infinity, we will use the Poincaré compactification in ℝ 3 and analyze the flow at infinity for the local charts U i and V i , i = 1, 2, 3. See [36] for more details on the Poincaré compactification of a polynomial vector field.
3.2.1. In the Local Chart U 1 and V 1 . Making the change of variables ðS, I, RÞ = ðz −1 3 , z 1 z −1 3 , z 2 z −1 3 Þ and the time rescaling dτ = z −1 3 dt, we obtain the Poincaré compactification of system (1) in the local chart U 1 The plane z 3 = 0 at infinity is invariant, which corresponds to the points on the sphere at infinity, and so system (1) restricted to z 3 = 0 becomes which has a line of singular points z 1 = 0 and an isolated singular point ðz 1 , z 2 Þ = ð−1, 0Þ. Moreover, system (40) has a first integral Φ = ðz 1 + 1Þ/z 2 . Using this first integral, the phase portrait on the local chart U 1 is described in Figure 3. The flow on the local chart V 1 is the same with the flow of U 1 by reversing the time since the compactified vector field in V 1 coincides with the vector field in U 1 multiplied by −1 [36].

Advances in Mathematical Physics
When z 3 = 0, system (41) becomes Clearly, system (42) coincides with system (40) by reversing the time. Then, the phase portrait on the local chart U 1 is described in Figure 4. Again, the flow in the local chart V 2 is the same as the flow in the local chart U 2 by reversing the time.

In the Local Chart
On the invariant plane z 3 = 0, system (43) is reduced to System (44) has two lines of singular points z 1 = 0 and z 2 = 0. It also has a first integral Φ = z 1 + z 2 , which helps us get its phase portrait as shown in Figure 5. The flow at infinity in the local chart V 3 is the same as the flow on the local chart U 3 by reversing the time.
To summarize the above analysis, one obtains a global picture of the dynamical behavior of system (1) on sphere S 2 at infinity: it has two closed curves of equilibria fx = 0, y 2 + z 2 = 1g ∪ fy = 0, x 2 + z 2 = 1g and two nodes (see Figure 1). Restricted to this invariant algebraic surface, system (1) becomes

Advances in Mathematical Physics
Next, we turn to study the infinite singularities by the Poincaré compactification in ℝ 2 . In the local chart U 1 , where S = z −1 2 , I = z 1 z −1 2 , and dτ = z −1 2 dt, we have It has two singularities ð0, 0Þ and ð−1, 0Þ on the invariant line z 2 = 0. The eigenvalues of the linear part at ð−1, 0Þ are −β and −β, which implies this point is a stable node. By Theorem 2.19 in [37], we see that ð0, 0Þ is a saddle-node point consisting of two hyperbolic sectors with one parabolic sector. Similarly, we take the change of variables S = z 1 z −1 2 and I = z −1 2 and the time rescaling dτ = z −1 2 dt to obtain the Poincaré compactification of system (45) in the local chart U 2 This system has an unstable node ð−1, 0Þ and a saddle node ð0, 0Þ.
Finally, we show that system (45) has no limit cycles in ℝ 2 . Since I = 0 is an invariant manifold, limit cycles (if exits) must be in I > 0 or I < 0. In the region fðS, IÞ: I > 0g, observing it follows from the Dulac Theorem that system (45) has no limit cycles. In the region fðS, IÞ: I < 0g, system (45) has either no singular points or a saddle, which implies system (45) has no limit cycles. Summarizing the above analysis, we obtain that the dynamics of the system (1) on F 1 = 0 is topologically described by one of Figure 2 on the Poincaré disk.

Conclusions
In this paper, by studying the invariant algebraic surfaces, we show that μ = 0 is the only value of the parameters for which the SIR epidemic model is integrable, and in this case, we provide its general solutions by implicit functions, two Lax formulations, and infinitely many Hamilton-Poisson realizations. When μ ≠ 0, the SIR model has no any algebraic first integral and we cannot get its exact solutions. However, based on the existence of the invariant algebraic surfaces, we characterize the topological structure of orbits for the SIR epidemic model with μ ≠ 0. Moreover, if the SIR model has a positive death/birth rate μ, the disease will ultimately approach either the disease-free steady state P 0 on F 1 = 0 or the endemic steady state P 1 on F 1 = 0, depending on the basic reproduction number R 0 . In case the R 0 is small, saying less than one, the SIR model has no positive equilibrium point and the disease approaches the disease-free steady state P 0 on F 1 = 0, which implies the disease will always fail to spread. In case the R 0 is big, saying larger than one, the SIR model admits a new positive equilibrium point through the transcritical bifurcation and the disease approaches the endemic steady state P 1 on F 1 = 0. In case R 0 is equal to one, all the solutions in ðℝ + Þ 3 tend to P 0 on F 1 = 0, which implies the disease is supposed to be controlled and the entire population tends to be healthy, but is susceptible to reinfection. These facts show that the basic reproduction number has played an important role in the spread of disease and the topological structure of orbits for the SIR model. The approach we used in this work may contribute to the understanding of the dynamics of the more general epidemic models.

Data Availability
All data included in this study are available upon request by contact with the corresponding author.

Conflicts of Interest
The authors declare that they have no conflicts of interest.