Multi-Parameter Reaction–Diffusion Systems with Quadratic Nonlinearity and Delays: New Exact Solutions in Elementary Functions

: The study considers a nonlinear multi-parameter reaction–diffusion system of two Lotka– Volterra-type equations with several delays. It treats both cases of different diffusion coefﬁcients and identical diffusion coefﬁcients. The study describes a few different techniques to solve the system of interest, including (i) reduction to a single second-order linear ODE without delay, (ii) reduction to a system of three second-order ODEs without delay, (iii) reduction to a system of three ﬁrst-order ODEs with delay, (iv) reduction to a system of two second-order ODEs without delay and a linear Schrödinger-type PDE, and (v) reduction to a system of two ﬁrst-order ODEs with delay and a linear heat-type PDE. The study presents many new exact solutions to a Lotka–Volterra-type reaction–diffusion system with several arbitrary delay times, including over 50 solutions in terms of elementary functions. All of these are generalized or incomplete separable solutions that involve several free parameters (constants of integration). A special case is studied where a solution contains inﬁnitely many free parameters. Along with that, some new exact solutions are obtained for a simpler nonlinear reaction–diffusion system of PDEs without delays that represents a special case of the original multi-parameter delay system. Several generalizations to systems with variable coefﬁcients, systems with more complex nonlinearities, and hyperbolic type systems with delay are discussed. The solutions obtained can be used to model delay processes in biology, ecology, biochemistry and medicine and test approximate analytical and numerical methods for reaction–diffusion and other nonlinear PDEs with delays.

