Adaptive optics observations of the flattened nuclear star cluster in the nearby edge-on spiral galaxy NGC 4244 using the Gemini Near-Infrared Integral Field Spectrograph (NIFS) have revealed clear rotation. Using these kinematics plus 2MASS photometry we construct a series of axisymmetric two-component particle dynamical models with our improved version of NMAGIC, a flexible Chi^2-made-to-measure code. The models consist of a nuclear cluster disc embedded within a spheroidal particle population. We find a mass for the nuclear star cluster of M=1.6^+0.5_-0.2 x 10^7 M_sun within ~42.4 pc (2). We also explore the presence of an intermediate mass black hole and show that models with a black hole as massive as M_bh = 5.0 x 10^5 M_sun are consistent with the available data. Regardless of whether a black hole is present or not, the nuclear cluster is vertically anisotropic (beta_z < 0), as was found with earlier two-integral models. We then use the models as initial conditions for N-body simulations. These simulations show that the nuclear star cluster is stable against non-axisymmetric perturbations. We also explore the effect of the nuclear cluster accreting star clusters at various inclinations. Accretion of a star cluster with mass 13% that of the nuclear cluster is already enough to destroy the vertical anisotropy, regardless of orbital inclination.