I Introduction
The electrocardiogram (ECG) is a quasiperiodic biomedical signal that provides information about cardiac muscle electrical activities. One cardiac cycle in a typical ECG signal is delineated by arrangements of P, the QRS complex, T waves, and PQ and ST segments. Correct Rpeak detection is the first and most critical step in almost all ECG analysis methods. The Rpeak is the highest and only positive peak within the QRS complex, reflecting the ventricular depolarization of the heart’s electrical activity. Precise detection of the Rpeak location plays a critical role in obtaining the morphology of the QRS complex and revealing the location of other ECG fiducial points. Furthermore, Rpeak localization serves as the basis for automated determination of the heart rate, which is a significant criterion for heart arrhythmia diagnoses such as premature atrial contraction, tachycardia, and bradycardia. Many other diseases can also be diagnosed in a noninvasive way using Rpeak detection due to the relationship between heart rate variability and several physiological systems (e.g., vasomotor, respiratory, central nervous, and thermoregulatory).
Various approaches have been proposed in the literature for detecting Rpeaks in an ECG signal [beraza2017comparative]. Typically, these methods consist of two main steps: preprocessing and detection. In the preprocessing step, the algorithm attempts to eliminate the noise and artifacts and to highlight the relevant sections of the ECG [castells2015simple, sharma2017qrs]
. In the second step, various methods are used to locate Rpeaks based on the result of the preprocessing step, and then other waves are detected by defining a set of heuristic rules
[hou2018real]. However, these approaches suffer from some critical drawbacks that limit their performance in practical applications. First, in realtime data processing and ambulatory care settings, where the collected data are highly noisy, preprocessingbased algorithms are less effective. Second, these algorithms can fail to detect Rpeaks in some determinant morphological patterns resulting from certain lifethreatening heart arrhythmias due to the timevarying morphology of the QRS complex. Incorrect detection of Rpeaks can affect the correct identification of subsequent waves.The Rpeak detection step can be generally accomplished either by implementing a thresholdbased technique or by employing an independent threshold technique. The amplitude of the peak and time duration between two consecutive R peaks (i.e., the RR interval) are typically used to determine a suitable threshold [song2015new]. A constant threshold is only efficient for detecting Rpeaks within records with normal morphological patterns. Therefore, recent studies have employed adaptive thresholds, for which there is no need to determine the threshold experimentally. In [nayak2019optimally] and [sahoo2016noising], the Hilbert transform with an adaptive thresholding technique was utilized to detect Rpeaks. Some thresholdbased techniques with other criteria have also been used to specify the threshold. In [hou2018real], an adaptive threshold concerning the geometric angle between two consecutive samples of the ECG signal was defined. The performance of the thresholdbased technique is highly dependent on the selection of initial parameters; hence, it can lead to a significantly higher number of false beats. Therefore, independent threshold techniques are more desirable than the thresholdbased technique.
Most of the stateoftheart methods for Rpeak detection are based on wavelet transform [park2017r, rakshit2017efficient, thirrunavukkarasu2019detection], simple mathematical operations [nayak2019optimally, PanTompkins1985, gutierrez2015novel]
, hidden Markov models, and machine learning. Wavelet transform is a suitable approach for considering the nonstationary behavior of the ECG signal. However, considering the various shapes of the QRS complex, it is difficult to select the optimal mother wavelet or find the required threshold in the detection step of the wavelet transform. Additionally, discrete wavelet transform fails to provide reliable results in a shortrecording duration. Mathematical operationbased algorithms have a low computational cost, which is more appropriate for realtime applications and large dataset analysis. However, achieving high performance when the signaltonoise ratio is high remains challenging for these algorithms. Hidden Markov models are also widely used in ECG segmentation because they are powerful tools for considering the temporal dependency among the waveforms
[akhbari2016ecg, monroy2020hidden, huque2019hmm]. The majority of the studies on machine learningbased methods have utilized sparse signal processing to represent an approximation of the nonlinear ECG signal using sparsity constraints [ning2013ecg, zifan2006automated, parviziomran2019optimization, tavakoli2013decoding, monti2018adaptive, tibshirani1996regression]. Some studies have also applied deep learning techniques to detect the ECG waveforms considering its high performance in various classification tasks
[xiang2018automatic, tavakoli2017learning]. However, the caveat with deep learningbased approaches is that they need largescale datasets for the training phase and often suffer from the imbalanced class problem [8683140, mousavi2020single].In this paper, we propose a new class of graphical models based on optimal changepoint detection models, named the graphconstrained changepoint detection (GCCD) model, to locate Rpeaks in the ECG signal. A changepoint detection model identifies abrupt changes in data when a property of the time series changes. In the nonstationary ECG signal, ECG waves can also be considered as abrupt up or down changes over time during the heart cycle. We exploit the model introduced by Hocking et al. [hocking2017log, hocking2018generalized], in which a graphbased optimal changepoint detection model was used for detecting abrupt changes in the genomics data. In their work, they propose a new class of functional pruning algorithms with loglinear time complexity in the amount of data, which is capable of handling the large datasets that are common to ECG analysis.
Only a few studies in the literature have applied changepoint detection models for cardiac analysis. Gold et al. [gold2018doubly]
adopted a changepoint detection method based on Bayesian inference to extract the onset of the QRS complex over a small time window containing just one QRS complex. In
[qi2014novel], a changepoint detection approach based on the Haar wavelet and KolmogorovSmirnov statistic was applied to find normal and abnormal ECG segments within the assembled ECG samples from different ECG datasets. Sinn et al. [sinn2012detecting] analyzed heart rate changes in ECG recordings by detecting abrupt changes in the ordinal pattern distributions, which are used to represent the order structure of a time series. Some studies have also applied changepoint detection models to investigate sleep problems by analyzing heart rate variability in the ECG signal during sleep [staudacher2005new, ma2020multiple].To the best of our knowledge, this is the first study in which changepoint detection models have been proposed to detect ECG fiducial points in long records of ECG signals. In this novel framework, prior biological knowledge about the expected sequence of changes can be specified in a constraint graph. Then, functional pruning dynamic programming algorithms can compute the globally optimal model (mean, changes, and hidden states) in fast loglinear time. We furthermore propose a new algorithm for learning the graph structure using labeled ECG data. Therefore, the main contributions of this study are:

