## 2 Materials and Methods
### Tomographic Kernel  
The kernel is divided into two sections: (i) Forward problem and (ii) Inverse problem. For the forward problem, analytical solutions to the Diffusion Equation under Extrapolated Boundary Conditions for a semi-infinite medium have been computed both for the homogeneous case and for the heterogeneous case. In the latter case, a perturbative solution under the Born approximation has been used over a voxelated volume [21]. Notwithstanding the large availability of more sophisticated models, such as finite-elements (FEM) [22] or Monte-Carlo [23; 24], the use of an analytical solution has the advantage to be quite faster than numerical methods, and to be independent of the temporal discretization. After the computation of the forward time-resolved curves, they are converted to photon-counts, accordingly to the desired count-rate, responsivity of the system, measurement time, source/detector areas, and source power. More in detail: the fluence is translated to actual photon counts, considering an average injection power of 1 mW at 800 nm. The acquisition time is set to 1 s, while the detectors are simulated with a 7 mm\({}^{2}\) active area and a 0.04 quantum efficiency. The simulated values represent a realistic estimate based on actual system setups. No optical losses are considered, since the detector is supposed to be in contact with the medium with full acceptance angle. These figures are used to estimate the maximum number of detected photons. Then, to account for the maximum count rate sustainable for the detection electronics, that number is clipped to an upper limit (typically, this value is \(10^{6}\) counts per second). If desired, an experimentally-derived Instrumental Response Function (IRF) can also be convolved to all the curves for simulating a specific detector. Afterwards, multiplicative Gaussian noise (0.1%) is added to take into account amplitude fluctuations and in general other instabilities of the system. Finally, Poisson noise is added to simulate the shot noise associated with photon detection. The two types of noise simulate the modification of an ideal distribution of time of flight due to the non-ideal experimental setup (Gaussian noise) and photon detection (Poisson noise). More details can be found in [5]. In particular, measurements with a gated detector can be also simulated [25]. Concerning the reconstruction of the absorption coefficient (\(\mu_{a}\)) in the volume, a Jacobian matrix is first computed by solving the following time convolution [21]:  
\[J(\mathbf{r}_{i},t_{k};\mathbf{r}_{s},\mathbf{r}_{d})=\int_{\mathbf{t}_{k}}^{t_{k}+\Delta t}G( \mathbf{r}_{s},\mathbf{r}_{i},t)\bigotimes G(\mathbf{r}_{i},\mathbf{r}_{d},t)\bigotimes h(t; \mathbf{r}_{s},\mathbf{r}_{d})dt, \tag{1}\]  
where \(\mathbf{r}_{i}\) is the position of the \(i^{\text{th}}\) voxel, \(t_{k}\) is the \(k^{\text{th}}\)\(\Delta t\) broad temporal window, \(\mathbf{r}_{s}\) and \(\mathbf{r}_{d}\) the source/detector position, respectively, \(G(\mathbf{r}_{l},\mathbf{r}_{m},t)\) is the [PERSON]'s Function for the semi-infinite medium calculated at the position \(\mathbf{r}_{m}\) for a source placed at \(\mathbf{r}_{l}\), and \(h(t;\mathbf{r}_{s},\mathbf{r}_{d})\) is the IRF related to the source/detector pairs placed at \(\mathbf{r}_{s}\) and \(\mathbf{r}_{d}\). The first convolution between the Green's functions has been analytically solved to speed up the calculation [21]. Figure 1: Schematic algorithm of the multipurpose simulation platform.