Exothermic reactions involving metastable neutrals and ions were recently proposed as sources of hot oxygen atoms in addition to the classical O₂⁺ and NO⁺ dissociative recombination. The Boltzmann equations for thermal and nonthermal populations of O atoms are solved with a Monte Carlo stochastic simulation method. It is shown that the calculated energy distribution functions of O atoms are significantly in nonequilibrium in the transition region between the thermosphere and the exosphere. It is found that the inclusion of additional sources leads to stronger disturbances of the energy distribution function and, as a consequence, increases the nonthermal fraction of hot O atoms. The variation of the vertical distribution of hot O between solar maximum and minimum conditions is also evaluated and shows good agreement with the available experimental evidence.