# On the Integrability of the SIR Epidemic Model with Vital Dynamics

Authors: Ding Chen
Journal: Advances in Mathematical Physics
 (2020)

---

## Abstract

In this paper, we study the SIR epidemic model with vital dynamicsṠ=−β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.


---

## Body

## 1. 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(1)Ṡ=−βSI+μN−S,İ=βSI−γ+μI,Ṙ=γI−μR,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–13] for instance) and ordinary differential equations (see [14–16] for instance). In particular, the existence of first integrals plays a crucial role in the integrability of ordinary differential equations [17–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–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–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 polynomialFS,I,R∈ℝS,I,R is called a Darboux polynomial of system (1) if it satisfies
(2)−βSI+μN−S∂F∂S+βSI−γ+μI∂F∂I+γI−μR∂F∂R=KF,for some polynomial KS,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=f1n1,⋯,fmnm is a Darboux polynomial iff each irreducible factor fi,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).Theorem 1.
All irreducible Darboux polynomials of system (1) consist of F1S,I,F=S+I+R−N with the cofactors K1S,I,R=−μ and F2S,I,F=I with the cofactor K2S,I,R=βS−γ−μ.In addition, Darboux polynomials can help us construct the algebraic (polynomial or rational) first integrals. By Theorem1, 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 integrableIt 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 GS,I,R=γlnS+β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=γlnS+βR can be built by the classical Jacobi last multiplier method, observing system (1) admits a Jacobi last multiplier M=1/SI. 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(3)S,I,R∣FS,I,R=c1,GS,I,R=c2,and its general solutions are given by the following implicit functions:
(4)R=−γβlnS+c2β,I=γβlnS−S−c2β+c1,∫1SβS−γlnS+c2−βc1dS=t+c3,where c1,c2,c3 are constants
(ii)
System (1) with μ=0 has two Lax formulations L̇i=Li,Ni,i=1,2, where the matrices Li and Ni are given by(5)L1≔0−2R2I2R0−2S−2I2S0,N1≔0−1βIS+2I2R10−βSI+2S−γI2R−βIS+2I2RβSI+2S−γI2R0,L2≔0−1i1+S+R10−1−S+iR−i1+S−R1+S−iR0,N2≔00IβS+iγ00IiβS−γ−IβS+iγ−IiβS−γ0(iii)
System (1) with μ=0 has infinitely many Hamilton-Poisson realizations parameterized by the group SL2,ℝ; that is, ℝ3,·,·ab,Hcd is a Hamilton-Poisson realization where the Poisson bracket ·,· reads(6)f,gab≔v∇Cab·∇f×∇g,for any functions f,g∈C∞ℝ+3,ℝ, v=1/SI, the Casimir function
(7)Cab=alnS+βR+bS+I+R,the Hamiltonian function
(8)Hcd=clnS+βR+dS+I+R,and the coefficients a,b,c,d∈ℝ such that ad−bc=1As 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 surfaceF1=0 and not at infinity are heteroclinic ones, which all positively approach the surface F1=0 and negatively go to infinity(b)
The dynamics of system (1) at the infinity S2 is topologically equivalent to the one described in Figure 1(c)
System (1) on the invariant algebraic surface F1=0 has four topologically different phase portraits, which are described in Figure 2Figure 1
Phase portrait at infinity on the Poincaré ball of system (1). Note that system (1) has two closed curves of equilibria x=0,y2+z2=1∪y=0,x2+z2=1 and two nodes q,p.Figure 2
Phase portraits of system (1) on the invariant planes F1=0.
(a)
R0<1(b)
R0=1(c)
R0>1,Δ≥0(d)
R0>1,Δ<0Theorem3 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 Section2, 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.


## 2. Proof of Theorem1 and Corollary 2

