Yes, for what I’m trying to do, I need two Neural networks.
The first, \vec{u}^{NN}(t), (represented by chain1) models the solution, \vec{u}(t) to the ODE system. It should have a 1-dimensional input (time) and a two-dimensional output (\vec{u}=(u_1,u_2)^T). The second, F^{NN}[\vec{u}], (represented by chain2) models the unknown ODE system dynamics given by (\beta u_1 u_2, \gamma u_1 u_2)^T. Its input should be 2-dimensional (\vec{u}), and its output should be 2-dimensional ((\beta u_1 u_2, \gamma u_1 u_2)^T).
I agree that it may simply not be possible to do this with NNODE, I just wanted to check here before attempting to build it from scratch.