We use the one-dimensional TRISTAN-MP particle-in-cell code to model the nonlinear evolution of the whistler heat flux instability that was proposed by Gary et al. (1999, 2000) to regulate the electron heat flux in the solar wind and astrophysical plasmas. The simulations are initialized with electron velocity distribution functions typical for the solar wind. We perform a set of simulations at various initial values of the electron heat flux and $beta_{e}$. The simulations show that parallel whistler waves produced by the whistler heat flux instability saturate at amplitudes consistent with the spacecraft measurements. The simulations also reproduce the correlations of the saturated whistler wave amplitude with the electron heat flux and $beta_{e}$ revealed in the spacecraft measurements. The major result is that parallel whistler waves produced by the whistler heat flux instability do not significantly suppress the electron heat flux. The presented simulations indicate that coherent parallel whistler waves observed in the solar wind are unlikely to regulate the heat flux of solar wind electrons.