A Bayesian Filtering Approach for Tracking Sympathetic Arousal and Cortisol-Related Energy From Marked Point Process and Continuous-Valued Observations

Multiple state variables governed by internal processes within the human body remain unobserved. On a number of occasions, these states are linked to point process bioelectric and biochemical observations coupled together with continuous-valued variables. These observations provide a means to estimate the latent states of interest. We develop a state-space method to estimate unobserved sympathetic arousal and energy production states from skin conductance and cortisol data respectively, comprising of a marked point process and a continuous-valued observation. The method involves Bayesian filtering applied within an expectation-maximization (EM) framework for state estimation and model parameter recovery. Results are evaluated on both simulated and experimental data. On experimental skin conductance data, high arousal levels are generally detected during cognitive stress periods and lower values are detected during relaxation. Results are also in conformity with general physiology for cortisol data. On separate experimental data, skin conductance-based estimates are validated/cross-checked with functional Near Infrared Spectroscopy features. Estimation is also performed on heart rate and skin conductance data to illustrate the method’s wider applicability. We also compare the method with earlier approaches. We show how it outperforms a previous method for cortisol-based energy estimation, and its superiority to earlier methods for estimating sympathetic arousal. The EM approach is thus able to estimate latent physiological states within the body from point process bioelectric and biochemical phenomena. The method could be applied in wearable monitoring and automated closed-loop therapy delivery for patients diagnosed with certain types of neuropsychiatric or hormone disorders.


I. INTRODUCTION
Several internal state variables maintained by physiological processes within the human body remain largely unobserved.Nevertheless, they are often tied to more easily-measured changes in biolelectric and biochemical phenomena that can be used for their estimation.For instance, energy metabolism and cognitive processes are difficult to observe directly.However, the changes they are linked to in hormones and other signals (e.g.heart rate, blood pressure, sweating) can be conveniently measured and used for their estimation.In a number of instances, the bioelectric and biochemical observations are pulsatile or impulse-like in nature.Variations in skin conductance are one such example.Bursts of neuroelectric activity to the sweat glands cause changes in the skin's conductivity and give skin conductance signals a ''spikey'' appearance.A skin conductance signal consists of a fast-varying phasic component superimposed on top of a slower-varying tonic level [1].Individual skin conductance responses (SCRs) make up the phasic component.Each SCR can be modeled as being generated by a single neural impulse to the sweat glands [2], [3].Since the sweat glands are innervated by sympathetic nerve fibers [4], a skin conductance signal becomes a sensitive index of sympathetic arousal [5], [6].Markers from both phasic and tonic components thus reflect changes in arousal.This has led to the use of skin conductance in a number of studies involving emotional arousal [7].The underlying pulsatile/impulse events can be extracted via deconvolution.The upper sub-panel depicts a skin conductance signal (blue), its tonic component (black) and the neural impulses underlying phasic variations (red).The lower sub-panel depicts blood cortisol measurements taken at 10 min intervals (blue dots), the reconstructed minute-wise blood cortisol concentrations (black) and the pulsatile secretions (red).
Cortisol is an example of a biochemical signal of similar nature.It is secreted in pulses and between 15-22 secretory events typically occur each day in a healthy adult [8], [9].A major function of cortisol is to raise blood glucose levels in response to external stressors [10], [11].Cortisol aids in energy production at the liver [12], [13].Deficiencies in glucocorticoid secretion are therefore frequently associated with symptoms of fatigue, weakness and tiredness [14].This can be the case, for instance, in conditions such as primary adrenal insufficiency.The number, timing and amplitudes of the pulses capture important information regarding the underlying energy production states that cortisol is linked to.Pulsatility, in particular, is a means of rapidly increasing hormone concentrations and providing signals to target cells [15].For both skin conductance and cortisol, a deconvolution procedure can be used to extract the underlying pulsatile or impulse events (Fig. 1).In the case of cortisol, the deconvolved pulsatile secretions are, in reality, abstractions of CRH (corticotropin-releasing hormone) pulses.It is these neuroendocrine CRH pulses that are ultimately responsible for cortisol secretion [16].
The constituent components of both skin conductance and cortisol are primarily in the form of a marked point process (MPP) coupled together with a continuous-valued variable (Fig. 1).For skin conductance, the neural impulse profile (responsible for phasic SCRs) and the tonic component form the MPP and continuous variable respectively.Similarly for cortisol, the pulses form the MPP and the blood cortisol concentrations make up the continuous-valued variable.

II. RELATED WORK
Physiological state-space models have often been used to estimate latent variables within the human body.These include states that are linked to point process phenomena.Methods in this family of point process estimators have found applications across the following fields: behavioral learning [17], [18], [19], sleep studies [20], movement/position decoding [21], [22], [23], [24], [25], [26], anesthesia and comatose state regulation [27], [28], [29], and heart rate variability analysis [30], [31].Bayesian filters are typically developed within an expectation-maximization (EM) framework in these methods.This is for both estimating the latent states and recovering unknown model parameters.
We primarily focus on the case of cortisol and skin conductance here.With cortisol, labeled data is scarce.Likewise, there is an accompanying scarcity of research into point process-based energy state estimation for cortisol (our previous work in [32] being one of the few).Unlike the case of cortisol (a biochemical signal), more work has been conducted on bioelectric skin conductance.This is especially so within the field of automated emotion recognition where important applications are present.Here, subject feedback (e.g. on rating scales) is typically used as labels for classifying features from a range of physiological signals including skin conductance.Initial work largely relied on earlier generations of machine learning methods including support vector machines, Gaussian naive Bayes classifiers etc. (e.g.[33], [34], [35]).Deep learning and transformer-based methods have more recently been applied to automated emotion recognition (e.g.[36], [37], [38]).

A. RELATED WORKS AND LIMITATIONS -MACHINE LEARNING METHODS
Due to the scarcity of research into energy state prediction based on cortisol measurements using machine learning, here we instead focus on emotion recognition using physiological signals in reviewing the relevant literature.We focus on emotion recognition since it is closely related to one of our own objectives, namely, estimating sympathetic arousal from skin conductance.Machine learning research in this field can broadly be divided into works focusing on continuous-valued prediction and discrete category prediction.The literature here is quite broad and we summarize noteworthy aspects considering representative examples under the two categories.
There are a few well-known datasets for automated emotion recognition in both categories mentioned above.The RECOLA [39] and CASE datasets [40] for instance, provide continuously-annotated metrics of emotion.The RECOLA dataset contains continuous emotion ratings from external observers along with physiological recordings of subjects including electrocardiography (EKG) signals and skin conductance.Here, observers annotated psychological valence and arousal continuously within a range of (-1) to 1 using a slider.Meanwhile in the CASE study, the subjects themselves provided annotations of valence and arousal using a joystick.
Subsequent machine learning works utilizing these datasets can also be divided based on the types of methods used for classification.While initial work made use of earlier machine learning algorithms, more recent methods utilize deep learning and transformer-based methods to decode emotion.For instance, Ringeval et al. [41] investigated the use of EKG, skin conductance and audiovisual modalities for emotion recognition using the RECOLA dataset with long short-term memory recurrent neural networks.They reported concordance correlation coefficients of 0.804 for arousal and 0.528 for valence in predicting emotion.Despite the continuous nature of the emotion labels in the RECOLA and CASE datasets, some works utilize these to only predict discrete categories.Zhang et al. [42] for instance, used deep multiple instance learning to classify labels in the CASE dataset into the low, neutral and high categories.Here, classification was based on features extracted from skin conductance, blood volume pulse (BVP), skin temperature and heart rate.The authors employed weakly-supervised training on post-stimuli labels to do so and reported classification accuracy values of 75.63 % and 79.73 % for valence and arousal respectively [42].Several other methods also utilize the RECOLA and CASE datasets for emotion prediction using machine learning methods (e.g.[43], [44], [45], [46]).
Other datasets for emotion recognition containing ratings on discrete-interval scales are also available.Among them are the K-EmoCon [47], MAHNOB-HCI [34] and DEAP datasets [33].The K-EmoCon dataset contains psychological arousal and valence ratings on scales of 1-5 along with physiological signals including skin conductance, BVP, heart rate and EKG.Emotion annotations from subjects, partners and external observers are also available [47].
In the MAHNOB-HCI study, subjects were presented with 20 emotional videos and were then required to fill-out selfassessments based on their evoked emotional responses.Valence and arousal had to be rated on a scale having nine options ranging from low to high [34].Likewise in the case of the DEAP study, participant emotional ratings on a scale of 1-9 along with electroencephalography (EEG) and peripheral physiological signals were collected [33].A number of machine learning works have applied thresholds to these subject-provided annotations to classify physiological signal features into discrete emotion categories (e.g.[45], [48], [49], [50], [51], [52]).Li et al. [48] for instance, recently performed multicategory emotion recognition by learning discriminative graph topologies in EEG brain networks.They classified emotion into four categories for both single and multiple subjects with the MAHNOB-HCI and DEAP datasets.Their model, which involved multiple emotion-related spatial network topology patterns, was able to reach an online experimental accuracy of 84.56 % for 14 participants [48].Another recent study also proposed a transformer-based selfsupervised framework to decode emotion using peripheral physiological signals such as skin conductance, BVP and skin temperature, and was tested on the CASE and K-EmoCon datasets for automated emotion recognition [45].
Supervised learning methods such as these, all require labels.Now acquiring labeled data can be resourceexpensive.Consequently, data can be scarce, especially in the clinical domain.For instance, the blood cortisol dataset we subsequently utilize does not contain any labels and prevents the use of traditional supervised learning techniques on it.Moreover, in cases where continuous-valued prediction is required, the labels cannot merely be in high vs. low (or other discrete) categories.Rather they must be on a more continuous-valued scale.Not only are labels required for supervised learning, they can also suffer from inter-subject variability [33].Assume, for instance, that ratings are acquired on a scale of 1-9 for emotional valence or arousal in an experiment from different subjects.A rating of ''6'' for one individual may not necessarily mean the same intensity/level for another.Such inconsistencies can exist even within the same subject.If multiple trials are present in an experiment for example, a subject may not provide ratings in exactly the same manner from trial to trial.Such factors can complicate subsequent model training due to inconsistent labeling.
A number of machine learning methods such as neural networks and support vector machines also have difficult interpretability.Given a set of features and a trained model, it is somewhat challenging to explain what exactly a model is doing to provide a prediction.This can especially be problematic when explaining a model to clinicians.Several supervised learning methods also have the general drawback of requiring a large amount training data.This is typically the case with neural networks which have gained much recent popularity.Such models consequently take longer to train.Moreover, a number of machine learning methods segment physiological data into small time windows and provide a prediction for each window independently.However, internal physiological states may not be completely independent of each other during consecutive windows.Rather, emotional states and feelings of energy/lethargy may evolve over time gradually and hence may not necessarily be independent of one another from one moment to the next.The Bayesian point process state-space method we propose here does not require labels.As such, it does not suffer from a number of the drawbacks just stated.Moreover, the method can be executed even on a single data record and does not require a large amount of training data.Time dependencies between states are also explicitly modeled and taken into account during estimation.

