There is a systematic way to treat this question not just in the above particular case but in that of any four of the around $40,000$ triangle centres documented in the "Encyclopedia of Triangle Centers" using the $p,q$ method. We suppose that the centers are described by the triangle functions $f_i$ and put $x_i=f_1(b,1,a)+pf_i(1,a,b)$, $y_i=qf_i(1,a,b)$. Then the four centres are cyclic if and only if the determinant of the $4×4$ matrix with rows $(x_i^2+y_i^2,x_i,y_i,1)$ is identically zero. This can sometimes be done by hand for simple centrrs but it is easy to mechanise it using, say, Mathematica to churn out as many of the $40,000^4$ potential results as desired. We are using the notation of the above Encyclopedia which is available online. We are also assuming, as we may, that the functions are homogeneous and sum symmetrically to $1$. In addition, we regard the above determinant as a function of $a$ and $b$,, using the relationshios $p=\frac 12(b^2÷1-a^2)$ and $q^2=b^2-p^2$.
It is perhaps of interest to point out the following thene. It is a surprise if three arbitrary points in the plane lie on a line Applied to triangle centres this means that usually three such are not collinear, unless the triangle is isosceles in which case all centres lie on the symmetry axis. However there are special cases where they are collinear for all triangles--think euler.line. A similar situation holds for four centres and cyclicity. There are certain foursomes for which it is valid for any triangle (as in the OP) andthere are, for a given foursome, certain special triangle shapes for which it holds. All of these questions can be attacked by the above method.