Time-domain diffuse optical imaging technique for monitoring rheumatoid arthritis disease activity: theoretical development and in silico validation

Objective. Effective early treatment—within 3–5 months of disease onset—significantly improves rheumatoid arthritis (RA) prognosis. Nevertheless, 1 in 3 patients experiences treatment failure which takes 3–6 months to detect with current monitoring techniques. The aim of this work is to develop a method for extracting quantitative features from data obtained with time-domain diffuse optical imaging (TD-DOI), and demonstrate their sensitivity to RA disease activity. Approach. 80 virtual phantoms of the proximal interphalangeal joint—obtained from a realistic tissue distribution derived from magnetic resonance imaging—were created to simulate RA-induced alterations in 5 physiological parameters. TD-DOI images were generated using Monte Carlo simulations, and Poisson noise was added to each image. Subsequently, each image was convolved with an instrument response function (IRF) to mimic experimental measurements. Various spatiotemporal features were extracted from the images (i.e. statistical moments, temporal Fourier components), corrected for IRF effects, and correlated with the disease index (DI) of each phantom. Main results. Spatiotemporal Fourier components of TD-DOI images were highly correlated with DI despite the confounding effects of noise and the IRF. Moreover, lower temporal frequency components ( ⩽ 0.4 GHz) demonstrated greater sensitivity to small changes in disease activity than previously published spatial features extracted from the same images. Significance. Spatiotemporal components of TD-DOI images may be more sensitive to small changes in RA disease activity than previously reported DOI features. TD-DOI may enable earlier detection of RA treatment failure, which would reduce the time needed to identify treatment failure and improve patient prognosis.


Introduction
Rheumatoid arthritis (RA) is a common chronic inflammatory disorder associated with fatigue, pain, depression, and reduced quality of life (Pollard et al 2006, Matcham et al 2014).While the disease has no cure, early treatment within the first 3-5 months of disease onset has been shown to significantly improve patient prognosis (Möttönen et al 2002, Monti et al 2015, van Nies et al 2015).Unfortunately, at least 1 in 3 RA patients fails their first treatment and with current monitoring methods-typically based on clinical examination, patient self-assessment, and laboratory testing (England et al 2019)-it can take 3-6 months to detect this failure (Goekoop-Ruiterman et al 2005, Moreland et al 2012).As a result, these patients not only incur unnecessary expenses and side effects, but also lose the benefits of effective early treatment and are at a higher risk of developing irreversible joint damage (Machold et al 2007).Therefore, there is a need for