B. RELATED WORKS AND LIMITATIONS -STATE-SPACE METHODS
We have already developed several such Bayesian point process-based state-space estimators for both skin conductance and cortisol previously [32], [53], [54], [55], [56], [57], [58].We thus provide a brief comparison with our earlier approaches motivating the need for the model we propose here.Of the methods based solely on skin conductance, the purely MPP-based estimation method in [58] outperformed two earlier ones in [53], [54], and [55].This method estimated arousal solely from the neural impulses underlying a skin conductance signal considering it as an MPP.Details of the comparisons can be found in [58].The first of the two earlier methods described in [53] and [54] discarded the mark amplitudes and only considered the discrete impulse (or SCR peak) occurrences.The second one in [55] used the same discrete events along with two continuous-valued signals.The first signal was the tonic component and the second was an interpolation over the SCR peak amplitudes.However, such a continuous signal interpolating over the peaks clearly does not match the inherent constituents of a skin conductance signal (as shown in Fig. 1).Despite its superior performance compared to [53], [54], [55], the purely MPP-based estimation approach in [58] too did not conform to this inherent nature.
We only developed one previous state-space estimator for cortisol [32].It too utilized two hypothetical continuous-valued signals derived from the upper and lower blood cortisol envelopes.The point process comprised of the cortisol pulses.Again, this too is unlike the inherent cortisol signal constituents in Fig. 1.Now biologically-active unbound cortisol in the bloodstream, which forms only a single continuous-valued variable, contributes to energy production [59].Consequently, the purely MPP-based method in [58] too cannot viably be applied to cortisol.Considering these factors as well as the similar dynamics of cortisol and skin conductance [2], [8], we require a different method-one which is capable of estimating an underlying latent state based on an MPP and a single continuous-valued signal.
We also note a further limitation of some of our earlier state-space estimators in their tendency to overfit to the continuous-valued variable(s).In our earlier methods in [32], [55], and [56], we used a form of early stopping to mitigate this effect.However, there was no means to control or slow-down the rate at which overfitting occurred.Consequently, it was somewhat difficult to prevent overfitting if it began to occur early.This factor also necessitates an approach whereby better fits can be obtained to the data, especially in the presence of continuous variables.As we show later on, our current method helps address this.
In a very recent work [60], we adapted the neural network approach for state-space estimation in [61] to cortisol and skin conductance.We used it to estimate cortisol-related energy production and sympathetic arousal from the MPP and continuous-valued signal formulation just described.This is a non-traditional use of neural networks and does not require labels in the typical sense.However, as we pointed out in [60], the neural network method of estimation does have certain drawbacks.For instance, different neural network configurations can result in inverted point process event probabilities and states being learned.This is contrary to physiological expectations.Additionally, the final estimates also had small noise fluctuations in them.This was likely due the large number of parameters the neural networks had to learn from limited data.Here we develop a simpler, EM equivalent of the state-space model in [60].As we point out in our results later on, this method does not have these drawbacks.Furthermore, the method does not suffer from the typical interpretability challenge that neural networks have.

C. PROPOSED APPROACH
The rate of SCR occurrence, their amplitudes and the tonic component are among the three most frequently used skin conductance markers of autonomic arousal in the literature [62].Recall that the SCR rates and amplitudes are related to the rates and amplitudes of their underlying neural impulses.Likewise for cortisol, both the pulsatile secretory events and the blood cortisol concentrations reflect changes in the underlying metabolic energy state process.Our previous EM-based point process state-space estimators for skin conductance and cortisol have not considered this inherent MPP coupled with together with continuous signal formulation [32], [53], [54], [55], [56], [57], [58].They thus either ignore certain types of information altogether or use interpolations that do not match the form of the input data.
Here we present a state-space method to estimate sympathetic arousal and energy production from skin conductance and cortisol respectively based on this formulation.The approach could find applications in wearable monitoring and automated closed-loop therapy delivery.The subsequent section describes our methods.We thereafter provide the results we obtained on both simulated and experimental data.A discussion section follows.We conclude making some final observations and note future directions of research.A list of abbreviations we use here is provided in Table 1.

A. DATA
We primarily use two datasets (skin conductance and cortisol) involving multiple subjects to evaluate our methods.Supplemental to this however, we also use data from two additional sources (single subjects).This is for the purpose of cross-checking/validation and to demonstrate the method's wider applicability.The four datasets are described below, the first two being the primary ones.

1) SKIN CONDUCTANCE DATA
We first use the Non-EEG Biosignals Dataset for Assessment and Visualization of Neurological Status [63] for the estimation of sympathetic arousal from skin conductance.Here we use data from the same set of subjects as in our previous work in [54] since a number of records are contaminated with noise.We also focus on the same portion of the experiment as in [54] consisting of a task instruction period, two cognitive stressors, relaxation and an emotional stressor.The two cognitive stressors consisted of a counting task and the Stroop color-word association test.A clip from a horror movie was used to create emotional stress.The cognitive stress period, relaxation and the emotional stress period each lasted for approximately 5 min [63].The instruction period just prior to the cognitive tasks was also categorized among the stressors by the authors [63].We downsampled the skin conductance data to 4 Hz, lowpass filtered at 0.5 Hz and separated out the tonic and phasic components using cvxEDA [64].We next deconvolved the phasic component using the method outlined in [54] to yield the underlying sequence of neural impulses.

2) CORTISOL DATA
We use blood cortisol measurements from the study described in [65] for energy estimation.In the study, blood samples were drawn intravenously at 10 min intervals over a 24 h period from both patients and matched controls.These were later processed to yield cortisol concentrations.Pednekar et al. [66] deconvolved the blood cortisol measurements from 31 of the 36 data record pairs that were available.Deconvolution yields the underlying pulsatile secretion timings and amplitudes at a 1 min resolution.It also yields the cortisol infusion and clearance rates necessary to reconstruct a minute-by-minute serum cortisol profile (to be taken as the continuous-valued variable).Thirteen of the 31 patients were diagnosed with chronic fatigue syndrome (CFS -we label this the CFS patient group).The remaining 18 subjects were diagnosed with fibromyalgia syndrome (FMS) or both FMS and CFS (which we label the FMS patient group).In the FMS patient group, patients 1 to 3 had only FMS while the other 15 had both FMS and CFS.

3) FNIRS AND SKIN CONDUCTANCE DATA
We conducted an experiment in which a subject had to perform a series of n-back tasks first while listening to calming music and secondly while listening to vexing music.In an n-back experiment, the subject is shown a sequence of characters on a screen and has to correctly recall whether the current character is the same as the n th previous one.There were eight 1-back trials and eight 3-back trials in each music session, and a brief relaxation period was provided between the two sessions.Each trial consisted of a 22letter sequence.The music was chosen by the subject.Functional Near Infrared Spectroscopy (fNIRS) signals and skin conductance were recorded.These were recorded using the NIRSport2 headset and Biopac sensors respectively.The study protocol was approved by the University of Houston Institutional Review Board.We processed the fNIRS data using the nirsLAB 2019.4software.We first bandpass filtered the fNIRS signals between 0.01 and 0.2 Hz (default settings in nirsLAB 2019.4) and then computed the hemodynamic states.We next extracted the total hemoglobin (hbT) signals over the pre-frontal cortex (two bad channels were discarded based on a default nirsLAB coefficient of variation criterion of 7.5%).We also used cvxEDA [64] for tonic-phasic separation of the skin conductance data followed by deconvolution [54].

4) HEART RATE AND SKIN CONDUCTANCE DATA
We finally use the Smart Reasoning for Well-being at Home and at Work (SWELL) dataset [67].We use this for estimating 137208 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
a quantity similar to sympathetic arousal from heart rate and skin conductance data.The subjects in the experiment had to perform office-like work (e.g.writing reports, preparing presentations).The tasks had to be first completed in the absence of any stressors and then in the presence of two different stressors-time pressure and e-mail interruptions.
Here we use the data from the same subject we did in our prior work in [57] (subject 16) during the first session devoid of time constraints or e-mails.We detected R-peaks using MATLAB's findpeaks() function and manually corrected erroneous beats.R-peaks accompany ventricular contractions in an EKG signal.We next re-sampled the RR-interval sequence at 4 Hz through interpolation [68].We then performed the continuous wavelet transform to extract a moment-by-moment low-to-high (LF to HF) frequency energy ratio.This type of LF to HF ratio of the RR-interval spectrogram has been used in a number of studies as an indicator of stress levels [69].The LF to HF ratio time series has a somewhat ''spikey'' appearance.We extract peaks in it (i.e., an MPP) with a simple threshold applied to avoid small values.LF and HF energies were calculated after re-sampling using MATLAB's cwt() function with the morse option.We again used cvxEDA [64] for extracting the tonic skin conductance component with data downsampled to 4 Hz.

