A comparative study of analytical models of diffuse reflectance in homogeneous biological tissues: Gelatin‐based phantoms and Monte Carlo experiments

Information about tissue oxygen saturation (StO2) and other related important physiological parameters can be extracted from diffuse reflectance spectra measured through non‐contact imaging. Three analytical optical reflectance models for homogeneous, semi‐infinite, tissue have been proposed (Modified Beer–Lambert, Jacques 1999, Yudovsky 2009) but these have not been directly compared for tissue parameter extraction purposes. We compare these analytical models using Monte Carlo (MC) simulated diffuse reflectance spectra and controlled gelatin‐based phantoms with measured diffuse reflectance spectra and known ground truth composition parameters. The Yudovsky model performed best against MC simulations and measured spectra of tissue phantoms in terms of goodness of fit and parameter extraction accuracy followed closely by Jacques' model. In this study, Yudovsky's model appeared most robust; however, our results demonstrated that both Yudovsky and Jacques models are suitable for modeling tissue that can be approximated as a single, homogeneous, semi‐infinite slab.


INTRODUCTION
Tissue oxygen saturation ( 2 ) is a key metric that could be used in surgery to determine tissue viability, tumour localisation, and vessel flow [1][2][3] .There are, however, no currently established, spatially-resolved, intra-operative methods of quantitatively determining  2 .In many clinical applications indocyanine green fluorescence angiography (ICG) is used to display perfusion as blood flow through vessels can be visualised by the fluorescence of the indocyanine dye given to the patient [4] .Whilst this method has shown clinical benefits [5,6] , it cannot quantify the  2 of the tissues being visualised leading to subjective interpretation.There has been increasing promise that diffuse reflectance spectroscopic techniques could be used to obtain  2 intra-operatively.Hyperspectral Imaging (HSI) is one such technique that obtains a diffuse reflectance spectrum for each pixel of an image [7,8] , enabling intra-operative, non-contact, spatially-resolved  2 extraction without the need for a contrast agent.
There have been many proposed methods to extract  2 parameters from tissue diffuse reflectance spectra ranging from the simplest ratiometric two-wavelength methods to more sophisticated analytical models [9] or deep-learning based approaches [10] .Using two or three wavelength models can provide good results but require data to be captured at these precise wavelengths placing large constraints on the measurement devices used and significant model assumptions.Analytical models can be applied to ranges of wavelengths for which they are developed and so these are favoured in this work.Tissue models tend to either be based on modifications of the Beer-Lambert model, or on the diffusion approximation particularly using the Kulbelka-Munk approximation method [9] .Three models that can be applied to homogeneous semi-infinite tissue in the visible wavelength range (450-650nm) include 1) a modified Beer-Lambert model [11] ; 2) a more elaborate Beer-Lambert-based model proposed by Jacques [12] ; and 3) a semi-empirical model based on the Kulbelka-Munk theory as described by Yudovsky [13] .Alternatively, Monte Carlo methods may be used to model tissues for a range of wavelengths.Whilst Monte Carlo is an established approach for modelling light-tissue interaction, it is computationally challenging to apply in an inverse problem setting making the estimation of optical properties from tissue spectra non trivial on the basis of Monte Carlo simulations alone.To address the need for improved inverse modelling utilising Monte Carlo, Inverse Adding Doubling (IAD) [14] has been proposed, however this is also slow iterative process and therefore not as efficient as approaches based on analytical models.A limitation of IAD is that, although a single reflectance measurement can be used to analyse semi-infinite samples with IAD, its primary focus is the extraction of optical properties from slabs of tissue with known thickness using both transmission and reflection measurements.In this work, IAD is used to provide information on optical properties from tissue phantom slabs but is not assessed in the same capacity as the other analytical models listed above.
Since there is no standard method of determining spatially resolved, ground truth, optical properties of tissues, simulations and physical phantoms must be used for testing purposes.Gelatin-based phantoms are commonly used to allow tunable absorption and scattering properties in solid phantoms [15,16] .In this work we compare three major analytical models against both Monte Carlo simulated spectra, and measured spectra from gelatin-based tissue phantoms with known ground-truth constituent quantities.These datasets cover a large parameter range mimicking that expected in biological tissue [17] .Gelatin phantoms are constructed to optically mimic biological tissue absorption and scattering ranges as closely as possible, whilst maintaining welldefined optical properties with known ground truth.We present a first direct comparison in the performance of these models against simulated and measured spectra, and aim to evaluate their potential for accurate  2 extraction.This evaluation is quantified in terms of the forward models (utilising ground truth values within the models to predict the diffuse reflectance spectra), and for the inverse problem solving (evaluating the fidelity of the parameter extraction when fitting to the simulated or measured data).Our results provide a comparison which can be used to inform the translation of these models into clinical applications.

