1.1. Trabajos afines

De acuerdo con Charlton et al. ^[10], los algoritmos PPG-RR pueden resumirse en cinco pasos:

Algunos de los primeros enfoques que estiman la RR a partir de PPG utilizan dos filtros Butterworth de segundo orden, uno pasa-alto y uno pasa-bajo con frecuencias de corte 0,1 Hz y 5 Hz, respectivamente, para eliminar el ruido de la PPG. Luego se estima la señal respiratoria aplicando un filtro pasa-bajo de 0,4 Hz a la señal PPG ^[11]. Madhav et al. ^[12] aplicaron una relación de valor singular (Singular Value Ratio, SVR) para extraer la periodicidad respiratoria de la PPG, y utilizaron análisis de los componentes principales (Principal Component Analysis, PCA) para estimar la actividad respiratoria del primer componente principal. En otro enfoque propuesto por Nam y colaboradores ^[13], se capturó una señal PPG utilizando la cámara de un teléfono inteligente, y luego se estimó la RR encontrando la frecuencia correspondiente al pico más alto del espectro generado por el periodograma de Welch ^[7].

Se compararon cinco metodologías diferentes en el trabajo realizado por Pimentel y colaboradores ^[4], incluyendo la desarrollada por los autores. La más antigua fue propuesta por el equipo encabezado por Nilsson ^[14], y consiste en aplicar un filtro Bessel pasabanda de orden 16 con frecuencias de corte 0,13 Hz y 0,48 Hz (7,8 a 28,8 respiraciones/min). La segunda propuesta ^[8], utiliza la transformada de Fourier sobre la señal PPG y obtiene la frecuencia con la mayor amplitud en el ancho de banda entre 0,08 Hz y 0,4 Hz. En la metodología sugerida por Tarassenko ^[5], los modelos autorregresivos (AR) estiman un filtro utilizando la señal PPG, y de acuerdo con ellos la frecuencia del polo de mayor magnitud dentro del intervalo del ancho de banda de la RR (0,08 Hz a 0,7 Hz) corresponde a la RR. En lugar de utilizar la propia PPG ^[6] se estiman tres características temporales relacionadas con la RR: variación de la intensidad inducida por la respiración (Respiratory-Induced Intensity Variation, RIIV), variación de la amplitud inducida por la respiración (Respiratory-Induced Amplitude Variation, RIAV) y variación de la frecuencia inducida por la respiración (Respiratory-Induced Frequency Variation, RIFV). Luego, se calcula el espectro de potencia de cada una y se estima la RR a partir de cada una tomando la frecuencia de mayor amplitud en la banda de frecuencias de la RR (0,067 Hz – 1,08 Hz) y calculando la media de las tres frecuencias. Finalmente, la metodología mostrada por Pimentel et al. ^[4] combina el modelado AR propuesto en Tarassenko ^[5], con la fusión de resultados de las tres características propuestas por Karlen y su equipo ^[6].

Los resultados presentados ^[4] muestran que para el conjunto de datos CapnoBase ^[15], el método que presentó el error más bajo ^[6], es decir, con un error absoluto de la mediana (Median Absolute Error, MdAE) de 0,8 respiraciones/min utilizando una ventana de 62 s; sin embargo, 46 % de los datos de PPG fueron excluidos debido al ruido. El segundo mejor en este conjunto de datos fue el citado por Tarassenko ^[5], que alcanzó un MdAE de 1,1 respiraciones/min y conservó 92 % de la data. El conjunto de datos BIDMC parece ser más desafiante, debido a que el mejor MdAE fue obtenido por el método propuesto por Shelley ^[8] (2,3 respiraciones/min), el cual también fue uno de los métodos que conservó la mayoría de la data (94 %). Un trabajo reciente ha obtenido resultados aún más precisos ^[16], en el que un método deriva múltiples ondas de una señal PPG, mide su calidad y luego la utiliza como ponderación para combinarlas utilizando un suavizador de Kalman. Con este enfoque obtuvieron un error mediano de 0,2 respiraciones/min