B. STATE-SPACE MODEL
Random walks and first-order autoregressive models have frequently been used in the estimation of latent physiological states underlying point process behavioral and bioelectric measurements [17], [18], [19], [20], [70].Here too we assume that our unobserved state x k evolves with time as where ε k ∼ N (0, σ 2 ε ) is process noise.We first consider just the occurrence of the point process events n k = {0, 1} in the MPP.Now n k alone is a Bernoulli-distributed random variable with mass function p n k k (1 − p k ) 1−n k where p k = P(n k = 1).Based on the theory of generalized linear models [71], we apply a logit transformation to relate Here, β 0 and β 1 are constant coefficients to be determined.When β 1 is positive, the probability of a point process event occurring rises with an increase in x k .Next we consider the amplitudes of the MPP observations r k .We assume a linear relationship between r k and x k .
where v k ∼ N (0, σ 2 v ) is sensor noise, and γ 0 and γ 1 are coefficients to be determined.We assumed an identical linear relationship between x k and r k in our previous work in [58] (MPP-based estimator) for skin conductance.Here we apply the same formulation for cortisol as well, assuming that the pulse amplitudes are linearly related to the energy state x k .Similar linear-type relationships have also been assumed in other physiological point process state estimation problems (e.g.[18], [19]).The joint density for the MPP observations is [58] We finally consider the continuous variable s k .We assume that it too is linearly related to x k .Therefore, where w k ∼ N (0, σ 2 w ), and δ 0 and δ 1 are similar to γ 0 and γ 1 .Again, this is identical to our formulation in [55].We also assumed similar linear relationships between x k and the upper and lower blood cortisol envelopes in [32].This was for estimating energy.While the upper and lower blood cortisol envelopes remain hypothetical in [32], here we assume a direct relationship with the actual blood cortisol level itself.Despite these linear relationships assumed between x k and some of the observations, the manner in which x k is related to p k causes the overall state-space model to be nonlinear.
In the case of skin conductance, the neural impulses form the MPP observation pairs (n k , r k ) and the tonic component makes up s k .In the case of cortisol, the pulsatile secretions make up (n k , r k ) and s k consists of the minute-by-minute blood cortisol concentrations.Fig. 2 graphically illustrates how x k evolves over time and the observations it gives rise to.

FIGURE 2. Graphical Representation of Latent States and Observations.
The latent state x k gives rise to the MPP observations (n k , r k ) and the continuous-valued observation s k as it evolves over time.
We estimate x k and determine the model parameters using an EM algorithm.We provide the derivations of the algorithm in Appendices A and B where we follow the general approach in [19] for obtaining the E-step and M-step updates.Taking y k to collectively denote the MPP and continuous-valued observations for convenience, we derive the E-step filter updates based on the posterior density shown below.

p(x k |y
Likewise, the M-step updates are derived based on the following likelihood term.
where represents the model parameters, K is the set of indices at which n k = 1 and K denotes the final time step index.

C. EM ALGORITHM
The E-step comprises of both a forward filter and a backward smoother.

1) E-STEP -FORWARD FILTER
The filter update equations switch between those given in [18] and [55] depending on n k .The filter equations for k = 2 : K are shown below.Predict: Update: If ) Now appears in both ( 12) and (15).Performing this substitution causes x k|k to appear on both sides of the equations and they are thus solved numerically using Newton's method.Note that the E-step equations switch between two possibilities based on the value of n k .Despite this, our model is unlike the conventional switching state-space models described in [72].Rather, it is an extension of our previous model in [58] where a similar switching occurs.Here we include s k as well.We do not assume that a set of discrete variables actually underlies x k , which then causes the state-space model itself to switch between possibilities over time.Instead, the model remains the same while the equation switching merely occurs due to the point process nature of the (n k , r k ) observations.

2) E-STEP -BACKWARD SMOOTHER
The backward smoother equations to obtain x k|K and σ 2 k|K are shown below [19].
3) M-STEP Making use of the state-space covariance algorithm in [73], we first define Therefore, the (m + 1) th M-step updates for the model parameters are Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Two alternative methods can be used to estimate β 0 and β 1 .First, they can be calculated at the M-step by solving the following two simultaneous equations.
A second approach is to use the equation set β 1 = 1, assume that x k ≈ 0 at the very beginning and calculate β 0 empirically as In this second approach, p 0 is approximated by the baseline probability of point process event occurrence p0 in the complete signal for each subject.Both these methods for estimating β 0 and β 1 have been used in the literature [17], [18], [19].We too have utilized both methods when estimating sympathetic arousal and cortisol-related energy in our previous work [32], [53], [54], [55], [56], [57], [58].
In our experience, estimating β 0 and β 1 at the M-step provides better estimates for p k in the case of cortisol since the pulses are more sparse.The alternate approach of calculating them empirically works well for skin conductance.Note that all the model parameters (including β 0 and β 1 ) are calculated for each subject separately.This ensures better personalization and the ability to handle inter-subject variability.
Similar to [55], we consider the EM algorithm to have converged once the mean difference of the model parameters in successive iterations does not exceed 10 −8 .

A. SIMULATED DATA
We first simulated two sets of data to check the ability of the EM algorithm to estimate an unobserved state x k and recover model parameters from a set of observations comprising of an MPP and a continuous-valued variable.Recall that β 0 and β 1 can be calculated in one of two different ways.Here we use the second approach similar to [55], [56], [58] to illustrate the performance of our method.In this approach, we set β 1 = 1 and calculate β 0 empirically based on the starting probability p 0 .Assuming that x k ≈ 0 at the very beginning, we have . While p 0 is unknown, it can be approximated by p0 which is the average probability of a point process event occurring in the data.Assume, for instance, an experiment comprising of both stress periods and relaxation intervals as in our first skin conductance dataset.In such an experiment, the average probability p0 may be used as an approximation for the subject's normal baseline probability of SCR occurrence in daily life [55], [56], [58].This would be equivalent to the corresponding neural impulse occurrence probability as well.We simulated two sets of data where we set p 0 = 0.05 so that β 0 = log[0.05×(1−0.05)−1 ], but where the approximation p0 was first slightly less than and secondly slightly greater than 0.05 in the resulting datasets [58].Fig. 3 shows the state estimation results.The model parameters used to generate the data and their final EM estimates are shown in Table 2.
In general, there is good agreement between the recovered model parameters, the state estimates and their true values.The quantile-quantile (QQ) plots also indicate a Gaussian distribution of the residual errors helping to demonstrate the accuracy of the algorithm.Note also that we have not placed any constraints on the sign of the marks (i.e., positive or negative).As such, the method can be used in more general state estimation applications involving an MPP and a continuous-valued variable.For instance, the algorithm could be used to determine a latent state underlying a set of financial time series observations.Here the observations could consist of point process amplitudes, either negative losses or positive gains, along with another continuous-valued financial variable.

B. EXPERIMENTAL DATA -SKIN CONDUCTANCE
Skin conductance-based sympathetic arousal estimation results for the study data in [63] are shown in Fig. 4  (participant 1) and in Appendix C (all participants).After the EM algorithm converges, we have x k ∼ N (x k|K , σ 2 k|K ) where x k|K and σ 2 k|K are the smoothened mean and variance estimates for x k .The results in the figures depict x k|K and its confidence intervals calculated using σ 2 k|K .We also provide a high arousal index (HAI) which is the probability that x k > x thresh similar to [54] in the figures.We take x thresh as the median state value [54].The HAI expresses the probability that x k is above a certain baseline and is similar to the certainty level of the ideal observer in [17] and [18].
On experimental data, the EM algorithm can result in x k overfitting to s k .Therefore, we apply a two-part modification at the M-step to prevent this from occurring.Firstly, we take θ = {δ 0 , δ 1 , σ 2 w } to be the model parameters governing x k 's fit to s k , and modify the (m + 1) th M-step update for θ to where 0 < λ ≤ 1. Setting λ = 1 is equivalent to performing the full update pred .Instead we choose λ < 1 taking a smaller step in the direction of the next value.Secondly, we initialize σ 2 w,init at a larger value and allow θ to only update until σ 2 w reaches a threshold σ 2 w,stop (σ 2 w is driven very low when overfitting begins to occur).This is a form of the early-stopping used in [55] but with a λ term now included.The hyper-parameters σ 2 w,init , λ and σ 2 w,stop can all be fine-tuned to achieve a desired level of overfitting control and an entire range of fits to s k .Additional results illustrating this are provided in Appendix D. Overfitting can also be controlled by completely halting the update of all the parameters when σ 2 w reaches σ 2 w,stop .In general, for most of the participants, the arousal estimates are high during the cognitive stressors and gradually diminish thereafter during relaxation.A slight increase is to be noted around the time of the emotional stressor.Minor variations between participants, with the exception of participant 4, can be seen.For participant 4, there is an unusually large increase in the tonic component during the relaxation period.This large increase occurs due to two sharp transitions at both the beginning and end of the relaxation period.These transitions may have been due to a baseline shift arising from sensor movement.The corresponding increase in arousal during the relaxation period for participant 4 may therefore have been due to artifact contamination.

