A reaction-diffusion system with cross-diffusion: Lie symmetry, exact solutions and their applications in the pandemic modeling
aa r X i v : . [ n li n . PS ] D ec A reaction-diffusion system with cross-diffusion:Lie symmetry, exact solutions and their applications inthe pandemic modeling
Roman Cherniha and Vasyl’ Davydovych Institute of Mathematics, National Academy of Sciences of Ukraine,3, Tereshchenkivs’ka Street, Kyiv 01004, Ukraine
Abstract
A nonlinear reaction-diffusion system with cross-diffusion describing the COVID-19outbreak is studied using the Lie symmetry method. A complete Lie symmetry classifica-tion is derived and it is shown that the system with correctly-specified parameters admitshighly nontrivial Lie symmetry operators, which do not occur for all known reaction-diffusion systems. The symmetries obtained are also applied for finding exact solutions ofthe system in the most interesting case from applicability point of view. It is shown thatthe exact solutions derived possess all necessary properties for describing the pandemicspread under 1D approximation in space and lead to the distributions, which qualitativelycorrespond to the measured data of the COVID-19 spread in Ukraine.
Keywords: reaction-diffusion system; cross-diffusion; Lie symmetry; exact solution; model-ing pandemic spread.
The outbreak of the coronavirus called COVID-19 in China has attracted extensive attentionof many mathematicians working in mathematical modeling. The first papers were alreadypublished in February–April 2020 (see, e.g., [1–7]. At the present time, the COVID-19 outbreakis already spread over the world as a pandemic. There were 65.5 mln. coronavirus cases andalmost 1.5 mln. deaths caused by this coronavirus up to date December 2 [8].Nowadays, there are many mathematical models used to describe epidemic processes andthey can be found in any book devoted to mathematical models in biology and medicine (see,e.g., [9–12] and papers cited therein). The paper [13] is one of the first papers in this direction.The authors created a model based on three ordinary differential equations (ODEs), whichnowadays is called the SIR model. There are several generalizations of the SIR model and theSEIR model (see the pioneering works [14, 15]), which involves four ODEs, is the most commonamong them. These two models are mostly used for numerical simulations in mathematicalmodeling the COVID-19 outbreak (see, e.g., [2, 4, 5]). Corresponding author. E-mail: [email protected] E-mail:[email protected]
1n the other hand, one may note that the spread of many epidemic processes, includingthe COVID-19 pandemic, is often highly non-homogenous in space. This fact can be takenin different ways but the most common approach consists in dividing the large domain (saya country) into many small sub-domains (regions of the country) and to apply the standardmodels based on ODEs to each sub-domain. However, there is another way – to use thereaction diffusion equations in order to model the spread of the infected population as a diffusionprocess [16, 17] (see also earlier papers cited therein). A possible model was also suggested inour previous work [18]. The model has the form u t = d ∆ u + u ( a − bu γ ) ,v t = d ∆ u + k ( t ) u, (1)where the lower subscript t means differentiation with respect to (w.r.t.) this variable, ∆ is theLaplace operator, u = u ( t, x, y ) and v = v ( t, x, y ) are two unknown functions, k ( t ) is the givensmooth positive function, d and d are diffusivities.The function u ( t, x, y ) describes the density (rate) of the infected persons (the number of theCOVID-19 cases) in a vicinity of the point ( x, y ), while v ( t, x, y ) means the density of the deathsfrom COVID-19. The diffusivity coefficients d and d describe the random movement of theinfected persons, which lead to increasing the pandemic spread. Formally speaking, one maytake d = d . However we believe that d > d because the movement of the infected personsleads firstly to higher rate of new COVID-19 cases but only some of them cause new deaths.Each coefficient in the reactions terms, a, b, γ and k ( t ), has the clear meaning described andverified in [18] (see Pages 2 and 3 therein).Of course, this model is an essential simplification because many factors causing the spreadof COVID-19 are not taking into account. In paper [16], for example, the authors constructthe diffusion model, which is essentially based on the SEIR model. As a result, their modelconsists of five PDEs, which can be analyzed only using numerical methods. Our idea was toconstruct a simpler model, which can be solved using analytical approaches, in particular, theLie symmetry method [19–21] and to show its applicability for the spread of the coronaviruspandemic. It is interesting to note that equations (3) and (5) from the model developed in [16]under the natural assumptions produce equations with a similar structure to those in (1). Infact, the density of infected population is proportional to that of the total living population andis proportional to the exposed population density (generally speaking, the relevant coefficientsare some functions but we keep constants). Having such assumptions, one arrives at the system(in our notations) u t = ∇ ( d u ∇ u ) + u ( a − bu ) ,v t = k ( t ) u, (2)The difference between (1) with γ = 1 and (2) consists only in the diffusion terms. In ourmodel, the diffusivity is taking to be a constant, while one is a linear function in paper [16].Notably, the diffusivity is a time-dependent function in paper [17].2he remainder of this paper is organized as follows. In Section 2, a complete Lie symmetryclassification (LSC) (group classification) of system (1) is derived. In particular, we have provedthat there are systems with correctly-specified parameters d , d , a and γ when the systemin question admits highly nontrivial Lie symmetry, which have no analogs for other knownreaction-diffusion systems. The results were obtained using the Lie–Ovsiannikov method [19],which is a combination of the classical Lie method and the technique for finding equivalencetransformations (ETs). The modern description of this method, its extension and applicationscan be found in [20] (Chapter 2).In Section 3, exact solutions of the specified system of the form (1) are constructed usingits Lie symmetry operators. In particular, the traveling wave type solution is derived and itsapplicability is extensively discussed. It is shown that this exact solution describes adequatelythe spread of the coronavirus pandemic provided the 1D approximation of space is assumed.Finally, we discuss the main results of the paper in the last section. In this section, it is identified that the basic system (1) for the pandemic modeling possesses avery reach Lie symmetry depending on the parameters γ, d and d and the function k ( t ). Firstof all, we note that for the LSC we need only the restrictions d + d = 0 (otherwise the systemin question degenerates into the ODE system, which was solved in [18]), k = 0 (otherwise thesystem in question seems to be useless for applications) and b = 0 and γ = 0 , − Theorem 1
The group of the continues ETs transforming system (1) to that with the samestructure, i.e., u t = d ∆ u + u ( a − bu γ ) ,v t = d ∆ u + k ( t ) u, (3) is the infinite-parameter Lie group generated by the transformations t = β t + α , x = β ( x cos β + y sin β ) + α , y = β ( y cos β − x sin β ) + α ,u = β u, v = β v + f ( x, y ) ,k = β β β k, d = β β d , d = β β β β d , a = aβ , b = bβ β γ , γ = γ. (4) Here α i ( i = 0 , , and β j ( j = 0 , . . . , are the real group parameters with the restrictions β β β = 0 , β > , and f ( x, y ) is an arbitrary smooth function. emark 1 In order to obtain system (3) with the nonnegative parameters (this is the biolog-ically motivated requirement explained above), the additional restrictions β > and β > should take place.Sketch of Proof of Theorem 1 is based on the known technique for constructing the groupof ETs. It is nothing else but a modification of the classical Lie method. In the case of system(1), one should start from the infinitesimal operator E = ξ ( t, x, u, v ) ∂ t + ξ ( t, x, u, v ) ∂ x + η ( t, x, u, v ) ∂ u + η ( t, x, u, v ) ∂ v + ζ ( t, x, u, v, k ) ∂ k + µ ∂ d + µ ∂ d + µ ∂ a + µ ∂ b + µ ∂ γ (5)being ξ , ξ , η , η and ζ to-be-determined functions, while µ i ( i = 1 , . . . ,
5) to-be-determinedconstants. The operator E involves the additional terms with the coefficients µ i and ζ , because d , d , a, b, γ and k ( t ) should be treated as a new variables.In order to find the operator E , we should apply Lie’s invariance criteria to the system ofequations consisting of (1) and a set of differential consequences of k ( t ) w.r.t. the variables t, x, u, v . Of course, each consequence is equal to zero, excepting ∂k∂t = k ′ ( t ) (the latter isnot useful because it is identity). As a result, we obtain a multicomponent system consistingequations from (1) and primitive equations like ∂k∂x = 0. Applying to this system Lie’s invariancecriteria, i.e., the second prolongation of the infinitesimal operator E , for deriving the system ofdetermining equations, the coefficients ξ , ξ , η , η , ζ and µ i ( i = 1 , . . . ,
5) were found. Theyhave the form ξ = C t + C , ξ = C x + C y + C , ξ = C y − C x + C ,η = C u, η = C v + h ( x, y ) , ζ = ( C − C − C ) k, µ = (2 C − C ) d ,µ = (2 C − C + C − C ) d , µ = − C a, µ = − ( C + C γ ) b, µ = 0 , (6)where C i , i = 1 , . . . , h ( x, y ) is an arbitrary smooth function. Theoperator (5) with the coefficients (6) generates the Lie group (4). The sketch of the proof is now completed.
In order to provide a complete LSC of system (1), one should to identify the principalalgebra of invariance (see definition, for example, in [20], page 23) from the very beginning.In fact, system (1) involves an arbitrary function k and several parameters (some of them canvanish). Thus, it should be considered as a class of systems of partial differential equations(PDEs), if one is going to provide a rigorous LSC. Theorem 2
The principal algebra of invariance of system (1) is infinite-dimensional Lie al-gebra generated by the operators: ∂ x , ∂ y , y ∂ x − x ∂ y , F ( x, y ) ∂ v , (7) where F ( x, y ) is an arbitrary smooth function. d = 0 and d = 0, leading to absolutelydifferent results. Theorem 3
System (1) with d = 0 admits the extension of the principal algebra (7) only infive cases. These cases and the corresponding Lie symmetry operators are as follows1) k ( t ) = 1 : ∂ t ;2) a = 0 , k ( t ) = t : t ∂ t + x ∂ x + y ∂ y − γ ( u∂ u + v∂ v ) ;3) a = 0 , d = 0 , k ( t ) = t p , p = − , : t∂ t + x ∂ x + y ∂ y − γ ( u∂ u + (1 − γp − γ ) v∂ v ) ; d = 0 , k ( t ) = e pt , p = 0 : ∂ t + pv∂ v ;5) d = 0 , a = 0 , k ( t ) = 1 : ∂ t , t∂ t + x ∂ x + y ∂ y − γ ( u∂ u + (1 − γ ) v∂ v ) .Here p is arbitrary constant.Any other system (1) with d = 0 admitting an extension of the principal algebra (7) isreduced by an ET from (4) to one of the listed in cases 1)–5). Remark 2
Using the simple ET from (4), one can set k k ( t ) ( k is an arbitrary constant)instead of k ( t ) in each case of Theorem 3 without any changes in Lie symmetry operators. Itis useful from the applicability point of view. Theorem 4
System (1) with d = 0 (then automatically d = 0 ) admits the extension of theprincipal algebra (7) only in four cases. In each case the additional operators have the structure X = ξ ( t, u ) ∂ t + ξ ( x, y ) ∂ x + ξ ( x, y ) ∂ y + η ( t, u ) ∂ u + ( G ( t, x, y, u ) + ( ξ t − ξ x + η u ) v ) ∂ v , (8) where the functions ξ , ξ form an arbitrary solution of the famous Cauchy-Riemann system ξ x = ξ y , ξ y = − ξ x , (9) and G is an arbitrary solution of the linear first-order PDE G t + u ( a − bu γ ) G u = k ( t ) (cid:0) η + 2 uξ x − uη u (cid:1) + uk ′ ( t ) ξ + u k ( t ) ( a − bu γ ) ξ u . (10) In the operator X , the functions ξ and η depending on the parameters γ and a have theforms : ) if γ = 1 and a = 0 then ξ = α + α e at au ( a − aγ + bγu γ ) ( a − bu γ ) − γ , η = α e at ( a − bu γ ) γ ;
2) if γ = 1 and a = 0 then ξ = α − γα t, η = α u ;
3) if γ = 1 and a = 0 then ξ = e at Z f (cid:0) a − buu e at (cid:1) u du + g ( t ) ,η = uae at (cid:16) f ( ω ) + ω Z e − at (cid:0) ω + be at (cid:1) f (cid:0) ωe − at (cid:1) dt (cid:17) , ω = a − buu e at ;
4) if γ = 1 and a = 0 then ξ = Z f (cid:0) u − bt (cid:1) u du + g ( t ) ,η = u (cid:18) f ( ω ) + b Z ( ω + 2 bt ) f ( ω + bt ) dt (cid:19) , ω = 1 u − bt. Here α and α are arbitrary parameters while the functions f , f and g are such that theidentity ξ t = − η u should take place. Proofs of Theorems 3 and 4 are based on the technique, which is a combination of theclassical Lie method and the group of ETs. This technique is often called the Lie–Ovsiannikovmethod because L.V. Ovsyannikov was the first who applied such technique for solving the LSCproblem (group classification problem) for a class of the nonlinear heat equations [19]. Here wepresent the proof of Theorems 4 because that is more complicated comparing with the proof ofTheorem 3.
Proof of Theorem 4.
System (1) with d = 0 , d = 0 can be rewritten as u t = u ( a − bu γ ) ,v t = ∆ u + k ( t ) u. (11)using an appropriate ET from (4). As usually, we start from the most general form of a Liesymmetry operator X = ξ ( t, x, y, u, v ) ∂ t + ξ ( t, x, y, u, v ) ∂ x + ξ ( t, x, y, u, v ) ∂ y + η ( t, x, y, u, v ) ∂ u + η ( t, x, y, u, v ) ∂ v . (12)6n order to find all Lie symmetry operators of the form (12) of system (11) one should applythe following invariance criterion : X ( u t − u ( a − bu γ )) (cid:12)(cid:12)(cid:12) M = 0 ,X ( v t − ∆ u − k ( t ) u ) (cid:12)(cid:12)(cid:12) M = 0 , (13)where operators X and X are the first and second prolongations of the operator X and themanifold M is defined by the system of equations u t = u ( a − bu γ ) , v t = ∆ u + k ( t ) u,u tt = ( a − b ( γ + 1) u γ ) u t , u tx = ( a − b ( γ + 1) u γ ) u x , u ty = ( a − b ( γ + 1) u γ ) u y . It should be stressed, that the manifold M involves not only the equations of the system inquestion but also the first-order consequences of the first equation of (11). These consequencesguarantee a complete solving the LSC problem. Notably, such peculiarity does not occur forscalar PDEs but one was noted for some systems of PDEs involving equations of different order(see, e.g., the relevant discussion in [21] Section 1.2.5).Having the correctly-defined manifold M , the invariance criterion (13) after rather standardcalculations leads to the system of determining equations as follows ξ x = ξ y = ξ v = 0 , ξ t = ξ u = ξ v = 0 , ξ t = ξ u = ξ v = 0 , (14) η x = η y = η v = 0 , η v = ξ t − ξ x + η u , (15) ξ x = ξ y , ξ y = − ξ x , (16) η uu = 2 ( a − b (1 + γ ) u γ ) ξ u + u ( a − bu γ ) ξ uu , (17) η t = ( a − b (1 + γ ) u γ ) η + u ( a − bu γ ) ξ u + u ( a − bu γ ) (cid:0) ξ t − η u (cid:1) , (18) η t + u ( a − bu γ ) η u = k ( t ) (cid:0) η + 2 uξ x − uη u (cid:1) + uk ′ ( t ) ξ + u k ( t ) ( a − bu γ ) ξ u . (19)Equations (14)–(15) can be easily integrated, hence the general form of the infinitesimaloperator (12) can be specified as (8). Now substituting the function η = G ( t, x, y, u ) + (cid:0) ξ t − ξ x + η u (cid:1) v into equation (19) and splitting the equation obtained w.r.t. the variable v , we arrive at thelinear equation (10) for G ( t, x, y, u ) and the equation u ( a − bu γ ) ξ tu + ξ tt + u ( a − bu γ ) η uu + η tu = 0 . (20)Obviously, equations (16) coincide with (9).Thus, we need only to solve the overdetermined system of equations (17), (18) and (20)w.r.t. the functions ξ and η . This is a nontrivial task because the function ξ depends on the7ependent variable u in contrast to the standard situation for the systems of reaction-diffusionequations (see [24, 25] and the papers cited therein). Since ξ u = 0, we used the differentialconsequence of equations (17) and (18) w.r.t. the variables t and u . Differentiating equation(17) w.r.t. t , taking the second-order consequence of equation (18) w.r.t. u and making therelevant calculations, we were able to derive the simple relation ξ t = (1 − γ ) η − uη u u . (21)Substituting the derivatives ξ tu and ξ tt derived from (21) into equation (20) one arrives at theclassification equation(1 − γ ) (cid:16) u ( a − bu γ ) η u − u ( a − bu γ ) ξ u + ( aγ − a + bu γ ) η (cid:17) = 0 . Thus, the following two cases must be examined separetely : γ = 1 and γ = 1 . In the case γ = 1, one immediately obtains η u = u ( a − bu γ ) ξ u + ( a − aγ − bu γ ) η u ( a − bu γ ) . (22)Differentiating equation (22) w.r.t. the variable u and substituting the expression obtained for η uu into equation (17), one arrives at the equation ξ u = a ( γ − η u ( a − bu γ ) . (23)Equations (22) and (23) can be easily integrated. The general solutions are ξ = g ( t ) + f ( t ) au ( a − aγ + bγu γ ) ( a − bu γ ) − γ , η = ( a − bu γ ) γ f ( t ) , (24)if a = 0 and ξ = g ( t ) , η = f ( t ) u, (25)if a = 0. Here f ( t ) and g ( t ) are arbitrary smooth functions at the moment. In order to find thefunctions f ( t ) and g ( t ), one needs to substitute (24) and (25) into equation (18). As a result,cases and of Theorem 4 were identified.In the case γ = 1, equations (17), (18) and (21) take the forms ξ t = − η u ,ξ tu + u ( a − bu ) ξ uu + 2( a − bu ) ξ u = 0 ,η t + 2 u ( a − bu ) η u = ( a − bu ) η + u ( a − bu ) ξ u , (26)8hile equation (20) is satisfied identically. Integrating the last two equations of system (26),we arrive exactly at Cases (for a = 0) and ( a = 0) of Theorem 4. Notably, we used thetransformations u ∗ = a − buu e at (in the case a = 0) and u ∗ = u − bt (in the case a = 0) for solvingthe second equation of system (26).As one can note, the functions ξ and η involve arbitrary functions f , f and g in cases and of Theorem 4. However, they should be specified from the equation ξ t = − η u . At thefinal stage, both functions, ξ and η , should be inserted into (10). The equation obtained isan integrable first-order PDE and its general solution is easily constructed in an explicit formprovided the function k ( t ) is given. Thus, all the coefficients of operator (8) are identified. The proof is now complete.
Let us present examples of highly nontrivial Lie symmetries in cases and (see The-orem 4). Setting f = 1 into expressions arising in cases and and using the equation ξ t = − η u , one can specify ξ and η as follows ξ = (cid:18) b a − u (cid:19) e at + α ba e at + α a e − at + α ,η = − ( a + 2 bu )( a − bu ) a u e at + α ( a − bu ) a e at + α ua e − at , and ξ = b t α bt − α t + α − u ,η = ( α − α bt − b t ) u + b t + α in cases and , respectively (here α , α , and α are arbitrary constants).Now one should use the functions ξ and η for finding the function G from equation (10). Inorder to avoid cumbersome formulae we additionally set α = α = α = 0 (just for simplicity),take the particular solution ξ = y, ξ = − x of the Cauchy-Riemann system (9) and fix thefunction k ( t ) : k ( t ) = e − at (case ) and k ( t ) = t − (case ). As a result, one arrives at theLie symmetry operators X = (cid:16) b a − u (cid:17) e at ∂ t + y∂ x − x∂ y − ( a +2 bu )( a − bu ) a u e at ∂ u + (cid:16) b a ln u + bu − a au + H (cid:0) x, y, a − buu e at (cid:1)(cid:17) ∂ v (27)and X = (cid:16) b t − u (cid:17) ∂ t + y∂ x − x∂ y + b t (1 − btu ) ∂ u + (cid:16) b ln( tu ) + btu − t u + H (cid:0) x, y, bt − u (cid:1)(cid:17) ∂ v (28)in cases and , respectively. Here H is an arbitrary smooth function.9 Exact solutions and their interpretation
Theorems 3 and 4 allows us to reduce the basic system (1) to that of lower dimensionality.In fact, using the Lie symmetry operators (or their linear combinations) listed in Theorems 3and 4 one can reduce (1) to the corresponding (1 + 1)-dimensional system and the latter to anODE system. Here we examine only two cases in order to show that those lead to useful exactsolutions.First of all, one may simplify the nonlinear system (1) using the ETs (4) with the correctly-specified parameters t = at, x = r ad x, y = r ad y,u = (cid:18) ba (cid:19) /γ u, v = (cid:18) ba (cid:19) /γ v, to the form u t = ∆ u + u (1 − u γ ) ,v t = D ∆ u + a k (cid:0) ta (cid:1) u, D = d d . (29)Here and in what follows we preserve the old notations for all the variables. Example 1 . Let us apply the operator ∂ y from the principal algebra (7) for reduction of thebasic system (29). Obviously, this operator produces the trivial ansatz u = u ( t, x ) , v = v ( t, x ),so that we arrive at the system u t = u xx + u (1 − u γ ) ,v t = Du xx + a k (cid:0) ta (cid:1) u, D = d d . (30)This system is nothing else but the initial system under assumption that the distributionof the infected persons is one-dimensional in space (i.e., the diffusion w.r.t. the axis y is verysmall). In this case, the distribution of the total number of deaths will be also one-dimensional.Making the further plausible assumption d ≫ d , i.e., the space diffusion of the infectedpersons leads mostly to increasing the total number of the COVID-19 cases and not so muchto new deaths, we may put D = 0. Following our previous paper [18], we specify the function k ( t ) = k e − αt (hereafter k > , α > u t = u xx + u (1 − u γ ) ,v t = k a exp (cid:0) − αta (cid:1) u. (31)Using Theorem 3 (see case therein), one notes that system (31) admits the Lie symmetry ∂ t − αa v∂ v . So, taking the linear combination of this operator and the operator ∂ x X = c∂ x + ∂ t − αa v∂ v , c ∈ R ,
10e obtain the ansatz u = φ ( ω ) , v = exp (cid:18) − αta (cid:19) ψ ( ω ) , ω = x − ct. (32)Substituting (32) into (31), one arrives at the ODE system φ ′′ + cφ ′ + φ (1 − φ γ ) = 0 ,cψ ′ + αa ψ = − k a φ. (33)The first equation in (33) is the known second-order ODE, which arises in many applications(e.g., for study the Fisher equation and its natural generalizations [26]). The general solutionof this ODE can be presented only in parametric form (see, e.g., [27]), which is not useful forfurther analysis. However, an exact solution in the explicit form can be constructed for thecorrectly-specified parameter c = γ +4 √ γ +2) . To the best of our knowledge, this parameter andthe relevant solution was established in [28] for the first time (see more references in [29]). Asa result, we obtain the traveling front solution of the first equation in system (31): u ( t, x ) = A exp γ p γ + 2) ω !! − /γ , ω = ± x − γ + 4 p γ + 2) t, A > . (34)Notably (34) with A < x → − x , therefore we may put ω = x − γ +4 √ γ +2) t in what follows.Having the function u in the explicit form (34), one easily derives the function v from thesecond equation of system (31): v ( t, x ) = k a Z exp (cid:18) − αta (cid:19) A exp γ p γ + 2) ω !! − /γ dt + g ( x ) , (35)where g ( x ) is an arbitrary smooth function. The integral in right hand side of (35) cannot beexpressed in the terms of elementary functions for arbitrary parameters a, α and γ , thereforewe study below a particular case.In order to avoid cumbersome formulae, let us set γ = 1 in system (31). In this case, theabove exact solution takes the form u ( t, x ) = (cid:16) A exp (cid:16) √ ω (cid:17)(cid:17) − , ω = x − √ t,v ( t, x ) = k a R exp (cid:0) − αta (cid:1) (cid:16) A exp (cid:16) √ ω (cid:17)(cid:17) − dt + g ( x ) . (36)11 emark 3 The expression for u ( t, x ) in (36) presents the well-known traveling front of thefamous Fisher equation u t = u xx + u (1 − u ) , which was firstly identified in [30]. The integral in right hand side of (36) can be expressed in the terms of elementary functionsfor several values of the parameter αa . Taking αa = , for example, we obtain u ( t, x ) = (cid:16) A exp (cid:16) √ ω (cid:17)(cid:17) − , ω = x − √ t,v ( t, x ) = g ( x ) − k a exp (cid:16) − x √ (cid:17) (cid:16) A + exp (cid:16) − √ ω (cid:17)(cid:17) − . (37)Now we turn to a possible interpretation of solution (37). First of all, the functions u and v should be nonnegative for any t > x ∈ I (here I ⊂ R ) because they represent thedensities. Obviously, the functions u is always positive. It is easily seen that each function g ( x )satisfying the inequality g ( x ) ≥ k a exp (cid:18) − x √ (cid:19) (cid:18) A exp (cid:18) − x √ (cid:19)(cid:19) − guarantees also nonnegativity of v . In particular, one may take the function g ( x ) = 6 k a exp (cid:18) − x √ (cid:19) (cid:18) A exp (cid:18) − x √ (cid:19)(cid:19) − , (38)which guarantees that the zero density of the deaths in the initial time t = 0, i.e., v (0 , x ) = 0.Examining the space interval I = [ x , x ] , x < x , we can calculate the total number of theCOVID-19 cases and deaths on this interval as follows U ( t ) = R x x u ( t, x ) dx,V ( t ) = R x x v ( t, x ) dx. (39)So, substituting solution (37) into (39), we arrive at the formulae U ( t ) = ( x − x ) − √ h(cid:16) A exp (cid:16) x − √ t √ (cid:17)(cid:17) − − (cid:16) A exp (cid:16) x − √ t √ (cid:17)(cid:17) − + ln (cid:16) A exp (cid:16) x − √ t √ (cid:17)(cid:17) − ln (cid:16) A exp (cid:16) x − √ t √ (cid:17)(cid:17)i ,V ( t ) = R x x g ( x ) dx − √ k a e − t h ln (cid:16) A exp (cid:16) √ t − x √ (cid:17)(cid:17) − ln (cid:16) A exp (cid:16) √ t − x √ (cid:17)i . (40)Obviously the functions U ( t ) and V ( t ) are increasing and bounded, because( U, V ) → (cid:16) ( x − x ) , Z x x g ( x ) dx (cid:17) as t → + ∞ . γ = 1. The function u ( t, x ) (leftsurface) describes the density of the of the COVID-19 cases while the function v ( t, x ) (rightsurface) describes the density of deaths. The parameters are: k = 0 . , a = 0 . , A = 1.Moreover, taking the appropriate function g ( x ), we can guarantee that U (0) = U ≥ , V (0) = V ≥ . Thus, one may claim that the exact solution (37) possesses all necessary properties for thedescription of the distribution of the COVID-19 cases and the deaths from this virus in time-space (under 1D approximation).Examples of the exact solution (37) with the specified parameters and the correspondingfunctions (40) are presented in Fig. 1 and Fig. 2. Notably, the parameters a and k were takenapproximately the same as in [18]. It follows from Fig. 1 that the spread of the COVID-19cases in space has the form of a traveling wave and this coincides (at least qualitatively) withthe real situation in many countries. In Ukraine, for example, the COVID-19 pandemic startedin the western part and then spread to the central and eastern parts of Ukraine (the majorexception was only the capital Kyiv, in which the total number of the COVID-19 cases was highfrom the very beginning). The distribution of deaths in space has more complicated behavior(see the right plot in Fig. 1). On the other hand, it is easily seen from Fig. 1 and Fig. 2 that v ( t, x ) ≪ u ( t, x ) what is in agreement with the measured data in many countries [8]. Notably,the behavior of the function v ( t, x ) can be essentially changed by the appropriate choice of thefunction g ( x ). Example 2 . Let us apply the operator y∂ x − x ∂ y from the principal algebra (7) forreduction of the basic system (29). Obviously, this operator produces the well-known ansatz13igure 2: The functions (40) with (38). The function U ( t ) (left curve) shows the time-evolutionof total number of the COVID-19 cases on the space interval [0 , V ( t )(right curve) shows the time-evolution of total deaths. The parameters k , a and A are thesame as in Fig. 1. u = u ( t, r ) , v = v ( t, r ) , r = x + y , i.e., we examine the radially-symmetric case. In this case,we arrive at the system u t = 1 r ( ru r ) r + u (1 − u γ ) ,v t = D r ( ru r ) r + 1 a k (cid:18) ta (cid:19) u, D = d d . Making the same assumptions k ( t ) = k e − αt and d ≫ d , i.e., D = 0, as in Example 1, weobtain the system u t = r ( ru r ) r + u (1 − u γ ) ,v t = k a exp (cid:0) − αta (cid:1) u. (41)We have proved that system (41) again admits the Lie symmetry ∂ t − αa v∂ v (however, theoperator of the space translation ∂ r is absent in this case). So, using this symmetry, one obtainsthe ansatz u = φ ( r ) , v = exp (cid:18) − αta (cid:19) ψ ( r ) , r = p x + y . (42)14ubstituting (42) into (41), one arrives at the system φ ′′ + 1 r φ ′ + φ (1 − φ γ ) = 0 , − αa ψ = k a φ. In contrast to the first equation in (33), the first ODE in the above system is much morecomplicated. To the best of our knowledge, exact solutions of this equation are unknown.It can be only noted that φ = 1 is the steady-state solution. As a result, we arrive at thespace-homogeneous solution of the nonlinear system (41) u = 1 , v = − k α exp (cid:18) − αta (cid:19) . Of course, the restriction α < v means the density. Onthe other hand, it means the exponential growing of the function v what is rather unrealisticbecause one obtains the total extinction of the population in question for a finite time. The main part of this paper is devoted to the LSC of the class of reaction-diffusion systemwith the cross-diffusion (1). The system in question was suggested in work [18] as the naturalgeneralization of the mathematical model for describing the COVID-19 outbreak.Firstly, we present a statement about the group of ETs of system (1) (see Theorem 1) inorder to establish possible relations between systems that admit equivalent invariance algebras.Secondly, we find the principal algebra of system (1), i.e., the maximal invariance algebra of thissystem with arbitrary coefficients (see Theorem 2). And lastly, we present two main theorems(Theorems 3 and 4) describing reaction-diffusion systems of the form (1) admitting nontrivialLie symmetry, i.e., present the LSC of system (1). In Section 3, we demonstrate that the Liesymmetries identified in Section 2 are useful for finding exact solutions, which can describe thespread of the COVID-19 pandemic.From the mathematical point of view, the most interesting Lie symmetry operators of system(1) occur when d = 0 and are presented in Theorem 4. One sees that the coefficient ξ ofthe infinitesimal operator X (see(8)) depends on the variable u (excepting case . Moreover,this dependence is nonlinear. To the best of our knowledge, it is the first example of suchdependence for systems of evolution equations, in particular, reaction-diffusion systems. Weassume that such unusual Lie symmetry of system (1) with d = 0 can be a consequence of itsintegrability. In fact, one may consider the first equation as an ODE with the variables x and y as parameters. Solving this ODE, one obtains u ( t, x, y ) = a /γ e at C ( x, y ) (cid:16) a + b C γ ( x, y ) (cid:0) e aγt − (cid:1) (cid:17) − /γ , C ( x, y ) is an arbitrary function. Substituting this expression for u into the secondequation of the system, one again obtains the integrable ODE to find the function v .Finally, it should be pointed out that Lie symmetries operators, which are nonlinear w.r.t.unknown functions were recently identified for a simplification of the Shigesada–Kawasaki–Teramoto system in [25] (see Section 3). Such peculiarity of Lie symmetry also occurs fora special Schr¨odinger type equation [31], which can be rewritten in the form of a reaction-diffusion system with the cross-diffusion. However, the coefficient ξ (see operator (12)) inall known Lie symmetries of a wide range of reaction-diffusion systems [25, 31–36] (see morereferences in Chapter 2 of [37]) do not depend on the unknown function(s) in contrast tothose in Theorem 4 and examples (27)–(28). Moreover, we may conclude that the well-known‘people theorem’ stating that the coefficient ξ in each Lie symmetry of an arbitrary scalarevolution PDE of the order two and higher can depend only the time variable (no dependence onspace variables and/or dependent variable!), cannot be generalized on the systems of evolutionequations without additional restrictions. The problem how to define these restrictions is anopen question.From the applicability point of view, the most interesting system of the form (1) admittingnontrivial Lie symmetry is presented in case of Theorem 3. Here the function k ( t ) of system(1) has the form that can be useful for describing the COVID-19 outbreak [18]. Moreover, thediffusivity d = 0 as it is stated in [16]. In Section 3, we demonstrate how the Lie symmetriesobtained can be applied for constructing of exact solutions. Furthermore we prove that anexact solution (with correctly-specified parameters) possesses all necessary properties for thedescription of the distribution of the COVID-19 cases and the deaths from this virus in timeand space. Although it was done under 1D space approximation, this solution can be usefulfor the prediction of the COVID-19 pandemic if its spread has a favorite direction (a typicalexample is Ukraine). Of course, one needs to identify all the parameters in system (1) in orderto calculate correct numbers of the COVID-19 cases and make a plausible prediction but thislays beyond the scope of this work. References [1] X. Luo et al., Analysis of potential risk of COVID-19 infections in China based on apairwise epidemic model. (2020) doi:10.20944/preprints202002.0398.v1.[2] L. Peng et al., Epidemic analysis of COVID-19 in China by dynamical modeling.
ArXiv (2020) arXiv:2002.06563.[3] N. Shao et al., Dynamic models for coronavirus disease 2019 and data analysis.
Math.Meth. Appl. Sci. (2020) 4943–4949.[4] J. Tian et al., Modeling analysis of COVID-19 based on morbidity data in Anhui, China. MBE (2020) 2842–2852. 165] D. Efimov, U. Ushirobira, On interval prediction of COVID-19 development based on aSEIR epidemic model. Research report. Inria Lille Nord Europe–Laboratoire CRIStAL–Universite de: Lille, France (2020).[6] W.C. Roda, M.B. Varugheseb, D. Han, M.Y. Li, Why is it difficult to accurately predictthe COVID-19 epidemic? Infectious Disease Modelling (2020) 271–281.[7] R. Cherniha, V. Davydovych, A mathematical model for the COVID-19 outbreak. ArXiv (2020) arXiv:2004.01487v2.[8] Available online: (accessed on 2 December2020).[9] F. Brauer, C. Castillo-Chavez,
Mathematical Models in Population Biology and Epidemi-ology . Springer, New York, 2012.[10] M.J. Keeling, P. Rohani,
Modeling Infectious Diseases in Humans and Animals . PrincetonUniversity Press, Princeton, 2008.[11] O. Diekmann, J. Heesterbeek,
Mathematical Epidemiology of Infectious Diseases: ModelBuilding, Analysis and Interpretation . Chichester, John Wiley, 2000.[12] J.D. Murray,
Mathematical Biology II: Spatial Models and Biomedical Applications .Springer, Berlin, 2003.[13] W.O. Kermack, A.G. McKendrick, A contribution to the mathematical theory of epi-demics.
Proc. Roy. Soc. A (1927) 700–721.[14] R.M. Anderson, R.M. May, Directly transmitted infectious diseases: control by vaccina-tion.
Science (1982) 1053–1060.[15] K. Dietz, The incidence of infectious diseases under the influence of seasonal fluctuations.
Lecture Notes in Biomathematics , Springer, Berlin, 1976, pp. 1–15.[16] A. Viguerie et al., Simulating the spread of COVID-19 via a spatially-resolved susceptible-exposed-infected-recovered-deceased (SEIRD) model with heterogeneous diffusion. Appl.Math. Lett. (2021) 106617, 9 pp.[17] Y. Mammeri, A reaction-diffusion system to better comprehend the unlockdown: appli-cation of SEIR-type model with diffusion to the spatial spread of COVID-19 in France.
Comput. Math. Biophys. (2020) 102–113.[18] R. Cherniha, V. Davydovych, A mathematical model for the COVID-19 outbreak and itsapplications. Symmetry (2020), 990, 12 pp.1719] L.V. Ovsiannikov, The Group Analysis of Differential Equations . Academic Press, NewYork, 1980.[20] R. Cherniha, M. Serov, O. Pliukhin,
Nonlinear Reaction-Diffusion-Convection Equations:Lie and Conditional Symmetry, Exact Solutions and their Applications . Chapman andHall/CRC, New York, 2018.[21] G.W. Bluman, A.F. Cheviakov, S.C. Anco,
Applications of Symmetry Methods to PartialDifferential Equations . Springer, New York, 2010.[22] I.S. Akhatov, R.K. Gazizov, N.H. Ibragimov, Nonlocal symmetries. Heuristic approach.
J.Sov. Math. (1991) 1401–1450.[23] N.H. Ibragimov, M. Torrisi, A. Valenti, Preliminary group classification of equations v tt = f ( x, v x ) v xx + g ( x, v x ). J. Math. Phys. (1991) 2988–2995.[24] R. Cherniha, J.R. King, Lie symmetries and conservation laws of nonlinear multidimen-sional reaction-diffusion systems with variable diffusivities. IMA J. Appl. Math. (2006)391–408.[25] R. Cherniha, V. Davydovych, L. Muzyka, Lie symmetries of the Shigesada–Kawasaki–Teramoto system. Commun. Nonlinear Sci. Numer. Simulat. (2017) 81–92.[26] J.D. Murray, Mathematical Biology . Springer, Berlin, 1989.[27] A.D. Polyanin, V.F. Zaitsev,
Handbook of Ordinary Differential Equations for Scientistsand Engineers . CRC Press Company, Boca Raton, 2018.[28] M.A. Abdelkader, Travelling wave solutions for a generalized Fisher equation.
J. Math.Anal. Appl. (1982) 287–290.[29] B. Gilding, R. Kersner, Travelling Waves in Nonlinear Diffusion-Convection Reaction .Birkh¨auser, Basel, 2004.[30] M. Ablowitz, A. Zeppetella, Explicit solutions of Fisher’s equation for a special wave speed.
Bull. Math. Biol. (1979) 835–840.[31] W. Fushchych, R. Cherniha, V. Chopyk, On unique symmetry of two nonlinear general-izations of the Schr¨odinger equation. J. Nonlinear Math. Phys . (1996) 296–301.[32] A.G. Nikitin, Group classification of systems of non-linear reaction-diffusion equations. Ukrainian Math. Bull. (2005) 153–204.[33] R. Cherniha, H. Wilhelmsson, Symmetry and exact solution of heat-mass transfer equa-tions in thermonuclear plasma. Ukrainian Math. J. (1996) 1434–1449.1834] M. Torrisi, R. Tracina, A. Valenti, A group analysis approach for a nonlinear differentialsystem arising in diffusion phenomena. J. Math. Phys. (1996) 4758–4767.[35] J.M. Stewart, P. Broadbridge, J.M. Goard, Symmetry analysis and numerical modellingof invasion by malignant tumour tissue. Nonlinear Dynamics (2002) 175–193.[36] M.I. Serov, T.O. Karpaliuk, O.G. Pliukhin, I.V. Rassokha, Systems of reaction-convection-diffusion equations invariant under Galilean algebras. J. Math. Anal. Appl. (2015)185–211.[37] R. Cherniha, V. Davydovych, Nonlinear Reaction-Diffusion Systems — Conditional Sym-metry, Exact Solutions and their Applications in Biology,
Lecture Notes in Mathematics2196