  The squared singular values of the product of $M$ complex Ginibre matrices
form a biorthogonal ensemble, and thus their distribution is fully determined
by a correlation kernel. The kernel permits a hard edge scaling to a form
specified in terms of certain Meijer G-functions, or equivalently
hypergeometric functions ${}_0 F_M$, also referred to as hyper-Bessel
functions. In the case $M=1$ it is well known that the corresponding gap
probability for no squared singular values in $(0,s)$ can be evaluated in terms
of a solution of a particular sigma form of the Painlev\'e III' system. One
approach to this result is a formalism due to Tracy and Widom, involving the
reduction of a certain integrable system. Strahov has generalised this
formalism to general $M \ge 1$, but has not exhibited its reduction. After
detailing the necessary working in the case $M=1$, we consider the problem of
reducing the 12 coupled differential equations in the case $M=2$ to a single
differential equation for the resolvent. An explicit 4-th order nonlinear is
found for general hard edge parameters. For a particular choice of parameters,
evidence is given that this simplifies to a much simpler third order nonlinear
equation. The small and large $s$ asymptotics of the 4-th order equation are
discussed, as is a possible relationship of the $M=2$ systems to so-called
4-dimensional Painlev\'e-type equations.