C. EXPERIMENTAL DATA -CORTISOL
Cortisol-related energy state estimation results are shown in Fig. 5 (CFS patient 1) and Appendix E (all subject pairs) for a selected set of σ 2 w,init , λ and σ 2 w,stop values.
137212 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Adjustments to these three hyper-parameters can provide varying degrees of fits to the data (Appendix D).The mean p k and x k estimates for all patients and controls are also shown in Appendix F. Now cortisol secretion is circadian [74].Serum cortisol levels typically start to increase early morning.Peak concentration values then tend to be reached during the 30-45 min period subsequent to awakening.The concentrations gradually decline thereafter and usually reach nadir values at about midnight during sleep [75], [76].This general pattern is visible for individual patients and controls.Care must be taken, however, when interpreting the group means.While blood cortisol concentrations generally lie between 0-25 µ g/dl for all subjects, the raw x k estimates can lie in slightly different ranges.A further discussion on cortisol differences in CFS and FMS is provided in the following section.

D. EXPERIMENTAL DATA -OTHER
While the EM algorithm provides reasonable sympathetic arousal and energy state estimates, the absence of a ground truth presents several difficulties.For instance, the question of how the estimates are finally validated is challenging to answer.Similar problems arise in biomedical applications elsewhere as well.For instance, modeling techniques that calculate central aortic blood pressure using only peripheral measurements have also faced validation challenges [77].
A possible option to address this is to cross-check/validate the estimates using another physiological signal(s)/feature(s).We illustrate this using experimental skin conductance and fNIRS data from a subject in the n-back experiment described earlier.Fig. 6 shows the sympathetic arousal estimates along with the mean and mean envelope of a 30 s sample-by-sample running average of the total hemoglobin (hbT) signal energy.The signal energies are calculated from optodes over the pre-frontal cortex (the fNIRS channel locations are shown in Appendix G).As seen in the figure, both sympathetic arousal and the fNIRS blood flow features are lower during calming music and are higher during vexing music.It is likely that the subject had to put in more cognitive effort while listening to vexing music to correctly perform the n-back tasks.Thus the use of other physiological signals to corroborate the state estimates is a possible alternative.This is especially since a ground truth is unavailable for the types of problems we address here.We finally illustrate the ability of the point process estimator to generalize to other types of physiological data with an additional example.We estimate a quantity akin to sympathetic arousal by using an MPP extracted from heart rate and the tonic skin conductance component (as s k ).The data we use is from a single subject in the experiment in [67].The estimation results are shown in Fig. 7 and illustrate how data from more than one physiological source can be utilized to estimate a latent internal variable.The estimates are high towards the beginning of the session and tend to diminish thereafter.Since this was the subject's first task session [67], it is possible that unfamiliarity with the experimental setup and the assigned work generated some stress.We also made a similar observation in our earlier work in [57] for this particular subject.There we noted a general decrease in stress state estimates over all three task sessions.Here we observe the same with our new estimator as well.

V. COMPARISON WITH PREVIOUS METHODS
We finally provide comparisons between the performance of the current method and our earlier state-space approaches.We begin with cortisol.As noted earlier, we only developed one cortisol-related energy state estimator [32].There we took the observations to consist of the pulses and the upper and lower cortisol envelopes (one binary and two continuous variables).The estimates for CFS patient 1 obtained using this method are shown in Fig. 8.As evident from a comparison between Figs. 5 and 8, none of the subtle, fine-grained energy state variations are visible with the earlier method.These variations are indeed visible with the current method, and serve as an illustration of how it outperforms the previous approach in [32].
We perform a similar comparison for skin conductance.We noted earlier that the purely MPP-based sympathetic arousal estimator in [58] outperformed earlier approaches in [53], [54], and [55].We thus compare our current method with this MPP approach.Since a ground truth is unavailable, the comparison is largely qualitative in nature (the comparisons shown in [58] were similarly qualitative).In [58], we illustrated one aspect of how the MPP method outperformed the others by not overly reacting to small impulses in an example.Namely, at the beginning of the horror movie clip for participant 1, the HAI of the MPP-based estimate was only moderately elevated.This was in keeping with the two small neural impulses seen at the beginning of this period.The HAI was elevated much more according to the estimates from the other approaches in [53], [54], and [55].We show in Appendix H how we can obtain a similar moderate elevation in the HAI with our current method by simply tuning σ 2 w,init , λ and σ 2 w,stop and using early stopping.We further show that similar HAIs can be obtained with different levels of fits to the tonic component s k .In each of these cases, the resulting x k estimates capture more subtle, FIGURE 8. Cortisol-based Energy State Estimation using the Method in [32].The sub-panels respectively depict: (a) the raw blood cortisol data z k ; (b) the deconvolved cortisol pulses (blue); (c), (d) the lower and upper envelopes of the reconstructed blood cortisol profile s k (orange) and the estimated fits (red); (e) the probability of pulse occurrence p k with its 95% confidence limits; (f) the energy production state x k with its 95% confidence limits.
fine-grained variations compared to the solely MPP-based estimate.The overfitting control technique we present here with λ included to govern the rate of fit to s k is also superior to that used in an earlier work in [55].In [55] there was no means to slow this rate.The state-space model we use here also conforms more naturally to the skin conductance and cortisol constituents.Overall, our current method is superior compared to the earlier sympathetic arousal estimators in [53], [54], [55], and [58].Note that a slightly different set of σ 2 w,init , λ and σ 2 w,stop values are used in Appendix H and in Fig. 4 and Appendix C.These are also the values we used for the subject in the n-back experiment and the one in the office work-like experiment to maintain consistency.However, as evident from the estimates, the method offers the capability to obtain an entire range of results based on the choice of these three hyper-parameters.
Recall that the supervised learning methods we described earlier all required labels.We do not have such labels in our multi-subject datasets.Hence, in order to compare our current estimator with machine learning-related methods in general, we use our recent neural network-based estimation method in [60] for comparison.Note that this represents a non-traditional use of machine learning as labels in the typical sense are not required.Rather, training relies on the maximization of a probability term relating observations to the latent state.Here, neural networks learn the state-space model and the estimator (a recurrent neural network is used for the estimator).Details of network architectures, sizes, hyper-parameters, weighted loss functions, training epochs, overfitting control, out-of-sample validation etc. can be found in [60].This method too utilized the same MPP and continuous-valued variables as observations that we use here.We provide the estimation results for CFS patient 1 and FMS patient 1 in Fig. 9 and Appendix H. Complete results are available in [60].As evident from the results, tiny noise fluctuations remain in the neural network estimates for x k and the fits to the observed data.A simple visual comparison of Fig. 5 (where estimates are from our current EM method) and Fig. 9 illustrates this.The neural networks have thousands of network weights to learn from only a little data.This likely contributes to the noise in the final estimates.In contrast, our EM-based method has only a handful of parameters to determine and the estimates are much smoother.The neural networks are also more difficult to control.As we pointed out in [60], the neural networks occasionally provide inverted estimates for x k and p k .This is implausible physiologically and would imply, for instance, that energy is low during the day and high at night for cortisol.There is also the challenge of interpretability.Since neural networks learn the state-space model, it is difficult to determine how exactly the observations are related to the state variable.Likewise, it is difficult to interpret how x k evolves over time as well since a neural network governs the equations.Finally, although minor, the neural networks take a considerably longer time to train.However, the time difference between the neural network method and the current EM method is not as large during inference.While it is difficult to perform a quantitative comparison since x k has no ground truth, nevertheless, our proposed EM method does have certain considerable advantages over the method in [60].
In summary, while we have developed several point process-based state estimators in the past, our current method has several advantages over them.As just described based on the results above, the current estimator offers the ability to obtain more fine-grained state estimates and has better overfitting control.It is also a new filter algorithm which conforms more naturally to the input data (i.e., x k is estimated from an MPP and a continuous-valued variable).Moreover, the resulting estimates do not have the tiny fluctuations present in [60], and the method also does not suffer from a difficulty in controlling as do the neural networks.

VI. DISCUSSION
The electrical and chemical changes caused by unobserved internal state variables within the human body provide an opportunity for estimating them.Here we develop a state-space method for estimating sympathetic arousal and energy production.The method is based on skin conductance and cortisol measurements respectively where the observations form an MPP and a continuous-valued variable.

A. SKIN CONDUCTANCE
With experimental skin conductance data, arousal was generally high during the cognitive stressors and only brief increases were seen during emotional stress.The cognitive tasks require concentration on the part of the subjects.In contrast, the horror movie clip only requires watching a screen.This is likely the reason for obtaining the arousal estimates that we did.The only exception was participant 4 for whom noise contamination could have been a major factor.

B. CORTISOL
For cortisol, the energy estimates are largely consistent with known physiology.As stated earlier, cortisol levels are known to rise towards morning as a part of the body's cortisol awakening response (CAR).Overall, energy estimates are consistent with the CAR across subjects.The nighttime drop in energy is also consistent with the known physiology of cortisol.We also investigated group-level differences.The cortisol data come from patients with CFS, FMS or both, and matched controls.CFS is a condition with unknown etiology.Patients suffering from it characteristically experience debilitating fatigue [78].CFS is also accompanied by symptoms such as difficulty concentrating, disturbed sleep, pain and post-exertional malaise [78].There has been a lack of consensus in studies investigating cortisol differences between CFS patients and healthy controls.A review of the data prior to 2003 can be found in [79], and a more recent review of the studies after 2003 in [80].Single sample studies, which have been criticized as being limited, have reported no change, increases and decreases between CFS patients and healthy controls [79], [80].Two of the most important serial sample studies in [65] and [81] found no significant cortisol differences between patients and controls during the day [80].Here, blood samples were drawn every 10-15 min over a 24 h period.While hypocortisolism may exist in at least some CFS patients, it is likely that no single cause for CFS exists [80].Moreover, it is also unclear whether changes in hypothalamic-pituitary-adrenal (HPA) axis function in CFS are a cause of it or are merely a consequence.FMS is similar to CFS with patients often experiencing similar symptoms [82].Again, studies investigating cortisol differences between FMS patients and healthy controls have yielded mixed results [83], [84].The serial sample study in [65] also found no significant differences in mean cortisol levels between patients (CFS, FMS, and FMS and CFS) and matched controls.We use the data from this same study here.The mean p k and x k estimates (shown in Appendix F) appear similar possibly due to this reason.

