We show that the average characteristic polynomial Pn(z) = E[det(zI−M)] of the random Hermitian matrix ensemble Z 1 n exp(−Tr(V (M) − AM))dM is characterized by multiple orthogonality conditions that depend on the eigenvalues of the external source A. For each eigenvalue aj of A, there is a weight and Pn has nj orthogonality conditions with respect to this weight, if nj is the multiplicity of aj. The eigenvalue correlation functions have determinantal form, as shown by Zinn-Justin. Here we give a different expression for the kernel. We derive a Christoffel-Darboux formula in case A has two distinct eigenvalues, which leads to a compact formula in terms of a Riemann-Hilbert problem that is satisfied by multiple orthogonal polynomials.