Lie symmetries of nonlinear parabolic-elliptic systems and their application to a tumour growth model

A generalisation of the Lie symmetry method is applied to classify a coupled system of reaction-diffusion equations wherein the nonlinearities involve arbitrary functions in the limit case in which one equation of the pair is quasi-steady but the other not. A complete Lie symmetry classification, including a number of the cases characterised being unlikely to be identified purely by intuition, is obtained. Notably, in addition to the symmetry analysis of the PDEs themselves, the approach is extended to allow the derivation of exact solutions to specific moving-boundary problems motivated by biological applications tumour growth). Graphical representations of the solutions are provided and biological interpretation addressed briefly. The results are generalised on multi-dimensional case under assumption of radially symmetrical shape of the tumour.


Introduction
It is well-known that nonlinear reaction-diffusion systems (i.e. systems of second-order parabolic PDEs with reaction terms) are used to describe various processes in physics, chemistry, biology, ecology, etc (see, e.g., books [1,21,22,24]). Indeed, there are many special cases of such systems describing real world processes, among the most famous being Lotka-Volterra type systems [20,29] and systems used for modelling the chemical basis of morphogenesis [28] (see, e.g., [21] for details in modern terminology).
The most common are two-component reaction-diffusion systems of the form 1 where U(t, x) and V (t, x) are unknown functions, the nonnegative smooth functions D 1 (U) and D 2 (V ) are coefficients of diffusivity (conductivity), and F and G are arbitrary smooth functions describing interaction between components U(t, x) and V (t, x). Hereafter the t and x subscripts denote differentiation with respect to these variables. It should be noted that the search of Lie symmetries of RD systems was initiated many years ago and probably paper [30] was the first in this direction. A complete Lie symmetry classification was done much later because this problem is much more complicated in comparison with the case of a single reaction-diffusion equation. In particular, it was necessary to treat separately the cases of constant diffusivities (i.e. semilinear equations) and nonconstant (quasilinear) ones. All possible Lie symmetries of (1) with constant diffusivities were completely described in [7][8][9]. In the case of non-constant diffusivities, it has been done in [10,19] (in [11], the result was extended to the multi-component RD systems).
It should be stressed that there are models used to describe real world processes that are based on singular limits of a RD system (1), namely with D 2 (V ) = 0 or with the V t term absent. The system (1) with D 2 (V ) = 0 has been used to describe interactions involving species (cells) V , which are not able to diffuse in space (typical examples include some models for plankton). Some particular results about Lie symmetries of such systems have been recently obtained in [27], however the problem of a complete Lie symmetry classification is still not solved.
Here, however, we examine the class of systems (1) with the V t term negligible so that V is governed by a quasi-steady equation, i.e.
Obviously that this system is reduced to the form by the Kirchhoff substitution for the component V → V * . Thus, the class of parabolic-elliptic systems (3) is the main object of this study. The main motivation for this follows from paper [5], in which a model for tumour growth was derived. Under some biologically motivated assumptions (see p.570 in [5] for details) the model was simplified to a boundary-value problem with governing equations of the form (3). System (2) is of course potentially of relevance to any two-component reaction-diffusion model in which one species diffuses much faster than the other; it is also of interest purely from the symmetry point of view, since its symmetries differ essentially from those of (1). The paper is organized as follows. Section 2 is devoted to search for form-preserving transformations for the class of systems (3) in order to establish possible relations between RD systems that admit equivalent Lie symmetry algebras. Section 3 is devoted to the identification of all possible Lie symmetries that any system of the form (3) can admit. The main results of this section are presented in the form of Tables 1 and 2. In Sections 4 and 5, we apply the Lie symmetries derived for reduction of a nonlinear boundary value problem (BVP) with a free boundary modeling tumour growth in order to construct its exact solution. Moreover, possible biological interpretations of the results derived is discussed. Finally, we briefly discuss the result obtained and present some conclusions in the last section.
2 Form-preserving transformations for the class of systems (3) First of all we simplify notations in systems (3) in a natural way: The class of systems (4) contains three arbitrary functions and the Lie symmetry of its different representatives depends essentially on the form the triplet (D, F, G). Thus, the problem of a complete description of all possible Lie symmetries (the so called group classification problem) arises. In order to solve this, we can apply the Lie-Ovsiannikov approach (the name of Ovsiannikov arises because he published a remarkable paper in this direction, [26]) of the Lie symmetry classification, which is based on the classical Lie scheme and a set of equivalence transformations of the differential equation in question. However, it is well-known that this approach leads to very long list of equations with non-trivial Lie symmetry provided the given equation (system) contains several arbitrary functions (system (4) involves three such).
During the last two decades new approaches for solving group classification problems were developed, which are important for obtaining the so called canonical list of inequivalent equations admitting non-trivial Lie symmetry algebras and allow the solution of this problem in a more efficient way than a formal application of the Lie-Ovsiannikov approach. Here we use the algorithm based on so called form-preserving transformations [17,18], which were used initially for finding locally-equivalent PDEs, especially those nonlinear PDEs that are linearizable by a point transformation. Interestingly such transformations were implicitly used much earlier in [23] in order to find all possible heat equations with nonlinear sources that admit the Lie symmetry either of the linear heat equation or the Burgers equation. In [16], these transformations are called 'admissible transformations' and they were used to classify the Lie symmetries of a class of variable coefficient Korteweg-de Vries equations.
It was noted later (probably, paper [13] was the first in this direction) that the formpreserving transformations (other terminology is 'additional equivalence transformations' or 'admissible transformations') allow an essential reduction of the number of cases obtained via the classical Lie-Ovsiannikov algorithm (see extensive discussions on this matter in [9,10,14]). For example, it was proved in [14] using a set of form-preserving transformations that the canonical list of inequivalent reaction-diffusion-convection equations consists of 15 equations only (not the 30 derived by the Lie-Ovsiannikov algorithm). Now we present the main result of this section, namely the theorem describing a general form of form-preserving transformations for the class of systems (3).
Theorem 1 An arbitrary system of the form (4) be reduced to another system of the same form by the local non-degenerate transformation if and only if the smooth functions a, b, ϕ and ψ have the form where the functions α(t), β(t, x), K(t, x), P (t, x), L(t, x) and Q(t, x) are such that the following equalities: hold.
The proof of Theorem 1 is quite similar to that for the class reaction-diffusion systems (1) with D k = constant (k = 1, 2) presented in [6].
The group of equivalence transformations can be easily extracted from Theorem 2 by assuming that D, F and G are arbitrary smooth functions. Solving the system (7) under such condition one arrive at the following statement.
One sees that (7) is a nonlinear system of functional-differential equations and its general solution seems to be impossible without further restrictions. For example, assuming a constant diffusivity D, one may derive that β(t, x) is linear with respect to (w.r.t.) the variable x; however, β(t, x) can be nonlinear for the non-constant diffusivity D (see Tables 3 and 4 below for examples).
While the form of the transformations (6) is still quite general, it will be shown in the next section that only particular cases play important roles in solving the Lie symmetry classification problem for the nonlinear parabolic-elliptic systems (3).

