The following derivation follows closely Vivas-Barber, A. et al. Dynamics of an "SAIQR" influenza model. 2014.
I have used Wolfram Mathematica for calculations as shown below. You should validate yourself that the calculations are correct.
First, we define the state vector $E=(S,A,I,Q,R)$ and the total population size $$N=S+A+I+Q+R.$$ Solving $\dot{E}=0$ for the trivial, disease-free equilibrium, we obtain:
$$N = \frac{\Lambda -\mu_1 Q}{\mu}.$$
Next, we calculate the Jacobian of the system at $E_0 = \left((\Lambda -\mu_1 Q)/\mu,0,0,0,0\right)$:
$$\left(
\begin{array}{ccccc}
-\mu & -\frac{\beta_1 \Lambda }{\mu } & -\frac{\beta_2 \Lambda }{\mu } & 0 & 0 \\
0 & -\alpha +\frac{\beta_1 \Lambda }{\mu }-\mu & \frac{\beta_2 \Lambda }{\mu } & 0 & 0 \\
0 & \alpha & -\theta -\mu & 0 & 0 \\
0 & 0 & \theta & -\gamma -\mu -\mu_1 & 0 \\
0 & 0 & 0 & \gamma & -\mu \\
\end{array}
\right)$$
and the corresponding eigenvalues:
$$
\begin{array}{l}
\lambda_1 = -\mu, \\
\lambda_2 = -\mu, \\
\lambda_3 = -\frac{1}{2 \mu }\left(\sqrt{2 \beta_1 \Lambda \mu (\theta -\alpha )+\mu \left(4 \alpha \beta_2 \Lambda +\mu (\alpha -\theta )^2\right)+\beta_1^2 \Lambda ^2}+\mu (\alpha +\theta +2 \mu )-\beta_1 \Lambda \right), \\
\lambda_4 = \frac{1}{2 \mu }\left(\sqrt{2 \beta_1 \Lambda \mu (\theta -\alpha )+\mu \left(4 \alpha \beta_2 \Lambda +\mu (\alpha -\theta )^2\right)+\beta_1^2 \Lambda ^2}-\mu (\alpha +\theta +2 \mu )+\beta_1 \Lambda \right), \\
\lambda_5 = -\gamma -\mu -\mu_1.
\end{array}$$
Rephrasing the conditions $\lambda_3<0$ and $\lambda_4<0$ as $\mathfrak{R}_0 < 1$ leads to:
$$\mathfrak{R}_0 = \frac{\alpha \beta_2 \Lambda }{(\theta +\mu ) (\mu (\alpha +\mu )-\beta_1 \Lambda )},$$
and a necessary condition on the parameters: $\beta_1 < (\alpha +\theta +2 \mu )\mu/\Lambda $.
Clear["Global`*"];
eqs = {
S'[t] == Λ - β1 S[t] A[t] - β2 S[t] Ι[t] - μ S[t],
A'[t] == β1 S[t] A[t] + β2 S[t] Ι[t] - (μ + α) A[t],
Ι'[t] == α A[t] - (μ + θ) Ι[t],
Q'[t] == θ Ι[t] - (μ + μ1 + γ) Q[t],
R'[t] == γ Q[t] - μ R[t]
};
vars = DeleteDuplicates@Cases[eqs, _Symbol[t], All];
params = DeleteDuplicates@Cases[eqs, _Symbol, All];
$Assumptions = Thread[params > 0];
Eliminate[{Total[Last /@ eqs] == 0, Ν == Total[vars]}, {R[t]}];
Ν0 = First@SolveValues[%, Ν];
J = D[Last /@ eqs, {vars}] /. S[t] -> Ν0 /. (vars -> 0 // Thread);
λ = Eigenvalues[J];
FullSimplify[Reduce[λ[[3]] < 0 && λ[[4]] < 0, Reals]]
(* β1 Λ < μ (α + θ + 2 μ) && β2 < ((θ + μ) (-β1 Λ + μ (α + μ)))/(α Λ) *)