A new class of graphical models based on optimal changepoint detection models to detect Rpeak positions in the ECG signal. The proposed method does not require any noise removal preprocessing step as it uses the sparsity of changepoints to detect abrupt changes.

A new algorithm to learn the graph structure and parameters using labeled ECG data. Thus, the model’s performance is no longer dependent on an expert to encode prior knowledge into the constraint graph.

Comparison of the learned graphs with the manually constructed graphs in terms of graph structure and detection accuracy. Results demonstrate that there can exist different optimum graph structures for one subject, and the proposed graph learning algorithm can find global optima depending on the initial graph structure.
The rest of the paper is organized as follows. In the next section, we describe the proposed model for Rpeak detection in the ECG signal. We explain the GCCD model in Section IIA and the constraint graph in Section IIB. Section IIB also defines the manual graph and the proposed graph learning algorithm. Section III provides a description of the dataset used in this study and a discussion of the results as well as a comparison between the performance of the manually defined graphs and learned graphs. Finally, Section IV summarizes this research work and its contributions.
Ii Methodology
The proposed method treats ECG wave detection as a changepoint detection problem for a nonstationary ECG signal. It extracts the Rpeaks in the raw ECG signal by representing the periodic nonstationary ECG signal as a piecewise locally stationary time series with constant mean values (i.e., each piece is the mean of one segment of datapoints). The model takes a raw ECG signal and a constraint graph as inputs and computes the onset/offset and the mean of desired segments (i.e., hidden states). Then, the center of each state is associated with the location of a peak. The constraint graph allows the incorporation of prior knowledge into the model and regularizes the model. Figure 1 illustrates an overview of the proposed algorithm in the detection of Rpeak positions in the ECG signal. It is worth reemphasizing that the model takes the raw ECG signal as the input, without applying any preprocessing step, as it leverages the sparsity of changepoints to denoise the signal and to detect abrupt changes.
The constraint graph, which encodes the expected sequence of changes in the ECG signal, can be defined manually by an expert or automatically from the data. In the following sections, we describe the details of various parts of the proposed model.
Iia GraphConstrained Changepoint Detection Model
ECG fiducial points detection can be defined as the problem of finding abrupt changes over one cardiac cycle caused by changes in statistical characteristics. From this point of view, a proper changepoint detection algorithm can be employed to detect ECG waves in a fast and effective way. We applied the optimal changepoint detection model introduced in [hocking2017log] to localize Rpeak positions in the ECG signal. In this model, prior biological knowledge about the expected sequence of changes is incorporated into the model as a graph constraint. Then, a dynamic programming algorithm using functional pruning computes the globally optimal model (mean, changes, and hidden states) in fast loglinear time.
We assumed a directed graph as the constraint graph, where the vertex set represents the hidden states/segments (not necessarily a waveform), and the edge set represents the expected changes between the states/segments. Each edge incorporates the following associated prior knowledge about the expected sequences of changes:

