Consider the parabolic system \begin{align} \begin{cases} u_t - \Delta\Big((a_1 + a_{11} u + a_{12} v) u\Big) = 0, & t >0, \ x \in \mathbb R^n \\ v_t - \Delta\Big((a_2 + a_{22} v + a_{21} u) v\Big) = 0, & t >0, \ x \in \mathbb R^n \end{cases} \end{align} where $a_1,a_2,a_{11},a_{22},a_{12},a_{21}$ are non-negative constants and the unknowns are $u: (0,\infty)\times \mathbb R^n \to \mathbb R$ and $v: (0,\infty)\times \mathbb R^n \to \mathbb R$.
Main Question: Formally, can we compute a (possibly radial) self-similar solution, that is a solution of this system the form $$\left(\frac{1}{t^{\gamma}} u(x/t^{\alpha}), \ \frac{1}{t^{\kappa}} v(x/t^{\beta})\right)$$ for suitable $\alpha$, $\beta$, $\gamma$, $\kappa$?
Subquestion: Partial results assuming something (reasonable) extra on the coefficients are also very welcome.
Remark 1. The computation that I'm looking for in the system above is very classical in the particular case of the heat equation $u_t -\Delta u = 0$ (i.e. $a_{11} = a_{12} = a_{2} = a_{22} = a_{21} = 0$). For your convenience, I'm typing it up below (following what can be found in Evans' book): Let us look for a solution of the form \begin{equation*} u(x, t)=\frac{1}{t^{\alpha}} v\left(\frac{x}{t^{\beta}}\right) \quad\left(x \in \mathbb{R}^{n}, t>0\right) \end{equation*} where the constants $\alpha, \beta$ and the function $v: \mathbb{R}^{n} \rightarrow \mathbb{R}$ must be found. Plugging this into the heat equation, we arrive at \begin{equation*} \alpha t^{-(\alpha+1)} v(y)+\beta t^{-(\alpha+1)} y \cdot D v(y)+t^{-(\alpha+2 \beta)} \Delta v(y)=0 \end{equation*} for $y:=t^{-\beta} x$. In order to transform this into an expression involving the variable $y$ alone, we take $\beta=\frac{1}{2}$, so that the equation reduces to \begin{equation*} \alpha v+\frac{1}{2} y \cdot D v+\Delta v=0 \end{equation*}
Assuming also that $v$ is radial -- that is, $v(y)=w(|y|)$ for some $w: \mathbb{R} \rightarrow \mathbb{R}$, we get \begin{equation*} \alpha w+\frac{1}{2} r w^{\prime}+w^{\prime \prime}+\frac{n-1}{r} w^{\prime}=0 \end{equation*} for $r=|y|, '=\frac{d}{d r}$. Setting $\alpha=\frac{n}{2}$, we get \begin{equation*} \left(r^{n-1} w^{\prime}\right)^{\prime}+\frac{1}{2}\left(r^{n} w\right)^{\prime}=0 \end{equation*} Thus \begin{equation*} r^{n-1} w^{\prime}+\frac{1}{2} r^{n} w=a \end{equation*} for some constant $a$. Assuming $\lim _{r \rightarrow \infty} w, w^{\prime}=0$, we conclude $a=0$; whence \begin{equation*} w^{\prime}=-\frac{1}{2} r w \end{equation*} But then for some constant $b$ \begin{equation} w=b e^{-\frac{r^{2}}{4}} \end{equation} In conclusion, we arrive at the expected Gaussian profile $$\frac{b}{t^{n / 2}} e^{-\frac{|x|}{4 t}}$$
Remark 2 For the same question for the Porous Medium Equation $u_t - \Delta u^m = 0$, see Section 4.4 Source-type solutions. Selfsimilarity from page 69 of https://verso.mat.uam.es/~juanluis.vazquez/BKPME2006six.pdf