Simulations
For reference, tissue-specific changes in physiological parameters associated with RA disease activity are summarized in table 1. Simulations were performed using virtual phantoms based on a realistic tissue geometry derived from a segmented MRI image of a healthy adult finger (Wang et al 2020)-details of the segmentation have been previously described (Ioussoufovitch and Diop 2021).The segmented finger was composed of 7 tissue types (skin, subcutis, tendon, synovial membrane, synovial fluid, cartilage, and bone) and was used as the base tissue geometry for a model with no RA disease activity (i.e.'None').Additional tissue geometries for 'Mild' , 'Moderate' , and 'Severe' models of disease activity were created by linearly transforming transverse slices of the None geometry to incrementally increase the volume of the synovial fluid and synovial membrane (Ioussoufovitch and Diop 2021).This approach was adopted from the work of Dolenec et al (2019), and is a straightforward method of simulating controlled expansion of the synovial fluid and synovial membrane within a desired region of an underlying tissue geometry.These four geometries served as 4 possible states of the 'Joint volume' parameter in table 1.
Optical properties for various RA disease activities were derived by considering changes to the remaining parameters in table 1.Our previous models of None and Severe RA disease activity (Ioussoufovitch and Diop 2021) were used to bound the range of possible values for each physiological parameter (table 2).Note that these values-and the resulting optical properties for each of the 7 tissue types-were amalgamated using a multitude of sources where data was obtained from real tissues.This process and the resulting properties are described extensively in our previous work (Ioussoufovitch and Diop 2021).The value of each parameter for Mild and Moderate disease activity were estimated by linearly interpolating between its values for None and Severe.Specifically, Mild and Moderate RA disease activity were defined as 33 and 66%, respectively, within the 0%-100% disease activity range bounded by None and Severe.
To create one virtual phantom, each of the five physiological parameters in table 1 is assigned a value that corresponds to one of the four disease activities.That is to say that a virtual phantom is created by converting a combination of 5 physiological parameters into a tissue geometry with tissue-specific optical properties (Ioussoufovitch and Diop 2021).To investigate the effects of changes in individual physiological parameters under different baseline disease activities, we assembled 20 sets of virtual phantoms.Here, a set of virtual phantoms is defined as a group of 4 phantoms for which 4 of the 5 physiological parameters are fixed to the same baseline disease activity, while the value of the remaining physiological parameter is sequentially set to one of the four disease activities for each individual phantom.For example, one set of virtual phantoms was composed of four virtual phantoms with different µ SF values ('None' , 'Mild' , 'Moderate' , and 'Severe') while all the other parameters were set to values that correspond to Mild disease activity.Therefore, investigating the effects of changes in the 5 physiological parameters under 4 possible baseline disease activities, resulted in 20 sets of virtual phantoms containing 4 virtual phantoms per set (80 virtual phantoms total).
MCXLAB was then used to simulate TD-DOI of the phantoms (Fang and Boas 2009).An 8 × 8 grid of light source positions was defined on the dorsal side of each phantom's PIP joint over an area (≈12 × 24 mm 2 ; L × W) that was centered around, and encompassed, the joint cavity.Next, by enabling one source at a time, the joint was raster scanned to create an image: for each source position, 10 8 photons were simulated and a single TPSF of all of the photons transmitted through the palmar side of the finger was recorded.Note that in this simple raster scanning approach, the TPSF detected for each source is simply assigned to the image pixel corresponding to the position of that source.As such, no additional reconstruction process is necessary to generate these images.TD-DOI was simulated at 800 nm for all the phantoms, which corresponds to one of the transmittance peaks in PIP joints (Dolenec et al 2019).A subset of simulations were repeated at 650, 750 and 830 nm to investigate the applicability of subsequent analyses to other wavelengths.Note that simulations were conducted with a Mild baseline disease activity since it is most likely to represent the state of early diagnosis prior to treatment (Aletaha and Smolen 2018).
The simulated TD-DOI data was organized as 3D data structures with spatial information on the x-y plane and TPSFs on the z-axis i.e. each image pixel represented a TPSF; hereafter, we refer to this structure as a TPSF image.Similar to Lighter et al's approach, we averaged each image in the transverse direction and normalized it by its total number of photons (Lighter et al 2019).As a result, each of the 80 virtual phantoms was linked to a transverse-averaged, normalized TPSF image (figure 1).
The large amount of multi-dimensional information within each image makes a straightforward comparison between images quite challenging.To address this, we used principal components analysis (PCA) to identify the largest sources of variance between all 80 images-this type of approach is often used to build classifiers for information-rich imaging data (Rodarmel and Shan 2002).PCA was applied using the pca function in MATLAB 2020b (The MathWorks Inc., USA) and by treating the value of each time bin in the TPSF image as a variable.The analysis revealed that the first principal component (PCA-1) accounted for 95% of the total explained variance between all 80 TPSF images.Further, the sensitivity of PCA-1 to changes in individual physiological parameters agreed with the anticipated sensitivity of TD-DOI to each physiological parameter (Ioussoufovitch and Diop 2023).To further quantify this sensitivity, the values of PCA-1 were z-score normalized and the mean change in these values for each set of virtual phantoms was calculated.Importantly, if changes to the disease activity of a parameter did not produce a consistent direction of change in PCA-1 values (e.g.consistent decrease in PCA-1 with worsening disease activity), we assumed that our implementation of TD-DOI was insensitive to RA-associated changes in that physiological parameter.It is noteworthy that in sets where PCA-1 changed consistently it always decreased as disease activity worsened.Given this consistency and the large amount of variance explained by PCA-1 alone, its values will be referred to as a 'disease index' (DI) hereafter, and a decrease in DI will be interpreted as an indicator of worsening RA disease activity.Further, DI will be used as a quantitative reference metric against which more intuitive and generalizable TD-DOI features (e.g.statistical moments, Fourier components) will be compared.
To better mimic experimental conditions, each TPSF image was convolved with a randomly-scaled IRF image acquired using an in-house TD-DOI system.Specifically, each pixel of a TPSF image was convolved with the IRF of the corresponding pixel in the IRF image: where each image has a spatiotemporal dependence (x, t), n is the index of a pixel located at x in the discrete image, t is the photon time-of-flight, and * is the convolution operator.Thereafter, we added Poisson noise (i.e.shot noise)-which is the most relevant type of noise for photon-counting applications such as TD-DOI-to each DTOF image (Becker 2023).Specifically, 10 versions of the resulting DTOF images were created by adding simulated Poisson noise 10 separate times using MATLAB's imnoise function, which resulted in a total of 800 DTOF images (10 for each virtual phantom).Note that, by definition, the amount of Poisson noise added to each time bin of a DTOF image is solely determined by the amount of photons within that time bin.Thus, as long as DTOF images contain an experimentally reasonable number of photons prior to noise addition, the imnoise function will simulate an experimentally reasonable level of Poisson noise.To this end, before applying the imnoise function, each DTOF image was scaled such that the total number of photons in the image approximately matched the number of photons typically seen in measurements acquired using the aforementioned TD-DOI system (10 5 photons).Note that no filtering was used to remove the effects of the noise prior to the data analysis described below.