C. CHALLENGES
In addition to the validation challenge noted earlier, the absence of a ground truth also creates a difficulty for feature selection.In the case of skin conductance, our choice of features largely relies on the medical/psychophysiology literature.SCR rates, their amplitudes and the tonic levels are commonly used indicators of sympathetic arousal in the literature [62].These form an MPP and a continuousvalued variable; hence the reason for including them in our state-space model.Feature selection and validation is more challenging for cortisol.Firstly, there is a lack of serial sample cortisol studies investigating what features correlate well with energy and fatigue.Moreover, a satisfactory, person-specific solution to determining the number of secretory events, and their timings and amplitudes from a series of blood cortisol measurements is comparatively recent [16].Here we take the pulsatile cortisol profile to estimate x k owing to its inherent secretory nature.We also choose to incorporate the continuous-valued cortisol concentration as a feature.This is because unbound blood cortisol is active biologically and has a contribution towards energy production [59].The final state estimates are generally consistent with the known physiology of cortisol (e.g. with the CAR).However, we do not possess additional subject data such as responses regarding how they felt (e.g.fatigue level) throughout the experiment that could be used for cross-checking/validation. Feature selection and validation thus remain a challenge with cortisol and we acknowledge the difficulty as a limitation of our study.
137216 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
We also provide some additional results with held-out data in Appendix I to illustrate the method's performance on data hitherto unseen.

D. STUDY LIMITATIONS
There are a few limitations of our study.As noted earlier, a ground truth for x k is absent.Hence, there is first a difficulty validating the continuous-valued state estimates.Selecting the best set of features with which to perform estimation can also be challenging.Consequently, we rely on physiological reasons for feature selection rather than numerical comparisons.Ideally, the best features for estimating sympathetic arousal or cortisol-related energy could have be selected in such a manner if a ground truth were available.Both of these challenges were discussed in sub-sections IV-D and VI-C above.Here we point out two additional, though related, limitations.Firstly, since a ground truth for x k unavailable, it is difficult to perform a quantitative comparison during evaluation.Rather we largely have to rely on qualitative methods.Our comparisons were mostly qualitative in nature in our earlier work in [55], [56], [58], and [60] as well.Nevertheless, this largely qualitative nature of the evaluations does remain an additional limitation in our current study also.
A final limitation is the complexity of the EM algorithm, and in particular, the growth in complexity with each new variable/observation added.This is largely due to the nature of the point process-based state estimators themselves.In our current work, we add just a single continuous-valued variable to our earlier model in [58].The addition of this one variable requires the estimation of δ 0 , δ 1 and σ 2 w .All these add to the number of equations that need to be solved at the M-step, and also add to the complexity of the E-step equations for x k|k and σ 2 k|k .Adding a few more variables to estimate x k would result in the complexity growing even further.This too is a limitation of our method as it is a hindrance to the inclusion of more observations to the state-space model.

E. MODEL PARAMETER SELECTION
There are two types of parameters in our proposed model.Firstly, there are the state-space model parameters that relate x k to the MPP and continuous-valued observations.These include γ 0 , γ 1 , σ 2 v etc.These state-space model parameters get selected during the EM algorithm iterations with the objective of maximizing the likelihood in (8).However, there are three additional hyper-parameters σ 2 w,init , λ and σ 2 w,stop in our model that need to be chosen as well.These help control the degree of overfitting to s k .We propose two methods whereby they can be chosen.Firstly, a human subjects experiment could be performed where skin conductance or cortisol levels are recorded along with some form of subject feedback of emotional arousal or energy levels respectively.The feedback would likely need to be on a rating scale and be recorded periodically.Next, state estimation could be performed using our proposed model.Since ratings are now available, the three hyper-parameters σ 2 w,init , λ and σ 2 w,stop could be tuned so that x k best matches them (e.g. based on a least squares criterion).This would provide one method of selecting the hyper-parameters.A second approach would be to embed the proposed state estimator within a larger control loop for regulating arousal or energy [85], [86], [87].An objective criterion could be set to maintain a patient's arousal levels or energy within a pre-defined range.The hyper-parameters could then be selected such that x k is best maintained within the required range.These two methods could be used to tune σ 2 w,init , λ and σ 2 w,stop .Note however, that these two options would require a separate human subjects experiment and the development of a control algorithm similar to [85], [86], [87] respectively.
EM algorithms can be sensitive to parameter initialization.In our case, we randomly initialize the model parameters.Point process-based state estimators such as the ones we have developed can have a tendency to be sensitive to the initialization of σ 2 ε .In some of our earlier testing work, we ran one of our state estimators for different parts of the data to obtain two possible extremes for σ 2 ε , i.e., we ran the estimator separately on the stress and relaxation portions of the skin conductance data in [63], and on the sleep and wake portions of the cortisol data in [65].After obtaining these separate noise variance estimates, we then re-initialized at those values and re-ran the EM algorithm.In general, the final values at which convergence occurred were close to those obtained by having the random initialization.Visual comparisons also showed negligible differences.A similar procedure could be used here as well to mitigate the sensitivity to the initialization of σ 2 ε .Other types of initialization strategies could also be utilized and the final estimates calculated based on the aggregated results.

F. EVALUATION METHODS
As evident from the results described above, we have largely used four methods for evaluating our point process Bayesian state estimator.These evaluation methods, along with reasons from their choice, are described below.

1) EVALUATION WITH SIMULATED DATA
We derive a new algorithm in our work-one which is capable of recovering a latent state x k giving rise to MPP-type and continuous-valued observations.As such, it is necessary to first verify that the algorithm itself is able to correctly recover x k adequately and is not faulty.Thus we first evaluate our method on simulated data, where the ground truth is known, to ensure that x k and the unknown model parameters can be recovered.This form of evaluation is identical to our earlier work in [55] and [58].The evaluation is twofold after the EM algorithm converges.Firstly, the recovered state-space model parameters need to be sufficiently close to the original ones and the estimated x k needs to closely follow the original x k as well.These results on simulated data are shown in Table 2 and Fig. 3 and demonstrate the accuracy of our algorithm.Moreover, the QQ-plots of the residuals also indicate that they are sufficiently close to a Gaussian distribution.A non-Gaussian distribution may indicate that the method is unable to correctly recover x k and is inadequate.

2) CORTISOL -QUALITATIVE EVALUATION BASED ON MEDICAL PHYSIOLOGY
We performed evaluations on experimental data next.Here it is necessary to verify that the estimates are indeed physiologically plausible and agree with medical knowledge.Disagreement here may indicate that the model may be insufficient to capture/explain the physiological reality we are attempting to investigate.Also note that, as pointed out earlier, ground truths for the subjects' internal physiological states are unavailable and we consequently need to rely on qualitative comparisons.Therefore, here too, we follow [60] for evaluations on cortisol data and rely on an agreement between the energy state estimates and known physiology.Cortisol secretion is circadian, and is known to exhibit an awakening response pattern and daytime vs. nighttime patterns [74], [75], [76].The estimates we obtain for x k are all consistent with these known variations.For instance, nighttime decreases and increases during morning awakening are all generally visible across subjects in the dataset.This form of evaluation helps corroborate the viability of the model.

3) SKIN CONDUCTANCE -QUALITATIVE EVALUATION BASED ON EXPERIMENTAL CONDITIONS AND MEDICAL PHYSIOLOGY
Evaluation of skin conductance-based sympathetic arousal estimates is similar.Here too, results from the EM algorithm need to agree with physiological expectations.Contradictory estimates may indicate a problem with the model itself or the estimator.Here again, we perform qualitative evaluations following [54], [55], [58], [60].Now stress activates the sympathetic nervous system.As such, higher levels of arousal are to be expected during the cognitive stressors and lower levels during relaxation for the experiment in [63].Our sympathetic arousal estimates are consistent with this expected physiology.In the experiment in [63], however, the emotional stressor may not have been sufficiently stressful to generate a response.Furthermore, we do point out that the overall results we obtain are generally similar to those obtained with our earlier methods as well.

4) SKIN CONDUCTANCE -QUALITATIVE COMPARISON WITH FNIRS-BASED BLOOD FLOW OVER THE PRE-FRONTAL CORTEX
We finally performed an additional evaluation with fNIRS data.fNIRS is a neuroimaging technique through which blood flow in the brain can be estimated.The method relies on the different levels of optical absorption by oxygenated and deoxygenated hemoglobin to perform measurements.Greater blood flow is correlated with increased neural activation in the cortex.Consequently, fNIRS permits the study of different cortical regions associated with different tasks [88].The n-back task in our experiment involves working memory and requires cognitive effort on the part of the subjects.fNIRS studies in n-back experiments have shown the importance of the pre-frontal cortex and have noted blood flow changes in this region associated with the mental workloads imposed thereof [89].Heavier loads (e.g.3-back task instead of 1-back), are also associated with increased cognitive demands, and blood flow changes have been observed accordingly [90].Furthermore, fNIRS studies have shown that increasing mental workloads are associated with higher hbT concentrations [91].
Thus, given the importance of the pre-frontal cortex in nback tasks and the increases in hbT associated with greater workloads, we chose to investigate an hbT signal energy feature here as well.As evident from Fig. 6, this hbT signal energy over the pre-frontal cortex is higher during vexing music than during calming music.This is consistent with expectations based on the literature.The subject in our experiment likely had to concentrate more in order to identify the correct answers in the presence of the vexing music disturbance.While a more comprehensive experiment further investigating hbT signal energy in a larger group of subjects could be performed, the results we obtained in Fig. 6 nevertheless agree with expectations.Thus, despite the absence of a ground truth, this fNIRS comparison helps corroborate the stress-related sympathetic arousal estimates we obtained from skin conductance.

