Multivariate Krawtchouk and Meixner polynomials are constructed explicitly as birth and death polynomials, which have the nearest neighbour interactions. They form the complete set of eigenpolynomials of birth and death processes with the birth and death rates at population \({\varvec{x}}=(x_1,\ldots ,x_n)\in \mathbb {N}_0^n\) are \(B_j({\varvec{x}})=\bigl (N-\sum _{i=1}^nx_i\bigr )\) , \(\bigl (\beta +\sum _{i=1}^nx_i\bigr )\) and \(D_j({\varvec{x}})=p_j^{-1}x_j\) , \(c_j^{-1}x_j\) , respectively with positive \( N, \beta , p_j, c_j\) , \(j=1,\ldots ,n\) . The corresponding stationary distributions are the multinomial distribution with the probabilities \(\{\eta _i\}\) , \(\eta _i{\mathop {=}\limits ^{\text {def}}}p_i/(1+\sum _{j=1}^np_j)\) and the negative multinomial distribution with the probabilities \(\{c_i\}\) , respectively. The polynomials, depending on \(n+1\) parameters ( \(\{p_i\}\) and N) and ( \(\{c_i\}\) and \(\beta \) ) satisfy the second order difference equations with the coefficients \(B_j({\varvec{x}})\) and \(D_j({\varvec{x}})\) \(j=1,\ldots ,n\) , which are the straightforward generalisation of the difference equations governing the single variable Krawtchouk and Meixner polynomials. The polynomials are terminating \((n+1,2n+2)\) hypergeometric functions of Aomoto-Gelfand. The bivariate Rahman polynomials are identified as the dual polynomials with a special parametrisation.