The Hybrid Excess and Decay (HED) model: an automated approach to characterising changes in the photoplethysmography pulse waveform [version 1; peer review: awaiting peer review]

Photoplethysmography offers a widely used, convenient and non-invasive approach to monitoring basic indices of cardiovascular function, such as heart rate and blood oxygenation. Systematic analysis of the shape of the waveform generated by photoplethysmography might be useful to extract estimates of several physiological and psychological factors influencing the waveform. Here, we developed a robust and automated method for such a systematic analysis across individuals and across different physiological and psychological contexts. We describe a psychophysiologically-relevant model, the Hybrid Excess and Decay (HED) model, which characterises pulse wave morphology in terms of three underlying pressure waves and a decay function. We present the theoretical and practical basis for the model and demonstrate its performance when applied to a pharmacological dataset of 105 participants receiving intravenous administrations of the sympathomimetic drug isoproterenol (isoprenaline). We show that these parameters capture photoplethysmography data conclude by discussing the possible value in using the HED model as a complement to standard measures of photoplethysmography signals.


Introduction
Photoplethysmography is an inexpensive and non-invasive means of measuring heart rate (HR) and peripheral oxygen saturation. It is ubiquitous in healthcare settings, and increasingly present in consumer-level wearable devices (Zhang et al., 2020). Photoplethysmography is also of growing interest in cognitive neuroscience research, where heart rate variability (HRV) can be derived as an indirect measure of autonomic arousal. Moreover, it is frequently employed to provide objective data on HR, against which individuals' subjective awareness of their own heartbeats can be compared i.e. as a measure of interoception (Khalsa et al., 2008;Legrand et al., 2022;Smith et al., 2021).
Here, we suggest that the fundamental waveform of the signal -the pulse waveform -is a rich and often overlooked source of information capable of complementing the more frequently-derived measures of HR and HRV. The haemodynamic events underlying the shape of the waveform are complex and, ultimately, offer the potential for a granular assessment of cardiovascular features which, in turn, are of relevance to physiological and psychological states. Thus, focusing on morphological features of the photoplethysmography signal could lead to new characterizations of psychophysiological states relevant to mental health and cognition.
The photoplethysmography signal is acquired using a paired light emitter and sensor (Jennings et al., 1980). The sensor measures the quantity of emitted light that returns after transmission or reflection through the tissue. Changes in the composition of the tissue -such as fluctuations in local peripheral blood volume -are therefore reflected in the signal. Typically, the signal is described in terms of its alternating current (AC) and direct current (DC) components (Akl et al., 2014;Allen, 2007). The AC component comprises the visible fluctuation in signal generated by arterial pulses, with each pulse giving rise to a distinct waveform spanning the length of the cardiac cycle. The DC component represents the light absorbed by venous blood and the remaining tissue. Fluctuations in the DC component are reflected in the signal baseline and are influenced by other physiological phenomena such as respiration, sympathetic nervous system activity, and thermoregulation (Allen, 2007).
The resulting waveform measured using photoplethysmography has a distinctive morphology (see Figure 1) related to that of the central aortic pressure wave (Imholz et al., 1998). From the field of arterial haemodynamics, there are two distinct but complementary models for understanding how the waveform is shaped: (i) The Windkessel model: Cardiac systole generates a pressure wave that is sustained by the elastic compliance of the arterial vessels, resulting in an initial surge followed by a slowly decaying pressure wave during diastole (Westerhof et al., 2009).
(ii) The Wave model: The initial systolic ejection is followed by one or more reflectance waves arising from distal sites in the arterial tree where impedance to flow is greater (van de Vosse & Stergiopulos, 2011). The precise nature and origin of these reflectance waves is not fully understood (Westerhof et al., 2008) but they have an appreciable effect on the overall shape of the wave during diastole.
Both models suggest that the waveform's shape is highly sensitive to cardiovascular state and therefore rich in information regarding it. It may also vary within and between individuals in ways that are highly informative. Indeed, existing work has demonstrated the value of the pulse waveform as a marker for physiologically meaningful information (Allen, 2007). Contour analysis, the identification of fiducial points on the waveform from higher-order derivatives (Elgendi, 2012;Millasseau et al., 2006), is a purely descriptive means of quantifying waveform morphology. It has identified a number of associations between the waveform and cardiovascular factors, including central arterial stiffness (Millasseau et al., 2002), systemic vascular resistance (Lee et al., 2011), blood lipid levels (Sattar et al., 2014), and response to vasoactive drugs (Takazawa et al., 1998). Importantly for our purposes, it has also demonstrated a sensitivity of the waveform to changes in psychological state due to associated autonomic and cardiovascular effects (Abbod et al., 2011;Charlton et al., 2018;Rinkevičius et al., 2019) (see Figure 1). Other general approaches to plethysmography waveform analysis include frequency domain (Nitzan et al., 1994), non-linear (Sviridova & Sakai, 2015, and machine learning methods (Zhang et al., 2021). Whilst not considered here, they hint further at the potential richness of the signal.
There are thus grounds for believing that the waveform offers a window onto processes highly relevant to brain-body interaction (i.e., interoception (Khalsa et al., 2018)). To speculate on the nature of this, one can imagine that, in addition to top-down drivers of autonomic activity, the waveform may also reflect bottom-up signals sent to the brain via vascular afferents, and this wave could thus be used by the brain to gain information about the peripheral body state. These could subsequently be useful in computing actions for goal-directed behaviour, or to augment exteroceptive input with interoceptive afferent information to inform adaptive responses. In short, the shape of the waveform could be a part of the array of interoceptive signals modulating brain function.
Given the existing interest in psychological and physiological factors influencing, and potentially influenced by cardiovascular state (Mulcahy et al., 2019;Penninx, 2017), a comprehensive characterisation of the features of this waveform could be a powerful complement to standard measures such as HR and HRV in characterizing and monitoring psychophysiological states. An initial goal, therefore, is to develop a nuanced means of characterising the waveform's shape, allowing quantitative assessments in its variability that can be related to physiological and psychological states. Whilst contour analysis has proven useful, the haemodynamic events underlying the waveform's shape (alluded to by the Wave and Windkessel models) are complex and may warrant a more data-driven analysis of the waveform (rather than one that is focused on pre-selected points of interest).
Pulse decomposition analysis (PDA) is a pulse wave analysis method which assumes that each waveform is composed of a set of underlying component waves. The method aims to infer these waves from the overall waveform, thereby moving beyond the purely descriptive approach of contour analysis. Once identified, the component waves help to parametrise changes in the overall waveform shape, and offer potential insights into the haemodynamic events underlying a change. The PDA approach has been applied previously to both simulated and experimentally-acquired photoplethysmography data. For camera-acquired data, it has been shown that PDA-derived measures correlate closely with contour-based measures whilst demonstrating greater robustness to noise (Fleischhauer et al., 2020). It remains to be determined, however, what the optimal number of composite waves is, and which kernels are used to represent them.
Here, we sought to develop and extend the PDA approach and to examine its performance in comparison to more standard, feature-based analyses of photoplethysmography data (HR, HRV and contour analysis). We present a prototype for characterising the waveform: the Hybrid Excess and Decay (HED) model. It is comparable to PDA insofar as it expresses the waveform and its variability in terms of component waves, but differs in that it involves the addition of three parameters to characterise the quasi-exponential decay observable in the pulse waveform during diastole. We describe and present the model along with an automated means of applying it to photoplethysmography data, now available as a free-to-use R package. We go on to demonstrate its utility by applying it to data from a study in which the bolus effects of intravenous isoproterenol, a potent peripheral beta-adrenergic agonist, were compared to those of saline in a sample of 118 participants (from which 105 datasets were suitable for study). We assessed the model's performance, in terms of goodness-of-fit, on this dataset as well as its ability to recapitulate the key fiducial points (O, S, N, and D, henceforth OSND) seen in Figure 1. Finally, we assessed its ability to classify between waveforms recorded during isoproterenol infusion versus waveforms during saline infusion, as compared to contour analysis, HR, and HRV. We conclude by considering the model's concurrence with existing concepts in haemodynamics, and its potential for optimising the usefulness of the signal in characterizing psychological and physiological states.

