Systemic Sclerosis (SSc, aka Scleroderma) is a complex, rare, connective tissue disease (CTD) involving the collagen, major organs, the blood vessels and the immune system (van den Hoogen et al., 2013; Di Battista et al., 2021), in which extensive fibrosis and vascular alterations take place. It has significant morbidity and mortality (van den Hoogen et al., 2013), and in the UK an estimated prevalence of 307 per million (95% CI: 290–323), with the highest occurrence in the 70–84 years age group (Royle et al., 2018). SSc (the two most common variants are limited cutaneous variant (lcSSc), and diffuse cutaneous (dcSSc)) is often associated with Raynaud’s phenomenon (RP), a condition in which recurrent, reversible vasospasm of the digital small arteries, arterioles, pre-capillary and post-capillary venules occurs on exposure to cold or emotional stress (Hughes and Herrick, 2016; Silva et al., 2016). RP is common and in the UK is reported to affect up to 10 million people (SRUK, 2023). About 1 in 16 women and 1 in 50 men with Raynaud’s develop SSc, usually between the ages of 25 and 55 (NHSinform, 2023; Haque, 2020; Belch, 2017). It is usually sub-categorised into: a) Primary RP (PRP) when no underlying cause condition is known (idiopathic); b) Secondary RP when RP is linked to an underlying disease such as SSc or dermatomyositis or to the intake of certain drugs. Secondary RP is typically seen in approximately 90%–96% of patents with SSc and often precedes the development of SSc by an average time of 10.4 years (Spencer-Green, 1998; Pauling et al., 2019). Clinical specialists differentiate secondary RP from PRP by checking for symptoms associated with secondary RP such as the age at onset (secondary RP is usually after 30 years of age), detecting abnormal immunology e.g., certain autoantibodies, observing nailfold capillaries, ulceration of digits, checking for fibrosis in the lungs or other organs, and skin thickening which is the hallmark of SSc. However, diagnosing SSc is not always easy as its symptoms resemble other conditions such as PRP and early symptoms of diseases such as systemic lupus erythematous (where 10% and 45% of patients show Raynaud’s phenomenon). Early detection and management of the disease is a must, to improve the morbidity and mortality in patients (Walker et al., 2007). This in turn requires a multi-disciplinary and collaborative effort involving clinical specialists and testing. It can take more than one visit to an expert Rheumatology specialist to diagnose the disease, especially in the early stages. Identification of internal organ involvement and its severity is also important.

These involve extensive and costly testing for autoantibodies and markers of organ involvement. Sometimes, it is difficult to distinguish between SSc and non-SSc cases as patients have overlap conditions. Nailfold capillaroscopy (NFC) is another key technique used to help diagnose SSc. NFC is a non-invasive, optical imaging technique (Allen and Howell, 2014; Eriksson et al., 2014) that is used by an expert operator to visually inspect the microcirculation in the nailfold capillaries of the distal papillae and hence assess pathological/morphological changes associated with SSc such as capillary “dilatation”, distribution and density (for “drop-out”), bushing (for “angiogenesis”) and microhaemorrhage (extravasation). Tests such as NFC, however, are usually currently performed in specialist hospitals and are not at all readily available to all patients.

There is huge scope to look for alternative, low-cost technologies to assess SSc. Photoplethysmography (PPG) is one such technique, which is non-invasive and optically assesses the circulation. The working principle of a PPG-based system uses a suitable light source such as infrared or near-infrared light to study the heart-synchronous changes in blood volume in the microvascular bed of tissues such as skin (Allen, 2007; Elgendi, 2012; Kyriacou and Allen, 2021). Additional key advantages of PPG are its portability and its versatility to be used in a range of settings such as measurement labs as well as ambulatory assessments (wearable sensors, Charlton et al., 2023). PPG is currently widely used in different clinical applications (Johnson et al., 2020) including for pulse oxygen saturation measurement (SpO2) (Ma et al., 2018), cardiovascular health (heart rate, blood pressure, blood vessel and arterial stiffness) monitoring (Castaneda et al., 2018), and for studying hypertension (Liang et al., 2018).