Linear and nonlinear differential equations, ordinary and partial, with delay are not uncommon for mathematical modeling of phenomena and processes in various areas of theoretical physics, mechanics, control theory, biology, biophysics, biochemistry, medicine, ecology, economics and technical applications.
Below are a few factors that lead to the need to introduce delay into mathematical models described by differential equations.In biology and biomechanics, the delay is due to the limited transmission rate of nerve and muscle reactions in living tissues.In medicine, in studying the spread of infectious diseases, the delay time is determined by the incubation period (the time from the moment of infection to the first signs when the disease begins to manifest itself).In the dynamics of populations, a delay arises because individuals participate in reproduction only after reaching a certain age.In control theory, delays are usually associated with the finite speed of signal propagation and the limited speed of technological processes.
The presence of delay in mathematical models and differential equations is a complicating factor, which, as a rule, leads to a narrowing of the stability region of the solutions obtained.The study and solution of ordinary differential equations (ODEs) with a delay are comparable in complexity to the study and solution of partial differential equations (PDE) without delays.Importantly, unlike simpler PDEs without delay, which often have self-similar solutions, PDEs with delay do not admit self-similar solutions.
The present article deals with a Lotka-Volterra-type reaction-diffusion system with two components and four constant delays described by two coupled PDEs with quadratic nonlinearities: where u = u(x, t) and v = v(x, t) are the unknown functions; ūi = u(x, t − τ i ), i = 1, 3; vj = v(x, t − τ j ), j = 2, 4; τ i ≥ 0 and τ j ≥ 0 are delay times; a 1 > 0 and a 2 > 0 are diffusion-type coefficients; and b 1 , b 2 , p 1 , p 2 , q 1 , and q 2 are some constants.
Systems of the form (1) are used in the mathematical modeling of various processes in biology, ecology, biochemistry, medicine, etc.For example, in population theory, the delays τ 1 and τ 4 characterize the mean reproductive age of species.At the same time, τ 2 and τ 3 are responsible for the time required for changes in the population size of one species to cause changes in the other.All the delays are non-negative and can be zero in some models.The terms with nonzero coefficients q 1 and p 2 make the model different from a single equation, while the coefficients are responsible for the interaction between individuals of the two populations.In the case of cooperation, when one species persists in the absence of the other or when species mutually increase each other's growth rate, the coefficients q 1 and p 2 are positive.In the case of competition, an increase in one population leads to a decrease in the other (for example, an increase in the number of predators leads to a decrease in the prey population), and the coefficients q 1 and q 2 are negative.Lotka-Volterra cooperative models with delays, described by PDEs (1), were treated in [1,2]; for competitive models, see [3][4][5].
The system of PDEs (1) without delays admits simple traveling wave solutions with the functions u(z) and v(z) described by a system of ODEs.The studies [20][21][22] obtained some exact solutions to this system, of the form (2), which are expressible in terms of elementary functions.
The system of delay PDEs (1) admits simple traveling wave solutions (2) (e.g., see [1,5,45]).At present, there is only one publication [46] where two other exact solutions to system (1) were obtained and some reductions of related nonlinear reaction-diffusion systems of PDEs were described.
1.2.The System of Delay PDEs in Question Exact Solutions and Generalized Separable Solutions 1 • .In general, the multi-parameter nonlinear Lotka-Volterra-type diffusion system of delay PDEs (1) contains 12 free parameters, four of which can be removed by scaling the dependent and independent variables.We will consider a special case of system (1) while imposing one constraint, p 1 /p 2 = q 1 /q 2 , on the free parameters.With this in mind, we set where ūi = u(x, t , τ i ≥ 0, and τ j ≥ 0 are free parameters; c 1 = 0 and c 2 = 0 are some constants to be determined later.The present study does not deal with the degenerate cases of k 1 = 0 or k 2 = 0, when system (3) becomes semi-coupled, meaning that one equation of the system is isolated, while the other one is linear in the unknown. 2 • .A solution is called exact if it satisfies the equation in question exactly (i.e., when substituted into the equation, the solution turns it into an identity).At the same time, any approximations or simplifications of the equation are not allowed, and no a priori assumptions are used.
The term 'exact solution' for nonlinear PDEs with delay is normally used in cases where the solution is expressed: (i) In terms of elementary functions; (ii) In terms of elementary and special functions and indefinite integrals; (iii) In terms of solutions to ODEs without or with delay or systems of such equations.
Combinations of solutions from items (i)-(iii) are allowed.In cases (i) and (ii), an exact solution can be represented in an explicit, implicit or parametric form.
Exact solutions play an essential role as standard reference results that can be used to verify the consistency and estimate errors of various numerical, asymptotic, and approximate analytical methods for solving nonlinear PDEs with delay.These solutions can also serve as a basis for perfecting and testing computer algebra software packages such as Maple and Mathematica.
The simplest and most preferred exact solutions for testing are solutions in terms of elementary functions (see Item (i)).Currently, the vast majority of the obtained exact solutions of such equations are more complex and relate to Item (iii).
3 • .We will be looking for exact generalized separable solutions to the system of delay PDEs (3) in the form or in an alternative form that can be obtained from (4) by replacing the functions ϕ(x) and ψ(x) with ϕ(t) and ψ(t), respectively.In some cases, we will look for more complicated solutions involving functions with two arguments.Special attention will be given to the construction of exact solutions expressible in terms of elementary functions.

Simplest Solutions (Stationary Points)
Constants (stationary points) are the simplest solutions to system (3): These constants are determined by the formulas The stationary points (6) will further be used to construct more complicated solutions to system (3).

Reductions and Exact Solutions of the System of Delay PDEs with Different
Diffusion Coefficients (a 1 = a 2 ) 2.1.Reduction of the System with Three Delays to a Single Second-Order Linear ODE without Delay The four delay times in system (3) are assumed to be related by the single constraint This suggests that any three of these delays can be set arbitrarily.Notably, relation (7) holds, for example, in the following three important special cases: We seek a generalized separable solution to system (3) in the form where u • and v • are the stationary points (5) of system (3).We assume that the function θ = θ(x) satisfies the second-order linear ODE The constants λ and µ appearing in (8) and (9) are to be determined in a subsequent analysis.
The functions (8) are selected so as to ensure that under condition (7), the following simple relations hold: Substituting ( 8) into (3) and taking into account relations (10) and Equation ( 9) followed by a simple rearrangement of terms, we arrive at the linear algebraic system It serves to determine the parameters λ and µ.If a 1 = a 2 , the solution of system (11) is expressed as For the first three stationary points in (6), the coefficients (12) are given by In case (d), the parameters degenerate, λ = µ = 0, resulting in solutions of little interest.
The general solution to the linear ODE (9) with µ = 0 reads where C 1 and C 2 are arbitrary constants.Formulas (8), ( 13) and ( 14) and the first three stationary points (6) determine six nondegenerate solutions, three for µ > 0 and µ < 0 each, to the nonlinear system (3) with a 1 = a 2 and four constant delays satisfying condition (7).The study [46] presents exact solutions to system (3) corresponding to the stationary point (b) in (6).These are written out below. Solution where C 1 and C 2 are arbitrary constants.
2.2.Reduction of the System of PDEs with Three Delays to a System of Two Second-Order ODEs without Delays We assume that the four delay times in system (3) are related by the single constraint (7).
Generalized separable solutions exponential in t.We seek a generalized separable solution to system (3) in the form with the functions θ = θ(x), ϕ = ϕ(x), and ψ = ψ(x) and parameter λ to be determined in the subsequent analysis.The functions (17) are selected so as to ensure that under condition (7), the following relations hold: Substituting ( 17) into (3) and taking into account (18), we obtain with the notation ρ = k 1 ϕ + k 2 ψ used for short.Relations (19) can be satisfied by setting The first two equations in (20) form a closed system of ODEs for determining the functions ϕ and ψ, while the other two equations make up an overdetermined system for the single function θ.On requiring that the last two equations in (20) coincide, we find the parameter λ and other constants: On substituting (21) into (20), we arrive at the following system of stationary secondorder ODEs: where be a solution of system (22).Then, the functions where A and B are arbitrary constants (B = 0) are also a solution of the system.
System (22) simplifies significantly if either ϕ or ψ is set to zero, in which case the system will only consist of two equations.
Properties and some transformations of the system of ODEs (22) with ψ = 0.The autonomous systems of ODEs (22) with ψ = 0 can be rewritten in the form where b = b 1 /a 1 , k = k 1 , and β = (b 2 − b 1 )/(a 2 − a 1 ).We assume that ϕ(x) ≡ const (the special case ϕ = u • = const was discussed above).Below are three linear transformations of the system of ODEs ( 24) and ( 25), which will be required for the construction of exact solutions.
1 • .The linear transformation takes system ( 24) and ( 25) to the same system with other determining parameters also takes system ( 24) and ( 25) to the same system with other parameters We will use (b, β) to refer to systems (24) and (25).Then, transformations G 1 and G 2 and their composition G 1 • G 2 connect this system with three other systems of the same form, which can schematically be depicted as These transformations allow one to multiply exact solutions to ( 24) and ( 25).In addition, the following statement holds true.

Statement 2. Suppose the functions
make up an exact solution to the system of ODEs ( 24) and (25).Then, exact solutions to the same system with other determining parameters b and β (k remains the same) can be obtained using three pairs of formulas specified in the last three rows of Table 1.
Table 1.Exact solutions to three related systems of ODEs of the form ( 24) and ( 25) that can be expressed in terms of solution (b, β) of the original system.

System of ODEs Function
reduces ( 24) and (25) to a system of the special form Exact solutions to the system of ODEs ( 24) and (25).The general solution to the autonomous ODE ( 24) can be represented in the implicit form where C 1 and C 2 are arbitrary constants.In general, the integral on the right-hand side of ( 34) is not expressible in terms of elementary functions.Solutions to ODE (24) are expressible in terms of the Weierstrass elliptic function ℘ = ℘(z, g 2 , g 3 ), which is defined explicitly via the elliptic integral [57]: where g 2 and g 3 are free parameters called invariants.The first ODE in (33) has a solution [57,58]), where C1 = g 3 is an arbitrary constant.By virtue of transformation H (32) and the first property in (23), Equation ( 24) with k = 0 admits the solutions where C1 and C2 are arbitrary constants.Notably, once ϕ(x) is known, the general solution to the linear homogeneous Equation ( 25) can be written as [58]: where θ 0 = θ 0 (x) is any nontrivial particular solution to Equation (25).
If β = b, Equation ( 25) with θ = ϕ coincides with Equation (24).Therefore, the function θ 0 = ϕ is a particular solution to Equation (25) in this case.On substituting θ 0 = ϕ into (36), we obtain the general solution to Equation ( 25) with β = b; the representation of the function ϕ in terms of the Weierstrass elliptic function (35) can also be used.