Feature extraction 2.2.1. Temporal feature recovery
For simplicity, we start by focusing on the relationship between the DTOF, IRF, and TPSF located at a single pixel.As shown by ( 1), it is well-known that the IRF causes temporal dispersion of the DTOF compared to the true TPSF (Liebert et al 2003a), and attempting to directly remove the effects of the IRF by performing a deconvolution is quite challenging since the problem is ill-posed (Hansen 2010, Diop and St Lawrence 2013).Nevertheless, it is possible to extract certain DTOF features from which the effects of the IRF can be removed in a straightforward manner.Common examples are the curves' statistical moments since these are known to have an additive relationship: and where m 1 is the first normalized moment, and V is the second centralized moment (Liebert et al 2003b).
Alternatively, temporal Fourier components can be recovered from each curve based on the well-known convolution property of the Fourier transform: where F t represents the Fourier transform operator in the temporal dimension.We note that in an experiment, it is a scaled version of the IRF that is typically measured (i.e.C 1 IRF where C 1 is an unknown positive scaling factor).In the case of ( 2) and ( 3), moment analysis inherently discards C 1 , and thus it is of no concern.In the case of ( 4), the effect of C 1 can be removed if each of the curves is first normalized by its area under the curve (AUC); thus, only the Fourier components of the normalized TPSF can be recovered once (4) is inverted.With this in mind, we computed m DTOF 1 , V DTOF , and F t [ DTOF]-where the hat symbol indicates a curve normalized by its AUC-from each pixel of the DTOF images.Next, we inverted the relations in (2)-(4) to estimate each pixel's m TPSF 1 , V TPSF , and F t [ TPSF].

