There has been growing research interest in the field of nanoscale thermal transport over the past two decades due its importance to a variety of fascinating applications, such as waste heat control, improved electronic functionality, and phononics building blocks. Much of this focus has been on solid-state systems for which advanced experimental characterizations and measurements are readily available. Molecule-based systems, which in principle exhibit no less structural richness than solid state systems and may show excellent energy transport capabilities, have been largely ignored until recently. This is mostly because of the difficulties associated with measuring heat transport on the molecular scale. However, a few recent experimental breakthroughs have brought molecular energy transport process into the spotlight, and at the same time established measurement techniques that can be tested, verified, and explained using theoretical tools. This dissertation examines and explores theoretical approaches for modeling heat transport in molecular systems. Specifically, we have developed a stochastic nonequilibrium molecular dynamics (MD) method which mimics the experimental setting of substrate-bridge-substrate structure, i.e., a molecular junction. We incorporate this approach, along with a quantum Landauer's formalism, into the open-source molecular simulation package--GROMACS, so that it can be applied to molecular systems with different topologies and thermal environments. Our simulations of heat conduction in hydrocarbon-based single molecule junctions yield excellent agreement with the recent state-of-the-art experimental data. Within the capacities of the new method,we have also investigated phononic interference effects in the heat conduction characteristics of benzendithiol molecules. Using the methods developed in this dissertation, we have mapped, for the first time, thermal fluxes down to the atomistic level. In the context of phononic energy transport, we develop a simulation method that integrates quantum effects into classical MD. This hybrid method, once fully implemented, will compensate for the disadvantages of classical approaches at low temperatures and for the difficulty in treating anharmonicites in Landauer-type quantum transport calculations. This method will improve the predictive power of classical heat conduction simulations. The second part of this dissertation explores an intriguing energy transport channel that has been newly discovered termed electron-transfer-induced heat transport (ETIHT), which is distinct from traditional heat transfer mechanisms that rely purely onmolecular vibrations. We construct a theoretical model that combines the two energy transport channels (ETIHT and phononic) into one general model and then we show analytically under certain parametric thresholds (e.g. reorganization energies) that ETIHT dominates while other conditions may magnify the phononic contributions. Although the work in this part of the thesis is currently purely theoretical, it may provide useful insights into future organic molecular thermoelectric devices.