### 2.1. Proof of Theorem1

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
(9)S=γ+μβx,I=γ+μβy,R=μβzand a time rescaling t=t~/μ+γ. Then, system (1) becomes
(10)x′=a−bx−xy,y′=−y+xy,z′=y−bz,where the prime ′ denotes a derivative with respect to the new time t~ 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.Proposition 4.
All irreducible Darboux polynomials of system (10) consist of fx,y,z=x+y+1−bz−a/b with the cofactors k1x,y,z=−b and f2x,y,z=y with the cofactor k2x,y,z=x−1.Supposefx,y,z is a Darboux polynomial of system (10) with the cofactor kx,y,z, that is,
(11)a−bx−xy∂f∂x+xy−y∂f∂y+y−bz∂f∂z=kf.Comparing the degree of both sides of (11), we get that the cofactor kx,y,z is a polynomial with degree less than two. Without loss of generality, we set kx,y,z=k0+k1x+k2y+k3z. Then, we claim k3=0. In fact, we write f as powers in the variable z, i.e., fx,y,z=∑j=0nfjx,yzj, where each fj is a polynomial in the variables x and y and fn≠0. Computing the coefficient in (11) of zn+1, we have 0=k3fnx,y which implies k3=0.We make the change of the variablesx=α−1X, y=α−1Y, z=Z, and t~=αT and rewrite (10) into
(12)dXdT=−XY−αbX+α2a,dYdT=XY−αY,dZdT=Y−αbZ.Letn be the highest weight degree in the weight homogeneous components of f in x,y,z with the weight exponent 1,1,0. Set
(13)FX,Y,Z=αnfα−1X,α−1Y,Z=F0X,Y,Z+αF1X,Y,Z+⋯+αnFnX,Y,Z,KX,Y,Z=αkα−1X,α−1Y,Z=k1X+k2Y+αk0,where Fi 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,
(14)−xy−αbx+α2a∑i=0nαi∂Fi∂x+xy−αy∑i=0nαi∂Fi∂y+y−αbz∑i=0nαi∂Fi∂z=k1x+k2y+αk0∑i=0nFi,where we still use x,y,z instead of X,Y,Z.Equating the terms withα0 in (14) yields
(15)LF0=k1x+k2yF0,where L stands for a linear partial differential operator
(16)L≔−xy∂∂x+xy∂∂y+y∂∂z.To solve (15), we introduce the change of variables
(17)u=x,v=x+y,w=xexpzand transform (15) into
(18)uu−v∂F0¯∂u=k1−k2u+k2vF0¯,where F0¯ is F0 written in terms of u,v,w. Solving (18) yields F0¯=G0v,wv−uk1u−k2 with G0 being an arbitrary smooth function in v,w. Clearly, in order for F0x,y,z=F0¯u,v,w=G0x+y,xexpzyk1x−k2 to be a weight homogeneous polynomial of degree n in x,y,z, we must have k1=m, k2=−s and
(19)F0=a0xsymx+yl,s+m+l=n,for some nonzero number c0∈ℝ/0 and nonnegative integers s,m,l∈ℕ∪0.Similarly, equating the terms withα in (14) leads to
(20)LF1=mxF1−syF1+bx∂xF0+y∂yF0+bz∂zF0+k0F0.SubstitutingF0=a0xsymx+yl into (20) yields
(21)LF1=mx−syF1+a0xsymx+y1−1bs+bl+m+k0x+bs+l+m+k0y.Proceeding as above, by (17) we rewrite (21) as
(22)uu−v∂F1¯∂u=mu+su−vF1¯+a0vl−1usv−umbs+bl+m+k0u+bs+l+m+k0v−u.Integrating the above equation with respect tou, we obtain
(23)F1=F1¯=a0usv−umvl−1bs+bl+m+k0lnu−v+bs+l+m+k0lnwu+G1v,w=a0xsymx+yl−1bs+bl+m+k0lny+bs+l+m+k0z+G1x+y,xexpz,with G1 being an arbitrary smooth function. In order for F1 to be a weight-homogeneous polynomial of weight degree n−1, we must have G1=c1∈ℝ which is a constant, k0=−bs−bl−m, and
(24)F1=a0xsymx+yl−1l1−bz+c1.Equating the terms withαj in (14) for j=2,3,⋯,n, we have
(25)LFj=mx−syFj+bx∂xFj−1+y∂yFj−1+bz∂zFj−1+k0Fj−1−a∂xFj−2.In particular, forj=2 we substitute F0,F1 into (25) and obtain
(26)LF2=mx−syF2+a0xsymx+yl−2l1−bz+c1−blx+y+bl−1x+l−1y+bl1−bx+yz−asx+y2x−alx+yor equivalently
(27)uu−v∂F2¯∂u=mu+su−vF2¯+a0usv−umvl−2−lnwu1−b2ll−1u−v−bc1+alu+c1bl+c1−c1l+alu−v−asv2u.Solving the above linear differential equation, we get(28)F2=F2¯=a0xsymx+yl−2ll−11−b2z22−bc1+al+asln∣y∣−c1bl+c1−c1l+al+asz−asx+yx+G1x+y,xexpz.In order forF2 to be a weight-homogeneous polynomial of weight degree n−1, we have bc1+al+as=0, G2=c2∈ℝ which implies
(29)F1=a0xsymx+yl−1l11−bz−l1ab−asb,F2=a0xsymx+yl−2l21−b2z2−l22a1−bzb+asl−1b−1zb−x+yx+c2.Proceeding as above, from (25) with j=3 and (29), we get
(30)LF3=mx−syF3+a0xsymx+yl−3bzx+yll−11−b2z−al+sl−11−bb+ll−11−b2z22+al+sl−1b−1zb+c2−2bx−bl+2−ly+as2b+2m+al+sb+l−1x+y2x+a2l−1l+sb−asl−1x+y−asx+y2l1−bzx−all−1x+y1−bz.Working in a similar way to solveF2, we can prove that
(31)F3=a0xsymx+yl−3M1z3+M2z2+M3+M4x+M4yz+M5+M6zlny+M7dilog−yx+M8x+yx+G3x+y,xexpz,where dilog· is the dilogarithm function defined by
(32)dilogx=∫1xlnt1−tdt,and the coefficients Mi,i=1,⋯,8 are given by
(33)M1=l31−b3,M2=−l33a1−b2b−as1−b2l−1,M3=as2b+2m+a2l−1+sb−c2bl+2c2−c2l−a2ll−1b,M4=as1−bl,M5=a2ll−1b−2bc2+asa2l+s−1b+2b+2m,M6=−as1−b2l−1,M7=−as1−b2l−1,M8=as2b+2m+al+sb+2l−1−bl.In order forF3 to be a weight-homogeneous polynomial of weight degree n−3, we have M5=M6=M7=0 and the function G3 is a constant, which implies
(34)s=0,c2=l2a2b2,F1=a0ymx+yl−1l11−bz−ab,F2=a0ymx+yl−2l21−bz−ab2,F3=a0ymx+yl−3l31−b3z3−3a1−b2z2b+3a21−bzb2+c3,with c3 being a constant. Using the same argument similar to that above, we solve F4 and get
(35)c3=−ab3,F4=a0ymx+yl−4l41−b4z4−4a1−b3z3b+6a21−b2z2b2−4a31−bzb3+c4,with c4∈ℝ. Furthermore, by recursive calculations, we can prove that
(36)Fj=a0ymx+yl−jlj1−bz−abj,j=0,⋯,n.Therefore, the Darboux polynomial of system (10)
(37)f=Fα=1=∑j=0na0ljymx+yl−j1−bz−abj=a0ym1−bz−abl,and the cofactor k=Kα=1=mx−b+m, which completes the proof.