Tissue models
Three analytical models are compared and evaluated in this work: modified Beer-Lambert [11] , Jacques 1999 [12] , and Yudovsky 2009 [13] .Due to the absorbing and scattering nature of biological tissue these models are used to analyse diffuse reflectance spectra which result from propagation of light via many optical paths through the tissue.

Common absorption and scattering model
All three forward models use the wavelength dependent absorption and reduced scattering coefficients (  () and  ′  ()) of the tissue as inputs to compute a diffuse reflectance spectrum.The reduced scattering coefficient comprises of the scattering coefficient (  ()) and tissue anisotropy () as follows: For biological tissue, the reduced scattering coefficient can be well approximated using Mie theory [17] : where  and  are Mie scattering coefficients that range between 8 −1 and 70 −1 , and 0.1 and 3.3 respectively depending on the tissue microstructure [17] .The absorption coefficient of most internal, homogeneous, single-layer tissues is dominated by haemoglobin in the visible region (450-650 nm) [18] and can be modelled using the following equation [13] : Here  2 refers to oxygen saturation,   refers to the total concentration of haemoglobin in whole blood commonly taken as 150 −1 [19] , with  2 denoting the fraction of this that is oxygenated ( 2 ) and the remainder is deoxygenated (), and ranges between 0-100% [13] .  is the volume fraction of tissue occupied by blood which has absorption coefficient  , () and is examined in the range of 0.2-7% [13] , with the remainder of tissue having background absorption  , () [13] .Finally,   2 () and   () denotes the wavelength-dependent extinction coefficients of the chromophores  2 and  which are found in the literature [19] .

Modified Beer-Lambert
The modified Beer-Lambert utilises these inputs to model absorption (()) and diffuse reflectance (()) as follows: Conventionally,   () and  ′  () are quoted in units of cm −1 , however in this model units of mm −1 are required for diffuse reflectance units of % hence a factor of 100 is included. describes a differential path length to account for the variety of photon path lengths through scattering media. is often simplified to be equal to 1 [11] and  ′  () is often modelled as a wavelengthindependent constant [11,20] .To allow for further flexibility in this model, we introduce linear scaling hyperparameters ( 1−3 ), as shown in Equation (5), that can be fitted to Monte Carlo simulations at the refractive index of interest.

Jacques 1999
The Jacques model is also based on the Beer-Lambert model, where an assumption is made that the ensemble of path lengths experienced by photons in a tissue can be approximated by a single wavelength-dependent path length () = ()() where () and () are defined in Equation 6.This results in the following model with the hyperparameters ( 1−3 ) which the authors fit to Adding Doubling simulations [12] .In our work, we refit these to Monte Carlo simulations since these are considered the gold standard optical simulation method and improve the model fitting results.

Yudovsky 2009
The original Yudovsky model presents a complex analytical derivation [13] with a simplified formulation subsequently presented in their Erratum [21] .The model takes the reduced albedo ( ′ ()) as input and returns the diffuse reflectance spectra (()).This provides an easily applicable model with hyperparameters ( 1−6 ) which are quoted for a refractive index () of 1.44.We found that these can be fitted to Monte Carlo spectra to allow for use with other refractive indices.

Fitting model hyperparameters
All analytical models considered use hyperparameters to account for the media's refractive index.A Monte Carlo dataset is generated for each refractive index as in Section 2.2.1.To fit the hyperparameters, the ground truth tissue parameters are inputted into an analytical model for each spectrum in the dataset and a single set of hyperparameters are fitted using a non-linear least squares fitting approach.

Generating reference datasets of diffuse reflectance spectra
In this work, we aim to compare the three models discussed in Section 2.1 against datasets with known ground truth.These datasets of diffuse reflectance spectra are generated either by simulation using Monte Carlo (Section 2.2.1) or measurement of controlled, gelatin-based, tissue phantoms (Section 2.2.2).

Monte Carlo
The models take   () and  ′  () as inputs.These are given for bulk tissue in Equations ( 3) and ( 2) and used to generate Monte Carlo simulated reference spectra.Monte Carlo takes   () as input, not  ′  (), so a random value of anisotropy () is chosen per spectrum between 0.7-0.9 [13]and used in Equation (1) for conversion to   ().The variable parameters are , ,  2 , and   which are bounded to 8-70cm −1 , 0.1-3.3,0-100%, 0.2-7% respectively to mimic biological tissue [13,17] ; a value of each of these variables is randomly selected within these bounds to generate each of the 100 simulated spectra, meaning 100 values of each parameter are sampled.When fitting the inverse models to these spectra for parameter extraction (as described in Section 2.3), the same bounds are imposed on the fitting routine.These simulated spectra are generated using CUDAMCML [22] , which is a GPU accelerated adaptation of the well-established MCML programme [23] .It has been shown that the semi-infinite approximation is valid for thicknesses of above 1cm [24] so 100 spectra were simulated for each refractive index of 1.33, 1.35, and 1.44, by propagating 100,000 photons through the semi-infinite slabs approximated by a thickness of 3cm.The refractive indices were chosen to represent the common phantom and tissue use-cases of these models: 1.33 for use with water-based phantoms, 1.35 for use with gelatin-based phantoms (as in this work), and 1.44 for use with biological tissue.