Spatial Fourier analysis
Inversion of the relations described in (2)-(4) allows for the recovery of TPSF features at each image pixel.However, previous work has demonstrated the utility of using changes in the spatial distribution of image information to track RA disease activity (Montejo et al 2013a, 2013b, Lighter et al 2019).Therefore, we applied the spatial Fourier decomposition proposed by Lighter et al (2019) for CW-DOI to the simulated TD-DOI data to assess how the spatial Fourier components of various features of the TPSF image vary with changes in RA disease activity.
A one-dimensional discrete FFT was implemented using the MATLAB function fft and applied to the m TPSF 1 and V TPSF images to extract their spatial Fourier components.All components were then normalized by the value of the images' DC component (i.e.first Fourier coefficient) except for the DC component itself.
This ensured that the extracted non-DC spatial Fourier coefficients only depended on the spatial distributions of the features within the m TPSF 1 and V TPSF images (Lighter et al 2019).We noted in section 2.2.1 that the process of creating F t [ TPSF] images requires that the curve at each pixel be normalized by its AUC.Although this process does not affect the relative relationship between the temporal Fourier components at each pixel, it does alter the relationship between the temporal Fourier components of the different pixels of an image i.e. the spatial variation that would have been present in F t [TPSF] images prior to pixelwise AUC normalization.For example, the aforementioned normalization results in the value of the DC temporal Fourier component at each pixel to be equal to 1. Thus, prior to proceeding with the spatial Fourier analysis of the F t [ TPSF] images, an additional processing step was required to re-establish the spatial variation that existed in the images prior to the pixelwise AUC normalization.
First, let us consider the temporal summation of TD images which contain T discrete time bins.In general, the pixelwise relationship between temporally-summed DTOF, IRF, and TPSF images is (5) This relation can be easily confirmed by considering the relationship-shown in (4)-between the DC temporal Fourier component of each curve, which is equal to the total number of photons in that curve.( 5) shows that to recover the spatial variation of interest-contained in ∑ T t=1 I TPSF (x n , t)-we must correct for the effects of the IRF ∑ T t=1 I IRF (x n , t) at each pixel.We also know that the measurement of the IRF must be scaled by some unknown constant C 1 (as discussed in section 2.2.1).Thus, we multiply (5) by C 1 and attempt to isolate As mentioned, C 1 is unknown and cannot be directly separated from the image described by the right side of ( 6).However, we note that this image is just a scaled version of our desired spatial variation.Further, we know that we are interested in TPSF images which have been normalized by their total number of photons (see section 2.1); in other words, where N is the total number of pixels in the image.Thus, if we simply normalize the image on the right side of (6) by its total number of photons- ∑ T t=1 I TPSF (x n , t)/C 1 -we recover our spatial variation of interest: Let ∑ T t=1 I TPSF (x n , t) be a column vector with N rows.Now, let Ψ be a matrix of N rows where each row contains the normalized temporal Fourier components of an image pixel (i.e.F t [ TPSF]).If we scale each row of Ψ using its associated row in ∑ T t=1 I TPSF (x n , t), we can recover the two-dimensional Fourier transform of the normalized TPSF image I TPSF (x, t) as where ⊙ denotes element-wise multiplication, F x is the one-dimensional spatial Fourier transform, and F is the two-dimensional spatiotemporal Fourier transform.It is important to note that even though the result of (8) suggests that Fourier inversion may be used to recover the original normalized TPSF image, this is typically infeasible in an experimental setting due to the distortion of many higher-order Fourier components by experimental noise.

Statistical analysis
Following extraction of each feature's spatial Fourier components from all the DTOF images, the ability of the components to track the expected variation between TPSF images was assessed.Specifically, we performed a correlation analysis between DI (see section 2.1) and the components.Note that the magnitude and phase of the Fourier components were assessed separately and, to avoid duplicate results, symmetric Fourier components were removed prior to analysis.
In addition, we quantified the sensitivity of each component to individual physiological changes.For each set of virtual phantoms, phantoms were treated as the between-subjects variable and repetitions created by addition of Poisson noise (see section 2.1) were treated as the within-subjects variable.Generalized eta squared (η 2 G ) (Olejnik and Algina 2003)-a measure of effect size-was then computed using the values of each component within each phantom set.