1.2. Redes neuronales profundas (Deep Neural Networks, DNN)

Para Goodfellow et al. ^[17], el aprendizaje profundo es un enfoque que puede manejar un paso importante en un problema de aprendizaje de máquina: la extracción de características. Su concepto consiste en concatenar múltiples capas de modelos simples, donde cada modelo «aprende» parte de los problemas-conceptos y, por tanto, un problema complejo puede dividirse en modelos más simples para ser resuelto de mejor manera. Un punto interesante de este enfoque es que realiza las tareas de aprendizaje automático (clasificación, regresión, agrupamiento, etc.) y también la extracción de características, en contraste con las técnicas poco profundas, en el que el paso de extracción de características debe hacerse de manera separada, en una etapa previa, y exige conocimiento del área. La limitación del enfoque de aprendizaje profundo es el conjunto de datos. Normalmente, los modelos de aprendizaje profundo tienen un gran número de valores a ser ajustados y, por tanto, necesitan grandes cantidades de datos para optimizarlos.

Existen algunos enfoques que instrumentaron técnicas de aprendizaje profundo para resolver problemas relacionados con data fisiológica. Biswas et al. ^[9] implementaron una DNN para estimar la frecuencia cardíaca (Heart Rate, HR), teniendo como data de entrada una señal PPG obtenida de un sensor colocado en la muñeca del usuario. El modelo fue probado en un conjunto de datos desafiante, en el cual la señal PPG fue capturada durante varias actividades físicas, las cuales insertan ruido y artefactos en ella. La arquitectura del modelo estaba constituida por dos capas convolucionales de una dimensión seguidas de dos capas de gran memoria de corto plazo (Long Short-Term Memory, LSTM) y, finalmente, una capa densa. Esta arquitectura alcanzó un error absoluto medio (Mean Absolute Error, MAE) de 1,47 ± 3,37 latidos por minuto (Beats Per Minute, BPM).

En el trabajo ^[1] se presenta un enfoque de aprendizaje profundo para extraer la señal RR a partir de la señal PPG. Su modelo consta de varias capas convolucionales conectadas a varias capas deconvolucionales que transforman la señal PPG en la señal respiratoria, y luego calculan a partir de ella la RR. Su método fue probado en dos conjuntos de datos independientes, CapnoBase y Vortal ^[18], y alcanzó un error cuadrático medio (Mean Squared Error, MSE) de 0,262 en el primer conjunto de datos y de 0,145 en el segundo, lo que corresponde a una raíz cuadrada del error cuadrático medio (Root Mean Squared Error, RMSE) de aproximadamente 3,1 respiraciones/min y 2,3 respiraciones/min, respectivamente.

2.1. Conjuntos de datos

En las pruebas se utilizaron dos conjuntos de datos: el benchmark CapnoBase ^[15] y el de respiración BIDMC PPG ^[19]. El primero consta de 42 registros ^[20], de 8 minutos cada uno, que contienen la PPG y la señal de dióxido de carbono (CO₂) inhalado y exhalado, ambos recolectados a una frecuencia de muestreo de 300 Hz. El conjunto BIDMC PPG ^[21] contiene 8 minutos de registros de señales PPG y respiración de CO₂ de 53 voluntarios, recolectados a 125 Hz.

2.2. Características respiratorias del fotopletismograma

A partir de la PPG se generaron tres señales relacionadas con la respiración ^[6]:

La Figura 1 muestra una representación gráfica de las características RIIV, RIAV y RIFV.

Figura 1. Características relacionadas a la RR extraídas de la PPG 

2.3. Índice de calidad de la señal (Signal Quality Index, SQI)

