To evaluate the proposed model, the BIDMC PPG, and Respiration dataset [37] from PhysioNet [38], containing signals extracted from the MIMIC II matched waveform database [39], was used. The dataset consisted of 53 recordings of PPG and impedance respiratory signals acquired from adult patients aged 19–90 years at the Beth Israel Deaconess Medical Center (Boston, MA, USA). Each recording was made for 8 min and sampled at 125 Hz. The dataset was collected by an analog-to-digital converter (ADC) with 16-bit precision.

In this study, reference RR was obtained from the annotations of the dataset, which were sampled at 1 Hz. To minimize the influence of other components, a bandpass filter with a cutoff frequency between 0.1–0.6 Hz (6–36 breaths per minute) was used in the extraction of respiratory information. Furthermore, to ensure sufficient data for training and validation, we applied a data augmentation strategy. We augmented the data by overlapping the PPG signals at 1 s intervals. For each 16 s PPG segment, we obtained the next PPG segment overlapped by shifting 125 data points. As a result, the number of segments was increased 15 times through the data augmentation.

Various biologically plausible spiking neuron models, such as the Hodgkin–Huxley (HH), Izhikevich, integrate-and-fire (IF), and leaky integrate-and-fire (LIF), have been proposed to transmit information converted into spike trains [40,41,42]. In particular, the IF and LIF neuron models have been employed in numerous studies to leverage the advantages of both biological plausibility and computational efficiency. In this section, we describe the two neuron models used in the proposed network: soft-reset IF and recurrent IF neuron models.

(1) Soft-reset IF neuron model: This neuron model is utilized in the spike encoder to minimize the information loss during spike conversion. Figure 2 shows the differences between the hard and soft-reset mechanisms. In the hard-reset approach, the membrane potential is reset to the reset voltage regardless of whether it exceeds the threshold. However, in the soft-reset approach, the membrane potential is initialized with a specific voltage that exceeds the threshold. The soft-reset IF neuron model is defined as follows:(1)si+1(t)=1,ifui+1(t)−Vth≥00,ifui+1(t)−Vth<0 (2)ui+1(t)=ui+1(t−1)(1−si(t))+(ui+1(t−1)−Vth)si(t), (3)ui+1(t+1)=ui+1(t)+Wi,i+1si(t+1).Equation (1) explains the Heaviside step function for spike activation. si+1(t) denotes the output spike train, ui+1(t) is the membrane potential at the tth time step in the (i+1)th layer, and Vth is a constant threshold voltage set as hyperparameter. If the membrane potential reaches the threshold, an output spike train is fired. Equation (2) expresses the soft-reset mechanism, where the membrane potential ui+1(t) depends on the Equation (1). If the spike is fired, the reset voltage is set to the difference between the previous membrane potential and the threshold; otherwise, it is maintained at the current value. Equation (3) expresses the IF neuron model, where Wi,i+1 denotes the trainable weights between the ith and (i+1)th layers.

(2) Recurrent IF neuron model: The conventional IF neuron model in Equation (3) accumulates the temporal information of the spike trains for sequential processing. However, there is no direct dependency among the time instances. In other words, updating the neuron from the current timestep is not influenced by the information from the previous timesteps. We utilized the feedback-based IF neuron model to incorporate information from the previous timestep and update of the current state of the neuron, whose mathematical model is defined in Equation (4):(4)ui+1(t+1)=ui+1(t)+Wi,i+1si(t+1)+Wrecsi+1(t), where Wrec denotes the recurrent weight of the (i+1)th layer. It leverages the relationship between the time instances by merging the information from the current spike trains in the (i+1)th layer with the previous spike trains in the ith layer.

Spike encoding is a crucial step in the processing of real-valued data using spike trains. The effective conversion of information into spike trains with minimal loss is crucial to the performance of the SNN model. Rate and temporal encoding methods have demonstrated good performance in classification tasks. Nevertheless, these traditional methods have limitations owing to information loss during the encoding process. In this study, we used a combination of direct encoding approaches and convolution to minimize the loss between the model prediction and its ground truth, enabling the direct conversion of PPG segments into spike trains without the need for additional processing steps.

