The challenge to fruitfully merge state-of-the-art techniques from mathematical finance and numerical analysis has inspired researchers to develop fast deterministic option pricing methods. As a result, highly efficient algorithms to compute option prices in Lévy models by solving partial integro differential equations have been developed. In order to provide a solid mathematical foundation for these methods, we derive a Feynman-Kac representation of variational solutions to partial integro differential equations that characterize conditional expectations of functionals of killed time-inhomogeneous Lévy processes. We allow for a wide range of underlying stochastic processes, comprising processes with Brownian part, and a broad class of pure jump processes such as generalized hyperbolic, multivariate normal inverse Gaussian, tempered stable, and $α$-semi stable Lévy processes. By virtue of our mild regularity assumptions as to the killing rate and the initial condition of the partial differential equation, our results provide a rigorous basis for numerous applications, not only in financial mathematics but also in probability theory and relativistic quantum mechanics.