El SQI cuantifica la contaminación por ruido y artefactos en la ventana PPG analizada. Este índice se utiliza para excluir las señales que no alcanzan una calidad mínima. En este trabajo, el SQI se calcula abordando dos puntos: «líneas planas» en la señal, y concordancia entre diferentes detectores de picos ^[4]. Las «líneas planas» se definen como las partes de la señal PPG en las que la diferencia entre muestras consecutivas es cercana a cero. Para detectarlas se realizó el siguiente método:

La concordancia entre dos detectores de picos fue realizada de la siguiente manera: si la diferencia entre la posición encontrada por los detectores en el mismo pico es menor que 150 ms, se considera que ambos detectores concuerdan en ese pico y, por lo tanto, su posición es correcta; de lo contrario, el pico se considera como un ruido ^[4]. Entonces se aplica un score F1 para medir la concordancia entre dos detectores de picos en el segmento analizado de la señal. Los métodos de detección de picos fueron proporcionados por dos librerías diferentes de Python: Neurokit2 ^[22], que implementa el algoritmo de detección de picos PPG descrito en ^[23] y la librería Heartpy ^[24], en la que los autores desarrollan su propia técnica de detección de picos.

Finalmente, el SQI se calcula como: SQI = F1×K, donde F1 es la concordancia entre picos y K es la proporción de «línea plana»

2.4. Arquitectura DNN

La arquitectura del modelo se inspiró es la CorNET ^[9], que tiene dos bloques convolucionales, cada uno de los cuales contiene una capa convolucional (32 filtros con un kernel de tamaño 40), una capa de normalización batch, una capa de activación con unidades lineales rectificadas (Rectified Linear Units, ReLU) y, finalmente, una capa maxpooling con un tamaño de pool de 4 (para el entrenamiento se utilizó una capa dropout con una tasa de eliminación (drop rate) de 0.1). Luego de los bloques convolucionales se concatenan dos capas LSTM, cada una con 128 unidades y función de activación tangente hiperbólica (tanh). Para la predicción se utilizó una neurona con activación lineal.

La arquitectura propuesta para predecir la RR es el resultado de la optimización mediante algoritmos genéticos, de los parámetros de la CorNET que se presentan en la Tabla 1.

Tabla 1. Parámetros probados en la arquitectura DNN 

Los parámetros mencionados arriba son optimizados, mientras que para los parámetros no mencionados se utiliza la definición por defecto en la librería Keras ^[25]. Considerando que el número de combinaciones es muy grande (2916), los parámetros se optimizan utilizando algoritmos genéticos. Esta estrategia comienza definiendo una población de L modelos con parámetros aleatorios, y luego todos los modelos son entrenados, probados y categorizados de acuerdo con el resultado. Los M modelos con mejor desempeño se conservan mientras que los otros se descartan. El conjunto de los modelos conservados se denomina padres, y luego se generan M − J nuevos modelos tomando valores de tres modelos diferentes (un padre para los parámetros convolucionales, uno para las capas LSTM, y otro para las capas densas). Adicionalmente, existe una probabilidad de H % de que los parámetros de los nuevos modelos cambien aleatoriamente (estos cambios se denominan mutaciones); estos nuevos modelos se denominan descendientes. Asimismo, existen J modelos nuevos generados de forma completamente aleatoria, y a este conjunto se le asigna el nombre de extranjeros. Luego, los conjuntos de padres, descendientes y extranjeros se unen comienza una nueva interacción.

El procedimiento se repite por E épocas (epochs). La idea es que los parámetros que disminuyen el error de estimación prevalezcan en el conjunto de padres, y su combinación pueda generar un modelo en el conjunto de descendientes que mejore el resultado de sus padres. Asimismo, nuevos parámetros aleatorios se insertan con las mutaciones y los extranjeros creados en cada iteración. Estas adiciones aleatorias establecen una posibilidad de que la población de modelos salga de un mínimo local subóptimo. Para este experimento, los parámetros de la optimización mediante algoritmos genéticos fueron los siguientes: L = 45, M = 20, J = 5, E = 20 y H = 20, lo cual redujo el número de modelos entrenados de 2916 a 520 (25 nuevos modelos/época × 20 épocas + los primeros 20 padres). El mejor modelo encontrado es una DNN con 2 bloques convolucionales con 32 filtros y un MaxPooling de una dimensión (tamaño de kernel 2) en cada capa convolucional, un tamaño de kernel de 21, y activación ReLU. Luego, una capa LSTM con 32 neuronas, y finalmente, 2 capas densas, una con 32 neuronas y activación ReLU y la última con una neurona y activación lineal. La Figura 2 presenta un resumen de la arquitectura DNN final.