G. COMPARATIVE ADVANTAGES
Despite the limitations of point process-based Bayesian estimators discussed above, they nevertheless do have certain advantages over traditional machine learning methods.While these were described in Section II above, we briefly review them here.Supervised learning methods require labels.These labels can often have variability between subjects, and even contain variations within the same subject over time.In some instances, labels are also provided by external individuals.Accurately assessing an internal physiological quantity such as cortisol-related energy or sympathetic arousal of a person by another external individual can, however, be challenging.Moreover, a number of machine learning methods and their associated models are difficult to interpret, and occasionally require large amounts of data and lengthy training times.Since our Bayesian estimators do not rely on labels, they suffer from none of these associated drawbacks.Moreover, the algorithms can even be run on a single data record.As such they do not require large amounts of data nor require long training times.

H. REAL-TIME ESTIMATION
Note that our method incorporates both a forward filter and a backward smoother.Consequently, real-time estimates of arousal or energy cannot be provided due to the backward pass through the data.We therefore suggest two approaches for the use of the current method in the real-world.Firstly, 137218 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
if real-time estimates are required, the outputs of the forward filter x k|k alone can be used instead of x k|K .If the real-time requirement is not present, a second approach would be to window the data.Smoothened estimates of x k|K could then be provided at the end of each window.With either of these two methods, however, occasional retraining will be necessary to ensure that model parameters are kept up to date.

I. EXTERNAL EFFECTS, ALTERNATE METHODS AND FUTURE WORK
External influences such as body temperature and hydration can affect skin conductance [54].However, in a controlled experimental setting such as an air conditioned room, it may be expected that these factors remain reasonably constant.Hence, the variations in sweat secretions affecting the skin's conductivity largely depend on the subjects' psychology and emotion.Future work would involve developing a more advanced state estimator for arousal where external factor-based corrections could be applied to the features.As in our previous work [55], external inputs may also be considered to influence x k based on a suitable modeling of them.The external inputs could also be considered for biochemical signals as well.
Pattern recognition methods have also been used to estimate physiological quantities (e.g.emotion) and we have used them in our previous work as well [92], [93].Unfortunately, these methods often depend on subject-provided labels and thus suffer from inter-subject variability.Moreover, the classifiers themselves also suffer from a difficulty with regard to interpretability.Our focus here is to develop a physiologically-explainable state-space model for the estimation of unobserved quantities within the human body.We have applied these to two major areas related to bioelectric and biochemical phenomena.

VII. CONCLUSION
Changes in internal states within the human body often give rise to point process-like bioelectric and biochemical phenomena that can be measured.Here we have developed methods for estimating sympathetic arousal and energy production from skin conductance and blood cortisol measurements respectively.The observations are in the form of an MPP accompanied by a continuous-valued variable.High levels of arousal were typically detected during stressful periods in experimental data as were low levels during relaxation.General energy estimates based on cortisol data also agreed with known circadian physiology.Furthermore, results on additional data helped corroborate the estimates and demonstrate the algorithm's wider applicability.A comparison with estimates from previous approaches also indicate a superiority of the current method.Future work would involve the simultaneous estimation of states linked to bioelectric and biochemical measurements.This would enable obtaining a more comprehensive view of the body as well as developing the closed-loop control framework for therapy delivery.

APPENDIX A DERIVATION OF E-STEP FILTER EQUATIONS
As noted in the main text, we follow the general approach in [19] for deriving the E-step filter update equations.In a recent work [58], we derived a similar EM algorithm to determine a latent state x k underlying a set of MPP observations (n k , r k ).
Here we extend the work to include another continuous variable s k .In the case of [58], we pointed out that the E-step filter update equations switched between those given in [17] and [18] based on an if-else condition depending on the value of n k .The EM algorithm in [17] estimated x k solely based on a sequence of binary observations n k .Meanwhile the algorithm in [18] did so based on a sequence of binary and continuous observations n k and r k .As shown below, when another continuous variable is added, the filter update equations switch between those in [18] and those in our previous work in [55].Here, x k was estimated based on a binary observation n k and two continuous observations r k and s k .

A. CASE WHEN POINT PROCESS EVENT IS ABSENT
First consider the case when n k = 0 and p Therefore, We make a Gaussian approximation to the posterior p(x k |y 1:k ) to derive the mean and variance [19].First however, note that We first take the partial derivative of q 1 with respect to x k and set it to 0 to solve for the mean. Therefore, Since n k = 0, we can write −β Solving for x k yields, which is the filter equation for the state update in [18].Taking the second derivative yields, This yields the variance update in [18]

B. CASE WHEN POINT PROCESS EVENT IS PRESENT
Next consider the case when n k = 1 and This yields the posterior Using q 2 to denote the exponent, the first partial derivative yields, Since n k = 1, we may replace (1 − p k ) with (n k − p k ).Setting the partial derivative to 0 and solving for x k yields which is the state update equation in [55].Similarly, the second derivative yields Therefore, the variance update is [55]

APPENDIX B DERIVATION OF M-STEP UPDATE EQUATIONS
The Bernoulli probability term in p(x 1:K , y 1:K | ) can be rewritten as and therefore the expected log-likelihood becomes We shall consider the derivations of the model parameters separately.When doing so, similar to [17], we take A. DERIVATION OF M-STEP UPDATE EQUATIONS FOR γ 0 , γ 1 , δ 0 AND δ 1 Let Q 1 denote the term in Q containing γ 0 and γ 1 .
137220 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Setting the partial derivatives of Q 1 with respect to γ 0 and γ 1 to 0 yield, The solutions to these two simultaneous equations provide γ 0 and γ 1 .
δ 0 and δ 1 may be obtained similarly by using s k instead of r k .In this case, the summation would also change to k = 1, 2, . . ., K (instead of just being over the subset K ).
We take the partial derivative of Q 2 with respect to σ 2 v and set it to 0.        (e) the energy production state x k with its 95% confidence limits.
137228 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
The M-step update for σ 2 w may be obtained similarly by summing over k = 1, 2, . . ., K using the corresponding s k , δ 0 and δ 1 terms.

C. DERIVATION OF M-STEP UPDATE EQUATION FOR
While it is possible to treat x 0 as a separate model parameter, here we set x 0 = x 1 following an option in [17] and [18] to permit some bias at the beginning.Therefore, We take the partial derivative with respect to σ 2 ε and set it to 0.
D. DERIVATION OF M-STEP UPDATE EQUATIONS FOR β 0 AND β 1 Let Q 4 denote the expectation term containing β 0 and β 1 .
We first perform a Taylor expansion of the log term around x k|K [19].
Taking the expected value on both sides, we have = log 1 + e β 0 +β 1 x k|K + 0 Therefore, Now, And similarly, 137240 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Taking the partial derivative with respect to β 0 , we have These two simultaneous equations for the partial derivatives of Q 4 taken with respect to β 0 and β 1 may be solved numerically to obtain the corresponding M-step updates.

APPENDIX C SKIN CONDUCTANCE -STATE ESTIMATION RESULTS
After the model parameters converge, we have for our state estimates x k ∼ N (x k|K , σ 2 k|K ).Since we can use the change-of-variables formula to determine the probability density function for p k as This can be used to calculate the confidence limits for p k .Skin conductance-based sympathetic arousal estimation

APPENDIX D REDUCTION OF OVERFITTING
The degree of x k 's overfitting to s k can be controlled by changing the initial guess for the sensor noise variance σ 2 w,init , λ which governs the step size and the threshold σ 2 w,stop at which to halt updating θ.Fig. 11 shows some examples where we have varied these three values and the different levels of overfitting control that can thus be achieved.In general, a larger σ 2 w,init initialization, a smaller λ and an earlier stopping threshold σ 2 w,stop provide more stringent overfitting control.
As noted in the main text, the absence of a ground truth poses several difficulties when performing estimation with experimental data.For instance, the question could be posed as to how much x k should be permitted to fit to s k .Since a ground truth is unavailable, selecting the degree of x k 's fit to s k presents a challenge.Our emphasis here is to develop the methods that permit x k 's potential to overfit to the continuous-valued s k to be controlled as desired.Tuning σ 2 w,init , λ and σ 2 w,stop allows us to do so.A final choice of these hyper-parameters could be made based on the application.
137242 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