We adopted a single-layer trainable 1D-CNN to encode information specifically related to respiration. The proposed method efficiently encoded only the respiratory information by extracting the temporal features of the PPG signal, which were obtained via trained convolution filters. Spike trains were generated through the accumulation and firing of respiratory information using a soft-reset IF neuron model that enables encoding with minimal loss.

The most critical challenge in training SNNs is the non-differential problem of the spike activation function [43]. Backpropagation learning, which relies on differentiation, is typically used in the standard learning procedure for DNNs. However, the Heaviside step function was used for spike activation, as described in Equation (1). It imposes constraints on backpropagation learning owing to its non-differentiability at the instance when ui+1(t) equals Vth. To address this limitation, a surrogate gradient learning method is proposed, which introduces a differentiable surrogate function that approximates the behavior of the discontinuous Heaviside step function. The surrogate function enables the utilization of optimization techniques based on gradients, thereby facilitating the training process for SNNs, as described in Equation (5) with the first derivative. (5)S(x)=xα1+x2+12, (6)S′(x)=1α(1+x2)32, where x denotes the ui+1(t)−Vth, and α is the scaling parameter to adjust the slope.

Figure 3a illustrates the Heaviside step function and approximate functions corresponding to various parameter values. The parameter α is experimentally chosen as 0.65, which best approximates the Heaviside step function. Figure 3b depicts the gradient functions corresponding to the functions shown in Figure 3a. By replacing the non-differentiable spike activation function with the proposed surrogate function, the network can be effectively trained in an approximated environment.

To overcome the limitations of SNN prediction, a CNN-SNN architecture is proposed that combines convolution operations with the SNN architecture. The proposed network consists of three layers: a spike-encoding layer, spike-hidden layer, and spike-decoding layer; the spike hidden layer was employed iteratively twice. The information regarding the input signal is continuously processed across the T time steps within the network to ensure an accurate interpretation of the input signal. Consequently, the average of the output values over the T time steps is calculated to predict the RR. The proposed network paradigm is illustrated in Figure 4.

Spike Encoding Layer: In traditional encoding methods, the input data are transformed into spike trains before being fed into neural networks. In contrast, this study employs a trainable machine learning-based encoding method that directly utilizes a neural network for spike conversion. The features related to respiration were extracted through convolution operations. Subsequently, the soft-reset IF neurons described in Equations (1)–(3) was employed to generate spike trains.

Spike Hidden Layer: This layer includes two convolution blocks Each convolution block comprises a 1-D convolution, batch normalization, and recurrent IF neurons. The number of hidden layers was experimentally determined since the increase of the number of hidden layers causes a loss of accuracy.

Spike Decoding Layer: This layer is composed of a pooling layer, recurrent IF neurons, and a fully connected layer. A pooling operation was applied to aggregate the information and simultaneously reduce the number of spatial features. The respiratory information converted into spike form from the hidden layer was analyzed to estimate the RR.

To assess the model performance, we adopted Pearson’s Correlation Coefficient (PCC) and Mean Absolute Error (MAE) in breath per minute, depending on the sizes of the PPG segment [19]:(7)PCC=∑(Y−Y¯)(Y′−Y′¯)∑(Y−Y¯)2∑(Y′−Y′¯)2, (8)MAE=1nW∑|Y−Y′|, where Y and Y′ denote the true and estimated RR, respectively. Y¯ and Y′¯ are the averages of the true and estimated RRs, respectively. W is the window size of the PPG segment.