The potential of PPG for detecting patients with SSc has already been explored using conventional optical pulse wave analysis techniques but the literature here appears to be limited to date. The largest SSc PPG study reported to date is by Rosato et al. (2010) with 105 SSc (compared to 96 PRP and 85 healthy controls) using a Termoflow type PPG instrument. The authors found that the mean amplitude of the PPG sphygmic wave was significantly lower in the PRP group than in the SSc group (mean ± standard deviation = 11 ± 10, given in arbitrary units (a.u.) vs. 24 ± 24 a.u. for SSc). The mean amplitude was also significantly lower in SSc than in HC (56 ± 19 a.u. for HC). A further study by Rosato et al. (2011) using bilateral PPG measurements found a homogeneous pattern (meaning uniformity of morphology and amplitude of sphygmic PPG wave across all 10 fingers) for 95% of the HCs and 93% of the PRPs but was only 28% for the SSc group. McKay et al. (2014) investigated multi-site PPG in 19 SSc, 19 PRP and 23 HC by studying measures of arterial, endothelial and peripheral autonomic dysfunction under a dynamic 3 phase testing protocol. The authors found that measures attributed to endothelial function were significantly impaired in SSc (p < 0.02), but with no difference between the HCs and PRPs. The authors reported that the Receiver Operating Characteristic (ROC) based classification accuracy was 81% (sensitivity 90%, specificity 74%) for separating SSc from HCs, and 82% (sensitivity 84%, specificity 79%) for separating SSc from PRPs. Mamontov et al. (2020) employed imaging PPG (iPPG) and studied 19 SSc and 21 HC participants (with age- and sex-matched balanced classes) and from the pulse arrival time (PAT) showed a significant increase in its variability in SSc patients as compared to HCs (52 ± 47 ms vs. 24 ± 13 ms, p = 0.01). These earlier works show the potential of PPG as a tool to help investigate SSc but such approaches from the literature have relied on extracting a range of pulse features based on domain knowledge, and then performing specific feature selection to try and improve classification performance. To the authors’ best knowledge, there are no papers on contact type PPG measurements using analysis based on deep learning (DL) for the assessment of SSc patients. The aim of this pilot study was to utilise DL in a novel way by applying it to automatically identify the presence (or absence) of SSc from multi-site PPG measurements collected under a dynamic 3-phase test protocol.

Consecutive patients were approached by specific autoimmune connective tissue disease rheumatologists or the connective tissue disease nurse specialist from the Rheumatology outpatient population attending Freeman Hospital, Newcastle upon Tyne. The SSc participants were each diagnosed by an expert consultant Rheumatologist at Freeman Hospital, using the 1980 American College of Rheumatology (ACR) preliminary criteria for the classification of systemic sclerosis (Subcommittee for Scleroderma Criteria, 1980). PRP patients were diagnosed by the same clinical team, as having vasospastic symptoms for >2 years, with routine practice assessments and no other underlying medical or mechanical cause. SSc and PRP participants were recruited from the Rheumatology outpatient clinic at Freeman Hospital, Newcastle upon Tyne. Healthy Control (HC) participants were recruited from volunteers - largely from the University of the Third Age (U3A, Wearside Branch), and from staff and students of Newcastle Hospitals and Newcastle University. The HC participants had no known underlying health condition (e.g., diabetes mellitus and hypertension). For each patient a range of clinical and demographic data were collected for the study, including participants’ pertinent medical history, and the tests done in the hospital to assess and diagnose SSc patients. All participants were permitted to continue their regular medication, including vasodilator treatment, at the time of their recruitment and informed consent was taken. Ethics approval for the original study data collection was granted by the National Research Ethics Service (NRES) Committee Northeast (County Durham and Tees Valley 1 REC, 07/H0905/72 2008). Ethics permission for re-analysis of the anonymised PPG waveform data for Sadaf Iqbal’s PhD studentship project was obtained from Newcastle University (Reference 7273/2018, with an extension to the study 17138/2021).

The PPG data set measurements were originally collected by Dr Neil McKay, Rheumatologist, between 2009 and 2011 at Freeman Hospital’s microvascular research facility, with 3-phase measurement protocol developed and physiological measurement training supported by expert PPG operator JA. Participants were firstly asked to lay comfortably in a symmetrical supine position for a period of at least 15 min whilst acclimatising in a warm normothermic temperature-controlled (23°C ± 1°C) clinical measurement room. Multi-site PPG waveforms were then collected simultaneously for 20 min from 6 symmetrical, body sites namely: right and left earlobes, index finger pads and great toe pads respectively, using optically and electronically matched amplifiers (bandwidth 0.5–20 Hz) and captured to computer at a sampling frequency (Fs) of 2000 Hz. The 3 phases of measurement were: subject resting supine (10 min, Baseline phase); an arm pressure cuff inflated at 300 mmHg to stop the arterial blood flow into the left arm (for 5 min, Occlusion phase); then at 15 min the cuff pressure was quickly but carefully released, and the degree of reactive hyperaemia monitored for a further 5 min (Flush phase). Figure 1 shows examples of 3-phase PPG beat-to-beat amplitude data for the left finger measurement site of a Control participant and a SSc participant.

