<p>Phototherapy using indocyanine green (ICG) has been shown to utilize photothermal and photodynamic effects, and its effectiveness is determined through complex interactions of light transport, heat generation, and reactive oxygen species (ROS) production. In this paper, we report a computational framework for simulating the pulsed 808-nm photoactivation of ICG in a three-dimensional tumor spheroid. Photon transport in the heterogeneous spheroid domain is simulated using GPU-accelerated Monte Carlo photon transport, and the absorbed optical energy is coupled with a transient bioheat equation, an oxygen diffusion–reaction equation, and ROS production equation. The computational domain was defined as a voxelized spheroid with dimensions of 100 × 100 × 100 voxels, which represented the hypoxic tumor core and proliferative tumor areas. The simulation revealed that the fluence distribution is heterogeneous with significant attenuation, with about 70% being reduced within the first 400&#xa0;µm. The maximum temperature ranged from 41.2 ± 0.6 to 53.6 ± 1.4&#xa0;°C based on the pulse duration. The oxygen-dependent ROS generation was found to be highest in the oxygenated proliferative shell, whereas the hypoxic tumor core presented reduced photodynamic effects. The coupled cytotoxicity simulation revealed apoptotic-dominant areas in the oxygen-enriched tumor and necrosis due to heat in the case of increased pulse energy. The pulse duration, oxygen diffusion coefficient, and ICG absorption strength were found to be the most significant parameters affecting the treatment outcome based on the sensitivity analysis. The proposed multiphysics approach can be used as a quantitative tool to analyze laser-induced tumor phototherapy and can be helpful in the computational optimization of the treatment approach. The model is physics-based, and it predicts relative response tendencies to therapy rather than actual biological responses.</p>