Results
Recall that-as described in section 2.1-virtual phantoms were organized into sets with the same baseline disease activity within which the value of only one physiological parameter was altered.Figure 2 quantifies the change in DI caused by these simulated physiological changes under various baseline disease activity conditions.For example, the number in the 2nd row and 4th column shows, within a set of phantoms with Mild baseline disease activity, the average change in DI in response to a change in µ SF for all four levels of disease activity.Note that changes in DI are presented on a z-score normalized scale for ease of interpretation.In 12 of the 20 sets of phantoms, changes in the physiological parameters were not associated with consistent changes in DI (in black, excluded).In particular, changes in SO 2 and H 2 O did not consistently affect DI regardless of the baseline disease activity.Overall, G caused the largest decrease in DI, and this decrease became slightly larger as baseline disease activity worsened.Similarly, the changes caused by µ SF became larger with worsening baseline disease activity; however, DI was not sensitive to changes in µ SF when baseline disease activity was None.BVF caused a relatively modest decrease in DI, but only at the None baseline disease activity (1st row, 3rd column).
Figure 3 shows the results of the correlation analysis between the spatial Fourier components of the TPSF image features and their respective DI.Higher correlation coefficients correspond to spatial Fourier components that accurately tracked DI despite the addition of simulated noise and IRF to the original TPSF data from which the DI was derived.Note that images of F t [TPSF] > 1.6 GHz were excluded from figure 3 since these images did not have a single component with |R| > 0.5.All components derived from images of V also had low correlations (|R| < 0.5) (Mukaka 2012).For m 1 images, the magnitude of the spatial DC and 0.06 mm −1 components had very high (|R| ⩾ 0.9) and high (0.7 ⩽ |R| < 0.9) negative correlations, respectively (Mukaka 2012), and the phase of the 0.06 mm −1 component also had a high positive correlation.Many F t [TPSF] image components were highly correlated with DI.Notably, for both magnitude and phase, the 0.06 mm −1 spatial frequency of most F t [TPSF] images <1.6 GHz had very high correlation coefficients.As well, the 0.13 mm −1 magnitude and DC phase of many components was highly correlated with DI.Considered together, 64 components were strongly correlated with DI (|R| > 0.8) and 7 components had extremely high correlation coefficients (|R| > 0.99); this clearly suggests that it is possible to extract individual features from DTOF images which can accurately track changes in DI.
To further clarify the utility of various components in a clinical setting, we focused on the components' ability to track individual physiological changes (see section 2.3).In particular, we chose to focus on components' ability to track changes in µ SF at a Mild baseline disease activity (figure 4): the sensitivity profile in figure 2 suggests that this is the smallest physiological change detectable with TD-DOI under this baseline disease activity (which is most likely to be clinically relevant) (Aletaha and Smolen 2018).Figure 4 quantifies this ability using the effect size-measured as generalized eta squared (η 2 G )-of changes in µ SF on each component.For clarity, only η 2 G ⩾ 0.26 values-which typically correspond to a large effect size (Cohen 1988, Bakeman 2005)-are shown.Most notably, the phase of the F t [TPSF] 0.06 mm −1 spatial Fourier component was associated with very large η 2 G ; for both the 0.24 and 0.33 GHz F t [TPSF] images, it had the largest effect size (η 2 G = 0.91).