### 2.2. Proof of Corollary2

#### 2.2.1. Proof of Statement (1) of Corollary2

It follows from Theorem1 and the fact that a function f is a polynomial first integral if and only if it is a Darboux polynomial with the zero cofactor.


#### 2.2.2. Proof of Statement (2) of Corollary2

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.


#### 2.2.3. Proof of Statement (3) of Corollary2

It follows from statements (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.


## 2.1. Proof of Theorem1

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
(9)S=γ+μβx,I=γ+μβy,R=μβzand a time rescaling t=t~/μ+γ. Then, system (1) becomes
(10)x′=a−bx−xy,y′=−y+xy,z′=y−bz,where the prime ′ denotes a derivative with respect to the new time t~ 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.Proposition 4.
All irreducible Darboux polynomials of system (10) consist of fx,y,z=x+y+1−bz−a/b with the cofactors k1x,y,z=−b and f2x,y,z=y with the cofactor k2x,y,z=x−1.Supposefx,y,z is a Darboux polynomial of system (10) with the cofactor kx,y,z, that is,
(11)a−bx−xy∂f∂x+xy−y∂f∂y+y−bz∂f∂z=kf.Comparing the degree of both sides of (11), we get that the cofactor kx,y,z is a polynomial with degree less than two. Without loss of generality, we set kx,y,z=k0+k1x+k2y+k3z. Then, we claim k3=0. In fact, we write f as powers in the variable z, i.e., fx,y,z=∑j=0nfjx,yzj, where each fj is a polynomial in the variables x and y and fn≠0. Computing the coefficient in (11) of zn+1, we have 0=k3fnx,y which implies k3=0.We make the change of the variablesx=α−1X, y=α−1Y, z=Z, and t~=αT and rewrite (10) into
(12)dXdT=−XY−αbX+α2a,dYdT=XY−αY,dZdT=Y−αbZ.Letn be the highest weight degree in the weight homogeneous components of f in x,y,z with the weight exponent 1,1,0. Set
(13)FX,Y,Z=αnfα−1X,α−1Y,Z=F0X,Y,Z+αF1X,Y,Z+⋯+αnFnX,Y,Z,KX,Y,Z=αkα−1X,α−1Y,Z=k1X+k2Y+αk0,where Fi 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,
(14)−xy−αbx+α2a∑i=0nαi∂Fi∂x+xy−αy∑i=0nαi∂Fi∂y+y−αbz∑i=0nαi∂Fi∂z=k1x+k2y+αk0∑i=0nFi,where we still use x,y,z instead of X,Y,Z.Equating the terms withα0 in (14) yields
(15)LF0=k1x+k2yF0,where L stands for a linear partial differential operator
(16)L≔−xy∂∂x+xy∂∂y+y∂∂z.To solve (15), we introduce the change of variables
(17)u=x,v=x+y,w=xexpzand transform (15) into
(18)uu−v∂F0¯∂u=k1−k2u+k2vF0¯,where F0¯ is F0 written in terms of u,v,w. Solving (18) yields F0¯=G0v,wv−uk1u−k2 with G0 being an arbitrary smooth function in v,w. Clearly, in order for F0x,y,z=F0¯u,v,w=G0x+y,xexpzyk1x−k2 to be a weight homogeneous polynomial of degree n in x,y,z, we must have k1=m, k2=−s and
(19)F0=a0xsymx+yl,s+m+l=n,for some nonzero number c0∈ℝ/0 and nonnegative integers s,m,l∈ℕ∪0.Similarly, equating the terms withα in (14) leads to
(20)LF1=mxF1−syF1+bx∂xF0+y∂yF0+bz∂zF0+k0F0.SubstitutingF0=a0xsymx+yl into (20) yields
(21)LF1=mx−syF1+a0xsymx+y1−1bs+bl+m+k0x+bs+l+m+k0y.Proceeding as above, by (17) we rewrite (21) as
(22)uu−v∂F1¯∂u=mu+su−vF1¯+a0vl−1usv−umbs+bl+m+k0u+bs+l+m+k0v−u.Integrating the above equation with respect tou, we obtain
(23)F1=F1¯=a0usv−umvl−1bs+bl+m+k0lnu−v+bs+l+m+k0lnwu+G1v,w=a0xsymx+yl−1bs+bl+m+k0lny+bs+l+m+k0z+G1x+y,xexpz,with G1 being an arbitrary smooth function. In order for F1 to be a weight-homogeneous polynomial of weight degree n−1, we must have G1=c1∈ℝ which is a constant, k0=−bs−bl−m, and
(24)F1=a0xsymx+yl−1l1−bz+c1.Equating the terms withαj in (14) for j=2,3,⋯,n, we have
(25)LFj=mx−syFj+bx∂xFj−1+y∂yFj−1+bz∂zFj−1+k0Fj−1−a∂xFj−2.In particular, forj=2 we substitute F0,F1 into (25) and obtain
(26)LF2=mx−syF2+a0xsymx+yl−2l1−bz+c1−blx+y+bl−1x+l−1y+bl1−bx+yz−asx+y2x−alx+yor equivalently
(27)uu−v∂F2¯∂u=mu+su−vF2¯+a0usv−umvl−2−lnwu1−b2ll−1u−v−bc1+alu+c1bl+c1−c1l+alu−v−asv2u.Solving the above linear differential equation, we get(28)F2=F2¯=a0xsymx+yl−2ll−11−b2z22−bc1+al+asln∣y∣−c1bl+c1−c1l+al+asz−asx+yx+G1x+y,xexpz.In order forF2 to be a weight-homogeneous polynomial of weight degree n−1, we have bc1+al+as=0, G2=c2∈ℝ which implies
(29)F1=a0xsymx+yl−1l11−bz−l1ab−asb,F2=a0xsymx+yl−2l21−b2z2−l22a1−bzb+asl−1b−1zb−x+yx+c2.Proceeding as above, from (25) with j=3 and (29), we get
(30)LF3=mx−syF3+a0xsymx+yl−3bzx+yll−11−b2z−al+sl−11−bb+ll−11−b2z22+al+sl−1b−1zb+c2−2bx−bl+2−ly+as2b+2m+al+sb+l−1x+y2x+a2l−1l+sb−asl−1x+y−asx+y2l1−bzx−all−1x+y1−bz.Working in a similar way to solveF2, we can prove that
(31)F3=a0xsymx+yl−3M1z3+M2z2+M3+M4x+M4yz+M5+M6zlny+M7dilog−yx+M8x+yx+G3x+y,xexpz,where dilog· is the dilogarithm function defined by
(32)dilogx=∫1xlnt1−tdt,and the coefficients Mi,i=1,⋯,8 are given by
(33)M1=l31−b3,M2=−l33a1−b2b−as1−b2l−1,M3=as2b+2m+a2l−1+sb−c2bl+2c2−c2l−a2ll−1b,M4=as1−bl,M5=a2ll−1b−2bc2+asa2l+s−1b+2b+2m,M6=−as1−b2l−1,M7=−as1−b2l−1,M8=as2b+2m+al+sb+2l−1−bl.In order forF3 to be a weight-homogeneous polynomial of weight degree n−3, we have M5=M6=M7=0 and the function G3 is a constant, which implies
(34)s=0,c2=l2a2b2,F1=a0ymx+yl−1l11−bz−ab,F2=a0ymx+yl−2l21−bz−ab2,F3=a0ymx+yl−3l31−b3z3−3a1−b2z2b+3a21−bzb2+c3,with c3 being a constant. Using the same argument similar to that above, we solve F4 and get
(35)c3=−ab3,F4=a0ymx+yl−4l41−b4z4−4a1−b3z3b+6a21−b2z2b2−4a31−bzb3+c4,with c4∈ℝ. Furthermore, by recursive calculations, we can prove that
(36)Fj=a0ymx+yl−jlj1−bz−abj,j=0,⋯,n.Therefore, the Darboux polynomial of system (10)
(37)f=Fα=1=∑j=0na0ljymx+yl−j1−bz−abj=a0ym1−bz−abl,and the cofactor k=Kα=1=mx−b+m, which completes the proof.


## 2.2. Proof of Corollary2

### 2.2.1. Proof of Statement (1) of Corollary2

It follows from Theorem1 and the fact that a function f is a polynomial first integral if and only if it is a Darboux polynomial with the zero cofactor.


### 2.2.2. Proof of Statement (2) of Corollary2

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.


### 2.2.3. Proof of Statement (3) of Corollary2

It follows from statements (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.


## 2.2.1. Proof of Statement (1) of Corollary2

It follows from Theorem1 and the fact that a function f is a polynomial first integral if and only if it is a Darboux polynomial with the zero cofactor.


## 2.2.2. Proof of Statement (2) of Corollary2

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.


## 2.2.3. Proof of Statement (3) of Corollary2

It follows from statements (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.


## 3. Proof of Theorem3

### 3.1. Proof of Statement (a) of Theorem3

According to Theorem1, system (1) has a Darboux polynomial F1=S+I+R−N with the cofactor −μ<0, which implies that IS,I,R,t=S+I+R−Nexpμt is a time-dependent first integral of system (1), that is, dI/dt=0 along the flow of (1). Let ϕt=St,It,Rt be an arbitrary orbit of system (1) with the initial value ϕ0=S0,I0,R0 satisfying F1S0,I0,R0≠0. By
(38)St+It+Rt−Nexpμt=S0+I0+R0−N∈ℝ/0,we see that ϕt will approach the invariant surface F1=0 as t⟶+∞, whereas it goes to infinity as t⟶−∞. Statement (a) follows.


### 3.2. Proof of Statement (b) of Theorem3

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 Ui and Vi, i=1,2,3. See [36] for more details on the Poincaré compactification of a polynomial vector field.


#### 3.2.1. In the Local ChartU1 and V1

Making the change of variablesS,I,R=z3−1,z1z3−1,z2z3−1 and the time rescaling dτ=z3−1dt, we obtain the Poincaré compactification of system (1) in the local chart U1(39)dz1dτ=−z1Nμz32−βz1+γz3−β,dz2dτ=−Nμz2z32+βz1z2+γz1z3,dz3dτ=−z3Nμz32−βz1−μz3.The planez3=0 at infinity is invariant, which corresponds to the points on the sphere at infinity, and so system (1) restricted to z3=0 becomes
(40)dz1dτ=βz1z1+1,dz2dτ=βz1z2,which has a line of singular points z1=0 and an isolated singular point z1,z2=−1,0. Moreover, system (40) has a first integral Φ=z1+1/z2. Using this first integral, the phase portrait on the local chart U1 is described in Figure 3. The flow on the local chart V1 is the same with the flow of U1 by reversing the time since the compactified vector field in V1 coincides with the vector field in U1 multiplied by −1 [36].Figure 3
The phase portrait in the local chartU1 at infinity.


#### 3.2.2. In the Local ChartU2 and V2

To get the Poincaré compactification of system (1) in the local chart U2, we take the change of variables S,I,R=z1z3−1,z3−1,z2z3−1 and the time rescaling dτ=z3−1dt and transform (1) into
(41)dz1dτ=Nμz32−βz12+γz1z3−βz1,dz2dτ=−βz1z2+γz2z3+γz3,dz3dτ=z3−βz1+γz3+μz3.Whenz3=0, system (41) becomes
(42)dz1dτ=−βz1z1+1,dz2dτ=−βz1z2.Clearly, system (42) coincides with system (40) by reversing the time. Then, the phase portrait on the local chart U1 is described in Figure 4. Again, the flow in the local chart V2 is the same as the flow in the local chart U2 by reversing the time.Figure 4
The phase portrait in the local chartU2 at infinity.


#### 3.2.3. In the Local ChartU3 and V3

Proceeding as above, making the change of variablesS,I,R=z1z3−1,z2z3−1,z3−1 and the time rescaling dτ=z3−1dt, system (1) becomes
(43)dz1dτ=Nμz32−γz1z2z3−βz1z2,dz2dτ=−z2γz2z3−βz1+γz3,dz3dτ=−z32γz2−μ.On the invariant planez3=0, system (43) is reduced to
(44)dz1dτ=−βz1z2,dz2dτ=βz1z2.System (44) has two lines of singular points z1=0 and z2=0. It also has a first integral Φ=z1+z2, which helps us get its phase portrait as shown in Figure 5. The flow at infinity in the local chart V3 is the same as the flow on the local chart U3 by reversing the time.Figure 5
The phase portrait in the local chartU3 at infinity.To summarize the above analysis, one obtains a global picture of the dynamical behavior of system (1) on sphere S2 at infinity: it has two closed curves of equilibria x=0,y2+z2=1∪y=0,x2+z2=1 and two nodes (see Figure 1).


### 3.3. Proof of Statement (c) of Theorem3

System (1) has the invariant algebraic surface S,I,R:S+I+R−N=0. Restricted to this invariant algebraic surface, system (1) becomes
(45)Ṡ=−βSI+μN−S≔PS,I,İ=βSI−γ+μI≔QS,I.We first study the finite singular points of system (45). Clearly, system (45) has two singular points
(46)P0=N,0,P1=γ+μβ,μNβ−γ−μβγ+μ.One can check that the eigenvalues of the linear part atP0 are −μ and Nβ−γ−μ. Set R0=Nβ/γ+μ, called the basic reproduction number. If R0<1, P0 is a stable node. If R0>1, P0 is a saddle. The eigenvalues of the linear part at P1 are
(47)−μβN+Δ2μ+γ,−μβN−Δ2μ+γ,with Δ=μβN−2γ+μ22−4γγ+μ3. If R0<1, P1 is a saddle. If R0>1 and Δ<0, P1 is a stable focus. If R0>1 and Δ≥0, P1 is a stable node. In particular, if R0=1, P1 coincides with P0, which is a semihyperbolic saddle-node point.Next, we turn to study the infinite singularities by the Poincaré compactification inℝ2. In the local chart U1, where S=z2−1, I=z1z2−1, and dτ=z2−1dt, we have
(48)dz1dτ=z1β−γz2+βz1−Nμz22,dz2dτ=z2βz1+μz2−Nμz22.It has two singularities0,0 and −1,0 on the invariant line z2=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=z1z2−1 and I=z2−1 and the time rescaling dτ=z2−1dt to obtain the Poincaré compactification of system (45) in the local chart U2(49)dz1dτ=Nμz22−βz12+γz1z2−βz1dz2dτ=z2μz2−βz1+γz2.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 S,I:I>0, observing
(50)∂P/I∂S+∂Q/I∂I=−β−μI<0,it follows from the Dulac Theorem that system (45) has no limit cycles. In the region S,I:I<0, 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 F1=0 is topologically described by one of Figure 2 on the Poincaré disk.


## 3.1. Proof of Statement (a) of Theorem3

According to Theorem1, system (1) has a Darboux polynomial F1=S+I+R−N with the cofactor −μ<0, which implies that IS,I,R,t=S+I+R−Nexpμt is a time-dependent first integral of system (1), that is, dI/dt=0 along the flow of (1). Let ϕt=St,It,Rt be an arbitrary orbit of system (1) with the initial value ϕ0=S0,I0,R0 satisfying F1S0,I0,R0≠0. By
(38)St+It+Rt−Nexpμt=S0+I0+R0−N∈ℝ/0,we see that ϕt will approach the invariant surface F1=0 as t⟶+∞, whereas it goes to infinity as t⟶−∞. Statement (a) follows.


## 3.2. Proof of Statement (b) of Theorem3

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 Ui and Vi, i=1,2,3. See [36] for more details on the Poincaré compactification of a polynomial vector field.


### 3.2.1. In the Local ChartU1 and V1

Making the change of variablesS,I,R=z3−1,z1z3−1,z2z3−1 and the time rescaling dτ=z3−1dt, we obtain the Poincaré compactification of system (1) in the local chart U1(39)dz1dτ=−z1Nμz32−βz1+γz3−β,dz2dτ=−Nμz2z32+βz1z2+γz1z3,dz3dτ=−z3Nμz32−βz1−μz3.The planez3=0 at infinity is invariant, which corresponds to the points on the sphere at infinity, and so system (1) restricted to z3=0 becomes
(40)dz1dτ=βz1z1+1,dz2dτ=βz1z2,which has a line of singular points z1=0 and an isolated singular point z1,z2=−1,0. Moreover, system (40) has a first integral Φ=z1+1/z2. Using this first integral, the phase portrait on the local chart U1 is described in Figure 3. The flow on the local chart V1 is the same with the flow of U1 by reversing the time since the compactified vector field in V1 coincides with the vector field in U1 multiplied by −1 [36].Figure 3
The phase portrait in the local chartU1 at infinity.


### 3.2.2. In the Local ChartU2 and V2

To get the Poincaré compactification of system (1) in the local chart U2, we take the change of variables S,I,R=z1z3−1,z3−1,z2z3−1 and the time rescaling dτ=z3−1dt and transform (1) into
(41)dz1dτ=Nμz32−βz12+γz1z3−βz1,dz2dτ=−βz1z2+γz2z3+γz3,dz3dτ=z3−βz1+γz3+μz3.Whenz3=0, system (41) becomes
(42)dz1dτ=−βz1z1+1,dz2dτ=−βz1z2.Clearly, system (42) coincides with system (40) by reversing the time. Then, the phase portrait on the local chart U1 is described in Figure 4. Again, the flow in the local chart V2 is the same as the flow in the local chart U2 by reversing the time.Figure 4
The phase portrait in the local chartU2 at infinity.


### 3.2.3. In the Local ChartU3 and V3

Proceeding as above, making the change of variablesS,I,R=z1z3−1,z2z3−1,z3−1 and the time rescaling dτ=z3−1dt, system (1) becomes
(43)dz1dτ=Nμz32−γz1z2z3−βz1z2,dz2dτ=−z2γz2z3−βz1+γz3,dz3dτ=−z32γz2−μ.On the invariant planez3=0, system (43) is reduced to
(44)dz1dτ=−βz1z2,dz2dτ=βz1z2.System (44) has two lines of singular points z1=0 and z2=0. It also has a first integral Φ=z1+z2, which helps us get its phase portrait as shown in Figure 5. The flow at infinity in the local chart V3 is the same as the flow on the local chart U3 by reversing the time.Figure 5
The phase portrait in the local chartU3 at infinity.To summarize the above analysis, one obtains a global picture of the dynamical behavior of system (1) on sphere S2 at infinity: it has two closed curves of equilibria x=0,y2+z2=1∪y=0,x2+z2=1 and two nodes (see Figure 1).


## 3.2.1. In the Local ChartU1 and V1

Making the change of variablesS,I,R=z3−1,z1z3−1,z2z3−1 and the time rescaling dτ=z3−1dt, we obtain the Poincaré compactification of system (1) in the local chart U1(39)dz1dτ=−z1Nμz32−βz1+γz3−β,dz2dτ=−Nμz2z32+βz1z2+γz1z3,dz3dτ=−z3Nμz32−βz1−μz3.The planez3=0 at infinity is invariant, which corresponds to the points on the sphere at infinity, and so system (1) restricted to z3=0 becomes
(40)dz1dτ=βz1z1+1,dz2dτ=βz1z2,which has a line of singular points z1=0 and an isolated singular point z1,z2=−1,0. Moreover, system (40) has a first integral Φ=z1+1/z2. Using this first integral, the phase portrait on the local chart U1 is described in Figure 3. The flow on the local chart V1 is the same with the flow of U1 by reversing the time since the compactified vector field in V1 coincides with the vector field in U1 multiplied by −1 [36].Figure 3
The phase portrait in the local chartU1 at infinity.


## 3.2.2. In the Local ChartU2 and V2

To get the Poincaré compactification of system (1) in the local chart U2, we take the change of variables S,I,R=z1z3−1,z3−1,z2z3−1 and the time rescaling dτ=z3−1dt and transform (1) into
(41)dz1dτ=Nμz32−βz12+γz1z3−βz1,dz2dτ=−βz1z2+γz2z3+γz3,dz3dτ=z3−βz1+γz3+μz3.Whenz3=0, system (41) becomes
(42)dz1dτ=−βz1z1+1,dz2dτ=−βz1z2.Clearly, system (42) coincides with system (40) by reversing the time. Then, the phase portrait on the local chart U1 is described in Figure 4. Again, the flow in the local chart V2 is the same as the flow in the local chart U2 by reversing the time.Figure 4
The phase portrait in the local chartU2 at infinity.


## 3.2.3. In the Local ChartU3 and V3

Proceeding as above, making the change of variablesS,I,R=z1z3−1,z2z3−1,z3−1 and the time rescaling dτ=z3−1dt, system (1) becomes
(43)dz1dτ=Nμz32−γz1z2z3−βz1z2,dz2dτ=−z2γz2z3−βz1+γz3,dz3dτ=−z32γz2−μ.On the invariant planez3=0, system (43) is reduced to
(44)dz1dτ=−βz1z2,dz2dτ=βz1z2.System (44) has two lines of singular points z1=0 and z2=0. It also has a first integral Φ=z1+z2, which helps us get its phase portrait as shown in Figure 5. The flow at infinity in the local chart V3 is the same as the flow on the local chart U3 by reversing the time.Figure 5
The phase portrait in the local chartU3 at infinity.To summarize the above analysis, one obtains a global picture of the dynamical behavior of system (1) on sphere S2 at infinity: it has two closed curves of equilibria x=0,y2+z2=1∪y=0,x2+z2=1 and two nodes (see Figure 1).


## 3.3. Proof of Statement (c) of Theorem3

System (1) has the invariant algebraic surface S,I,R:S+I+R−N=0. Restricted to this invariant algebraic surface, system (1) becomes
(45)Ṡ=−βSI+μN−S≔PS,I,İ=βSI−γ+μI≔QS,I.We first study the finite singular points of system (45). Clearly, system (45) has two singular points
(46)P0=N,0,P1=γ+μβ,μNβ−γ−μβγ+μ.One can check that the eigenvalues of the linear part atP0 are −μ and Nβ−γ−μ. Set R0=Nβ/γ+μ, called the basic reproduction number. If R0<1, P0 is a stable node. If R0>1, P0 is a saddle. The eigenvalues of the linear part at P1 are
(47)−μβN+Δ2μ+γ,−μβN−Δ2μ+γ,with Δ=μβN−2γ+μ22−4γγ+μ3. If R0<1, P1 is a saddle. If R0>1 and Δ<0, P1 is a stable focus. If R0>1 and Δ≥0, P1 is a stable node. In particular, if R0=1, P1 coincides with P0, which is a semihyperbolic saddle-node point.Next, we turn to study the infinite singularities by the Poincaré compactification inℝ2. In the local chart U1, where S=z2−1, I=z1z2−1, and dτ=z2−1dt, we have
(48)dz1dτ=z1β−γz2+βz1−Nμz22,dz2dτ=z2βz1+μz2−Nμz22.It has two singularities0,0 and −1,0 on the invariant line z2=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=z1z2−1 and I=z2−1 and the time rescaling dτ=z2−1dt to obtain the Poincaré compactification of system (45) in the local chart U2(49)dz1dτ=Nμz22−βz12+γz1z2−βz1dz2dτ=z2μz2−βz1+γz2.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 S,I:I>0, observing
(50)∂P/I∂S+∂Q/I∂I=−β−μI<0,it follows from the Dulac Theorem that system (45) has no limit cycles. In the region S,I:I<0, 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 F1=0 is topologically described by one of Figure 2 on the Poincaré disk.


## 4. 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 P0 on F1=0 or the endemic steady state P1 on F1=0, depending on the basic reproduction number R0. In case the R0 is small, saying less than one, the SIR model has no positive equilibrium point and the disease approaches the disease-free steady state P0 on F1=0, which implies the disease will always fail to spread. In case the R0 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 P1 on F1=0. In case R0 is equal to one, all the solutions in ℝ+3 tend to P0 on F1=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.



---

*Source: 5869275-2020-07-06.xml*