Exact Solutions in Terms of Elementary Functions to the Reduced System of Two Second-Order ODEs
Consider a few special cases where solutions to the system of ODEs ( 24) and ( 25) are expressible in terms of elementary functions.1 • .For C 1 = 0 and b = 0, the function ϕ in ( 34) can be represented in explicit form as Substituting ( 37) into (25) and multiplying by (x + C 2 ) 2 , we obtain the equation Up to renaming, it coincides with a special case of an ODE discussed in the book [58] (see Equation ( 132) with a = 0 and n = 2 on page 535).
For β = 0, the general solution of Equation ( 38) is where C 3 and C 4 are arbitrary constants.For β = 0, the general solution of Equation ( 38) can be represented as where C For C 1 = 0 and b = 0, the integral in (34) can be expressed in terms of elementary functions, and the function ϕ can be written explicitly as where √ ∓b C 2 is an arbitrary constant.For β = 0, the general solutions of Equation ( 25) corresponding to solutions (40) are and A 2 and A 3 are arbitrary constants.
For β = 0 and b < 0, Equation ( 25) admits particular solutions of the form [58]: Substituting θ 0 into Equation (25) and in view of the first solution (40), we find the following allowed values of the parameters: As a result, taking into account (36), we obtain four general solutions for different β: For β = 0 and b > 0, Equation (25) admits the following particular solutions [58]: In view of (36), we obtain four general solutions for different β: In this case, the integral can be expressed in terms of elementary functions, and the function ϕ can be represented in explicit form as For β = 0, the general solutions of Equation ( 25) corresponding to solutions (43) are where For β = 0 and b < 0, Equation ( 25) admits the following exact solutions [58]: The respective general solutions of Equation ( 25) are where ξ = 1 2 √ −b x + A 1 .For β = 0 and b > 0, Equation (25) has the following exact solutions [58]: The respective general solutions of Equation ( 25) are where 4 • .Apart from the above exact solutions for ϕ(x), Equation ( 24) admits the solutions The following replacements should be made in the respective solutions for θ(x): {sin, cos, tan, arctan} with {cos, sin, cot, arccot} and {sinh, cosh, tanh, artanh} with {cosh, sinh, coth, arcoth}, respectively.
In particular, besides the first solution in (40) with b < 0 and the first solution in (41) with β = b, there are solutions Moreover, apart from the first solution in (43) with b < 0 and the first solution in (44) with β = − 5 4 b, there are solutions Table 2 summarizes the exact solutions in terms of elementary functions to the system of ODEs (24) and (25) for various values of the parameters b and β.The space variable x can everywhere be replaced with x + A 1 , where A 1 is an arbitrary constant.