Controlled gelatin-based tissue phantoms
Here we construct controlled, optical tissue phantoms with well-characterised components.By measuring these phantoms we constructed a dataset of measured spectra with well-defined ground truth which can be modelled using the analytical models.

Phantom composition and synthesis
The dyes acid red 1 (AR1, 210633, Merck, Germany) and acid red 14 (AR14, B22328, Fisher Scientific, England) are chosen to mimic the extinction coefficients of oxygenated and deoxygenated haemoglobin respectively, with a third dye of crystal violet (CV, C6158, Merck, Germany) chosen to investigate the effect of including further chromophores.As in the modelled tissue   () (Equation (3)), the total dye concentration is modelled independently of the relative ratios of each dye.To ensure the absorbance impact of each dye is approximately equal, a factor of 5 3 is included for AR14 and 1 2 for CV computationally by combining with the extinction coefficients to create effective extinction coefficients (   ) that have approximately equal impact.To reflect this experimentally, the concentrations of these dyes are modified by these factors.These effective extinction coefficients calculated from literature values [25] , alongside those of haemoglobin [26] , can be seen in Figure 1 a.
Intralipid is used to modulate the scattering coefficient of these phantoms which can be modelled with a Mie scattering function (Equation (2)) using parameters within the range of tissues.5 volume fractions of intralipid are chosen between 1% and 6% to cover a range of scattering parameters within those seen in biological tissue [17] .
Each phantom consists of 6% gelatin (Type A approximately 175g bloom, G2625, Merck, Germany) by mass and 0.5% of 4% formaldehyde (J60401, Fisher Scientific, England) to increase the melting point of the phantoms and allow their use at room temperature [15] .The remainder of each phantom consists of the dye solutions in a variety of ratios combined with intralipid at a variety of concentrations.Phantoms are constructed with the ratios given in Table 1 , where each configuration is constructed with each intralipid concentration.
For each dye, stock solutions are made at double the required concentrations: 2x, 20x, and 40x in phosphate buffered saline (PBS, P4417, Merck, Germany) to allow for 1:1 dilution with intralipid as the final step.The stock concentrations are multiplied by a factor of 5 3 for Acid Red 14 and 1 2 for Crystal Violet to ensure similar absorbance impact.For each intralipid concentration, the stock solutions are made at double the required volume fractions: 2% -12% to allow for 1:1 dilution with the dye solutions as the final step.The dye solutions combined in the correct ratios are heated to 45-50  C with 12% gelatin (i.e.double the required amount) until solvation.The solution is cooled to below 40  C where 1% formaldehyde (ie.double the quantity) is added and the TABLE 1 Table displaying the ratios of dyes [acid red 1 (AR1), acid red 14 (AR14), crystal violet (CV)] a used for each dye configuration alongside the total dye scaled concentration in arbitrary units.

Total dye scaled concentration
Dye ratio (arbitrary units) solution is combined in equal quantities with the solution with double the desired intralipid concentration.This final solution now has the intended concentration of all constituents and is poured into a silicon mould in an ice bath, as seen in Figure 1 b, to ensure homogeneity in the final phantom by rapidly setting prior to any density separation.Once cooled to below 10  C these are placed in a fridge at 4  C to fully set for at least 7 days before measurement.Each dye solution is also combined with gelatin with a 0% intralipid solution for absorbance measurements.Additionally, each intralipid solution is combined with a gelatin solution with no dyes to allow analysis of scattering due to intralipid concentration.
Moulds with a depth of 1 cm are used to allow measurement of diffuse reflectance within a semi-infinite regime [24] .Monte Carlo simulations were used to confirm that this thickness was sufficient to be within the semi-infinite regime.Some additional phantoms are constructed with a depth of 5 mm to allow for sufficient transmission for total transmittance measurements required for inverse adding doubling (IAD) analysis [14] .

Spectral measurements
Measurements of these phantoms are taken using a PerkinElmer Lambda 750s spectrophotometer.This dual-beam integrating sphere spectrophotometer allows for precise and accurate measurements of diffuse reflectance, total reflectance, total transmittance, and absorbance.Absorbance for each dye solution and pure dye gelatin phantom is measured with either PBS or pure gelatin as a reference respectively.The pure gelatin absorbance is also measured.All 1cm phantoms are used for diffuse reflectance measurements.A subset of these are also made into 5mm phantoms alongside pure intralipid gelatin 5mm phantoms, which are measured for total reflectance and total transmittance.All configurations are depicted for our system in Figure 2 .