For this study involving advanced analysis of the data, a visual analysis of the PPG data was carried out by operator SI to check for the presence of unexpected artifacts (such as that caused by unreliable PPG probe skin contact or when a study did not follow protocol) or a distorted flush response was evident. Some of the 92 participants originally entered into the study had to be excluded: original PPG data collected using too high a manual gain setting causing a high flush which saturated i.e., electronically clipping and making unusable a PPG trace for the left study arm (N = 3); participants not following the 20-min protocol (N = 7). In total, PPG and ECG data from a total of 20 SSc, 22 PRP and 29 HC participants were included and analysed using the techniques described in the next section. Figure 2 shows a participant flowchart summary for the included and excluded subjects in our pilot study.

Healthy controls and patient controls (i.e., PRP) participants were combined in a single control group as they both belong to clinical class of non-life-threatening conditions as opposed to SSc group which needs clinical attention and early identification and management. Also, noting here that although pulse amplitudes can be different between the groups and from the same subject if tested on different days, the literature suggests that PPG pulse measures are similar across the two groups. The study by Rosato et al. (2011) had shown that PPG shows a homogeneous pattern, i.e., a uniformity in morphology as well as amplitude of PPG sphygmic waves across the fingers, in 95% of the healthy control subjects and 93% of the PRP patients as opposed to this being present in only 28% of the SSc patients. McKay et al. (2014) also had found no differences between healthy controls and PRP for the case of dynamic physiological testing i.e., the derived PPG median flush response slope and flush response value.

In this study a PPG signal was pre-processed with only a normalisation stage, with each of the 6 PPG signals divided by their respective PPG amplifier channel gain setting. There was no additional signal filtering performed by the computer.

This data analysis was carried out on a 64-bit Windows 10 and 14 Core PC fitted with a single NVIDIA GeForce RTX 2080 Ti GPU card which used MATLAB software. Here, the PPG time series was converted into scalogram images for the study subjects using continuous wavelet transform (CWT). Scalograms give a time-frequency (T-F) image representing the percentage of signal energy contained in different frequency bands (Gandhi et al., 2014), with the x-axis representing time and y-axis representing frequency and a signal varying in colour intensity over the T-F plane.

Each PPG channel’s gain normalised time series was divided into consecutive 30 s non-overlapping windows and for each 1D time window, a 2D scalogram image along with its label was then generated using continuous wavelet transform (CWT) methods. The Morse mother wavelet was used in CWT as this type of analytic wavelet is very useful to analyse signals with time-varying amplitude and frequency (PPG here) (Wachowiak et al., 2018) and Voices per Octave was selected as 12 to keep the computational complexity low (Wachowiak et al., 2018). The 30 s period for plotting scalograms was selected based on initial exploratory scalogram plotting to select meaningful T-F resolution and hence extract meaningful features from the scalogram images. Noting that with the uncertainty principle, the greater length of time window means that frequency resolution is higher but time resolution is lower, and vice versa. Hence, for the data under consideration a 30 s epoch provided a reasonable trade-off between the time and frequency resolutions. Figure 3 shows a 30 s sample scalogram with frequencies ranging from 0–20 Hz (noting each PPG was already bandpass filtered in this range using analogue electronics during physiological data acquisition). The scalogram was calculated using Eqs 1, 2 below: S = c o e f . * c o e f w h e r e   c o e f   a r e   t h e   C W T   c o e f f i c i e n t s   o v e r   a   t i m e   w i n d o w (1) S c a l o g r a m = 100 × S / ∑ S (2)

In the first stage, the results were calculated using the number of images that were classified into SSc versus Control (‘image-based performance’). Two pretrained convolutional neural networks (CNN) from MATLAB were used (GoogLeNet (Szegedy et al., 2015) and EfficientNetB0 (Tan et al., 2019)) to learn the PPG T-F features and then perform the classification. These networks have already been previously trained on millions of high-resolution images from the ImageNet database (Krizhevsky et al., 2017). GoogLeNet (structure comprising 144 layers) was the winner of 2014 ImageNet competition and had least number of parameters (∼6.9 million) as compared to other pretrained models available in MATLAB at the time of study. EfficientNetB0 (290 layers) which was introduced in 2019 (parameters ∼5.3 million) represents a newer generation of CNNs, based on ResNet design, and had been designed to work better and faster than the hitherto available CNNs (Tan and Le, 2019), giving 2 key types of CNN to implement and explore respective performances. Figure 4 shows the basic building blocks for these two types of CNN network. Figures 5, 6 respectively show the network architectures of GoogLeNet (Szegedy et al., 2015) and EfficientNetB0 (Tan et al., 2019; Tan and Le, 2019).