Discussion
Novel monitoring techniques that are sensitive to subtle, subclinical changes in RA disease activity could substantially reduce current delays in detecting RA treatment failure, and help ensure that patients reap the benefits of effective early treatment while avoiding irreversible joint damage.To this end, various DOI techniques have been proposed (Scheel et al 2002, Schwaighofer et al 2003, Hielscher et al 2011, Montejo et al 2013a, 2013b, Lighter et al 2019); however, these techniques have typically focused on making the binary distinction between healthy and inflamed joints.TD-DOI is more information-rich than previously used DOI methods without the practical drawbacks of DOT methods.This work highlights how the additional information content of TD-DOI may be leveraged to track subtle changes in RA disease activity.Importantly, this investigation was conducted in silico to enable systematic assessment-using realistic tissue geometries and physiological properties-of the sensitivity of TD-DOI to RA disease activity with greater precision and flexibility than could be achieved in a phantom or human study.A key advantage of using computational  modelling is the ability to control for the confounding effects of physiological states and geometric variability between patients.This type of approach has been previously used not only to study joints (Milanic et al 2015), but also the brain (Barker et al 2014) and breast (Hassan et al 2023).Starting with a tissue geometry based on an adult PIP joint, we individually simulated 5 major physiological changes associated with RA (table 1) at 4 different levels of disease activity to create 80 virtual finger phantoms i.e. 80 states of disease activity.These phantoms were further grouped into sets of virtual phantoms, each containing 4 phantoms among which only one physiological parameter was changing.Employing a DOI imaging geometry which was previously tested for early RA diagnosis (Lighter et al 2019), we simulated TD-DOI of each virtual phantom, and used PCA to compare the resulting images.PCA-1 accounted for 95% of the variance between the TD-DOI images of the 80 simulated states of disease activity, and its value consistently decreased as disease activity worsened.As such, we selected PCA-1 as a disease index (DI).Then, we used it to estimate the sensitivity of TD-DOI to changes in each of the RA-associated physiological parameters by calculating the mean z-score normalized change in DI within each set of virtual phantoms (figure 2).This analysis revealed that TD-DOI was mostly sensitive to changes in joint volume, and had lower or negligible sensitivity to changes in other physiological parameters.TD-DOI sensitivity tended to increase with worsened baseline disease activity; this was likely caused by increases to the volumes of affected tissues which occur as baseline disease activity worsens.Clinically, starting treatment early-within the first 3-6 months of disease onset-is well known to significantly improve patient prognosis (Möttönen et al 2002, Monti et al 2015), and it is widely accepted that a therapeutic window of opportunity exists during which RA is much more susceptible to treatment.For example, based on data from clinical trials which followed 738 and 533 RA patients, van Nies et al (2015) found that patients who initiated treatment before 15-20 weeks of symptom appearance were disproportionately more likely to achieve remission than those who started treatment later.RA treatment monitoring is typically initiated at a baseline disease activity that is sufficient for diagnosis but still considered to be 'early RA'; thus, TD-DOI sensitivity to changes occurring at 'Mild' baseline disease activity is of particular interest.
While changes to joint volume and synovial fluid produced noticeable changes in TD-DOI data, changes in [H 2 O] only affected joint cartilage so-given the relatively small volume of the affected tissue-it is unsurprising that their impact was difficult to detect with TD-DOI.Physiological changes affecting the optical properties of the synovial membrane (SO 2 , BVF) and synovial fluid (µ SF ) may have been obscured for similar reasons.Nevertheless, TD-DOI still demonstrated moderate sensitivity to changes in µ SF ; this was likely because the optical properties of the synovial fluid change more as a function of RA disease activity than those of the synovial membrane (Ioussoufovitch et al 2020).Though changes in SO 2 were also simulated within the subcutis-which makes up much of the total finger volume-the low blood volume fraction in this tissue (0.8%) resulted in negligible changes to its optical properties as a function of SO 2 (Ioussoufovitch et al 2020).In addition to the tissue-specific factors mentioned above, low sensitivity to SO 2 could have been caused by the choice of the imaging wavelength (800 nm), which is close to the isosbestic point of hemoglobin.To test the generalizability of these findings to other wavelengths, a subset of simulations were repeated at 650, 750 and 830 nm, to sample the spectral range that is typically used for tissue optical imaging.Using the same methodology as in section 2.1, PCA analysis was conducted on the normalized TPSF images at these additional wavelengths; the DI for 650, 750, and 830 nm explained 90%, 92%, and 92% of the total variation in their datasets, respectively.We then performed a linear regression between the DI values derived from the 800 nm data and the DI values derived from 650, 750, and 830 nm data, and found a strong linear relationship (R 2 = 0.99) in all cases with slopes of 1.23, 1.06, and 0.99, respectively.Overall, this suggests that the aforementioned trends are not strongly influenced by the choice of wavelength in 650-830 nm spectral range.Further, as previously suggested by Lighter et al (2019), these results also highlight the potential of lower wavelengths to have slightly higher sensitivity to RA-associated physiological changes.
Though the DI provides an excellent way to distinguish between the simulated images, it is an abstract metric that is challenging to interpret and difficult to generalize to TD-DOI images of different finger geometries or captured under different experimental conditions.Therefore, we focused on extracting more intuitive and generalizable TD-DOI features, and used the DI as a quantitative reference metric to assess their sensitivity to RA disease activity.Further, to confirm the future utility of these features in an experimental setting, we added Poisson noise and convolved each TPSF image with a randomly-scaled IRF prior to feature extraction.Given the significant influence of the IRF on TD data (Hansen 2010, Diop and St Lawrence 2013), we focused on extracting image features from which the effects of the IRF could be easily removed: statistical moments and temporal Fourier components (see section 2.2.1).Notably, since the proposed feature extraction is based on an analytical approach, the computational time associated with our analysis is negligible even when performed with a consumer-grade laptop.Motivated by previous work-which showed notable changes in the spatial distribution of image information in response to changes in RA disease activity (Hielscher et al 2011, Lighter et al 2019)-we also performed spatial Fourier decomposition on extracted TD-DOI features.Despite the addition of noise and the effects of the IRF, the spatial Fourier components of many extracted features remained highly correlated with DI (figure 3).Further, the 0.06 mm −1 phase of lower temporal frequency Fourier components (⩽0.4 GHz) was most sensitive to small, likely subclinical, physiological changes (figure 4).
Using CW-DOI, Lighter et al (2019) found that changes in 0.05 mm −1 and 0.1 mm −1 spatial frequencies were the most accurate at detecting early inflammation, and that the length scales corresponding to these spatial frequencies approximately matched the size of dips observed in the intensity profiles of inflamed PIP joints.These findings are consistent with the results of our simulations: both the 0.06 mm −1 and 0.13 mm −1 spatial frequencies of many features were most sensitive to small changes in simulated RA disease activity (figure 4).The spatial frequency corresponding to the width of the joint area in the current work was 0.09 mm −1 which is between 0.06 mm −1 and 0.13 mm −1 .As such, our results suggest that the sensitivity of a given spatial frequency component to RA disease activity is related to how closely it matches the width of the joint cavity.Further, it appears that we can expect changes in RA disease activity to produce changes across a band of spatial frequencies.
Overall, the TD-DOI features revealed that lower frequency temporal Fourier components had substantially higher sensitivity than any of the extracted statistical moments (figures 3 and 4).One potential explanation for this observation is the well-known, considerable sensitivity of statistical moments to noise (Liebert et al 2003b); this is especially true when calculating higher-order statistical moments, and is consistent with the poorer sensitivity of V to RA disease activity compared to m 1 .Indeed, without the addition of Poisson noise, the 0.06 mm −1 and 0.13 mm −1 spatial frequencies of V were found to have a high correlation with DI and, in the case of m 1 , these frequencies also became more strongly correlated than originally observed.As for the F t [TPSF] features, Poisson noise is known to have a greater effect on higher frequency temporal Fourier components (Gu et al 2007); this effect can clearly be seen in the tendency of R and η 2 G values to decrease with increased F t [TPSF] frequency (figures 3 and 4).It is important to note that temporal Fourier decomposition provides insight into the relative sensitivity that is expected when using other DOI modes of acquisition.Specifically, the DC temporal component of F t [TPSF] (i.e. the total number of photons) is equivalent to a typical CW-DOI measurement, while the remaining F t [TPSF] measurements are equivalent to FD-DOI measurements taken at various modulation frequencies.From this perspective, our results show that the sensitivity of both CW and FD to RA disease activity is enhanced when looking at changes in spatial frequency components aside from the DC spatial component.Nevertheless, the DC spatial components of some F t [TPSF] frequencies still demonstrated relatively high correlations with DI (figure 3) and sensitivity to small DI changes (figure 4).FD instruments for diffuse optical applications typically operate in the 0.1-1 GHz range (Durduran et al 2010); for small tissue volumes (e.g.human finger), numerical simulations have identified 0.4 GHz as an optimal modulation frequency (Gu et al 2007).Hielscher et al (2011) acquired FD-DOT images of PIP joints at DC, 0.3 GHz, and 0.6 GHz, and found that DC images had considerably lower Youden indices (a statistical measure of diagnostic performance) than their counterparts.A similar trend can be noted in the phase of various spatial frequencies in figure 4: the largest η 2 G values generally occurred at F t [TPSF] frequencies slightly greater than the DC frequency, though the exact frequency varies depending on the spatial frequency component under consideration.
Overall, the results of this study confirm the previously demonstrated benefits of utilizing spatial frequency analysis to enhance sensitivity to changes in RA disease activity (Montejo et al 2013b, Lighter et al 2019).Further, we have observed that utilizing this analysis on low temporal frequency F t [TPSF] features extracted from TD-DOI images appears to yield the greatest sensitivity, and that this sensitivity is primarily related to increases in joint volume and changes in the optical properties of the synovial fluid.In particular, our results suggest that FD-DOI and TD-DOI may enable us to track the latter, likely subclinical, changes with greater sensitivity than previously tested CW-DOI (Lighter et al 2019).Additionally, since TD-DOI is able to simultaneously assess changes across a band of F t [TPSF] frequencies, it may offer a more robust confirmation of potential changes in disease activity compared to a single-frequency FD system.Of course, we note that this same information could be obtained by using an FD-DOI system to acquire multiple images at several modulation frequencies.
While encouraging, we note that these results were obtained in silico using a finger model which may not fully capture the true tissue geometry and physiological changes present in an in vivo human finger.For simplicity, this model could not account for localized changes in physiology and relied on assumptions, sourced from the literature, to define average changes in tissue geometry and physiological properties as a function of simulated RA disease activity.Notably, the use of linear transformations to model different tissue geometries is clearly a simplification of the heterogeneity in synovitis and synovial effusion observed in individual patients.Further, we limited our work to studying the individual effects of 5 RA-associated parameters whose effects were limited to an area encompassing the PIP joint.Therefore, we anticipate that the high sensitivity of the TD-DOI features presented herein to RA disease activity may be more pronounced than they would be in an in vivo setting, and the findings of this work inform about the average sensitivity of TD-DOI to various physiological changes, as opposed to TD-DOI's ability to assess the precise level of physiological changes within an individual patient.Future work will aim to experimentally validate these findings in solid tissue-mimicking phantoms using an in-house TD-DOI system.

