We propose novel numerical schemes based on the Boris method in curved spacetime, incorporating both hadronic and radiative interactions for the first time. Once the proton has lost significant energy due to radiative and hadronic losses, and its gyroradius has decreased below typical scales on which the electromagnetic field varies, we apply a guiding center approximation (GCA). We fundamentally simulate collision processes either with a Monte-Carlo method or, where applicable, as a continuous energy loss, contingent on the local optical depth. To test our algorithm for the first time combining the effects of electromagnetic, gravitational, and radiation fields including hadronic interactions, we simulate highly relativistic protons traveling through various electromagnetic fields and proton backgrounds. We provide unit tests in various spatially dependent electromagnetic and gravitational fields and background photon and proton distributions, comparing the trajectory against analytic results. We propose that our method can be used to analyze hadronic interactions in black hole accretion disks, jets, and coronae to study the neutrino abundance from active galactic nuclei.