Methods
Our approach comprised two main parts: Model development: This entailed (i) a theoretical analysis of the most likely constituents of the photoplethysmography waveform and (ii) refinement of the model and its assumptions with the use of a small pilot dataset.
Model validation: An analysis pipeline was developed to test the model on a larger dataset acquired on different hardware from that used in the pilot dataset. This dataset came from a study in which cardiovascular parameters were systematically manipulated using isoproterenol infusion. We aimed to examine (i) how well the model performed across both the drug and placebo infusions and (ii) its sensitivity to the pharmacological challenge (isoproterenol versus saline) in comparison to other standard measures, including HR, derived from photoplethysmography.

Model development Theoretical foundation and model components.
A single prevailing view on the nature of arterial waves and their reflections has not yet been reached, largely due to the inherent complexity of the arterial tree (Segers et al., 2017). Nonetheless it is generally accepted that backward travelling reflectance waves exist and are likely to be measurable as composites from multiple reflection sites, irrespective of origin. Hitherto, three component waves of the pulse waveform have been described: -A systolic wave caused by the increase in pressure, and therefore volume, arising from left ventricular contraction and ejection of oxygenated blood into the arterial system.
PDA models have taken a largely data-driven approach, using combinations of component waves of varying number and shape. Comparisons of these suggest the optimal number of component waves for capturing the waveform's morphology is three (Tigges et al., 2017). From this starting point, we modelled the waveform as a composite of three component waves (systolic, first reflectance wave (R1), and second reflectance wave (R2) as above). Nine parameters characterised the timing, width and amplitude of each component wave.
The HED model, however, differs from existing PDA models in its approach to the diastolic decay. We noted a tendency for waves of prolonged duration to decay below baseline, an established phenomenon in aortic pressure waveform studies (Parker, 2017). To our knowledge, this decay is not explicitly parameterised by current PDA models. To account for this, and to better elucidate the underlying factors driving waveform morphology, we incorporated three additional parameters: -Decay rate: The rate at which the signal decays exponentially to baseline in the absence of component wave influence (Wang et al., 2003).
-Baseline 1: The baseline towards which the signal decays during the systolic portion of the waveform.
-Baseline 2: The baseline towards which the signal decays during the diastolic portion of the waveform.
In summary, we modelled the waveform as a composite of signal due to the initial systolic pressure wave, an exponential decay, and two reflectance waves. For simplicity, we refer to the systolic and reflectance waves as the "excess" elements and the exponential decay as the decay element (see Figure 2), referring to the overall model as the HED model. The output of the model is a 12-parameter vector, which can be used to construct a fitted wave to the original photoplethysmography data (see Figure 2).