Inverse adding doubling (IAD)
IAD is used to obtain   () and  ′  () from total reflectance and total transmittance measurements [14] .We use this with the dual beam spectrophotometer setting and an incidence angle of 8  as per the experimental set-up in Figure 2 d, with  fixed to 0.8 for all calculations since it was found not to change the results.Since the outputted optical properties can contain a significant number of wavelengths without convergence, a Mie scattering curve is fitted to the output  ′  () and a second stage of IAD is run with the  ′  () fixed to this spectrum.This returns outputs which are very similar but removing noise or errors.Further details on this two-stage IAD fitting approach can be found in previous work [27] .
In order to obtain accurate optical properties from IAD a highly accurate sample depth must be provided.To measure these accurately for the highly compressible gelatin phantoms, a CT scan is taken of each phantom using a Small Animal Radiotherapy System (SmART+, precision X-ray).A mean of 10 digital measurements is used for each phantom and an example of these measurements is shown in Figure 1 c.Finally, a commercial calibrated phantom (BioPixS) with validated ground truth optical properties at three wavelengths is used to confirm that our output   () and  ′  () are correct.We highlight that only the single-stage IAD algorithm can be used here as a Mie scattering curve does not accurately represent the scattering of this sample.The BioPixS phantom is measured on three days to demonstrate the accuracy and reproducibility of IAD-outputted optical parameters.

Modelling tissue phantom optical properties
The   () of the two-dye configuration phantoms can be modelled as in Equation (8) where the background spectrum is determined by measuring the absorption of pure gelatin solution, whereas the scattering is modelled based on the intralipid concentration according to trends determined using IAD.
In this equation   is the total concentration of dye independent of dye identity, 1 and 14 are the relative ratios of each dye,  1,  () and  14,  () are the effective extinction coefficients of each dye, and  , () is the measured background absorbance from gelatin.This can be adapted to a three-dye configuration as follows: Where  is the relative ratio of crystal violet and   ,  () is the effective extinction coefficient of this dye.

Evaluation of model performance
In this work, the above models are evaluated against spectra with known ground truth.These reference spectra are either found by Monte Carlo simulation or measurement of controlled phantoms, however the evaluation of each is broadly the same.The forward models are evaluated in terms of the fit of their predicted spectra, and the inverse problem solutions (parametric fits) are investigated in terms of the quality of parameter extraction.Diffuse reflectance spectra can be predicted by inputting ground truth values into each forward tissue model and compared to corresponding reference spectra.The Normalised Root Mean Squared Error (), defined in Equation (10), is calculated to quantitatively evaluate the similarity between the reflectance predicted by the forward model and the associated reference reflectances.In the context of this work the wavelength ranges considered for this metric correspond to those ranges used for fitting the inverse models.
Here the modelled spectrum   at any given wavelength  is evaluated against the intensity at the same wavelength in the reference spectrum   using the normalized root mean squared error () calculated across all Λ wavelengths.These forward model reflectance spectra are also plotted against the reference reflectance spectra and a regression line calculated.This is evaluated using the Pearson correlation coefficient () and the p-value () for a hypothesis test with the null hypothesis that there is no correlation.Using a 95% confidence interval, significant  is less than 0.05, whereas a strong correlation is demonstrated by  close to 1.This gives an indication of similarity in shape, while disregarding offsets.
The inverse problem solutions can also be evaluated by determining how well they recover the ground truth input parameters using a non-linear least-squares fitting approach (using SciPy v1.10.0 scipy.optimize.least_squaresfunction).The cost function for this is shown in Equation (11) where () is the reference spectrum and  is the forward model function.This equation is modified to fit the parameters used to calculate   () and  ′  () directly.
All models are fitted by only considering wavelengths up to 600nm as this is the region the Yudovsky model has been shown to be effective by the authors [28] .When fitting to measured phantom spectra, only wavelengths up to 575nm are considered due to the forward model predicted spectra being highly unreliable after this point for all three models as discussed in Section 3.2.The quality of parameter recovery is examined by calculating the correlation between the fitted and ground truth parameters (using SciPy v1.10.0 scipy.stats.linregressfunction).This is evaluated with the Pearson correlation coefficient () and the p-value () as above.The quality of parameter recovery is also evaluated using absolute percentage errors ( ) as given by Equation (12) where  is the extracted parameter and  is the ground truth parameter.The median and inter-quartile range of these parameters are presented for each dataset.
These evaluation methods can be done for quantitative or relative spectra.Here relative spectra are defined as mean normalised spectra.This allows parameters to be extracted considering the shape of the spectrum but without considering absolute intensity.This is investigated as it is simpler to capture relative data in clinical environments [29] .