Reduction of a System of PDEs with Three Delays to a System of Three Second-Order ODEs without Delays
System (3) with four delays that satisfy the single relation ( 7) also admits generalized separable solutions linear in time t: An analysis similar to that above results in the following parameters of Equation ( 3): where σ is an arbitrary constant.In this case, the functions θ = θ(x), ϕ = ϕ(x), and ψ = ψ(x) are described by the stationary system of ODEs Let us look at the special case a 1 = a 2 = a.By adding up the first two equations of (49) multiplied by k 1 and k 2 , we obtain an isolated ODE for ρ = k 1 ϕ + k 2 ψ.We rewrite this equation, together with the third and first ODEs of ( 49) and the algebraic relation for ρ, as the mixed algebraic-differential system of equations which contains three second-order ODEs.

Exact Solutions in Terms of Elementary Functions to the Reduced System of Three Second-Order ODEs
First, let us describe two simple classes of exact solutions to system (50) with ρ = const.1 • .Exact solution to system (50) with ρ = −σ: where C 1 , . . ., C 4 are arbitrary constants.2 • .Exact solution to system (50) with ρ = 0: Table 3.The functions ρ(x) and ω(x) determining exact solutions (51) to the system of ODEs (50).

No. Parameter σ Function ρ
Function ω Function ξ The first two ODEs of system (50) coincide, up to renaming, with system (24) and ( 25) at k = 1 and β = b, whose exact solutions were described above.For known ρ = ρ(x), the second equation in (50) for θ is a second-order linear homogeneous ODE, which has a particular solution θ 0 = ρ(x).Hence, the general solution of this equation is given by Formula (36).Once θ is found, the third equation in ( 50) is a second-order linear nonhomogeneous ODE for ϕ = ϕ(x), with the particular solution ϕ 0 = ρ(x) of the homogeneous ODE already known.Considering the aforesaid and using relevant formulas from [58], one can find the general solution to the ODE for ϕ.As a result, we obtain the solution to system (50) in the form where C 1 , . . ., C 4 are arbitrary constants, and ϕ p (x) is a particular solution to the third equation of system (50), which is defined as Here, the function W(x) = ρω x − ωρ x is the Wronskian determinant.Simple computations show that W(x) = 1.Exact solutions to system (50) can be obtained by substituting the functions ρ(x) and ω(x) from Table 3 into formulas ( 51) and ( 52). 4 • .For arbitrary σ, the function ρ(x) in solution ( 51) can be taken to be the Weierstrass elliptic function: where C 5 is an arbitrary constant.

