Two methodologies for performing decay heat analysis with Monte Carlo simulations were developed and implemented on a representative nuclear thermal propulsion (NTP) system. This paper presents the underlying theory, discusses the methodology, and states the key results. This work investigated the importance of utilizing a time-dependent Q-value for fission in NTP systems due to their short burn time. Two approaches for deriving the Q-value were taken: one based on deconvolving the fission rate from the reactor power to yield the rate of fission energy deposition, and the other based on the convergence of the fission product decay power during a long burn. The fission product decay power method is hypothesized to be the more accurate representation of an NTP system as it captures more of the underlying physics occurring during burnup, such as fission product transmutation. The calculated Q-values were employed to derive decay power profiles that were compared to the current state-of-the-art NTP decay power model. According to these new models, it is shown that the cooling requirements for decay heat removal calculated with the state-of-the-art model differ from the developed methods by as much as 23.3%. There exists a need to experimentally validate, and by extension improve, the proposed methods to better understand the nature of decay heat production in NTP systems.