  We consider mixed type multiple orthogonal polynomials associated with a
system of weight functions consisting of two vectors. One vector is defined in
terms of scaled modified Bessel function of the first kind $I_\mu$ and
$I_{\mu+1}$, the other vector is defined in terms of scaled modified Bessel
function of the second kind $K_\nu$ and $K_{\nu+1}$. We show that the
corresponding mixed type multiple orthogonal polynomials exist. For the special
case that each multi-index is on or close to the diagonal, basic properties of
the polynomials and their linear forms are investigated, which include explicit
formulas, integral representations, differential properties, limiting forms and
recurrence relations. It comes out that, for specified parameters, the linear
forms of these mixed type multiple orthogonal polynomials can be interpreted as
biorthogonal functions encountering in recent studies of products of two
coupled random matrices. This particularly implies a Riemann-Hilbert
characterization of the correlation kernel, which provides an alternative way
for further asymptotic analysis.