Monte Carlo
Each model has hyperparameters which are fitted to Monte Carlo datasets for refractive indices 1.33, 1.35, and 1.44.These are listed in Appendix Table 1 alongside any literature hyperparameters [12,13] .It should be noted that the literature hyperparameters for Jacques could not be directly replicated by fitting to our Monte Carlo simulations.Refitting these hyperparameters leads to the mean (± standard deviation)  of the  = 1.33 dataset improving from 0.080(±0.056) to 0.021(±0.028).In contrast, Yudovsky's literature hyperparameters are similar to our refitting.For Yudovsky, despite the expected equivalence, we observed a small discrepancy in fitting quality between the "extensive model" [13] and the "simplified model" described in their Erratum [21] .This can be seen for a refractive index of 1.44 in Appendix Figure 1 , which matches the results quoted in their Erratum [21] .The "simplified model" improves the mean (± standard deviation)  of the  = 1.44 dataset from 0.050(±0.014) to 0.010(±0.003).For this reason the Erratum model is used for this work.
Example data comparing the analytical models to Monte Carlo simulations can be seen in Figure 3 .For each model, spectra are generated with each input parameter set from the Monte Carlo dataset and  is calculated per spectrum.The mean (± standard deviation) of these  per model per refractive index are shown in Table 2 .Only wavelengths until 600nm are considered for these metrics as in Section 2.3.An example spectrum with each of the forward models and modelled spectra using ground truth parameters are shown in Figure 3 a for a refractive index of 1.44.A regression line calculated between each forward model dataset and the Monte Carlo simulated dataset for each refractive index and the correlation coefficient () can be seen in Table 2 .Finally, the inverse problems with our fixed hyperparameters are fitted as in 2.3 to this same Monte Carlo dataset to recover the  2 ,   , , and  tissue parameters.A linear regression is fitted between these retrieved values and the ground truth counterparts.The Pearson  values of this, and the median (inter-quartile range) absolute percentage errors between the fitted and ground truth tissue parameters are shown in Table 3 for a refractive index of 1.44 with further parameters and refractive TABLE 2 Mean (± standard deviation)  (3.d.p.) between each forwards spectrum from each model and each of 100 Monte Carlo simulated spectra using the same ground truth variable parameters for each refractive index dataset and each analytical model.This is presented with the Pearson  (bold if Pearson  < 0.05) for the linear regression between all forwards spectra against Monte Carlo simulated spectra for each refractive index dataset and each analytical model.All metrics are evaluated for the wavelength region of 450-600nm.

Model
Refractive index   indices shown in Appendix Table 2 .An example of the fitted parameters compared to the ground truth parameters is shown for  2 at a refractive index of 1.44 in Figure 3 b.