APPENDIX E CORTISOL -STATE ESTIMATION RESULTS
Energy state estimates for all 31 subject pairs (patients and matched controls) are shown in Fig. 12. Again, although overfitting control could be customized based on each participant, we have chosen σ 2 w,init = v s , λ = 0.01 and σ 2 w,stop = 0.75v s for uniformity.A further point is to be noted when using our method on cortisol data.In certain cases, β 1 and γ 1 can converge to values that are negative.However, this would be physiologically implausible.For instance, a negative γ 1 would imply that energy reduces with larger pulse secretions, and a model with a negative β 1 would mean that energy again reduces when there is a greater tendency for pulses to be secreted.In such cases, we choose to halt the EM iterations prematurely at the step just before γ 1 or β 1 goes into a physiologically implausible region.A second issue arises if this transition into a physiologically implausible region happens very early on (e.g. at the first or second iteration).In such cases, we suggest changing the initialization of γ 0,init (we usually initialize this at the median value of r k in the data).Changing γ 0,init is also an option if the final estimates yield a largely negative estimate for r k .These are two points worth noting when performing estimation on cortisol data.variations in x k can also be captured with the current method.Moreover, different levels of fits to s k can be achieved as well due to our overfitting control technique and serve to illustrate the methodological superiority of the current approach.
Cortisol-based energy states estimates for two subjects (CFS patient 1 and FMS patient 1) using our earlier neural network method in [60] are shown in the main text in Fig. 9 and here in Fig. 17.As stated in the main text, these estimates are noisier compared to those from our current EM method.Moreover, other drawbacks are also present.These include the difficult interpretability, the longer training times and a difficulty in controlling the output of the neural networks.

FIGURE 1 .
FIGURE 1. Deconvolved Skin Conductance and Blood Cortisol Data.Skin conductance signals and blood cortisol profiles have an inherent ''spikey'' appearance owing to the mechanisms responsible for their generation.The underlying pulsatile/impulse events can be extracted via deconvolution.The upper sub-panel depicts a skin conductance signal (blue), its tonic component (black) and the neural impulses underlying phasic variations (red).The lower sub-panel depicts blood cortisol measurements taken at 10 min intervals (blue dots), the reconstructed minute-wise blood cortisol concentrations (black) and the pulsatile secretions (red).

FIGURE 3 .
FIGURE 3. State Estimation with Simulated Data [Smaller Approximation ( p0 < 0.05) and Larger Approximation ( p0 > 0.05)].Within each sub-figure, the sub-panels respectively depict: (a) the MPP observations (blue) and the estimated fit to r k (red); (b) the point process event occurrence probabilities p k (blue) and their estimate (red); (c) the continuous-valued variable s k (blue) and its estimate (red); (d) the state x k (blue) and its estimate (red); (e) the QQ plot for the residual error of x k .

FIGURE 4 .
FIGURE 4. Skin Conductance-based Sympathetic Arousal Estimation.The sub-panels respectively depict: (a) the skin conductance signal; (b) the deconvolved neural impulses (blue) and the estimated fit to r k (red); (c) the tonic component s k (orange) and its estimated fit (red); (d) the probability of impulse occurrence p k with its 95% confidence limits;(e) the sympathetic arousal state x k with its 95% confidence limits; (f) the HAI (the regions above 90% and below 10% are shaded in red and green respectively).The shaded backgrounds correspond to the period of instruction (red), the backward counting task (green), the Stroop color-word test (orange), relaxation (blue) and emotional stress (yellow).The cognitive stress portion consisted of the counting task and the Stroop test.

FIGURE 5 .
FIGURE 5. Cortisol-based Energy State Estimation.The sub-panels respectively depict: (a) the raw blood cortisol data; (b) the deconvolved cortisol pulses (blue) and the estimated fit to r k (red); (c) the reconstructed blood cortisol profile s k (orange) and its estimated fit (red); (d) the probability of pulse occurrence p k with its 95% confidence limits; (e) the energy production state x k with its 95% confidence limits.

FIGURE 6 .
FIGURE 6. Skin Conductance-based Sympathetic Arousal Estimation with fNIRS Validation/Cross-checking. The sub-panels respectively depict: (a) the skin conductance signal; (b) the deconvolved neural impulses (blue) and the estimated fit to r k (red); (c) the tonic component s k (orange) and its estimated fit (red); (d) the probability of impulse occurrence p k with its 95% confidence limits; (e) the sympathetic arousal state x k with its 95% confidence limits; (f) the HAI calculated based on the median state (the regions above 90% and below 10% are shaded in red and green respectively); (g) the mean (blue) and mean envelope (red) of the hbT energies.The background colors correspond to calming music (yellow), relaxation (green) and vexing music (orange).

FIGURE 7 .
FIGURE 7. Heart Rate and Skin Conductance-based Estimation.The sub-panels respectively depict: (a) the skin conductance signal; (b) the MPP extracted from the LF to HF peaks (blue) and the estimated fit to r k (red); (c) the tonic component s k (orange) and its estimated fit (red); (d) the probability of point process event occurrence p k with its 95% confidence limits; (e) the estimated state x k with its 95% confidence limits; (f) the HAI (the regions above 90% and below 10% are shaded in red and green respectively).

FIGURE 9 .
FIGURE 9.Cortisol-based Energy State Estimation using the Machine Learning Method in[60].The sub-panels respectively depict: (a) the raw blood cortisol data z k ; (b) the deconvolved cortisol pulses (purple) and the estimated fit to r k (blue); (c) the reconstructed blood cortisol profile s k (green) and its estimated fit (orange); (d) the probability of pulse occurrence p k ; (e) the energy production state x k (maroon) with its 95% confidence limits and the circadian term fit to x k (red).

FIGURE 10 .
FIGURE 10.Skin Conductance-based Sympathetic Arousal Estimation.Within each sub-figure, the sub-panels respectively depict: (a) the skin conductance signal z k ; (b) the deconvolved neural impulses (blue) and the estimated fit to r k (red); (c) the tonic component s k (orange) and its estimated fit (red); (d) the probability of impulse occurrence p k with its 95% confidence limits; (e) the sympathetic arousal state x k with its 95% confidence limits; (f) the HAI (the regions above 90% and below 10% are shaded in red and green respectively).The shaded backgrounds correspond to the period of instruction (red), the backward counting task (green), the Stroop color-word test (orange), relaxation (blue) and emotional stress (yellow).

FIGURE 10 .
FIGURE 10. (Continued.)Skin Conductance-based Sympathetic Arousal Estimation.Within each sub-figure, the sub-panels respectively depict: (a) the skin conductance signal z k ; (b) the deconvolved neural impulses (blue) and the estimated fit to r k (red); (c) the tonic component s k (orange) and its estimated fit (red); (d) the probability of impulse occurrence p k with its 95% confidence limits; (e) the sympathetic arousal state x k with its 95% confidence limits; (f) the HAI (the regions above 90% and below 10% are shaded in red and green respectively).The shaded backgrounds correspond to the period of instruction (red), the backward counting task (green), the Stroop color-word test (orange), relaxation (blue) and emotional stress (yellow).

FIGURE 11 .
FIGURE 11.Controlling the Level of Overfitting.Overfitting control can be applied by changing the initial guess for the sensor noise variance σ 2 w ,init , the step size parameter λ and the threshold σ 2 w ,stop .The sub-panels respectively depict: (a) the raw blood cortisol data z k ; (b) the deconvolved cortisol pulses (blue) and the estimated fit to r k (red); (c) the reconstructed blood cortisol profile s k (orange) and its estimated fit (red); (d) the probability of pulse occurrence p k with its 95% confidence limits; (e) the energy production state x k with its 95% confidence limits (v s denotes the variance of the s k values in the figure titles).

FIGURE 12 .
FIGURE 12. Cortisol-based Energy State Estimation for All Patients and Matched Controls.Within each sub-figure, the sub-panels respectively depict: (a) the raw blood cortisol data z k ; (b) the deconvolved cortisol pulses (blue) and the estimated fit to r k (red); (c) the reconstructed blood cortisol profile s k (orange) and its estimated fit (red); (d) the probability of pulse occurrence p k with its 95% confidence limits; (e) the energy production state x k with its 95% confidence limits.

FIGURE 12 .
FIGURE 12. (Continued.)Cortisol-based Energy State Estimation for All Patients and Matched Controls.Within each sub-figure, the sub-panels respectively depict: (a) the raw blood cortisol data z k ; (b) the deconvolved cortisol pulses (blue) and the estimated fit to r k (red); (c) the reconstructed blood cortisol profile s k (orange) and its estimated fit (red); (d) the probability of pulse occurrence p k with its 95% confidence limits; (e) the energy production state x k with its 95% confidence limits.

FIGURE 12 .
FIGURE 12. (Continued.)Cortisol-based Energy State Estimation for All Patients and Matched Controls.Within each sub-figure, the sub-panels depict: (a) the raw blood cortisol data z k ; (b) the deconvolved cortisol pulses (blue) and the estimated fit to r k (red); (c) the reconstructed blood cortisol profile s k (orange) and its estimated fit (red); (d) the probability of pulse occurrence p k with its 95% confidence limits; (e) the energy production state x k with its 95% confidence limits.

FIGURE 12 .
FIGURE 12. (Continued.)Cortisol-based Energy State Estimation for All Patients and Matched Controls.Within each sub-figure, the sub-panels respectively depict: (a) the raw blood cortisol data z k ; (b) the deconvolved cortisol pulses (blue) and the estimated fit to r k (red); (c) the reconstructed blood cortisol profile s k (orange) and its estimated fit (red); (d) the probability of pulse occurrence p k with its 95% confidence limits;(e) the energy production state x k with its 95% confidence limits.

FIGURE 12 .
FIGURE 12. (Continued.)Cortisol-based Energy State Estimation for All Patients and Matched Controls.Within each sub-figure, the sub-panels depict: (a) the raw blood cortisol data z k ; (b) the deconvolved cortisol pulses (blue) and the estimated fit to r k (red); (c) the reconstructed blood cortisol profile s k (orange) and its estimated fit (red); (d) the probability of pulse occurrence p k with its 95% confidence limits;(e) the energy production state x k with its 95% confidence limits.