Conclusion
This work presents a systematic in silico study that demonstrates the potential of TD-DOI to track changes in RA disease activity.Specifically, we investigated TD-DOI's sensitivity to 5 RA-associated physiological changes under 4 baseline disease activities.In total, TD-DOI images of 80 disease activity states-each modelled using a virtual finger phantom based on an adult PIP joint-were compared both with and without the addition of Poisson noise and the effects of an IRF.A DI based on the variance of the former images revealed that TD-DOI is primarily sensitive to changes in joint volume, but also has sensitivity to changes in the optical properties of the synovial fluid.In contrast, low to negligible sensitivity was observed to changes in oxygen saturation, blood volume fraction, and water concentration within the joint.Furthermore, many spatial Fourier components of features extracted from TD-DOI images which incorporated the influence of Poisson noise and the IRF remained highly correlated with DI.In particular, the phase of the 0.06 mm −1 spatial frequency of TD-DOI images' low frequency temporal Fourier components (⩽0.4 GHz) was both strongly correlated with DI and highly sensitive to small physiological changes.

Figure 1 .
Figure 1.Schematic of the creation of a normalized 2-dimensional TPSF image (right) from raw 3-dimensional TD-DOI data (left).

Figure 2 .
Figure 2. Heatmap of mean intra-set change in z-score normalized DI values for various physiological parameters (see table 1) at different states of baseline disease activity.Sets in which DI values did not maintain a consistent direction of change were excluded from the analysis.