Refinement of the model and assumptions.
Once established conceptually, the model was refined using a set of test data samples acquired from members of the experimental team at rest using a Bioradio device (Great Lakes Neurotechnologies, Ohio, USA), a multichannel portable device designed for simultaneous acquisition of multiple physiological signals (glneurotech, 2020). The accompanying sensors are transmission-based finger clips made of either hard or soft plastic. Both types were used at the index finger, with a sampling rate of 75 Hertz (Hz). This small test dataset served as the basis for assessment (see Figure 3) and correction of aberrant model behaviour, as described in the following paragraphs. A number of assumptions were made in developing the model, as is common in PDA (Fleischhauer et al., 2020). For the excess components, sine waves were used as they exhibit a number of suitable properties: first, a sine wave smoothly interpolates to zero which fits with the appearance of photoplethysmography waves, insofar as the wave onset is not sudden; second, we found it useful to look at the first derivative of the waveform in order to identify beats and this derivative is well-fitted by a sine wave. For the decay element, a second baseline spanning the diastolic portion was included to allow for more flexible modelling of varying dicrotic notch (see Figure 1) height. Additionally, we assumed that (i) component waves cannot have the same timing and must occur in order, (ii) R1 plays a relatively minor role in shaping the waveform and must therefore have the smallest amplitude and width, and (iii) the number of component waves need not be three if a more parsimonious fit can be achieved with fewer. Constraints of the model reflect these assumptions and are detailed in full in the PulseWaveform GitHub project.
Further, we noted that six of the 12 model parameters tended to converge on constant values across multiple resting heartbeats. These were (i) width of systolic wave (ii) width of R1 (ii) width of R2 (iv) timing of R1 (v) timing of R2 and (vi) decay rate. The option to fix these parameters across multiple beats during parameter optimisation was added. Doing so greatly reduced the number of free parameters across time series (in brief, the number of free parameters was reduced by (n/2 +6) where n is the number of free parameters (12) multiplied by the number of beats in the dataset) as well as allowing greater consistency in how adjacent beats were fitted. There was also a corresponding reduction in computational time required. We therefore maintain this function for the analyses herein but stress the optionality of this feature in our shared code such that waves can still be modelled individually if desired.

Model validation
Formal testing of the model took place on a dataset from an ongoing functional magnetic resonance imaging (fMRI) project investigating the effect of adrenergic receptor stimulation on physiological, subjective, and neural measures of interoception (ClinicalTrials.gov identifier: NCT02615119). All study procedures were approved by the Western Institutional Review Board, and all participants provided written informed consent before participation. This is hereafter referred to as the ISO study given its focus on the effects of isoproterenol (isoprenaline). For details regarding the development and implementation of this pharmacological perturbation approach, see (Hassanpour et al., 2018;Khalsa et al., 2009b;Khalsa et al., 2016;Teed et al., 2022). Participants in the study were predominantly female (93.3%; n=98) with a mean age of 24.2 (SD=5.7; range=18-39).
The sample was clinically heterogenous combining both healthy (n=45) and individuals with various psychiatric conditions (generalized anxiety disorder [n=29], major depressive disorder [n=7], and anorexia nervosa [n=24]). The experimental procedure entailed the randomized administration of two bolus infusions of saline and two bolus doses of isoproterenol (two micrograms (mcg) and 0.5 mcg, each dose repeated once) during fMRI scanning, though analyses undertaken here excluded data from the 0.5 mcg time series (rationale described in the "Model sensitivity to pharmacological perturbation" section).
Participants were blinded to the order of administration. Importantly, photoplethysmography data were collected using different hardware from that used to develop the model: GE MR750 MRI scanner equipment (i.e., hard plastic finger clip) affixed to the index finger of the nondominant hand and recorded at a sampling rate of 40 Hz while in the supine position.

Summary of the model pipeline.
The processing pipeline employed for the ISO study analysis is laid out in Figure 4 and comprised: (i) Pre-processing (ii) Beat segmentation (iii) Morphological feature extraction (contour analysis) (iv) Pulse decomposition modelling.
All beats underwent (iii) and (iv), with the former applied in order to generate a more "standard" analysis against which the performance of the model could be assessed. The output of the decomposition step was also assessed in terms of goodness of fit to the dataset as a whole (i.e., isoproterenol plus saline) (i) Pre-processing: Many forms of photoplethysmography hardware automatically detrend raw data to maintain a constant baseline. This can confound contour analysis by altering the shape of the waveform. We employed an optional step to reverse the effect of detrending in two stages: (a) Adjusting the fixed amplitude drop with each successive datapoint such that the central portion of the waveform did not drop below baseline, and (b) Correcting for the change in gradient of the entire time series resulting from step (a).