Figura 2. Arquitectura de red neuronal seleccionada 

2.5. Procedimiento de entrenamiento y prueba

Se utilizaron dos conjuntos de datos para entrenar y evaluar los modelos BIDMC ^[19] y CapnoBase ^[15]. Se empleó el método de validación cruzada con cuatro grupos: para cada base de datos, sus archivos se dividieron en cuatro grupos, cada uno de los cuales contenía el 25 % de los datos en cada conjunto. Luego, se realizó el mismo procedimiento en cada grupo: primero, se encontraron los picos de la PPG utilizando la librería Neurokit2 ^[22]. Posteriormente, utilizando una ventana deslizante de 64 s con un paso de 4 s, se calcularon las tres señales relacionadas a la RR (RIIV, RIAV y RIFV), que luego se interpolaron a una frecuencia de muestreo artificial de 4 Hz y se extrajeron las primeras 240 muestras (60 s) de la ventana, generando así una muestra con 240 líneas (pasos de tiempo) y tres columnas (correspondientes a RIIV, RIAV y RIFV). El paso final consiste en normalizar la data, de manera que cada canal de la nueva muestra tenga media 0 y varianza 1: Y = (X − µ(X))/σ(X), donde Y es la serie normalizada, X la serie de entrada, µ(.) es la operación que calcula la media y σ(.) es la operación que calcula la desviación estándar.

Para obtener el valor de referencia de la RR, se localiza el pico de CO₂ correspondiente a cada muestra utilizando el módulo Respiratorio de Neurokit2, y luego se calcula la mediana de la RR de la muestra de 64 s, tomando la inversa del valor de la mediana en el intervalo de tiempo entre picos consecutivos.

Posteriormente se dividen los cuatro grupos de muestras normalizadas en conjuntos de entrenamiento y prueba, donde el conjunto de entrenamiento tiene tres grupos y el de prueba tiene uno. Los procedimientos de entrenamiento y prueba se repitieron cuatro veces, para que cada grupo pudiera ser evaluado una vez en el conjunto de prueba. Cada uno de los 4 modelos (uno para cada grupo de prueba) fue entrenado por 1400 épocas utilizando el optimizador Adam con parámetros estándar de Keras ^[25]. El procesamiento se realizó en un computador Ryzen 5 3500X con 16 GB de memoria RAM y un NVIDIA GeForce RTX 2060 SUPER con 8 GB de memoria dedicada. Se utilizó la métrica del error cuadrático medio para calcular la pérdida (loss) en cada época. La Figura 3 presenta la curva de pérdida resultante.

Figura 3. Curva Ts para las 4 instancias de la red 

Los grupos de validación 1 y 4 tuvieron una pérdida decreciente hasta la época 1400, el grupo 2 tuvo su pérdida mínima alrededor de la época 400 y el grupo de validación 3 tuvo su mínimo alrededor de la época 200, y luego decreció nuevamente después de la época 800. El gráfico se utilizó para determinar el número de épocas de entrenamiento de las DNN, dado que la mitad de los grupos de validación se mantienen decreciendo hasta la época 1400 y la capa dropout parece evitar el sobreajuste en los otros grupos, debido a que la pérdida no se incrementa mucho luego de los mínimos. El proceso se repitió revolviendo la data y los grupos, y entrenando cuatro DNN (una para cada grupo de entrenamiento-prueba) durante 1400 épocas para realizar la validación cruzada.