The source and target are vertices/states for a changepoint from to .

A nonnegative penalty constant is the cost of changepoint .

A constraint function defines the possible mean values before and after each changepoint . If is the mean before the changepoint and is the mean after the changepoint, then the constraint is . These functions can be used to constrain the direction (up or down) and/or the magnitude of the change (greater/less than a certain amount).
Mathematically, given the input signal and the directed graph , the problem of finding changepoints , segment means , and hidden states can be solved using the following optimization problem:
(1)  
s. t.  (2)  
(3) 
The changepoints can be assigned to any of the predefined edges (). Consequently, indicates no change with zero cost, . Function (1) consists of a datafitting term and a model complexity term [amini2019iterative, amini2019iterative2]. represents the negative loglikelihood of each datapoint, and is a nonnegative penalty on each changepoint. In other words, regularizes the number of predicted changepoints/segments by the model so that a larger
reduces the number of changepoints by estimating a more sparse changepoint vector. The constraint function
also encodes the expected up/down change and the least amplitude gap between the mean of two states. When there is no change , Constraint (2) forces the model to stay in the current state with no change in mean . However, when there is a change , Constraint (3) imposes a change in the mean implied by the constraint function as well as a change in the state. An opensource implementation of the Generalized Functional Prunining Optimal Partitioning (GFPOP) algorithm is available in C++ code inside an R package named GFPOP on GitHub
[GFPOP].IiB Constraint Graph
The constraint graph in the optimization problem of Equation (1) encodes prior biological knowledge about the expected sequences of changes within one cardiac cycle. It can be designed manually by an expert or be learned from the data by the model. The two following subsections detail both the manual and learningbased designs.
IiB1 Manual Graph Definition
To manually define the constraint graph , we took into account the possible morphological categories for the ECG waves (i.e., P, QRS, and T waves) and the overall morphological properties of the signal in each record [fotoohinasab2020graph]. An expected hidden state/segment in the signal is characterized as a node in the constraint graph, and the required conditions for transition between states are encoded in the edges. The required conditions are determined based on the expected minimum amplitude difference of two successive states and the polarity of each transition (i.e., up/down).
The caveat with the manual definition of the constraint graph is that it can be inefficient for ECG signal analysis considering the various morphological patterns for each waveform. Furthermore, the model’s performance depends on the expert knowledge encoded into the constraint graph. In the next subsection, we explain the proposed graph learning algorithm for learning the constraint graph using the Rpeak labels provided by the gold standard.
IiB2 Constraint Graph Learning
To automate the Rpeak detection task, we modified the previous model by learning the constraint graph from the data (see the dashed part in Figure 1). In this new framework, the proposed model takes the raw signal and an initial graph structure as inputs and yields the desired outputs, including the onset/offset and the mean of segments specified in the nodes of the learned constraint graph [fotoohinasab2021graph]. Here, the model architecture is comprised of two stages: training and detection. The training step tries to heuristically find an optimum graph structure by which the label errors in the training set are minimized (the block named “Graph Learning Algorithm” in Figure 1). The detection step then extracts the Rpeaks in the raw ECG record constrained to the graph learned in the previous step (the block named “Changepoint Detection Model” in Figure 1).
The novelty of this new structure lies in the training step, which is comparable to the previous model in Section IIB1. The main idea of the training step is to automatically discover the desired topology of the constraint graph and the information about the edges from the data. As described in Section IIA, each edge contains the following information: (1) the expected up/down change in the segment means, (2) the least amplitude gap between the means of two states, and (3) a nonnegative penalty imposed by the edge transition. Suppose that the initial graph for each record is denoted as , where and are the corresponding graph node and edge sets, respectively. Each node in the set represents initial hidden states in the model. Each edge in the set represents a transition between two consecutive hidden states (i.e., a changepoint from the source to the target in section IIA) and also contains initial values for parameters of , , and , which are the initial type, the initial gap between two states, and the initial penalty, respectively. Figure 2a shows the simple initial graph used for the optimization process. It should be noted that the initial edge information was chosen based on the overall results obtained from the manual definition of the constraint graph.
A sketch of the proposed graph learning algorithm is summarised as Algorithm 1. The greedy graph search algorithm starts with the initial graph and iteratively optimizes the graph structure and edge parameters to find a graph that maximizes the accuracy regarding the provided labels. At the th iteration, the function finds the graph candidate set using the editing candidates for each edge of the output graph from the previous iteration . In this study, the algorithm considers 11 editing candidates per edge to optimize the graph topology and the three edge parameters. For example, in the iteration , if the parent graph (i.e., ) has two edges, the graph candidate set will have no more than 22 members . These editing candidates include three types of adding a node, two types of deleting a node, one type of adding two nodes, changing the type of the abrupt change, and increasing or decreasing the penalty and gap corresponding to an edge. We believe all morphological patterns of the ECG waves can be constructed using these editing candidates. Figure 2 illustrates the graph editing candidates related to the edge with an up change.
IiC Computational Complexity
As can be seen in Algorithm 1, the time complexity of the GCCD algorithm is theoretically proportional to the number of graph candidates at each iteration (Line 10) and the number of required iterations to achieve an optimum graph with minimum label errors (Line 5). There are three main aspects that characterize the time complexity of the algorithm:

Given a record with data samples and a graph candidate with vertices and edges, the time complexity to detect Rpeaks (Lines 11–15) is in the worst case (pathological simulated data) and in the average case (typical in real data). Also note that since we consider only graphs with a single circular path, , and the time complexity is further reduced to (for average case/nonpathological data).

Considering graph edit candidates in the iteration , the time complexity to compute all the models is (where is the time complexity of solving for optimal model parameters given a single graph). It should be noted that the number of graph candidates in the iteration depends on , which is the graph from the previous iteration (Line 10). The time complexity to compute the label error given labels is , which can be effectively ignored from the overall time complexity as this task is fast.

Finally, iterating over iterations to obtain the graph with the minimum label error (Line 5) causes the overall time complexity of the algorithm to be , where is the time to solve for a single graph, and is the number of edit candidates considered in each iteration.
Iii Experimental Studies
Iiia Dataset
We applied the wellknown MITBIH Arrhythmia (MITBIHAR) database to evaluate the GCCD model. This database contains 48 ECG recordings taken from 47 subjects. Each record’s duration is 30 min, and each recording is sampled at 360 Hz with a resolution of 200 samples over a 10 mV range [moody2001impact, goldberger2000physiobank]. Each recording consists of two ambulatory ECG channels from the modified lead II (MLII) and one of the leads V1, V2, V4, or V5. In this study, all 48 records with one MLII or V5 lead were used to evaluate the algorithm. The database has been annotated with both RR intervals and heartbeat class information by two or more expert cardiologists independently.
IiiB Results and Discussion
This section presents a comprehensive discussion of the results obtained by the proposed model and a detailed comparison between the manually defined graphs and the learned graphs. We also provide some suggestions for the future development of the GCCD model.
Figure 3 illustrates an example of the model’s performance with a manually defined constraint graph in the Rpeak detection task for a window of Record 230 of the MITBIHAR dataset. However, as mentioned in Section IIB1, the performance of the model using manually defined graphs depends on an expert with prior knowledge. Furthermore, manual annotation by an expert is time consuming and expensive. To address this issue, we proposed a new graph learning algorithm that searches for a locally optimal constraint graph using a greedy scheme on the labeled ECG data. Regarding the various morphological patterns for the ECG signal, the proposed graph learning algorithm can relax the model from the manual definition of the constraint graph for each record.
We adopted the intrapatient paradigm to train the constraint graph to address the intrapatient variation in ECG morphologies. Thus, the training and testing sets were generated by randomly splitting the intrasamples for each record with an approximate ratio of . We used a fold crossvalidation approach to evaluate the model performance with a size of . More specifically, we divided the intrasample data into folds so that each trial used four folds to train the model and one fold for validation.
Figures 4 and 5 show representative examples of the Rpeak detection task performed by the model integrated with the graph learning algorithm for two records from the MITBIHAR database. These figures illustrate how the proposed graph learning algorithm iteratively edits the graph structure to yield a model with maximum accuracy in detecting Rpeaks. We initialized the constraint graphs using the graph structure in Figure 2a with the initial values of and for Record 107 and and for Record 219. It should be noted that the initial edge information was assigned based on the overall results derived from the manually defined graphs in all experiments. However, graph candidates 2f and 2g can adjust the parameters and for the optimum values. For these two examples, we chose the initial edge information so that all the training steps could be completely displayed. Label errors are omitted from Figures 5a–5c to reduce clutter in the figures. The red part of the graph structure in each iteration presents the chosen editing candidate in the current iteration over the graph in the previous iteration. More interestingly, Figure 5 demonstrates the model’s capability to detect Rpeaks in the presence of a baseline wandering artifact, which is a typical artifact in the ECG signal. Baseline wandering can change the shape of the QRS complex and thereby causes incorrect detection of the Rpeak. The performance of the Pan and Tompkins [PanTompkins1985] algorithm, algorithms derived from the ECG signal slope [sabzevari2018ultra], and methods based on wavelet transform are highly dependent on the removal of this artifact. Figure 6 shows the test result for this record over two different time windows of data. Figure 7 illustrates the training progress for these two records, where the Yaxis shows the sum of false negative and false positive error rates. Indeed, the training progress curve reflects the number of label errors produced by the model in each iteration given the provided labels for the training and validation sets. It is worth mentioning that the proposed graph learning algorithm avoids possible overfitting issues as it tries to extract the morphology of the ECG signal that contains multiple various morphological patterns.
The proposed graph learning algorithm employs a greedy search scheme to select the best performing graph in terms of detection accuracy (see Section IIB2). Therefore, the performance of the algorithm depends heavily on the initial graph structure and will likely lead to local optima. Figure 8 compares the training progress for Record 106 of the MITBIHAR database initialized with the two simple (see Fig. 2a) and complex graph structures (i.e., a graph with eight nodes representing the morphology of a normal ECG signal). Figure 9 also presents a comparison of the final selected graphs and their performances for a window of this record. As these figures show, the model initialized with the complex graph structure can achieve higher accuracy (i.e., a lower number of label errors) in a lower number of iterations than the model initialized with the simple graph structure.
The investigation of the experimental results shows that the greedy graph search algorithm can achieve optimal performance for the model trained with the manually defined graphs, although its performance is affected by the initial graph. We noticed that for most of the records from the MITBIHAR database, the learned graphs could reach the performance of the manually defined graphs but with different graph structures. This means that the GCCD model can obtain global optima using various initialization structures, which will likely lead to different final graph structures. Figure 10 compares the constraint graph structures defined manually vs. those learned automatically using the initial graph structure in Figure 2a for Record 100 of the MITBIHAR dataset. As shown in this figure, the manually defined graph and the learned graph both achieved the optimal performance but with different graph structures. We also noticed that for some records from the MITBIHAR database, the graph learning algorithm chose the same structure as the manually defined graph structure. Figure 11 shows the model performance using the graph learning algorithm for Record 232 from the MITBIHAR dataset, for which the manually defined constraint graph and the learned graph had the same structures.
Different metrics were adopted to evaluate the performance of the proposed model with both the manual and learningbased graph designs. These metrics included sensitivity (), positive predictivity rate (), and detection error rate (), which are calculated by:
(4)  
(5)  
(6) 
where is true positives, is false positives, is false negatives, and is true negatives. Table I presents the performance of the proposed model regarding both the manually defined and learning graphs against other stateoftheart methods for Rpeak detection (QRS complex). As shown in the table, the proposed algorithm achieved Sen = %99.76, PPR = %99.68, and DER = 0.55 for the manual definition of the constraint graph and Sen = %99.64, PPR = %99.71, and DER = 0.19 for the learning constraint graph using the MITBIHAR database. Note that the model constrained to the manually defined graphs outperformed the model combined with the graph learning algorithm because in the latter, the model’s performance was largely dependent on the initial graph structure.
Method  Sen (%)  PPR (%)  DER (%) 