Reduction of the System of PDEs with a Single Delay to a Unsteady System of Two First-Order Delay ODEs
Let the four delays in system (3) be identical: We look for generalized separable solutions to system (3) in the form with the functions θ = θ(x), ξ = ξ(t), ϕ = ϕ(t), and ψ = ψ(t) to be determined in the subsequent analysis.The functions (54) are selected to ensure that the composite arguments of f (. . . ) are only dependent on t.We impose the following additional condition (a linear differential constraint) on the function θ = θ(x): where the constants µ and ε are to be found in the subsequent analysis.Notably, without loss of generality, one can set ε = 0 in (55) in the nondegenerate case of µ = 0, since the translation of θ by a constant only leads, by virtue of representation (54), to new definitions of ϕ(t) and ψ(t).Substituting ( 54) into ( 3) and taking into account (55), we get where the notations f = k 1 φ + k 2 ψ, φ = ϕ(t − τ), and ψ = ψ(t − τ) are used for short.Relation ( 56) can be satisfied if we set The system of ODEs ( 57) is overdetermined because it consists of four equations for three functions ξ, ϕ, and ψ.On requiring that the last two equations of system (57) coincide, we find the parameter µ and other constants: For b 1 = b 2 , it follows from ( 58) that µ = 0, and hence, we can set ε = 0 in Equations ( 55) and (57).The general solution of Equation ( 55) with ε = 0 is given by Formula ( 14), while the first two equations of (57) are independent of ξ = ξ(t) and make up a closed nonlinear system of first-order delay ODEs for ϕ = ϕ(t) and ψ = ψ(t): where φ = ϕ(t − τ) and ψ = ψ(t − τ).Integrating the last equation in Equation ( 57) and using relations (58), we find that the function ξ = ξ(t) is expressed in terms of ϕ = ϕ(t) and ψ = ψ(t) as where A is an arbitrary constant.
It is noteworthy that the system of delay ODEs (59) generally admits one-component solutions ϕ = ϕ(t), ψ = 0 and ϕ = 0, ψ = ψ(t).In the absence of delay, these solutions can be expressed in terms of elementary functions.