Gelatin-based tissue phantoms
Each dye absorbance is measured in both an aqueous and gelatin based solution with effective extinction coefficients calculated in each case.Whilst the aqueous measurements closely match the expected    () from literature [25] , the gelatin measurements show a shifting of peaks, as seen in Figure 4 a.This is likely due to interaction of the dyes with gelatin altering their interactions with solvent as has been noted with other dyes [30] .Since these peak shifts are seen in all subsequent data, these gelatin-based    () are used in the model analysis.A pure gelatin solution absorption is measured at the concentration used in each phantom.This is used as a background   (), as seen in Appendix Figure 2 , to account for any absorption not due to dyes.We further assume that intralipid is a purely scattering medium.This is considered a reasonable assumption since the IAD returned   for the purely intralipid phantoms are similar to the pure gelatin background   seen in Appendix Figure 2 .IAD [14] is used to analyse phantoms of 5mm depth with or without dyes in all intralipid concentrations.The  and  Mie coefficients fitted to the scattering in these cases are plotted against intralipid concentration.It was found that  had no significant trend with intralipid concentration and so a median value of 0.98 was used in this work.However, a clear trend was identified for  as can be seen in Appendix Figure 3 .This trend is described in Equation ( 13), where  is the intralipid concentration and  0 is the reference concentration of 1%v/v: A Pearson correlation coefficient () of 0.99 and p of 0.00 was observed, therefore this is used throughout this work to model scattering based on Intralipid concentration.The calculated   () using    () displayed in Figure 4 a is compared to the IAD outputs of phantoms including dyes with each intralipid concentration.An example of this is seen in Figure 4 b.This shows the overall spectrum is accurate, however the peak for CV (approx 600nm) appears shifted in wavelength, likely due to interactions with intralipid which cannot be captured in absorbance measurements.
Using a refractive index of 1.35 [15] diffuse reflectance spectra are generated from each model using   () and  ′  () from Equations ( 8) and ( 2) using the trend in Equation ( 13) and the median  value for each 2-dye configuration and intralipid concentration.The forward Yudovsky, Jacques, and Modified Beer-Lambert models in Section 3.1 are compared to diffuse reflectance measurements for each phantom.An example of this can be seen in Figure 5 a for quantitative spectra or Figure 5 c for mean normalised relative data.The data is normalised using the mean of the wavelength range 450-575nm as beyond  4 a shows    () calculated for each dye [acid red 1 (AR1), acid red 14 (AR14), crystal violet (CV)] measured in gelatin (dashed) or PBS (dotted) demonstrating the shift in peaks.Figure 4 b shows a calculated (black solid)   () for one 3-dye configuration compared to the IAD output   () (coloured dashed) from measurements of phantoms in this configuration with a range of intralipid concentrations.
this region the models appear to fit poorly as seen in Figure 5 b.The mean  (± standard deviation) between each spectrum and the measured spectrum, in the region 450-575nm, for this dataset can be seen in Table 4 .The Pearson correlation coefficient  between the forwards spectra compared to the measured spectra for each model is shown in Table 4 and further parameters in Appendix Table 3 .Whilst the overall magnitude of these is higher, the trend in quality of fit is similar as to the Monte Carlo simulations.
The inverse problems are solved by non-linear least squares approaches using measured spectrum as reference observations.The recovered parameters are correlated to the ground truth parameters.The associated errors and correlation coefficients can be seen in Table 4 (and further in Appendix Table 3 ) and an example of the fitted parameters compared to ground truth parameters is shown for 1 in Figure 5 d for fits to quantitative spectra or Figure 5 e for fits to relative spectra.Finally, some 3-dye configurations are investigated with the addition of CV.Some examples of both quantitative and relative spectral reconstruction using the ground truth parameters in the forward Yudovsky model for a single 3-dye configuration at a range of intralipid concentrations can be seen in Appendix Figure 4 .Despite considering the shift in CV peak due to gelatin, there appears to be a further shift, likely due to interactions with intralipid.This can be seen as an offset between the measured and reconstructed spectra in Appendix Figure 4 and appears to change with intralipid concentration.Due to the limited ground truth parameter range for these 3-dye configurations, only median percentage errors are presented for the inverse problem solution performance in Appendix Table 4 alongside their inter-quartile range.These show a dramatic loss in accuracy of 1 and 14 recovery, despite good  recovery, likely due to the imperfect  absorbance modelling distorting the spectrum.

