Photo-acoustic Tomography (PAT) is a hybrid imaging method which simultaneously takes advantage of a high spatial resolution attributed to ultrasound and a rich contrast provided by optics. In PAT, the tissue is irradiated by electromagnetic waves in the visible or near-infrared ranges. A portion of energy from the emitted pulses is absorbed, and induces local pressures which propagate outwards as acoustic waves, and are measured in time by ultrasound detectors. The objective is a reconstruction of a distribution of optical absorption coefficient from the measured data. This inverse problem involves two steps, the first of which is a reconstruction of the initial pressure distribution from the measured data (PAT), and the second step is a reconstruction of optical absorption coefficient from the initial pressure distribution (quantitative PAT). Variational approaches are one of the robust methods for PAT, and are often solved via an iterative implementation of a pair of forward and adjoint operators. In chapter 1, we present a detailed introduction for PAT. In chapter 2, a line search multi-grid method is developed for improving the speed of image reconstruction using variational approaches. We define the forward operator as a linear system of first order wave equations that can be adapted to heterogeneous media, and account for an acoustic attenuation following a frequency power law. We derive the adjoint of this operator using an adjoint-then-discretise method. We numerically approximated the forward and adjoint operators using a k-space pseudo-spectral method. In chapter 3, a continuous adjoint was derived for PAT of the brain. We considered a forward operator that describes the propagation of acoustic waves for linear isotropic heterogeneous and lossy elastic media with an acoustic attenuation following a frequency power law. Using a k-space pseudospectral method for an approximation of the forward and adjoint operators, we analytically show that the derived continuous adjoint matches an associated algebraic adjoint. In chapter 4, we solve a single-stage problem of quantitative PAT for realistic acoustic media. The optical portion of the forward operator is defined using a diffusion approximation (DA) model, and the acoustic portion is described in the same way as chapter 2. We developed two inexact Newton approaches for direct QPAT. Using these approaches, an associated nonlinear objective function is minimised as a sequence of linearised problems using matrix-free Jacobian-based methods.