Without loss of generality, all the $a_i$'s and $b_i$'s are nonzero.
Let $\tilde d$ denote the difference between the left- and right-hand sides of the conjectured inequality $(*)$, which then of course can be rewritten as $\tilde d\ge0$.
In the previous version of my answer, I rewrote $\tilde d$ in new variables, $x_i$ and $y_i$, after which the inequality $\tilde d\ge0$ could be (rigorously) verified with Mathematica (in about 22 min).
Here that expression for $\tilde d$ is further rewritten -- in new, "more-macro", variables -- so that the resulting expression can be rather easily analyzed, to prove the inequality $(*)$. Indeed, let
$p_i:=(x_i-y_i)y_i$, $x_i:=a_1 a_2 a_3/a_i$, $y_i:=b_1 b_2 b_3/b_i$,
\begin{equation}
c_1:=p_2^2 + p_2 p_3 + p_3^2\ge0,\quad c_2:=p_1^2 + p_1 p_3 + p_3^2\ge0,\quad c_3:=p_2^2 + p_2 p_1 + p_1^2\ge0, \tag{0}
\end{equation}
and $z_i:=y_i^2\ge0$.
Note that $x_1 x_2 x_3=(a_1 a_2 a_3)^2>0$ and
$y_1 y_2 y_3=(b_1 b_2 b_3)^2>0$; moreover,
\begin{equation}
(p_1+z_1)(p_2+z_2)(p_3+z_3)\ge0. \tag{1}
\end{equation}
The crucial identity is
$$
\tilde d\,y_1 y_2 y_3=d:=
p_1 p_2 p_3+c_1 z_1+c_2 z_2+c_3 z_3.
$$
Since $y_1 y_2 y_3>0$, $\tilde d$ equals $d$ in sign. So, it suffices to show that $d\ge0$ -- for any real $p_i$'s, the $c_i$'s as in $(0)$, and any nonnegative $z_i$'s satisfying $(1)$.
Note here that without loss of generality $p_1 p_2 p_3<0$ -- otherwise, $d\ge0$ immediately follows because the $c_i$'s and $z_i$'s are nonnegative. So, we may assume that the $p_i$'s are are all nonzero and hence the $c_i$'s are all strictly positive.
Take any nonzero real $p_i$'s and any nonnegative $z_i$'s such that $(1)$ holds. Let us then fix those $z_1$ and $z_2$, and let $z_3$ be decreasing as long as $z_3$ remains nonnegative and $(1)$ holds; clearly, this process can stop only when the value of $z_3$ becomes either $0$ or $-p_3$, and in the latter case we must have $-p_3>0$.
Moreover, since $c_i>0$ for all $i$, the value of $d$ will not increase after this process is complete.
We can then proceed similarly by decreasing $z_2$ (instead of $z_3$), and then by decreasing $z_1$.
Let now $(z_1,z_2,z_3)$ be any minimizer of $d$. Then it follows from the above reasoning that $z_i\in\{0,-p_i\}$ for each $i=1,2,3$; moreover, if at that $z_i=-p_i$ for some $i$, then we must have $-p_i>0$.
So, by the symmetry with respect to permutations of the indices, it is enough to consider the following four cases:
(i) $z_1=-p_1>0$, $z_2=-p_2>0$, $z_3=-p_3>0$;
(ii) $z_1=-p_1>0$, $z_2=-p_2>0$, $z_3=0$;
(iii) $z_1=-p_1>0$, $z_2=0$, $z_3=0$;
(iv) $z_1=0$, $z_2=0$, $z_3=0$, so that $(1)$ becomes $p_1 p_2 p_3\ge0$.
In case (i), $\min_{z_1,z_2,z_3}d=-(p_1 + p_2) (p_1 + p_3) (p_2 + p_3)>0$.
In case (ii), $\min_{z_1,z_2,z_3}d=-p_1 p_2 (p_1 + p_2) - p_1 p_2 p_3 + (-p_1 - p_2) p_3^2$, which is a convex quadratic polynomial in $p_3$, with discriminant $-p_1 p_2 (4 p_1^2 + 7 p_1 p_2 + 4 p_2^2)<0$, whence again $\min_{z_1,z_2,z_3}d>0$.
In case (iii), $\min_{z_1,z_2,z_3}d=-p_1 (p_2^2 + p_3^2)>0$.
In case (iv), $\min_{z_1,z_2,z_3}d=p_1 p_2 p_3\ge0$.
Thus, $\min_{z_1,z_2,z_3}d\ge0$ in all cases, and the inequality in question is proved.