(ii) Beat segmentation:
A custom peak detection algorithm identified peaks in the first derivative of the time series (Elgendi et al., 2018) (sensitivity 0.97, positive predictive value 0.99, further details are provided in Extended data (S1), Williamson, 2022). Identification of the first peak of this derivative (the so-called "w peak") provided a marker of systolic waves and, thereby, a convenient means of calculating inter-beat intervals, and thereby, HR and HRV. The minima preceding w, denoted O, gave an indication of the wandering baseline of the signal, and through spline interpolation were used to remove the baseline drift (or DC component) from the time series. This ensured that only the arterial signal was assessed for any given waveform. O points also defined the start and end of each beat and were used to derive individual beat segments.
(iii) Morphological feature extraction: Each waveform was assessed in turn to derive fiducial points OSND (see Figure 1), with use of the first derivative. In the absence of a clear notch, the inflection point closest to zero in the first derivative was instead taken as N (Lee et al., 2011).
To increase robustness against misidentification of points due to noise, OSND values for the average wave were identified first and informed identification on individual waves. Contour analysis features were then derived from OSND directly or with additional information from the segmented waveform (e.g., area under the wave). All waves were scaled for amplitude. Further derived measures are outlined in Extended data (S2, Williamson, 2022).
(iv) Pulse decomposition modelling: Initial parameter estimation: Initial parameter estimates were calculated per beat for the entire time series and were first estimated from the derived excess of a given waveform. In the excess, the systolic peak was tallest, followed by R2 and then R1. The nine peak parameters could therefore be estimated sequentially. The three decay parameters were initially given fixed values.
Parameter optimisation: Initial parameter estimates were optimised for each waveform (and across waveforms -see the "Refinement of the model and assumptions" section). To achieve this we used the downhill simplex method (Nelder & Mead, 1965), which is a straightforward and robust means of finding the minimum of a function in multi-dimensional space (Press et al., 2007). In this case the function to optimise was goodness of fit, defined for each waveform as the reduced chi-squared (χ²v): Where RSS = The sum of squared (weighted) residuals, V = number of datapoints − number of parameters. Residuals are weighted such that a fitting of the central morphology of the waveform (from w to D) is prioritised.
With no convention for the number of parameters used in PDA models, the model took an adaptive approach. Χ²v, by definition, penalizes non-parsimonious use of parameters. Since this function drives the optimisation process, each waveform was fitted by the minimum number of parameters necessary to generate an optimal fit. In practice, this drove model behaviour such that, for a given waveform, two component waves were selected over three if goodness of fit was not significantly worsened by doing so. This approach is comparable to some existing PDA models (Wang et al., 2013).
Finally, the downhill simplex routine was restarted multiple times to avoid convergence on local minima. Without an established means of assessing how many restarts are sufficient to converge on a global minimum, three restarts (four runs total) were chosen for simplicity and computational efficiency. Optimised parameters from the fourth run were taken as model outputs.

Assessment of model performance.
We aimed to demonstrate the external validity of our model by testing its performance on data from the ISO study (see above). Importantly, the hardware used to collect the photoplethysmography signal for that study was different from the hardware we used to develop the model in the pilot study, allowing an assessment of the generalisability of the model and associated processing pipeline. For each modelled waveform, goodness of fit was assessed. In the absence of a standard measure for goodness of fit in PDA, we used the normalized root mean square error (NRMSE) as reported by Sorelli et al., among others (see Extended data -S3 [Williamson, 2022]). A value of 1 indicates a perfect fit to the data, whilst a value of 0 indicates no improved performance over a null model (Sorelli et al., 2018) .
In a further assessment of the model's performance, we used it to attempt to recapitulate standard morphological points modelled from each waveform in the form of OSND values. Model-generated waves were created and these values were extracted and compared to the corresponding values taken from the real data. This allowed an assessment of concordance between the model and measures of proven clinical relevance.