Furthermore, we adopted the following methods to measure the energy efficiency of the proposed SNN and DNN model [44]:(9)E(DNN)=∑ℓ=2LFLℓ∗EMAC, (10)E(SNN)=FLℓ∗EMAC+∑ℓ=2LFLℓ∗EAC, where FLℓ denotes the number of floating-point operations at layer ℓ, EMAC is the energy consumption used for the multiply-–accumulate (multiplication and addition) operations, and EAC is the energy consumption for the accumulated (addition) operations. To count the number of floating point operations for each layer, we utilized the ptflops library. Furthermore, we assumed that the energy consumption for the addition process was 0.1 pJ and the multiplication process was 3.1 pJ, referring to [45].

We divided the training and test data at the subject level. From the BIDMC respiratory dataset, 40 subjects were randomly selected for the training process and 13 of them were used for the test process. Furthermore, a five-fold cross-validation method was applied during the training process. For the benchmark test, a CNN-LSTM model with three convolution layers and one LSTM layer [46], a CNN-RNN model with three convolution layers and one RNN layer [47], and a VGG-8 model were chosen [34]. The detailed parameter settings are presented in Table 1.

Table 2 lists the PCC and MAE performance corresponding to three different window sizes of the PPG segment. Experiments were conducted with 4, 8, and 16 timesteps of the spike trains. The proposed model demonstrated outstanding performance despite an increase in the number of time steps. The optimal performance was achieved when the time step T was set to 8. Generally, the mean firing rate plays a more critical role than the patterns of neuronal firing in SNN studies of classification problems [48,49,50]. Therefore, it was demonstrated that the performance improved with an increase in the timesteps during decision making. However, in the proposed model, longer timesteps of the spike trains did not improve performance.

Table 3 lists the PCC and MAE performances of the proposed model compared with those of other DNN approaches. The proposed model outperformed the CNN-RNN and VGG-8 models with MAE values of 1.37 ± 0.04 and 1.15 ± 0.07 bpm when the window sizes were 16 and 64, respectively. It also yielded better performance than the VGG-8 model, with an MAE of 1.23 ± 0.03 bpm at a window size of 32. The model of Osathitporn et al. [19] exhibited the best performance for the 16 and 32 window sizes with MAE of 1.34 ± 0.01 and 1.11 ± 0.01 bpm. The CNN-LSTM showed the best performance with MAE of 1.11 ± 0.03 bpm in 64 window size. However, the overall results of the proposed model achieved comparable performance with the CNN-LSTM, CNN-RNN and Osathitporn et al. [19] models.

Figure 5a–c display the training and validation loss curves of the proposed model and the other DNN models with window sizes of 16, 32, and 64, respectively. The x-axis represents the training and validation epochs, and the y-axis represents the mean squared error (MSE) losses. The validation losses converge in all cases, and these curves validate the reliability of the results presented in Table 2. In particular, Figure 5c displays the optimal convergence compared to Figure 5a,b.

In addition to the PCC in Table 4, we performed visualization to evaluate the reliability of each estimated RR. Figure 6 shows the Bland–Altman graphs between the estimated RR and ground-truth RR. Figure 6a,b illustrate the results for window sizes of 16 and 32 s using the best performing model, the Osathitporn et al. [19] network as indicated in Table 4. Figure 6c shows the results for a 64-second window size, using the CNN-LSTM network. Figure 6d–f visualize the results of the proposed model. The x-axis represents the average of two measurements, and the y-axis represents difference between the two measurements. For both the best performing models and the proposed model, the majority of the data points contain within the 95% confidence interval, demonstrating the reliability of the models.

Table 4 presents the floating point operations per seconds (FLOPs) and energy costs for the proposed SNN model and other DNN models. DNNs utilize MAC operations as metrics for FLOPs, whereas the proposed model employs synaptic operations [51]. The proposed model showed comparable performance to other DNN models in terms of MAE performance. Furthermore, compared to the best performing model, it demonstrated 18.6, 18.7, and 64.6 times higher energy efficiency, and 1.05, 1.1, and 3.97 times lower floating-point operation counts for the window sizes of 16, 32, and 64, respectively.