FIGURE 12 .
FIGURE 12. (Continued.)Cortisol-based Energy State Estimation for All Patients and Matched Controls.Within each sub-figure, the sub-panels respectively depict: (a) the raw blood cortisol data z k ; (b) the deconvolved cortisol pulses (blue) and the estimated fit to r k (red); (c) the reconstructed blood cortisol profile s k (orange) and its estimated fit (red); (d) the probability of pulse occurrence p k with its 95% confidence limits; (e) the energy production state x k with its 95% confidence limits.

FIGURE 12 .
FIGURE 12. (Continued.)Cortisol-based Energy State Estimation for All Patients and Matched Controls.Within each sub-figure, the sub-panels respectively depict: (a) the raw blood cortisol data z k ; (b) the deconvolved cortisol pulses (blue) and the estimated fit to r k (red); (c) the reconstructed blood cortisol profile s k (orange) and its estimated fit (red); (d) the probability of pulse occurrence p k with its 95% confidence limits; (e) the energy production state x k with its 95% confidence limits.

FIGURE 12 .
FIGURE 12. (Continued.)Cortisol-based Energy State Estimation for All Patients and Matched Controls.Within each sub-figure, the sub-panels respectively depict: (a) the raw blood cortisol data z k ; (b) the deconvolved cortisol pulses (blue) and the estimated fit to r k (red); (c) the reconstructed blood cortisol profile s k (orange) and its estimated fit (red); (d) the probability of pulse occurrence p k with its 95% confidence limits; (e) the energy production state x k with its 95% confidence limits.

FIGURE 12 .
FIGURE 12. (Continued.)Cortisol-based Energy State Estimation for All Patients and Matched Controls.Within each sub-figure, the sub-panels respectively depict: (a) the raw blood cortisol data z k ; (b) the deconvolved cortisol pulses (blue) and the estimated fit to r k (red); (c) the reconstructed blood cortisol profile s k (orange) and its estimated fit (red); (d) the probability of pulse occurrence p k with its 95% confidence limits; (e) the energy production state x k with its 95% confidence limits.

FIGURE 12 .
FIGURE 12. (Continued.)Cortisol-based Energy State Estimation for All Patients and Matched Controls.Within each sub-figure, the sub-panels respectively depict: (a) the raw blood cortisol data z k ; (b) the deconvolved cortisol pulses (blue) and the estimated fit to r k (red); (c) the reconstructed blood cortisol profile s k (orange) and its estimated fit (red); (d) the probability of pulse occurrence p k with its 95% confidence limits; (e) the energy production state x k with its 95% confidence limits.

FIGURE 12 .
FIGURE 12. (Continued.)Cortisol-based Energy State Estimation for All Patients and Matched Controls.Within each sub-figure, the sub-panels respectively depict: (a) the raw blood cortisol data z k ; (b) the deconvolved cortisol pulses (blue) and the estimated fit to r k (red); (c) the reconstructed blood cortisol profile s k (orange) and its estimated fit (red); (d) the probability of pulse occurrence p k with its 95% confidence limits; (e) the energy production state x k with its 95% confidence limits.

FIGURE 12 .
FIGURE 12. (Continued.)Cortisol-based Energy State Estimation for All Patients and Matched Controls.Within each sub-figure, the sub-panels respectively depict: (a) the raw blood cortisol data z k ; (b) the deconvolved cortisol pulses (blue) and the estimated fit to r k (red); (c) the reconstructed blood cortisol profile s k (orange) and its estimated fit (red); (d) the probability of pulse occurrence p k with its 95% confidence limits; (e) the energy production state x k with its 95% confidence limits.

FIGURE 12 .
FIGURE 12. (Continued.)Cortisol-based Energy State Estimation for All Patients and Matched Controls.Within each sub-figure, the sub-panels respectively depict: (a) the raw blood cortisol data z k ; (b) the deconvolved cortisol pulses (blue) and the estimated fit to r k (red); (c) the reconstructed blood cortisol profile s k (orange) and its estimated fit (red); (d) the probability of pulse occurrence p k with its 95% confidence limits; (e) the energy production state x k with its 95% confidence limits.

FIGURE 12 .
FIGURE 12. (Continued.)Cortisol-based Energy State Estimation for All Patients and Matched Controls.Within each sub-figure, the sub-panels respectively depict: (a) the raw blood cortisol data z k ; (b) the deconvolved cortisol pulses (blue) and the estimated fit to r k (red); (c) the reconstructed blood cortisol profile s k (orange) and its estimated fit (red); (d) the probability of pulse occurrence p k with its 95% confidence limits; (e) the energy production state x k with its 95% confidence limits.

FIGURE 12 .
FIGURE 12. (Continued.)Cortisol-based Energy State Estimation for All Patients and Matched Controls.Within each sub-figure, the sub-panels respectively depict: (a) the raw blood cortisol data z k ; (b) the deconvolved cortisol pulses (blue) and the estimated fit to r k (red); (c) the reconstructed blood cortisol profile s k (orange) and its estimated fit (red); (d) the probability of pulse occurrence p k with its 95% confidence limits; (e) the energy production state x k with its 95% confidence limits.

FIGURE 12 .
FIGURE 12. (Continued.)Cortisol-based Energy State Estimation for All Patients and Matched Controls.Within each sub-figure, the sub-panels respectively depict: (a) the raw blood cortisol data z k ; (b) the deconvolved cortisol pulses (blue) and the estimated fit to r k (red); (c) the reconstructed blood cortisol profile s k (orange) and its estimated fit (red); (d) the probability of pulse occurrence p k with its 95% confidence limits; (e) the energy production state x k with its 95% confidence limits.

FIGURE 13 .
FIGURE 13.Mean Pulse Probabilities and Energy States for Patients and Controls.Sub-figures (A) and (B) depict the group-wise mean pulse probabilities and state estimates for patients and controls.They correspond to: (A) CFS group; (B) FMS group.Within each sub-figure, the sub-panels respectively depict: (a) the mean probabilities for patients (red) and controls (blue); (b) the mean energy states for patients (red) and controls (blue).

FIGURE 14 .
FIGURE 14. fNIRS Channel Locations in n-back Experiment with Calming and Vexing Music.The red and blue dots in the figure denote the fNIRS sources and detectors respectively.

FIGURE 15 .
FIGURE 15.Comparison of Skin Conductance-based Sympathetic Arousal Estimates using the Method in[58] (Left) and the Current Method (Right).HAI values comparable to that in[58] be obtained with the current method.The sub-panels of the right sub-figure respectively depict: (a) the skin conductance signal z k ; (b) the deconvolved neural impulses (blue) and the estimated fit to r k (red); (c) the tonic component s k (orange) and its estimated fit (red); (d) the probability of impulse occurrence p k with its 95% confidence limits; (e) the sympathetic arousal state x k with its 95% confidence limits; (f) the HAI (the regions above 90% and below 10% are shaded in red and green respectively).The left sub-figure does not have the sub-panel corresponding to s k .However, the remaining sub-panels remain consistent with the right sub-figure.The shaded backgrounds correspond to the period of instruction (red), the backward counting task (green), the Stroop color-word test (orange), relaxation (blue) and emotional stress (yellow).

FIGURE 16 .
FIGURE 16.Skin Conductance-based Sympathetic Arousal Estimation Results Comparable to the Method in [58] obtained with the Current Method but with Different Levels of Overfitting Control.The sub-panels respectively depict: (a) the skin conductance signal z k ; (b) the deconvolved neural impulses (blue) and the estimated fit to r k (red); (c) the tonic component s k (orange) and its estimated fit (red); (d) the probability of impulse occurrence p k with its 95% confidence limits; (e) the sympathetic arousal statex k with its 95% confidence limits; (f) the HAI (the regions above 90% and below 10% are shaded in red and green respectively).The shaded backgrounds correspond to the period of instruction (red), the backward counting task (green), the Stroop color-word test (orange), relaxation (blue) and emotional stress (yellow).

FIGURE 17 .
FIGURE 17. Cortisol-based Energy State Estimation using the MachineLearning Method in[60].The sub-panels respectively depict: (a) the raw blood cortisol data z k ; (b) the deconvolved cortisol pulses (purple) and the estimated fit to r k (blue); (c) the reconstructed blood cortisol profile s k (green) and its estimated fit (orange); (d) the probability of pulse occurrence p k ; (e) the energy production state x k (maroon) with its 95% confidence limits and the circadian term fit to x k (red).resultsare shown in Fig.10.While overfitting control can be tailored to each subject, we have selected σ 2 w,init = 10v s ,

FIGURE 18 .
FIGURE 18. Cortisol-based Energy State Estimation with Held-out Data.The sub-figures depict results when FMS patient 1's model parameters are used to estimate energy for CFS patient 1 (left) and vice versa (right).Within each sub-figure, the sub-panels respectively depict: (a) the raw blood cortisol data z k ; (b) the deconvolved cortisol pulses (blue) and the estimated fit to r k (red); (c) the reconstructed blood cortisol profile s k (orange) and its estimated fit (red); (d) the probability of pulse occurrence p k with its 95% confidence limits; (e) the energy production state x k with its 95% confidence limits.

FIGURE 19 .
FIGURE 19.Cortisol-based Energy State Estimation with Values Evaluated based on the First Half of Data.The sub-figures depict results when model parameters are estimated only using the first half of the data for CFS patient 1 (left) and FMS patient 1 (right).Within each sub-figure, the sub-panels respectively depict: (a) the raw blood cortisol data z k ; (b) the deconvolved cortisol pulses (blue) and the estimated fit to r k (red); (c) the reconstructed blood cortisol profile s k (orange) and its estimated fit (red); (d) the probability of pulse occurrence p k with its 95% confidence limits;(e) the energy production state x k with its 95% confidence limits.

TABLE 2 .
Results of model parameter estimation on simulated data.