Figure 3 .
Figure3.Heatmap of correlation (R) between the spatial Fourier components of features extracted from normalized DTOF images and the reference PCA-1 (i.e.DI) value of their respective TPSF images.As described in section 2.2.1, three types of features are considered: the first normalized moment (m1), the second centralized moment (V), and the temporal Fourier coefficients (GHz).Note that the magnitude and phase of the components were correlated separately, and components from images of Ft[TPSF] > 1.6 GHz were excluded.

Figure 4 .
Figure 4. Heatmap of generalized eta squared (η 2 G ) for spatial Fourier components of features extracted from normalized DTOF images of four virtual phantoms with different µSF values and Mild baseline disease activity.Only components with η 2 G ⩾ 0.26 are shown.

Table 1 .
Changes in tissue-specific physiological parameters associated with worsening RA disease activity.

Table 2 .
Tissue parameter values (%) for None and Severe disease activities.Detailed references for parameter values are provided in our previous work (Ioussoufovitch and Diop 2021).S: oxygen saturation; W: water volume fraction; Bl: blood volume fraction; F: fat volume fraction; cF: fat concentration; M: melanin volume fraction; Bn: bone volume fraction; D: disease activity fraction.Reproduced with permission from Ioussoufovitch and Diop (2021).© (2021) Society of Photo-Optical Instrumentation Engineers (SPIE).
a Skin and bone parameters were set to the same values for all models.b Disease activity fraction represents the weighting given to optical properties obtained from arthritic synovial fluid (Milanic et al 2014) compared to the optical properties of water in the calculation of a particular model's synovial fluid optical properties.