Model sensitivity to a pharmacological perturbation.
During the ISO study, effects of bolus infusions of either saline or isoproterenol during photoplethysmography measurement were examined. Isoproterenol is a rapid-acting beta-adrenoceptor agonist similar to adrenaline that has been used for decades in psychophysiological research (Cameron & Minoshima, 2002;Cleaveland et al., 1972;Contrada et al., 1991;Khalsa et al., 2009a). It primarily targets peripheral adrenergic β1 receptors, exerting a positive inotropic and chronotropic effect on the heart, and equivalently targets peripheral adrenergic β2 receptors, inducing bronchodilation and peripheral vasodilation. To a much lesser extent it also targets α1 and β1 noradrenergic and α1 adrenergic receptors.
Much work has focused on the effect of cardiovascular drugs on the pulse waveform. Feinberg and Lax discovered that infusions of adrenaline and noradrenaline produced different waveform shapes, despite both agents' overall hypertensive effects (Feinberg & Lax, 1958). Adrenaline (i.e., epinephrine), which is pharmacologically similar to isoproterenol, resulted in generalised vasodilation and a lowered diastolic portion relative to the systolic wave. Noradrenaline (i.e., norepinephrine), which conversely increases peripheral resistance, resulted in a relatively heightened diastolic portion. Other studies have demonstrated similar reductions in relative diastolic amplitude in response to the vasodilatory action of systemic glyceryl trinitrate (GTN) infusion (Chowienczyk et al., 1999;Millasseau et al., 2003). Despite a high level of nuance to these pharmacological effects, the general picture is one of reduced peripheral resistance resulting in a lower relative amplitude of the waveform's reflectance waves.
The effect of isoproterenol-induced β2 agonism on peripheral blood flow has also been directly assessed. Cohen and Coffman confirmed with venous occlusion photoplethysmography that isoproterenol's vasodilatory effect increased blood flow to the vascular bed of the finger itself (Cohen & Coffman, 1981), while Shen et al. also noted a reduction in total peripheral resistance (measured using impedance cardiography) with isoproterenol in the supine position, specifically (Shen et al., 1999). Overall, then, it could be reasonably anticipated that the ISO study data is representative of a typical isoproterenol perturbation, and that this would result in a lowered diastolic portion of the waveform similar to the change seen with adrenaline. We aimed to assess the model's ability to classify waveforms acquired during isoproterenol versus saline infusion, and therefore its ability to effectively characterise a change in waveform shape through parameterization.
After subjecting all data (across both isoproterenol and saline) to contour and pulse decomposition analyses (see 'summary of model pipeline'), we assessed the performance of features derived from both methods as classifiers between dose levels. We aimed to determine (i) if traditional contour features are more sensitive to the pharmacological manipulation than HR or HRV, and (ii) if model parameters are more sensitive to the manipulation than either of the above. For simplicity of comparison and to ensure a reasonable cardiovascular effect size for characterization, photoplethysmography data corresponding to the 0.5 mcg dose level were not included.
The duration and subsequent effects of bolus isoproterenol infusions were short. Accordingly, we analysed the subset of the data comprising only beats from the infusion period of each time series and a matched sample from the saline infusion (for a summary of this process see extended data -S4). Given the risk of sampling bias in time series with few beats to analyse (Allen, 2007), sensitivity analyses were conducted excluding (i) no time series (ii) time series with fewer than 30 identified beats and (iii) time series with fewer than 60 identified beats.

Results
We present our results in two parts:

General model performance: Performance on the ISO test dataset as assessed by (i) goodness of fit (NRMSE) (ii) Recapitulation of important fiducial points (OSND).
Characterization of a pharmacological manipulation: The assessment of model-derived and standard photoplethysmography metrics as classifiers between dose levels.
The ISO study dataset comprised 118 human participants. Of these, six were excluded due to obtrusive clipping of the photoplethysmography signal, one for missing data, and three for suspected mislabelling. An additional three were excluded on running the pipeline due to excessively noisy signals. Thus, a total of 105 participants were included in the final analysis. Each participant had two time series available at each dose level, both of which were included. The mean number of beats sampled was 67 (+/-25), ranging from 12 to 161. Results presented below were consistent across sensitivity analyses.

General model performance ISO data (general model fit).
Waves in the dataset were robustly fitted with an average (median) NRMSE value of 0.9 obtained across participants ( Figure 5), comparable to previously reported models (Sorelli et al., 2018;Wang et al., 2013). See Extended data (Williamson, 2022) figure S1 for alternative NRMSE values). Figure 6 shows the accuracy of the model in regenerating fiducial points OSND. In order to illustrate this, a representative waveform and corresponding OSND values were plotted. Around each point model error (in the form of 65% and 95% confidence boundaries) is represented by green and red ellipses, respectively. A tendency for the model to displace both O and S vertically is apparent, as well as N and D diagonally. Figure 7 shows the average waves for each of the 105 participants' time series (from the subset of data related to the relevant periods of the isoproterenol or saline infusion) as well as an overall average across participants. These are presented separately for the isoproterenol and saline infusions. There appeared to be greater variability in morphology between participants at the saline dose. We therefore used the Levene test to compare homogeneity of variance between dose levels. This was carried out at each time point between w (corresponding to the first peak of the first derivative marking the systolic rise) and the end of the wave. This region of analysis lasted 0.865 seconds and consisted of 346 time points. For an alpha level of 0.05, the Bonferroni corrected p-value was 0.00014. The variance significantly differed across isoproterenol and saline for a continuous period of 0.425 seconds (from 0.233 seconds to 0.658 seconds after w, which corresponds approximately to the peak of R1 and the end of R2).

Characterization of a pharmacological manipulation Ability of HED model to characterise the effects of isoproterenol. For illustrative purposes,
Aside from this difference in variance within conditions, importantly, there were striking morphological differences between drug and saline (summarized in Figure 7c). The isoproterenol (2 mcg) waveform showed a significantly lowered diastolic portion, with a steepened systolic decay and shortened diastolic decay. The overall length of the isoproterenol waveforms also differed, as would be expected with an adrenergicallymediated increase in HR.
Differences in morphology between dose levels were first quantified using standard measures including fiducial point measures relating to contour analyses as well as HR and HRV. SDNN (The standard deviation of the inter-beat interval of normal beats) was selected as the measure of HRV since a number of  Figure S1, Extended data [Williamson, 2022]).