DISCUSSION AND FUTURE WORK
When compared to Monte Carlo simulations, the semi-empirical Yudovsky model performs best in terms of fit and parameter extraction, followed closely by the Jacques model.The Modified Beer-Lambert model, however, performs much more poorly.It is also noted that the Jacques model fits Monte Carlo simulations significantly better after refitting the hyperparameters, however the Yudovsky model performs well using literature values provided that the simplified erratum model is used.
When modelling tissue phantom spectra, it is important to have some prior knowledge of the phantoms.This allows accurate   () and  ′  () calculation, for example utilising the shift in dye peaks from their literature data, background gelatin absorbance, and the trend in Mie scattering parameters with intralipid concentration.Some variation is seen in IAD which may be due to inaccuracies in sample depth measurement or variations in spectral measurements.This suggests that it gives an indication of the optical properties of a sample, however it may not be precise.
When modelling the measured tissue phantom diffuse reflectance spectra using the forward models, Yudovsky performs best followed closely by Jacques, with Modified Beer-Lambert performing significantly worse.The  calculated for relative data is an improvement on those calculated for quantitative data, suggesting that these models are reproducing the shape of the spectra well even when relatively constant offsets are present.When introducing a third dye ( ) the absorbance peak of this dye is noticeably wavelength shifted in experimental measurements compared to the theoretical spectra which then appears distorted.This shift appears to change with varying concentrations of intralipid suggesting some interaction there.
When fitting the models to experimental data for parameter recovery, Yudovsky and Jacques produce excellent parameter recovery for 1 and 14 in both quantitative and relative regimes.For   and  no models were able to recover the parameters well, whereas all parameters were recovered well by Yudovsky and Jacques when fitted to Monte Carlo data.This is suggestive of some overlap in the effects of   and  on the spectra.In the three-dye configuration  is recovered reasonably, TABLE 5 The Pearson  (bold if  < 0.05) of the linear regression line between the fitted parameters and their ground truth displayed with their median (inter-quartile range) absolute percentage errors for each variable when extracted by fitting Yudovsky 2009 (Y), Jacques 1999 (J), or Modified Beer-Lambert (BL) for both quantitative and relative data.All metrics calculated for the wavelength range 450-575nm and presented to 3s.f.despite the peak being outside of the region considered for fitting.The shift in the  peak, however, distorts the rest of the spectrum leading to very poor recovery of all other parameters.This demonstrates the importance of prior understanding of the chromophores for any model to recover parameters accurately.This work is associated with some limitations.Firstly, whilst every effort was made to synthesise optical phantoms resembling biological tissue, measurements of true biological tissue could not be used for this work due to a lack of reliable ground truth parameters in tissue.Similarly, haemoglobin chromophores could not be utilised in the phantoms due to their oxygen sensitivity which could not be controlled within the spectrophotometer measurement set-up.This led us to rely on stable synthetic dyes instead.Secondly, the shift in absorbance of the third chromophore made it difficult to model thereby limiting the investigation of the impact of additional chromophores.Thirdly, whilst the IAD method produces results similar to our ground truth, there remain some differences between IAD output and the expected ground truth.This could impact the quality of the intralipid scattering model used in this work.Fourthly, these models assume a planar surface for measurement, however the surface may deviate from this.In these cases, Spatial Frequency Domain methods may be an alternative to obtain accurate   and  ′  values for tissues [31] .Finally, the precise measurement of extinction coefficients of haemoglobin in biological tissue is not possible due to the scattering nature of these media.However, it is shown in our work that the medium chosen for the dye can have an impact on the extinction coefficients and therefore the quality of the models.For this reason, the optical properties of tissue must be well defined to mimic the quality of results in this work.
Future work should include investigation of these models for use with hyperspectral imaging cameras which may not have as densely sampled spectra as these have been shown to be likely methods of obtaining these diffuse reflectance spectra intraoperatively. [32].There could also be greater investigation of including further chromophores whose absorbance behaviour is well modelled and optimisation of the fitting region.Finally, two-layer models should also be considered for tissue structures which exist in layered configurations, such as skin.
In conclusion, the Yudovsky single layer model works well for modelling of tissue that can be approximated by a semi-infinite, homogeneous slab.Jacques is also able to well approximate this with a simpler model.All models require appropriate prior knowledge of the   and  ′  properties of the tissues to allow for a good quality of parameter recovery.

A IAD VALIDATION
To demonstrate that optical properties output by IAD are reproducible and accurate, the output properties from IAD analysis of measurements of the same BioPixS phantom on three days are shown in Figure A5 a and Figure A5 b.Whilst there is some variability, the results show that the optical properties returned by this method accurately reflect the ground truth properties.

FIGURE 1
FIGURE 1 Figure regarding phantom development.Effective extinction coefficients of acid red 1 (AR1), acid red 14 (AR14), and crystal violet (CV) displayed with the extinction coefficients of oxygenated (HbO 2 ) and deoxygenated haemoglobin (Hb) (1 a).Image of 1cm depth (left) and 5mm depth (right) phantoms setting in ice bath (1 b).An example of a CT scan of a tissue phantom and 10 digital measurements taken of the phantom thickness calculated to be 6.11mm (±0.11mm) (1 c).

FIGURE 2
FIGURE 2 Depictions of measurement set-up for spectrophotometer measurements of gelatin based tissue phantoms of absorbance (2 a), total transmittance (2 b), diffuse reflectance (2 c), and total reflectance (2 d) where an 8  wedge is used to ensure the specular reflectance is also included.

FIGURE 3
FIGURE 3 Figure depicting evaluation of forward models (3 a) and inverse problem solutions (3 b). Figure 3 a depicts an example of the predicted spectra from each forward analytical model: Yudovsky 2009 (green dotted), Jacques 1999 (orange dashed), and Modified Beer-Lambert (purple dot-dashed), using ground truth variables for a refractive index of 1.44 compared to that predicted by Monte Carlo (blue solid).Figure 3 b shows an example of difference in quality of parameter recovery by fitting each inverse analytical model to Monte Carlo simulations in the wavelength range of 450-600nm: Yudovsky 2009 (green +), Jacques 1999 (orange ×), and Modified Beer-Lambert (purple * ), and their associated trend lines for a refractive index of 1.44 for  2 .
FIGURE 3 Figure depicting evaluation of forward models (3 a) and inverse problem solutions (3 b). Figure 3 a depicts an example of the predicted spectra from each forward analytical model: Yudovsky 2009 (green dotted), Jacques 1999 (orange dashed), and Modified Beer-Lambert (purple dot-dashed), using ground truth variables for a refractive index of 1.44 compared to that predicted by Monte Carlo (blue solid).Figure 3 b shows an example of difference in quality of parameter recovery by fitting each inverse analytical model to Monte Carlo simulations in the wavelength range of 450-600nm: Yudovsky 2009 (green +), Jacques 1999 (orange ×), and Modified Beer-Lambert (purple * ), and their associated trend lines for a refractive index of 1.44 for  2 .