The First Integral and Exact Solutions of the Reduced System of Two First-Order Delay ODEs
On eliminating the combination k 1 φ + k 2 ψ from the system of delay ODEs (59), we get Integrating yields the first integral where C 3 is an arbitrary constant, which may depend on τ.Substituting (61) into the first equation of system (59), we arrive at one first-order delay ODE Example 1. ODE (62) for the Lotka-Volterra system of PDEs without delay (3), where τ = 0 is a Bernoulli equation.The integration of it gives, in view of (61), an exact solution to system (59): where C 4 is an arbitrary constant.The respective function ξ(t) is defined by Formula (60) with τ = 0.It can be expressed in terms of elementary functions if, for example, In particular, by setting τ = 0 and C 3 = 0 in (60) and (63), we get Formulas ( 14), ( 54), ( 58) and (64) describe a new solution to the Lotka-Volterra system of PDEs without delay (3).
Similarly, if b 1 = b 2 = b and τ = 0, we get In general, the function ϕ for the delay ODE (62) must be defined on the interval [−τ, 0]: The Cauchy-type problem for the delay ODE (62) subject to the initial condition (66) can be solved using the method of steps [59,60].To this end, we split time t into segments of length τ and denote where t m = mτ and m = 0, 1, 2, . . .Then, on integrating Equation (62) from t m−1 to t on the interval [t m−1 , t m ], we get The function ϕ m (t) on the left-hand side of Formula (68) is sought on the interval [t m−1 , t m ], whereas ϕ m−1 (t) on the right-hand side is defined on the preceding interval [t m−2 , t m−1 ].The computations are carried out successively starting from m = 1, when the right-hand side is a known function defined on the initial interval (66).This step results in ϕ 1 (t).Then, one sets m = 2 and finds the function ϕ 2 (t) using (68) via the already known function ϕ 1 (t).The procedure continues further in a similar fashion.
For ε = 0, in order to determine ξ = ξ(t), ϕ = ϕ(t), and ψ = ψ(t), one has to solve the nonlinear system consisting of the last three equations in (57) with b 1 = b 2 and µ = 0. We assume that the four delay times in system (3) are linked by the single relation (7).On setting a 1 = a 2 = a, b 1 = b 2 = b, and c 1 = c 2 = 1, we look for solutions to system (3) in the form

Reductions and Exact
where λ is a free parameter, while θ = θ(x, t), ϕ = ϕ(x), and ψ = ψ(x) are functions to be determined in the subsequent analysis.Solutions of the form (69), where θ is a function of two independent variables and λ is an arbitrary parameter, generalize substantially solutions (17), in which θ is independent of t, and λ is expressed in terms of the system constants.We will refer to such solutions as incomplete separable solutions.
We assume the function θ = θ(x, t) to satisfy the additional periodicity condition One can verify that relations (18) remain valid under conditions ( 7) and (70).Considering the above and substituting (69) into (3), we obtain a closed system for ϕ and ψ consisting of two second-order ODEs without delay: The function θ = θ(x, t) satisfies the nonstationary Schrödinger equation and the periodic condition (70).Importantly, the coefficients of Equation (73) depend on the space variable x alone.We multiply ODEs (71) and (72) by k 1 and k 2 , respectively, and add up to obtain an ODE for ρ = k 1 ϕ + k 2 ψ.We write this equation in conjunction with ODE (71) and PDE (70) as the following mixed algebraic-differential system of equations: Systems ( 74)-( 77) can be solved sequentially from one equation to another.We start from the isolated stationary ODE (74).Divided by a, it coincides, up to obvious renaming, with Equation (24).PDE (77) is subjected to the additional conditions ( 7) and (70), which hold true, for example, in the stationary case θ = θ(x) and the nonstationary case of τ 2 = τ 1 and τ 4 = τ 3 .
(i) The isolated subsystem consisting of two ODEs (74) and (75) for ρ and ϕ coincides, up to notation, with the system of ODEs ( 24) and ( 25) with β = b.It follows that exact solutions to Equations ( 74) and (75) can be obtained using Table 2 and formulas with β = b, in which the functions and determining parameters must be renamed as follows: ϕ ⇒ ρ, θ ⇒ ϕ, b ⇒ b/a, β ⇒ b/a, and k ⇒ 1/a.(ii) The function ψ is determined by substituting the functions ρ and ϕ obtained in step (i) into (76), which results in (iii) The isolated subsystem of two ODEs, (74) and (77), for ρ and θ with θ t = 0 coincides, up to notation, with systems ( 24) and (25).Hence, exact solutions to Equations ( 74) and (77) can be obtained using the formulas from Table 2, where ϕ and the determining parameters must be renamed as follows: , some nonstationary exact solutions to systems (74)-(77) with ρ = 0 are defined by Formulas (78) and (79) for ϕ and ψ with any expression of the function θ = θ(x, t) specified below: where C and µ are arbitrary constants.Some other nonstationary exact solutions to systems (74)-( 77) with ρ = 0 can be obtained using Formulas (78) and (79) and the expression θ = e (b−λ)t ξ, where ξ = ξ(x, t) is any solution of the standard linear heat equation ξ t = aξ xx .5 • .If τ 2 = τ 1 and τ 2 = τ 4 , there are nonstationary exact solutions to systems (74)-(77) with ρ = −b that are defined by Formula (81) for ϕ and ψ and any expression of θ from (83) with b = 0. 6 • .For τ 1 = τ 2 , one can seek exact solutions to PDE (73) that satisfy the periodicity condition (70) and relation (7) in the form where m = 0 corresponds to a stationary solution.Substituting (84) into (71) yields the following linear stationary system of ODEs for ξ m = ξ m (x) and η m = η m (x): Since PDE (73) is linear in θ, an arbitrary linear combination of exact solutions (84), where α m are arbitrary constants, is also an exact solution to this equation.