Figure 7. a:
Averaged waveforms (black) during pharmacological manipulation (one per subject/time series) aligned by systolic peak. The mean of these (red) is overlaid. b: Averaged waveforms during infusion of saline (one per subject/time series) aligned by systolic peak. The mean is overlaid in red. c: The mean waves derived from a and b, overlaid to display the morphological change between saline and 2 mcg isoproterenol infusion. beats had to be rejected due to noise. Consequently, measures such as RMSSD (The root mean square of successive differences between normal heartbeats), which are dependent on distances between successive beats, were considered suboptimal. Those that characterized the difference most effectively are included in Figure 8. Consistent with Figure 7c, these  Table 1 for accompanying statistics.

Table 1. Statistical comparisons across isoproterenol and placebo conditions from standard and Hybrid Excess and Decay (HED) model-derived measures.
Standard measures (see Figure 8) HED Model (see Figure 9)   Table 1 for accompanying statistics. (Note that the truncated distribution (especially notable for R1 width) is a result of a constraints applied to parameters to prevent over-fitting. (The algorithm was tending to use parameters for R1 to optimise overall fit).

Parameter t (df) P value Cohen's d Parameter t (df) P value Cohen's d
measures were directly related to either height of the diastolic portion (N amplitude, D amplitude), HR (wavelength) or both (area under the curve). Figure 9 shows the model parameter outputs which most effectively characterized the difference in dose level, namely, via a combination of both excess and decay parameters (for remaining parameters see Extended Data S2 [Williamson, 2022]). The systolic and R2 were modelled as wider at the 2-mcg level, whilst the R1 was variable in width and earlier in time. The first baseline was lower at 2 mcg, and the decay rate was increased. Both contributed to the modelling of a much steeper decay at 2 mcg than at saline.
In Table 1, we summarise the statistics and effects sizes from each of the parameters shown in Figure 8 and Figure 9.
Receiver operating characteristic (ROC) curves were then generated to assess classifying performance of derived features, shown in Figure 10. The five features with the highest Area Under the Curve (AUC) values in 10a appeared to be similar in discriminative ability, which may indicate a redundancy between measures. Wavelength, a measure directly related to heart rate, performed best numerically, suggesting that for the purposes of classifying between stimulation with isoproterenol versus saline, traditional morphological assessment of the waveform was of no added value. SDNN, the standard deviation of inter-beat intervals across the selected region of the time series, performed less well. This may have been due to the limited number of inter-beat intervals recorded over short infusion periods of less than 5 minutes. Figure 10b shows the best performing model parameters.
Decay rate showed the highest AUC value, with, notably, a higher numerical performance in distinguishing isoproterenol from saline than any of the more standard measures (shown in Figure 10a) including wavelength and HR. These findings suggest that the added granularity offered by the proposed model was useful and informative in this instance. This, however, is a qualitative observation. Formal testing using DeLong's test for two correlated ROC curves was used to compare the "best-performing" measure from the HED model (decay rate) and that from more standard approaches (wavelength). These were not significantly different (p = 0.495; Z = 0.682).

Discussion
We present the HED model, with an accompanying automated pipeline for analysing morphological characteristics of the photoplethysmography waveform. While standard approaches focus primarily on basic features extracted from inter-beat interval data, such as HR and HRV, we argue that the waveform itself contains valuable information about underlying physiological states relevant to psychological and human neuroscientific research. Previous analyses of the waveform have used particular landmarks such as the onset, systolic peak, dicrotic notch and diastolic peak (OSND) to identify features of relevance to individual state and trait characteristics, including anthropometrics (Broyd et al., 2005), cardiovascular health (Girón et al., 2019) and current level of arousal (Arza et al., 2019). The HED model was motivated by this early promise and by the belief that a comprehensive, more nuanced, and automated approach to characterising the waveform would help to maximize the full potential of this signal. Below, we consider the HED model in relation to previous approaches to the waveform, its biological plausibility, and its sensitivity to a relevant pharmacological perturbation. We acknowledge that this model must be seen as early in development and consider its limitations before highlighting areas where we feel that the current analysis approach will be of value to neuroscience and psychiatry.
By formally modelling a decay component, in addition to its component waves, along with shifting baselines across the cardiac cycle, the model extends previous PDA approaches (Fleischhauer et al., 2020). A key test of the model lay in its application to an independent dataset. The average goodness of fit achieved by the model is comparable to that of other PDA models (Sorelli et al., 2018) though a standardized means of comparison has not been established. Additionally, fits were equivalent across both isoproterenol and saline infusions. In classifying between dose levels, one model parameter -the decay rate -showed numerical superiority over HR, HRV, and contour analysis measures in terms of ROC and effect size measures (see Table 1). Though it is important to note that this particular parameter did not show statistical superiority, our findings suggest that thoroughly characterizing the peripheral pulse waveform can yield useful information that are not captured by typical photoplethysmography-derived measures of HR and HRV. It is also notable that the dose of isoprenaline used for this manipulation (2 mcg) produces a marked increase in HR (by an average of 26 beats per minute in the current dataset). We might speculate that when assessing more subtle experimental manipulations, such as the imposition of psychological stress, which produce much smaller -if any -differences in heart rate, the parameterised morphology of the pulse wave may become increasingly valuable. Initial studies of waveform morphology in distinguishing physiological changes induced by cognitive demand would suggest so (Pelaez et al., 2019).
It is important to note that our model was informed by the existing literature on physiological factors shaping the waveform, and was only subsequently modified to best capture the actual data. We can only speculate therefore on how the parameters of the model, including the newly added decay element relate to specific physiological indices. Indeed, it seems likely that the parameters that emerged in the development of the model are influenced by several interacting factors which, themselves, are influenced by the pharmacological challenge that we sought to capture. Nonetheless, it is worth briefly considering the parallels between the HED model and existing concepts in physiology.
The HED model produces a waveform characterisation that is very similar to an existing model from the field of arterial haemodynamics. The hybrid reservoir-wave model (Wang et al., 2003) aims to account for aortic pressure waveform morphology by decomposing it into 'reservoir' and 'excess' pressures. Doing so enables both Windkessel and wave reflectance properties of the arterial tree to be captured. According to the hybrid model, the reservoir pressure accounts for arterial wall expansion properties and is constant throughout the arterial tree, while excess pressure accounts for local changes in pressure due to propagating waves at the measurement site, including the initial systolic and ensuing reflectance waves. Though lacking widespread acceptance (Mynard & Smolich, 2017;Segers et al., 2017), there is growing evidence that the hybrid model's parameters are of cardiovascular prognostic value (Hughes & Parker, 2020). Our model shares features with the reservoir-wave model, but it is important to note a key difference: the HED model relates to a peripheral (i.e., finger) measurement. Despite evidence of some equivalence between peripheral arterial and aortic measures (Peng et al., 2018), the photoplethysmography signal measures blood volume, rather than pressure, and at low wavelengths of light may incorporate time-delayed signals from different tissue depths (Liu et al., 2019). Nonetheless, the conceptual similarity to a model that has proven useful in central arterial studies lends credence to the one presented here.
The two-element Windkessel model (Westerhof et al., 2009) provides a useful framework for contextualizing our finding of the numerical superiority obtained by including the decay rate parameter. It models the time duration of exponential pressure decay in diastole (RC time) as the product of arterial compliance (C) and peripheral resistance (R). The model can therefore characterize changes in either R or C in terms of RC time. Since exponential decay rate is inversely related to exponential decay time, it is plausible that our modelled decay rate approximates RC time (insofar as photoplethysmography decay mimics central pressure decay). With this framework, we can link the reduced peripheral resistance induced by isoproterenol (Green, 1977;Smith et al., 1967) to the change in waveform morphology effectively characterized by decay rate. Furthermore, if decay rate is a parameter sensitive to key features of cardiovascular state (R and C), then it may in turn be sensitive to measures related to these (e.g., cardiac output, pulse pressure, stroke volume). Thus, decay rate has plausible links to peripheral cardiovascular features beyond those directly pertaining to the heart, which may explain its superior performance in distinguishing isoproterenol from saline infusion.
The model we propose, then, could be interpreted as following a simple underlying principle: An initial pressure wave, corresponding to cardiac systole, decays exponentially. The nature of this decay is determined both by arterial compliance and by ensuing reflectance waves which are, themselves, influenced by compliance and a number of other features relating to an individual's cardiovascular state. While we have examined the impact of a pharmacological perturbation on a within-subjects basis, it will be interesting in the future to determine whether the model is sensitive to other variations across individuals and groups.
There are several limitations that should be considered. The data used in the development of the proposed model were obtained from two different hardware set-ups, both of which were transmission-based. Differing placement of photoplethysmography devices, and the use of reflectance rather than transmission-based measures will have an impact on the measured pulse waveforms. Whilst our model has demonstrated robustness to different morphologies of waveform, we cannot confirm its validity on sources such as these, from which different signals may arise. Related to this, the population of participants in the ISO dataset were predominantly younger and female, and a range of psychiatric conditions were also represented. As with more widely used measures such as HR, characteristic waveforms will likely differ across sex and age (which themselves may influence the effects of isoproterenol (Stratton, 1999)) as well as other factors such as height, and it remains to be demonstrated that the model pipeline will generalise across these factors.  Lee et al., 2011;Millasseau et al., 2002;Takazawa et al., 1998) and mental (Kontaxis et al., 2021) health suggesting that this characterisation of the waveform across groups may prove informative.
In addition, it is likely that the performance of our model could be optimised further. When increasing degrees of freedom with additional parameters, the risk of overfitting should be considered. Standard deviations of waves fitted by the model indicated an increased within-participant variability in the R1 Time and R1 width parameters at the 2-mcg dose level (see extended material -S5). On visual inspection of these time series, a marked variability in the beat-to-beat fitting of the first reflectance wave was noted in some cases, indicating a tendency to optimise fit at the expense of consistency. This behaviour only becomes apparent when the diastolic portion of the waveform has dropped below a certain threshold, when the position of a first reflectance wave in the data becomes visually ambiguous. This limitation could be overcome by an initial calibration period at baseline when the parameters of R1 can be ascertained more reliably. Stricter constraints could then be applied to any periods of deviation from baseline. We suggest any future implementations of this model take such an approach to avoid similar instances of over-flexibility.
Despite these limitations, we argue that there are several good reasons for pursuing a robust, automated approach to characterising the morphology of the pulse waveform generated by standard photoplethysmography hardware. Given that waveform morphology reflects vascular characteristics which are the subject of several controlling factors, its study opens potentially important avenues of research. As an example, there has been increasing interest in the influence of the intra-cellular Adhesion Molecule-1 (ICAM) in a range of psychiatric disorders (Müller, 2019) and given too that this molecule is an important regulator of vascular elasticity, it may be that the pulse waveform offers itself as a valuable marker for alterations in ICAM activity in certain psychiatric illnesses.
The possibility that the waveform may signal subtle cardiovascular changes that occur in association with longer-term changes in psychological health has already been examined in the context of, for example, depression (Dregan et al., 2020) and observations that other conditions, such as anorexia nervosa are associated with altered arterial stiffness (Jenkins et al., 2021;Tonhajzerova et al., 2020) hint at its potential usefulness in psychiatry. Indeed, features present in the waveform have been related to aspects of mental health in the case of depression, but cardiovascular changes are also well recognized hallmarks of other long term mental illnesses such as schizophrenia and AN.
With the increasing prevalence of photoplethysmography-based wearables, the pulse waveform may offer a means of remotely monitoring physical and mental health (Myin-Germeys et al., 2018), its ease of use making it applicable at a large scale (Natarajan et al., 2020).
Related to the above, the growing field of interoceptive neuroscience and its focus on measures that characterize subjective and objective elements of brain-body interaction would appear to be complementary to the work presented here (Allen, 2020;Allen et al., 2020;Critchley & Garfinkel, 2015;Khalsa et al., 2018). Heartbeat detection tasks have been used frequently to measure interoceptive perceptual ability (Brener & Ring, 2016) since heartbeats represent a convenient interoceptive signal, and subjective perception of heartbeats reported by participants can be compared with objective measurement of these beats. However, photoplethysmography measures that go beyond simple indices of systole and diastole in relating interoceptive abilities to psychological states may be of value in developing our understanding of links between cardiovascular states and interoception (Schandry et al., 1993), including studies that employ peripheral perturbations to modulate pulsatility of the vascular tree in healthy (Khalsa et al., 2009b), psychiatric (Khalsa et al., 2015) and neurological populations (Khalsa et al., 2009a;(Khalsa et al., 2016), and those examining the accuracy and precision of metacognitive beliefs (Garfinkel et al., 2016;Legrand et al., 2022), possibly in combination. While the subjective experience that forms the basis for successful cardiac interoception (i.e., 'cardioception') is not clear, it seems feasible that the pulsatility of the vascular tree contributes to it. As such, we should consider the possibility that individual variation in cardioceptive abilities does not rely solely on higher-order sensitivity but also on peripheral influences on the signal itself, with the nature of pulse wave morphology potentially contributing to this. Characterising wave differences across individuals will thus be important when relating cardiac interoceptive abilities to psychological states, traits, and psychiatric disorders (Domschke et al., 2010). Moreover, relating characteristics of the waveform to the nature and magnitude of heartbeat evoked potentials (Babo-Rebelo et al., 2016;Coll et al., 2021;Park & Blanke, 2019;Pollatos & Schandry, 2010) may prove useful in understanding generators and modulators of these potentials. Crucially, too, if the waveform differs in individuals with psychiatric illness, it may be that their altered interoceptive capacity is secondary to an underlying physiological alteration, which might explain the limited benefits of extant attempts at cardioceptive biofeedback training (Rominger et al., 2021). Future work should focus on validating the model with additional datasets and assessing its ability to characterize a range of psychophysiological states across individuals. Obtaining reliable signals of sufficient quality will be important for this. Reassuringly, this in itself is a field of great interest and innovation (McDuff et al., 2020). Moreover, existing early observations that features of the waveform may outperform more standard measures like HRV (Pelaez et al., 2019) are promising in this regard. Thus, the more fully we understand the photoplethysmography signal, as well as how it varies across individuals, and across physiological and psychological states, the more realistically we might be able to understand interoceptive abilities and perhaps even develop a clearer idea of the causal chain between brain and body feedback interactions.
Overall, we argue that the shape of the pulse waveform, which has tended to be overlooked (Murray & Foster, 1996), contains information that will prove useful in a number of settings, including assessments of psychological states of arousal.
We present the HED model, a robust, automated and freely available model for characterising this waveform and show that its parameters compare favourably with more established ones such as HR and HRV in distinguishing the effects of a sympathomimetic drug from a saline control. In this respect, we show that the waveform conveys valuable information on acute changes in cardiovascular state, whilst extending existing approaches using pulse decomposition. How the parameters of the waveform perform in distinguishing states that differ more subtly than the drug effects seen here remains to be seen. The fact that features of the wave are influenced by a number of relevant factors give some grounds for supposing that it may provide a rich and sensitive objective index of physiological or psychological arousal.

Data availability
Underlying data The data in this paper were obtained from the Laureate Institute of Brain Research Database supported by grants from the National Institute of Mental Health (grant K23MH112949 to Dr Khalsa); National Institute of General Medical Sciences (grant P20GM121312 to Drs. Khalsa and Paulus). They were accessed by a Data Transfer Agreement with Dr Martin Paulus and Dr Sahib Khalsa, contactable on skhalsa@laureateinstitute.org.

Extended data
Data that is representative of the analysed dataset (along with all the analysis code) can be used to apply the methodology herein described and is available at: https://github.com/stw32/PulseWaveform (Williamson, 2022) Data are available under the terms of the BSD-3-Clause license.