Figure 3 b
FIGURE 3 Figure depicting evaluation of forward models (3 a) and inverse problem solutions (3 b). Figure 3 a depicts an example of the predicted spectra from each forward analytical model: Yudovsky 2009 (green dotted), Jacques 1999 (orange dashed), and Modified Beer-Lambert (purple dot-dashed), using ground truth variables for a refractive index of 1.44 compared to that predicted by Monte Carlo (blue solid).Figure 3 b shows an example of difference in quality of parameter recovery by fitting each inverse analytical model to Monte Carlo simulations in the wavelength range of 450-600nm: Yudovsky 2009 (green +), Jacques 1999 (orange ×), and Modified Beer-Lambert (purple * ), and their associated trend lines for a refractive index of 1.44 for  2 .

FIGURE 5
FIGURE5 Example of measured (blue) spectrum for one configuration at an intralipid concentration of 3.5% compared to the respective spectra generated from the Yudovsky 2009 (green dotted), Jacques 1999 (orange dashed), and Modified Beer-Lambert (purple dot-dashed) models using the ground truth parameters for quantitative (5 a) or relative (5 c) spectra (mean normalised in the range 450-575nm).The residuals between the modelled quantitative spectra and the measured spectra are shown in 5 b.Examples of the parameter recovery from the measured phantom spectra for the AR1 proportion by Yudovsky 2009 (green +), Jacques 1999 (orange ×), and Modified Beer-Lambert (purple * ) models compared with the ground truth parameters (blue solid) for fits to quantitative (5 d) or relative (5 e) spectra in the wavelength range 450-575nm.

1 FIGURE 3
FIGURE 3 Trend in  with intralipid concentration alongside linear regression.The trend in  with intralipid concentration is also plotted alongside the median  value.

FIGURE 4
FIGURE 4 Examples of Measured (blue) spectrum for one 3-dye configuration at a variety of intralipid concentrations compared to the respective spectra generated from the Yudovsky 2009 (green dotted) model using the ground truth parameters to generate quantitative (4 a) or relative (4 b) spectra.

Figure
Figure A5 a shows the IAD outputted   and Figure A5 b shows the IAD outputted  ′  compared to the well-characterised ground truth values for a BioPixS optical phantom measured on three separate days.

TABLE 3
The Pearson  (bold if  < 0.05) of the linear regression line between the fitted tissue parameters and their ground truth displayed with their median (inter-quartile range) absolute percentage errors ( ).This is shown for each variable and for a refractive index of 1.44 when extracted by fitting Yudovsky 2009 (Y), Jacques 1999 (J), or Modified Beer-Lambert (BL) to the Monte-Carlo dataset in the wavelength range 450-600nm.All presented to 3s.f.

TABLE 4
mean (± standard deviation) (3.d.p.)between the modelled spectrum using the ground truth parameters and the measured spectrum of a tissue phantom for each analytical model for either quantitative or relative spectra.Alongside this is the Pearson  (bold if  < 0.05) for each model for all forwards spectra compared to the measured spectra.All metrics evaluated for the wavelength range 450-575nm

TABLE 1
Table displaying single layer model hyperparameters fitted to Monte Carlo datasets of each refractive index (black) with available literature values (blue).Example of difference in forwards model similarity to Monte Carlo simulation (blue solid) between the literature extensive Yudovsky 2009 (red dashed) and simplified Yudovsky 2009 (green dotted) models using ground truth variables for a refractive index of 1.44.

TABLE 2
The gradient , offset , Pearson , and  values of the linear regression line between the fitted tissue parameters and their ground truth displayed with their median (inter-quartile range) absolute percentage errors.This is shown for each variable and for each refractive index dataset when extracted by fitting Yudovsky 2009 (Y), Jacques 1999 (J), or Modified Beer-Lambert (BL) to the Monte-Carlo dataset.All presented to 3s.f.Background   using absorption of pure gelatin solution compared to the IAD returned   of the phantoms containing intralipid but no dyes.

TABLE 3
The regression gradient , offset , Pearson , and  values, and the median (inter-quartile range) absolute percentage errors for each variable when extracted by fitting Yudovsky 2009 (Y), Jacques 1999 (J), or Modified Beer-Lambert (BL) to measured tissue phantom spectra for both quantitative and relative data.All presented to 3s.f.

TABLE 4
The median (inter-quartile range) absolute percentage errors for each variable when extracted by fitting Yudovsky 2009 (Y), Jacques 1999 (J), or Modified Beer-Lambert (BL) to measured 3-dye tissue phantom spectra for both quantitative and relative data.All presented to 3s.f.