If circles of radius $r_1$, $r_2$, $r_3$, $r_4$ are mutually tangent, and also circles of radius (say) $r_2$, $r_3$, $r_4$, $r_5$, then we have
$$\begin{align}
\left( k_1 + k_2 + k_3 + k_4 \right)^2 &= 2 \left( k_1^2 + k_2^2 + k_3^2 + k_4^2 \right) \\
\left( k_2 + k_3 + k_4 + k_5 \right)^2 &= 2 \left( k_2^2 + k_3^2 + k_4^2 + k_5^2 \right)
\end{align}$$
with $k_i = \pm 1/r_i$ as desired.
If you need the fourth circle's radius, then, of course, you can solve the equations in stages: get $k_4$ from the first, and use that in the second to get $k_5$.
If you don't really care about the fourth radius, we can eliminate $k_4$ from the equations. (For instance, subtract one from the other to get rid of the $k_4^2$ terms, and solve for $k_4$ in the resulting linear equation. Then substitute this $k_4$ into either of the original equations.) Assuming $k_3 \neq k_5$, we arrive at this relation:
$$16 k_1^2 + 16 k_2^2 +
k_3^2 + k_5^2 + 16 k_1 k_2 - 8 k_1 k_3 - 8 k_1 k_5 - 8 k_2 k_3 - 8 k_2 k_5 - 2 k_3 k_5 = 0$$
Thus,
$$k_5 = 4 k_1 + 4 k_2 + k_3 \pm 4 \sqrt{k_1 k_2 + k_2 k_3 +k_3 k_1}$$
Let $P := (a \cos\phi, b \sin\phi)$ on an origin-centered ellipse with radii $a$ and $b$; define $c := \sqrt{a^2-b^2}$, so that the ellipse's eccentricity is $e := c/a$. The line through $P$, normal to the ellipse —that is, in direction $(b\cos\phi,a\sin\phi)$— meets the $x$-axis at $K:= (k,0)$, where $k:= c^2/a \cos\phi$. So, $K$ is the center of a circle internally tangent to the ellipse at $P$, and its radius, $r$, is given by
$$r^2 = |PK|^2 = \frac{b^2(a^2-c^2\cos^2\phi)}{a^2} = \frac{b^2(c^2-k^2)}{c^2} \tag{1}$$
so that
$$\frac{r^2}{b^2}+\frac{k^2}{c^2}=1 \tag{2}$$
This allow us to write, for some $\theta$,
$$r = b\sin\theta \qquad k = c \cos\theta \tag{3}$$
Now, suppose $\bigcirc K_0$ and $\bigcirc K_1$ are circles internally tangent to the ellipse, with respective centers and radii given by $(3)$ for $\theta = \theta_0$ and $\theta=\theta_1$. If these circles are tangent to each other (with $K_1$ "on the right" of $K_0$), then
$$\begin{align}
k_0 + r_0 &= k_1 - r_1 \\[4pt]
\to\quad -2 c \sin\frac{\theta_0 + \theta_1}{2} \sin\frac{\theta_0 - \theta_1}{2} &= -2 b \sin\frac{\theta_0 + \theta_1}{2} \cos\frac{\theta_0 - \theta_1}{2} \\[6pt]
\to\quad \tan\frac{\theta_0 - \theta_1}{2} &= \frac{b}{c} \\[6pt]
\to\quad \theta_1 &= \theta_0 - 2\arctan\frac{b}{c} \\[6pt]
&= \theta_0 - 2\arccos e \tag{4}
\end{align}$$
More generally, if circles $\bigcirc K_i$, defined by $\theta = \theta_i$ in $(3)$, form a tangent chain, then
$$\theta_i = \theta_0 - 2 i \arccos e \tag{5}$$
where index $i$ is subject to certain viability conditions (eg, $\theta_i \geq 0$) that we'll assume hold. Thus, defining $\varepsilon := 2\arccos e$, we have
$$\begin{align}
\frac{r_{i+j} + r_{i-j}}{r_i} &= \frac{b\sin(\theta_0-(i+j)\psi)+b\sin(\theta_0-(i-j)\varepsilon)}{b \sin(\theta_0-i\varepsilon)} \\[6pt]
&= 2\cos j\varepsilon = 2\cos( 2j \arccos e ) \\[4pt]
&= 2\,T_{2j}(e) \tag{6}
\end{align}$$
where $T_{2j}$ is the $2j$-th Chebyshev polynomial of the first kind. Notably, the value of $(6)$ is independent of $i$. In particular, if we take $j=3$ and both $i=4$ and $i=7$, we can write
$$\frac{r_{4-3}+r_{4+3}}{r_4} = 2\;T_{2\cdot 3}(e) =\frac{r_{7-3}+r_{7+3}}{r_7} \tag{7}$$
which gives the result. $\square$
Addendum. In this follow-up question, @g.kov asks when an ellipse allows a "perfect packing" of $n$ tangent circles along its axis. It seems reasonable to append here a justification of the condition given there.
In a perfect packing, the first and last circles in a chain are tangent to the ellipse at the endpoints of the axis, so that their radii match the ellipse's radius of curvature (namely, $b^2/a$) at those points. Thus, we have
$$r_0 = r_{n-1} = \frac{b^2}{a} \quad\to\quad \sin\theta_0 = \sin\theta_{n-1} = \frac{b}{a} \quad\to\quad \cos\theta_0 = \cos\theta_{n-1} = e \tag{8}$$
We can say that $\theta_0 = \pi - \arccos e$ and $\theta_{n-1} = \arccos e$. By $(5)$, this implies
$$\arccos e = \theta_{n-1} = \theta_0 - 2(n-1)\arccos e = (\pi - \arccos e) - 2(n-1)\arccos e \tag{9}$$
so that
$$\pi = 2n\arccos e \qquad\to\qquad \cos \frac{\pi}{2n} = e \tag{10}$$
This is equivalent to @g.kov's condition for a perfectly-packable ellipse. $\square$
Best Answer
I'm using the same notation as in the link mentioned by David Mitra. By similarity we have $$ \frac{r}{OA}=\frac{r_a}{OA-r-r_a}=\sin\left(\frac{\alpha}{2}\right) $$ where $\alpha=\angle A$. Then $$ \frac{r_a}{r} = \frac{1-\sin(\alpha/2)}{1+\sin(\alpha/2)} = \frac{1-\cos((\beta+\gamma)/2)}{1+\cos((\beta+\gamma)/2)} = \tan\left(\frac{\beta+\gamma}{4}\right)^2 $$ where we used that $\alpha+\beta+\gamma=\pi$. Similarly $$ \frac{r_b}{r} = \tan\left(\frac{\alpha+\gamma}{4}\right)^2 \qquad\qquad \frac{r_c}{r} = \tan\left(\frac{\alpha+\beta}{4}\right)^2 $$ Therefore, $\sqrt{r_ar_b}+\sqrt{r_ar_c}+\sqrt{r_br_c}$ can be written as $$ \begin{split} \frac{\sqrt{r_ar_b}+\sqrt{r_ar_c}+\sqrt{r_br_c}}{r} = \frac{\sin(\rho)\sin(\sigma)\cos(\tau)+\sin(\rho)\cos(\sigma)\sin(\tau) +\cos(\rho)\sin(\sigma)\sin(\tau)}{\cos(\rho)\cos(\sigma)\cos(\tau)} \end{split} $$ where $\rho=(\beta+\gamma)/4$, $\sigma=(\alpha+\gamma)/4$, and $\tau=(\alpha+\beta)/4$. Note that $\rho+\sigma+\gamma=\pi/2$. Consider $x$, $y$, $z$ such that $x+y+z=\pi/2$. Then $$ \begin{split} \sin(x)\sin(y)\cos(z) &= \frac{1}{2}(\cos(x-y)-\cos(x+y))\cos(z) \\&= \frac{1}{4}(\cos(x-y+z)+\cos(x-y-z) -\cos(x+y+z)-\cos(x+y-z)) \\&= \frac{1}{4}( \cos\left(\frac{\pi}{2}-2y\right) + \cos\left(-\frac{\pi}{2}+2x\right) - \cos\left(\frac{\pi}{2}\right)-\cos\left(\frac{\pi}{2}-2z\right) \\&=\frac{1}{4}\left( \sin(2y)+\sin(2x)-\sin(2z)\right) \end{split} $$ And similarly $$ \begin{split} \cos(x)\cos(y)\cos(z) &= \frac{1}{2}(\cos(x-y)+\cos(x+y))\cos(z) \\&= \frac{1}{4}(\cos(x-y+z)+\cos(x-y-z) +\cos(x+y+z)+\cos(x+y-z)) \\&=\frac{1}{4}\left( \sin(2y)+\sin(2x)+\sin(2z)\right) \end{split} $$ This means that $(\sqrt{r_ar_b}+\sqrt{r_ar_c}+\sqrt{r_br_c})/r$ is equal to $$ \frac{\sqrt{r_ar_b}+\sqrt{r_ar_c}+\sqrt{r_br_c}}{r} = \frac{\sin(2\rho)+\sin(2\sigma)+\sin(2\tau)}{\sin(2\rho)+\sin(2\sigma)+\sin(2\tau)}=1 $$