Example 2.
For any of the four simplest solutions (5), the general solution of Equation (73) satisfying the periodicity condition (70) can be represented as where The constants A m , B m , C m , and D m , otherwise arbitrary, must ensure the convergence of the series (86) and derivatives θ t and θ xx ; the convergence is certainly ensured if, for example, we set A m = B m = C m = D m = 0 for all m > M, where M is an arbitrary positive integer.
On substituting (87) into the system of PDEs (3), we arrive at a nonlinear system of delay ODEs for ϕ and ψ, and a linear parabolic PDE with variable coefficients for θ: The system of delay ODEs (88) coincides with system (59) at b 1 = b 2 = b.Hence, it can be integrated as described in Section 2.7, where some exact solutions can also be found.With the substitution Equation (89) can be reduced to the standard linear heat equation whose exact solutions can be found, for example, in [61].

Reductions and Exact Solutions of a Multidimensional System with Different Diffusion Coefficients and Delays
Now, let us look at the following Lotka-Volterra-type nonlinear diffusion system of PDEs with n space variables: where is the Laplace operator and x 1 , . . ., x n are Cartesian coordinates, while the other notations are the same as in system (3).As above, we assume that relation (7) holds.System (90) admits reductions to simpler equations, which are similar to the reductions of system (3).We will illustrate this with specific examples and, in addition, present some exact solutions to (90) as well as solutions expressible in terms of solutions to linear PDEs. 1 • .There is an exact solution to system (90) with a 2 = a 1 that generalizes solutions ( 15) and ( 16): where x = (x 1 , . . ., x n ), and θ = θ(x) is a function satisfying the Helmholtz equation ∆θ = µθ. (92)
In the three-dimensional case, a few exact solutions to the linear PDEs (104) can be obtained using Formula (93) in which one must set µ = −b/a and µ = −(b − λ)/a, respectively.For other solutions to the linear PDEs (104), see [61].