For the PPG dataset, 16,188 labelled scalogram images were generated for the 71 participants covering all 6 PPG sites, i.e., channels of the multi-site PPG system. The scalogram images were rescaled to 224 × 224 × 3 size to match input dimensions of GoogLeNet and EfficientNetB0 using a standard MATLAB function “imresize” which applies a simple scale transformation to the original image using bicubic interpolation which is a standard algorithm used in image resizing (Hashemi et al., 2015). Data were divided into 10 mini batches (batch size 1620 images) to reduce computational time. The learning rate was 0.005. Optimisation utilised the stochastic gradient descent method (Tian et al., 2023) as this is computationally faster and can converge quicker than other optimisation algorithms. The loss function used was cross-entropy loss. Ten-fold stratified cross validation was carried out on 71 participants wherein 9/10 of the participants’ images were used in training and 1 mutually exclusive (1/10 of the participants’ images) were used for its testing, this process was repeated 10 times. The image wise combined confusion matrix obtained after 10-fold CV was used to calculate measures of classification performance. The two networks train on the input training data and adapt their weights, to learn the features of the data. The approximate training time in any fold of training was about 240–300 min. The approximate testing time for any test fold was between 20–30 min.

The SSc versus Control classification was performed based on number of images and this is called as image-based classification throughout. A post-processing step was also applied which calculated the number of images classified in each category per each individual test participant, to provide SSc versus Control classification based on number of participants. In this case, the majority class of images in each test participant was considered as the class of the output label. This was compared against the ground truth label (i.e., SSc diagnosis) which was clinically determined beforehand. The steps and methodology of the DL analysis are summarised in the block diagram in Figure 7.

This is a comparator ML method using the discrete wavelet transform (DWT) to compare with the performance of the T-F based DL classification. DWT method has advantages over other T-F methods such as CWT such in being computationally faster and the ability to analyse the input signal into desired constituent frequency bands. The same participant-wise partitions of the training/test data sets were used as for the DL classification work to allow a fair inter-comparison between methods.

Ten-level DWT decomposition (Mohamed and Deriche, 2014) was carried out per channel PPG using Daubechies 4 (db4) mother wavelet thereby producing 10 levels of detailed coefficients (d1, d2, ……. , d10) and approximate coefficients (shortened as ap in this work) for each channel. The mother wavelet db4 was chosen as it matches the shape of the PPG pulse more than any other wavelet, and hence is the best choice to calculate the wavelet transform of the signal. Since the PPGs had already been bandpass filtered (i.e., 0.5–20 Hz) at data acquisition then only the detailed coefficients containing these relevant frequencies namely, d6 (frequency range 15.63–31.25 Hz), d7 (7.81–15.63 Hz), d8 (3.91–7.81 Hz), d9 (1.95–3.91 Hz), d10 (0.98–1.95 Hz) and ap (0.49–0.98 Hz) were selected for further analysis. To make the computational complexity less and thus the algorithm faster, 4 features were extracted from each of the 6 frequency bands for each channel, thereby giving 6*6*4 = 144 features per participant. PPG features extracted were Energy, Entropy, Mean absolute value and Skewness, as defined in Eqs 3–5.

Energy: Energy of the PPG DWT coefficients is calculated as the sum of the squares of all the sampled amplitudes for a single channel and is defined by Eq. 3: E n e r g y = ∑ i = 0 N x i 2 (3)

Where x_i, the i^th data instance and N are the total number of sample values.

Entropy: Entropy describes the irregularity, complexity, or unpredictability characteristics of a signal. In this work entropy was used as a feature to quantify the irregularity of the PPG time series and hence Shannon entropy was calculated as given by Eq. 4: E n t r o p y = ∑ i = 1 N x i   ∗   log 2 1 / x i (4)

Mean absolute value: The absolute value of the average of each channel within the data matrix.

Skewness: Skewness can be used as a measure to describe the asymmetry around the mean of the data sample. If the value of skewness is less than zero, the data has more spread around the left-hand side of the mean, if greater than zero, it is more towards the right-hand side of the mean and if equals to zero, the data can be considered as symmetrically distributed. For a dataset, skewness can be described by Eq. 5: S k e w n e s s = E x i − μ 3 σ 3 (5) Where µ is the mean of x, σ is the standard deviation of x, and E(t) represents the expected value of the quantity t.

The extracted features were then fed into the linear discriminant analysis (LDA) and into the K-Nearest Neighbour (KNN, set at K = 9) analysis to classify SSc versus Control.

Demographic data were expressed using mean (±standard deviation, SD) values. Since the subjects are grouped into 2 distinct SSc and Control classes with each group having independent and different participants, hence unpaired Student’s t-test was used to study the mean values of SSc and Control groups. A p-value <0.05 was the level of statistical significance. A crosshair plot illustrating a meta-analysis overview of classifier performance (test sensitivity versus False Positive Rate i.e., 1-specificity) was produced using the mada command from R (RStudio version 1.4.1106).