Park et al. [park2017r]  99.93  99.91  0.163 
Farashi [farashi2016multiresolution]  99.75  99.85  0.40 
Sharma and Sunkaria [sharma2017qrs]  99.50  99.56  0.93 
CastellsRufas and Carrabina [castells2015simple]  99.43  99.67  0.88 
GCCD Model with Manual Definition of the Constraint Graph  99.76  99.68  0.55 
GCCD Model with Learning of the Constraint Graph  99.64  99.71  0.19 
ECG recordings in the MITBIHAR database were chosen to challenge the Rpeak detection task because they represent a wide variety of QRS morphologies with realworld variability. Our proposed model yielded outstanding results when detecting Rpeaks in these tricky records. Records 103, 104, 105, 108, 111, 112, 116, 200, 201, 203, 205, 208, 210, 217, 219, 222, and 228 are comprised of abrupt changes in ECG morphology, and they are severely affected by noise and artifacts. Figure 5 shows the capability of the model to detect Rpeaks in the presence of baseline wandering noise. We reemphasize that these comparable results were obtained without applying any preprocessing operations, as opposed to other methods in the literature. Records 108, 113, 117, 201, 202, 203, 213, 219, 222, 223, 231, and 232 contain many peaks with unusual amplitudes. Smallamplitude Rpeaks or highamplitude P and Tpeaks embedded in highamplitude QRS complexes can lead to high FN and FP errors in the Rpeak detection task. As a representative example, Figure 4 illustrates the efficiency of the GCCD model in Rpeak detection for Record 117, which contains many beats with highamplitude Tpeaks.
The experimental results obtained using the proposed model justify changepoint detection models as a potential approach to extract ECG fiducial points. In this study, we demonstrated the capability of the GCCD model in locating Rpeaks within various morphological patterns of ECG. The proposed greedy graph search algorithm can potentially detect ECG waves other than the R wave (i.e., P, Q, S, and T waves) by considering corresponding prior knowledge of the graph editing candidates. We noticed that in Records 114, 200, 203, 207, and 210, the Sen and PPR values were less than 99%. These records contain multiple different morphological patterns, including negative QRS complexes, and Records 200 and 203 have several QRS complexes with ventricular arrhythmias. The constraint graph for these records involves learning a graph with more than one optimum graph path. Learning a multipath constraint graph is also required to detect all ECG waves due to the various morphological patterns of each wave incorporated into the graph. The other point that should be considered here is that the GCCD model estimates the ECG signal using a Gaussian function. A modified model with a multiGaussian fitting method can drastically improve the ECGrelated changepoint detection task.
Future work should focus on developing the proposed model with a multiGaussian fitting and a multipath graph learning algorithm. Incorporating these modifications into the proposed model could provide a promising platform for evolving new graphbased tools to detect and classify heart arrhythmias. A multipath graph learning algorithm could reveal the morphology of the ECG signal (time duration, amplitude, and direction of each wave) in each cardiac cycle. Subsequently, new graphbased features could be extracted from the constraint graph path for an ECG cycle to classify heartbeats.
Iv Conclusion
The accurate delineation of Rpeaks in the ECG signal plays a crucial role in most automated ECG analysis tools. This paper proposed a novel graphical model based on changepoint detection techniques for detecting Rpeaks within a nonstationary ECG signal. The proposed model was highly successful at detecting Rpeaks in noisy ECG data without applying any preprocessing steps. To our knowledge, this is the first time that a changepoint detection model has been applied for ECG fiducial points detection. In this new framework, prior biological knowledge about the expected sequences of changes was incorporated into the model using a graph. We defined the constraint graph manually and automatically using a proposed greedy graph search algorithm. Using the proposed graph learning algorithm, the initial graph structure can develop into a structure containing edge parameters with maximum detection accuracy for a record. The experimental results provided in this paper demonstrate that the GCCD model can be a promising approach for detecting ECG waves and developing new graphbased tools for further ECG analysis. The proposed graphical model approach can be advanced by learning a multipath constraint graph and fitting a multiGaussian curve model to the ECG signal, which should be considered in future studies.
Acknowledgment
This material is based on work supported by the National Science Foundation under Grant Number 1657260. Research reported in this publication was supported by the National Institute on Minority Health and Health Disparities of the National Institutes of Health under Award Number U54MD012388.
Comments
There are no comments yet.