Further Generalizations and Modifications
1 • .The nonlinear system of delay PDEs (3) can be generalized to the anisotropic case (when the diffusion coefficients depend on the spatial coordinate): where σ(x) is a given function.
Exact solutions of the linear ODE (107) for various σ(x) can be found in the handbook [58].In particular, for σ(x) = x k and σ(x) = e kx , the general solutions of ODE (107) can be expressed in terms of Bessel functions or modified Bessel functions.
Under conditions (7) and (21), system (106) admits more complicated generalized separable solutions of the form (17). 2 • .A further generalization of the system of delay PDEs (3) is Here, F(x, ρ) is an arbitrary function of two arguments, while L[u] is an arbitrary linear differential operator with respect to the spatial variable of the form where σ j (x) are some given functions and m is a positive integer.Under conditions ( 7) and ( 21), systems (108) and (109) admit a generalized separable solution of the form (17), where the functions ϕ = ϕ(x), ψ = ψ(x), and θ = θ(x) satisfy the system of ODEs where ρ = k 1 ϕ + k 2 ψ.
If the function F(x, ρ) is explicitly independent of x, so that F(x, ρ) = f (ρ), then system (110) admits particular solutions of the form In this case, if all the coefficients σ j of the differential operator (109) are independent of x, then the last equation in (110) for θ will be a linear ODE with constant coefficients.
3 • .Nonlinear Klein-Gordon-type systems with delay, which are obtained from the reaction-diffusion delay systems (3), ( 106) and (108) by formally replacing the first time derivative u t with the second derivative u tt , also admit exact generalized separable solutions of the form (4).

Brief Conclusions
We have described several reductions of a Lotka-Volterra-type nonlinear multi-parameter reaction-diffusion system with several delay times to simpler systems of ODEs with or without delay, as well as reductions to mixed systems consisting of ODEs and a single linear PDE without delay.
We have found many new exact solutions to the reaction-diffusion delay system in question, including over 50 solutions expressible in terms of elementary functions.Moreover, setting the delay times in these solutions to zero can result in many new exact solutions in terms of elementary functions to a simpler nonlinear Lotka-Volterra-type system without delays.
The exact solutions presented in this paper all involve several free constants of integration.Such solutions may be suitable for testing approximate analytical and numerical methods for solving reaction-diffusion equations as well as other nonlinear PDEs with delays.

1 .
Preliminary Remarks.The Lotka-Volterra-Type Diffusion System with Delays 1.2.The System of Delay PDEs in Question.Exact Solutions and Generalized Separable Solutions 1.3.Simplest Solutions (Stationary Points) 2. Reductions and Exact Solutions of the System of Delay PDEs with Different Diffusion Coefficients (a 1 = a 2 ) 2.1.Reduction of the System with Three Delays to a Single Second-Order Linear ODE without Delay 2.2.Reduction of the System of PDEs with Three Delays to a System of Two Second-Order ODEs without Delays 1. Introduction 1.1.Preliminary Remarks The Lotka-Volterra-Type Diffusion System with Delays 3 and C 4 are arbitrary constants, J 5/2 (ζ) and Y 5/2 (ζ) are Bessel functions of the first and second kind, and I 5/2 (ζ) and K 5/2 (ζ) are modified Bessel functions of the first and second kind.These can be expressed in terms of elementary functions as [58]:

a 1 and ψ = 0 .Remark 3 .
The value β = b in Table ) for various values of the parameters b and β.Notation: J 5/2 (ξ) and Y 5/2 (ξ) are the Bessel functions of the first and second kind, I 5/2 (ξ) and K 5/2 (ξ) are the modified Bessel functions of the first and second kind, and A 1 , A 2 , A 3 , C 2 , C 3 , and C 4 are arbitrary constants.No. Parameters b and β Function ϕ, Equation (24) Function θ, Equation (25) Function ξ = ξ(x) Solutions of the System of Delay PDEs with Identical Diffusion Coefficients (a 1 = a 2 ) 3.1.Reduction of the System of PDEs with Three Delays to a System of Two Second-Order ODEs without Delay and One Linear PDE

3. 3 .
Reduction of the System of PDEs with a Single Delay to a System of Two First-Order Delay ODEs and One Linear PDE We look for incomplete separable solutions to system (3) with a 1 = a 2 = a, b 1 = b 2 = b, and c 1 = c 2 = 1 and a common delay time (53) in the form u System (103) generalizes systems (74)-(77).Presented below are two simple stationary solutions to system (103) with ρ = const.1 • .Solution with ρ = 0 expressible in terms of solutions to two independent Helmholtz equations: ∆ϕ + (b/a)ϕ = 0, ∆θ + [(b − λ)/a]θ = 0,

Table 2 .
Solutions of systems (