Lie symmetries of a class of parabolic-elliptic systems
In order to find Lie symmetry operators of a system of the form (4) via the Lie method [2-4, 15, 25], one needs to consider manifold (S 1 , S 2 ) in the space of the following variables: The system (4) is invariant under the transformations generated by the infinitesimal operator when the following invariance conditions are satisfied: The operator X 2 is the second prolongation of the operator X, i.e.
where the coefficients ρ and σ with relevant subscripts are calculated by well-known formulae (see, e.g., [4,15,25]). Substituting (12) into (11) and eliminating the derivatives U t and V xx using (4), we can obtain the following system of determining equations (DEs): Obviously, the general solution of above system of DEs essentially depends on the form of the triplet (D, F, G). Thus, the problem of a complete description of all possible Lie symmetries (the so called group classification problem) arises. First of all we find the so called trivial algebra (other terminology used for this algebra is the 'principal algebra ' and the 'kernel of maximal invariance algebras') of the class in question, i.e. the maximal invariance algebra (MAI) admitted by each equation of the form (4). Under the natural assumption that D, F and G are arbitrary smooth functions the above system of DEs can be easily solved and the two-dimensional Lie algebra generated by the basic operators is obtained. In order to solve the problem of Lie symmetry classification for system (4), i.e. to find all possible systems admitting three-and higher-dimensional MAI, we use the algorithm based on form-preserving transformations, which was briefly discussed in Section 2. Moreover, it will be shown below that, similarly to the case of standard nonlinear reaction-diffusion systems [10], such an approach is more efficient also for parabolic-elliptic systems (compared to with the classical Lie-Ovsiannikov algorithm).
Let us formulate a theorem which gives complete information on the classical symmetry of the nonlinear system (4).
Theorem 2 All possible maximal algebras of invariance (up to equivalent representations generated by transformations of the form (22)) of the system (4) for any fixed nonconstant function D and nonconstant function vectors (F, G) are presented in Tables 1 and 2. Any other system of the form (4) with non-trivial Lie symmetry is reduced by a local substitution of the form to one of those given in Tables 1 and 2 (the constants C with subscripts are determined by the form of the system in question, some of them necessarily being zero in any given case).
The sketch of the proof. Consider the determining equations (13)- (20) under the assumption ∂D ∂U = 0 (we remind the reader that the case ∂D ∂U = 0 was examined in [7,8]). Obviously, equations (13)-(16) can be easily solved and the formulae are obtained (here a, b, r and p are arbitrary smooth functions and the dot means the time derivative). Now we can consider equations (17)- (20) as classification equations to find the function D and the pairs of the functions (F, G) for which the system (4) has a non-trivial Lie symmetry, i.e. its MAI is larger than the trivial algebra (21).
It follows from (17)-(18) that there is a special case 2b x = a t when the diffusivity D can be an arbitrary smooth function. Indeed equation (17) (20) have the form: The general solution of (24) depends essentially on the functions r and p (r 2 + p 2 = 0, otherwise a trivial Lie symmetry is obtained). In order to derive a complete result, one needs to examine two subcases. Let us consider the first subcase r = 0, p = 0, i.e. the first equation in (24) takes the form pF V = −2b 1 F . A simple analysis says that the general solution of (24) leads to the functions F, G and coefficients (23), which produce the system and MAI arising in case 1 of Table 1 provided b 1 = 0. If b 1 = 0 then case 3 of Table 2 is derived.
The second subcase r = 0 leads to case 2 of Table 1 and cases 1 and 2 of Table 2.
The condition 2b x − a t = 0 is a generic inequality and produces all other cases listed in Tables 1 and 2. Using this condition and (23), one easily solves the equation (18) and obtains where d = 0, C 1 , C 2 and k = 0 are arbitrary constants. Moreover, substituting (25) into (17), One notes that the system (4) with diffusivities (25) and (26) can be simplified to one with d = 1, C 1 = 0, C 2 = 1 by the appropriate equivalence transformations from (8). Now one needs to examine separately other determining equations with power law diffusivity, exponential diffusivity and the special power law coefficient k = −4/3. We present here only the last case, i.e. the power low diffusivity D = U −4/3 , which is the most complicated. Let us consider system (4) with D = U −4/3 . It turns out that we may assume that ξ 1 xx = 0 (otherwise only particular cases of (4) with D = U k will be derived). The classification equations (19)- (20) for D = U −4/3 take the form (because the function b depends only on x we use the standard notations for its derivatives and b (k) = d k b dx k , k = 3, 4 . . . ). Now we construct differential consequences of (27)-(28) w.r.t. x: Taking into account that b ′′ = 0 one easily finds the general solution of (29)-(30), where f and g are arbitrary smooth functions. The further analysis essentially uses the fact that the functions F and G depend only on U and V (not explicitly on t and x!). This means that relevant constraints on the functions f, g, b and p must take place.
Let us consider the most interesting case when the functions f and g are still arbitrary. This happens if and only if the constraints (here α 1 , α 2 and α 3 are some constants) hold in (31)-(32). Because the functions (31)-(32) must satisfy equations (27)-(28) for arbitrary functions f and g, the additional condition p = γr, springs up. Direct calculations show that the constant γ arising in the system obtained can be reduced to γ = 0 by equivalence transformation V + γ → V . Thus, we arrive at system (4) with D = U −4/3 and i.e. the parabolic-elliptic system where f and g are arbitrary smooth functions, α is an arbitrary constant. It turns out that MAI of system (35) essentially depends on the parameter α. Setting α = 0 and substituting (33)-(34) into (27)-(28), one obtains the systeṁ in order to find MAI of (35). The linear ODE system (35) can be easily solved, and as a result case 8 of Table 1 is obtained.
x reduces (35) to the system listed in case 3 of Table 3. If α > 0, then τ = α 4 t, y = √ α 2 x produces the system from case 4 of Table 3. It turns out that both systems can be reduced to one listed in case 8 of Table 1. The relevant transformations are presented in the 3rd column of Table 3 and they have been derived by applying Theorem 1 to the nonlinear system (35).
Thus, the most general system with the diffusivity D from (26) and its MAI are presented in case 8 of Table 1 and any other system with such a diffusivity admitting a nontrivial Lie symmetry is reduced to one listed in case 8 of Table 1. It should be noted that the corresponding substitutions are highly non-trivial and were found using a special case of form-preserving transformations (not equivalence transformations!) derived in Theorem 1.
Now we need to find all possible correctly-specified functions f and g leading to extensions of MAI of system (35). Such functions have been found by analysis of the differential consequences F xV = F tV = 0, G xU = G tU = 0, where the functions F and G are of the form (31)-(32). Finally, we have established only two special cases leading to the above mentioned extensions. ( where α 1 , α 2 and γ are arbitrary constants (α 2 1 +α 2 2 = 0). The systems (35) with these f and g are presented in case 9 of Table 1 and cases 5-6 of Table 3.
(ii) f = const, g = const. The constant f = 0 leads to case 12 of Table 2 and cases 6-7 of Table 4. The zero constant f leads to case 13 of Table 2 and cases 8-9 of Table 4.
At the final stage, we note that the form-preserving transformations derived above are still working, so that only case 9 of Table 1 and cases 12-13 of Table 2 Table 1 and cases 4-11, 14-21 in Table 2. All other parabolic-elliptic systems of the current class with non-trivial Lie symmetry can be reduced to those in Tables 1-2 by the appropriate equivalence transformations from (8) and/or form-preserving transformations and these are listed in Tables 3-4.
The sketch of the proof is now completed. ✷ Table 1 contains semi-coupled systems (i.e. those with F V G U = 0) as special subcases. A system specified in such a way has the same MAI provided it is not listed in Table 2 as a special case. For example, the system arising in case 1 with an arbitrary function f (U) and g(U) = constant is a semi-coupled system; however, this is not listed in Table 2, so that its MAI is still a 3-dimensional Lie algebra. Tables 1 and 2 We would like to conclude this section by the following observation. As it was noted in Section 2, the form-preserving transformations are used in order to solve the problem of the Lie symmetry classification for the class of systems (4). We have proved that there are 35 inequivalent systems admitting non-trivial Lie algebras (i.e. their MAI is three-and higher-dimensional) Table 1: Lie symmetries of system (4) in the case F V G U = 0

Continuation of
u τ = (e u u y ) y + α 1 e u + λ U = u − λτ, t = 1 λ e λτ , 16 u τ = (e u u y ) y + λ U = u − λτ, t = 1 λ e λτ , 21 0 = v yy + e βu + γ V = v + γ 2 y 2 e −βλτ and they are listed in Tables 1 and 2. The systems are inequivalent up to point transformations of the form (22), which are particular cases of form-preserving transformations described in Theorem 1. It is a standard routine to show that there are no other point transformations that allow the reduction of a system from Table 1 or 2 to another system from these tables. Thus, a canonical list of nonlinear parabolic-elliptic systems of the form (4) admitting a non-trivial Lie symmetry consists of the 35 inequivalent systems listed in Tables 1 and 2. A natural question arises: How many systems of the form (4) with non-trivial Lie symmetry can be derived using the Lie-Ovsiannikov algorithm? This means that we should take into account only the group of equivalence transformations (8). A formal answer is very simple: one should extend the list of above 35 systems by those from Tables 3 and 4, therefore 57 systems with the three-and higher-dimensional MAI will be derived (in the sense that the relevant MAI are inequivalent up to the equivalence transformations (8)).

Boundary value problems for a one-dimensional tumour growth model with negligible cell viscosity
In this section we apply the results of Lie symmetry classification derived in Section 2 for exact solving a real world model proposed in [5] (see p. 570). Consider the pressure difference with r = 1 [5] (see p. 569 therein) the equation (see (23) in [5]) for the tumour cell concentration, α(t, x), takes the form In order to simplify the boundary conditions we set (wherein c corresponds to the level of nutrient in [5], but might instead correspond to that of a chemotherapeutic drug) and arrive at the BVP where the concentrations U(t, x) and V (t, x) and the moving boundary location R(t) are unknown functions, while the functions S and Q are given smooth functions and m, α * and c ∞ are nonnegative parameters. One may note that the governing equations of this BVP with S = f (U)V −β and Q = g(U)V −β+1 are a particular case of case 2 in Table 1, so that they admit the Lie symmetry operator It turns out that the following statement can be easily proved, using the definition of invariance for BVPs with moving boundaries [12].

Theorem 3
The nonlinear BVP (39) with S = f (U)V −β and Q = g(U)V −β+1 is invariant with respect to the two-dimensional MAI generated by the Lie symmetry operators (40) and ∂ t .
Clearly the operator ∂ t cannot help to find any realistic solutions of BVP (39), while the Lie symmetry (40) guarantees a reduction of the problem, which leads to some interesting results. In fact, the operator generates the ansatz which immediately specifies the function R(t) = ω 0 √ t with the constant ω 0 > 0 to be found. The reduced ODE problem is The nonlinear BVP (42) is still a complicated problem and we were unable to solve it analytically in general. Let us consider a special case, namely (42) with β = 1 : In this case the exact solution can be found for arbitrary m in the form provided the functions f (ϕ) and g(ϕ) have the form Remark 4 A similar result can be obtained for arbitrary diffusivity in the governing equation (37), i.e. the drag coefficient k(α), however, the structure of the function f (ϕ) will essentially depend on the diffusivity. Now we present the final formulae for the original concentrations of tumour cells α(t, x) and nutrient (or drug) c(t, x) using formulae (38), ansatz (41) and solution (44): and the boundary conditions where the function f is defined by (45), in particular one obtains the cubic polynomial for m = 0: hence This solution for some specified parameters is presented in Fig. 1, Fig. 2 and Fig. 3. We note that the positive parameter α * < 1 can be interpreted as a natural concentration of tumour cells, while c ∞ is the concentrations of nutrient (drug) in the medium surrounding the tumour. The positive parameter ω 0 may be found from additional biologically motivated conditions. Notably the exact solution (46) can be easily generalised to one with two arbitrary parameters by the time translation t → t + t 0 .
It should be noted that the functions f (ϕ) and g(ϕ) in the cases m = 0 (constant diffusion) and m = 1 (as in the porous medium equation) are positive for the solution (44), meaning that the corresponding functions S(α, c) and Q(α, c) (see p.570 in [5]) are positive. The cases m = 0, 1 lead to a much more complicated structure for f (ϕ) and need to be examined separately for correctly-specified values of m. For example, the functions S(α, c) and Q(α, c) are positive on the solution for m = 1/2 in the case of the parameters α * and ω 0 used in Fig. 2.
We now turn to regimes in which the results may be susceptible to biological application, largely returning to the notation of [5]. The effective diffusivity of the cells is given by Figure 4: Surfaces representing the function S(α(t, x), c(t, x)) (left) and Q(α(t, x), c(t, x)) (right) on the known concentrations given by (46) (other parameters as in Fig. 3). That cell-cell repulsion can be expected to increase at high densities implies that but this increasing contribution to (51) may be offset by the (1 − α) 2 prefactor (which itself will be less significant because k(1) = 0 typically holds) that in part reflects the overall mass balance between cell and water velocities, v c and v w , respectively (so when 1 − α is small so typically is the cell velocity). The upshot is nevertheless that cases with m ≥ 0 are those of most relevance, with increasing m corresponding to stronger cell-cell repulsion at high densities. S(α, c) will be positive when c represents a nutrient, with S an increasing function of c in both cases (corresponding to β > 0 in the above). In this case, S will be an increasing function of α for sufficiently small α − α * ( and linear in α) and S = 0 at α = α * . One may note that the curves in Fig. 5 satisfy these properties. In fact, all the curves are increasing functions, the green and red curves in Fig. 5 behave as linear functions, while the blue curve is linear for sufficiently small ϕ = α − α * (the curves are built on the interval [0, ϕ max ] with ϕ max prescribed by the first formula in (44)). Moreover, f (ϕ) = 0 for ϕ = 0 because the free parameter ω 0 was taken to be 1 (see α 1 in (49)). Notably, S may decrease again for larger α, reflecting the limited supply of water to make new cells. The function S defined by (50) possesses such property provided α > α max ∈ (α * , α 2 ) , where α max is a local maximum of f from (49) and can be easily calculated.
In the drug case, S would typically be taken positive for small c and negative for large c (with S a decreasing function of c) -this behaviour cannot be captured by the above nonlinearities, so in interpretation of the results, c should be regarded as everywhere sufficiently large that S is negative, again with β > 0. Various dependencies on α can then be identified, and even S = 0 at α = 0 need not hold (though additional constraints need imposing if α reaches zero).
Q should be positive in both interpretations (absorption, rather than generation, though the generation of toxic materials in the drug case could allow Q < 0 to be contemplated) and an increasing function of c and α. According to Theorem 3 Q = g(U)V −β+1 (with U and V given in (38)), therefore one needs β > 1 and an increasing function g in order to satisfy the above requirement. Q > 0 implies by the maximum principle that c will be an increasing function of x, having a minimum at x = 0. Such behaviour of c is pictured in Fig. 1-3 for different values of m. The qualitative form of α will depend on the initial data and the form of S; in the case S > 0, the nutrient supply is highest at the edge of the tumour, but α is held at α * there by the stress balance (more water is taken up to prevent α exceeding α * due to cell division), so, perhaps counterintuitively, it is quite possible for α to be maximal at x = 0 if S remains sufficiently large, despite the decrease in c. However, such behaviour is more natural for S < 0 (the drug case) since the cell-kill rate can be expected to be highest at the edge of the tumour.

5
Higher-dimensional tumour growth model without cell viscosity In Section 4, we examined the one-dimensional model (39) describing the tumour growth. Here we generalise the model to higher dimensions under assumption that tumour cells create a radially symmetrical shape (a circular cylinder for n = 1 and a sphere for n = 2). In this case, the nonlinear BVP (39) takes the following form: where n = 1, 2 and the same notation are used as in Section 4, however, the space variable x is replaced by the variable r = x 2 1 + · · · + x 2 n + x 2 n+1 . One may check that the governing equations of this BVP with S = f (U)V −β and Q = g(U)V −β+1 are invariant w.r.t. the Lie symmetry operator which has the same structure as (40).

Theorem 4
The nonlinear BVP (54) with S = f (U)V −β and Q = g(U)V −β+1 is invariant with respect to the two-dimensional MAI generated by the Lie symmetry operators (55) and ∂ t .
Thus, the Lie symmetry (55) generates the ansatz which immediately specifies the function R(t) = ω 0 √ t and reduces the two-dimensional BVP (54) (with S = f (U)V −β and Q = g(U)V −β+1 ) to the one-dimensional problem The nonlinear problem (57) is still not integrable for arbitrary given coefficients, however, one with β = −1 is again a solvable case. Direct calculations show that formulae (44) and (45) can be generalised in this case as follows and g = q 0 ϕ, q 0 > 0, In order to obtain the final formulae for the original concentrations of tumour cells α(t, x) and nutrient (drug) c(t, x), we use formulae (38), ansatz (56) and solution (58). The formulae obtained coincide with those presented in (46) for α(t, x) and R(t) while for the nutrient (drug) concentration we have the solution The biomedical interpretation of the various regimes coincides with that in Section 4 and is omitted here.

Conclusions
The main part of this paper is devoted to the Lie symmetry classification of the class of parabolic-elliptic systems (3). First of all, the structure of form-preserving (admissible) transformations for the systems of the form (3) is established (see Theorem 1) in order to establish possible relations between systems that admit equivalent MAIs. Theorem 2 gives a complete list of inequivalent parabolic-elliptic systems consisting of 35 systems, which are invariant under non-trivial MAIs, i.e. under the three-and higher-dimensional Lie algebras. Moreover, we have established all possible subclasses of parabolic-elliptic systems (see Tables 3 and 4) with non-trivial Lie symmetry, which are reducible to those listed in Tables 1 and 2 by the relevant form-preserving transformations. Notably, our result is a new confirmation of importance of finding form-preserving transformations because they allow to reduce essentially (in contrary to the group of equivalence transformations) the number of inequivalent systems. In fact, the traditional Lie-Ovsiannikov algorithm, which is based on the group of equivalence transformations, would lead to 57 parabolic-elliptic systems instead of 35 systems listed in Tables 1 and 2  because 22 systems from Tables 3 and 4 should be added. From the applicability point of view, the most interesting parabolic-elliptic systems occur in Table 1 because the systems from Table 2 are semi-coupled. It is worth to note that MAIs for the systems arising in all the cases excepting 7 and 14 are some analogs of those for single reaction-diffusion equations. In fact, fixing the diffusion coefficient in the first equation of the corresponding system one may easily identify its generic analog in Table 4 [14] (see cases 3-5 and 7 therein). For example, MAIs of the systems from cases 8 and 9 contain the well-known algebra sl(2, R) as a subalgebra (see the second, third and fourth operators in the last column). The same subalgebra occurs for the diffusion equation with the same diffusivity (see case 3 in Table 4 [14]). By the way, Ovsiannikov was the first who identified this special case for nonlinear diffusion equations [26]. Of course, the systems in Table 1 have more general structures (many of them contain arbitrary functions as coefficients) in contrast to their analogs listed in Table  4 [14]. However, it happens on the regular basis if one considers systems instead of single equations (see, e.g., case 10 in Table 1 [10] involving the same diffusivity).
Cases 7 and 14 in Table 1, have no analogs among single reaction-diffusion equations. Moreover, the relevant MAIs generate infinite-dimensional Lie algebras because they contain the operators X ∞ involving an arbitrary smooth function ϕ(t) (see the last column of Table 1). Interestingly that the both MAIs from Cases 7 and 14 contain again the well-known algebra sl(2, R) as a subalgebra. In fact, if one sets ϕ(t) = t and ϕ(t) = t 2 then the operators ∂ t , X ∞ with ϕ(t) = t and X ∞ with ϕ(t) = t 2 produce nothing else but sl(2, R) with new representations (it is a simple task to check the commutator relations), which do not happen in the case of single reaction-diffusion equations.
In order to demonstrate applicability of our pure theoretical result, we apply the Lie symmetry corresponding to one-parameter group of scale transformation for solving a (1+1)dimensional BVP with a free boundary modeling tumour growth [5]. We reduce the given nonlinear BVP to that governed by ODEs. As a result, its exact solution was constructed under additional restrictions on the coefficients arising in the governing equations. Moreover, possible biological interpretations of the solution is discussed. In particular, it was shown that the results obtained allow plausible interpretation in the case when the model describes growth of tumour cells consuming a nutrient. Finally, the results for the above (1+1)-dimensional BVP were generalised on multi-dimensional case under assumption that tumour cells create a radially symmetrical shape.