We present the results of numerical studies of the whistler wave parametric decay instability in the system with the suppressed Landau damping of ion acoustic waves (IAWs) based on the self-consistent Darwin particle-in-cell (PIC) model. It has been demonstrated that a monochromatic whistler wave launched along the background magnetic field couples to a counter-propagating whistler mode and co-propagating ion acoustic mode. The coupling of the electromagnetic mode to the electrostatic mode is guided by a ponderomotive force that forms spatio-temporal beat patterns in the longitudinal electric field generated by the counter-propagating whistler and the pump whistler wave. The threshold amplitude for the instability is determined to be δB_w/B₀ = 0.028 and agrees with a prediction for the ion decay instability: δB_w/B₀ = 0.042 based on the linear kinetic damping rates, and δB_w/B₀ = 0.030 based on the simulation derived damping rates. Increasing the amplitude of the pump whistler wave, the secondary and tertiary decay thresholds are reached, and cascading parametric decay from the daughter whistler modes is observed. At the largest amplitude (δB_w/B₀ ∼ 0.1) the primary IAW evolves into a short-lived and highly nonlinear structure. The observed dependence of the IAW growth rate on the pump wave amplitude agrees with the expected trend; however, quantitatively, the growth rate of the IAW is larger than expected from theoretical predictions. We discuss the relevant space regimes where the instability could be observed and extensions to the parametric coupling of whistler waves with the electron acoustic wave (EAW).