Vida útil acelerada
Maximiliano Kniazev
2025-03-26
1 Introducción
El objetivo básico de los ensayos acelerados de vida útil es realizar estudios a temperaturas elevadas, donde el deterioro sucede a mayor velocidad, para luego predecir el deterioro a temperaturas de almacenamiento menores.
2 Cinética de reacción
La cinética de la pérdida de la calidad de un alimento se representa de la siguiente forma:
\[ -\frac{dA}{dt} = k A^n \] Donde:
- \(A\) es el valor del factor seleccionado para determinar la calidad (como función del tiempo),
- \(t\) es el tiempo,
- \(k\) es una constante que depende de la temperatura,
- \(n\) es el exponente indicativo del orden de reacción.
Para \(n = 0\), se dice que la reacción es de orden cero, y se tiene que
\[ -\frac{dA}{dt} = k A^n\Big|_{n = 0} \overset{\text{Integración}}{\Rightarrow} A = A_o - kt \]
Algunos parámetros de calidad cuya pérdida en función del tiempo es orden cero son la degradación enzimática en frutas frescas y vegetales, alimentos congelados y pastas refrigeradas; pardeamiento no enzimático en cereales y en productos lácteos deshidratados; y oxidación de lípidos en alimentos congelados y deshidratados.
Para \(n = 1\), se dice que la reacción es de orden uno, y se tiene que
\[ -\frac{dA}{dt} = k A^n\Big|_{n = 1} \overset{\text{Integración}}{\Rightarrow} A = A_o e^{-kt} \overset{\text{Linealización}}{\Rightarrow}ln(A) = ln(A_o) - kt \] Algunos parámetros de calidad cuya pérdida en función del tiempo es orden uno es la rancidez en aceites o alimentos deshidratados; el crecimiento de microorganismos y sus defectos (aparición de mucílagos y sabores extraños); pérdida de vitaminas en alimentos enlatados y deshidratados; o pérdida en la calidad de proteínas de alimentos deshidratados.
Uno de los aspectos más importantes a recalcar es que cuando hay menos del 50% de conversión de reacción (el deterioro no ha concluído), es difícil distinguir entre orden cero y primer orden (Figura 1), se requiere un mínimo de 6 puntos (tiempos) para determinar el modelo.
Las reacciones que poseen tiempo de demora o incubación (log phase) pueden inducir a un error en el cálculo.
#Gráfico que simula qué sucede en orden cero y orden 1 para parámetros cualesquiera.
# Parámetros
dias <- seq(0, 300, by = 1)
A0 <- 100
k0 <- 0.77
k1 <- 0.01
# Crear dataframe con ambas curvas
datos <- tibble(
dias = dias,
orden0 = pmax(0, A0 - k0 * dias),
orden1 = A0 * exp(-k1 * dias)
) %>%
pivot_longer(cols = c(orden0, orden1), names_to = "orden", values_to = "calidad") %>%
filter(calidad > 0)
# Renombrar para etiquetas más limpias
datos$orden <- recode(datos$orden,
orden0 = "Orden cero",
orden1 = "Primer orden")
# Graficar
ggplot(datos, aes(x = dias, y = calidad, color = orden, linetype = orden)) +
geom_line(size = 0.5) +
geom_hline(yintercept = 50, color = "red", linetype = "dotted",
linewidth = 0.5, alpha = 0.7) +
scale_color_manual(values = c("black", "black")) +
# scale_linetype_manual(values = c("solid", "dashed")) +
labs(
title = "Desplome de A vs tiempo de almacenamiento",
caption = "Figura 1: Se muestra cómo se comporta la curva de cinética de orden 0 y 1 para datos arbitrarios. \n
Notar que por encima de la linea roja (conversión del 50%) \n no tenemos forma de decidir a qué modelo ajustan nuestros datos.",
x = "Almacenamiento (días)",
y = "A (% de calidad)",
color = NULL,
linetype = NULL
) +
theme_minimal(base_size = 14) +
theme(legend.position = "top") +
coord_cartesian(ylim = c(0, 100))
3 Ecuación de Arrhenius
El parámetro \(k\) de la cinética de reacción es una función de la temperatura (\(k = f(T) \sim k(T)\)), en particular, en la gran mayoría de las reacciones de pérdida de calidad alimentaria, el valor de \(k\) puede modelarse según la ecuación de Arrhenius:
\[ k(T) = k_{\text{ref}} \ e ^{-\frac{E_A}{R}(\frac{1}{T} - \frac{1}{T_{\text{ref}}})} \overset{\text{Linealización}}{\Rightarrow} ln(k) = ln(k_{ref}) - \frac{E_A}{R}(\frac{1}{T} - \frac{1}{T_{\text{ref}}}) \]
Donde:
- \(k\) es la constante de velocidad de reacción a una temperatura fija \(T\).
- \(k_{ref}\) es la constante de velocidad de reacción a una temperatura de referencia.
- \(E_A\) es la energía de activación (\(\text{cal}.\text{mol}^{-1}\)).
- \(R\) es la constante general de los gases (\(\text{cal}.\text{mol}^{-1}.\)ºK\(^{-1}\))
- \(T\) es la temperatura (ºK)
- \(T_{ref}\) es una temperatura de referencia en (ºK)
La temperatura de referencia se debe seleccionar en función del intervalo de temperaturas de trabajo. Por ejemplo, si se está realizando un estudio acelerado en el intervalo 15ºC - 40ºC, se puede seleccionar como referencia 27ºC (300ºK).
4 Ejemplo delucidado
Para poder inferir sobre las temperaturas más bajas, será necesario aplicar la ecuación de Arrhenius, cuyo parámetro clave es la energía de activación.
En este ejemplo, se trabaja sobre un caso donde existen datos a temperatura ambiente, lo cual no es usual, pero permite entender la metodología de forma inductiva.
## Cargar datos
df <- data.frame(
Temperatura = c(rep(20, 7), rep(35, 7), rep(45, 6)),
Dias = c(122,145,164,183,201,224,245,
11,20,30,39,48,56,61,
7,14,21,28,35,42),
Sabor_oxidado = c(4.0, 8.1, 6.6,16.9,19.3,21.2,28.2,
8.1,13.9,19,23,31.6,33.2,37.2,
2.3,7.5,15.3,24.4,36.9,43.4)
)Los datos cargados tienen este formato:
Tabla 1: Desarrollo del sabor oxidado de una muestra comercial de mayonesa en función del tiempo, a tres temperaturas de almacenamiento.
| 20°C | 35°C | 45°C | |||
|---|---|---|---|---|---|
| Días | Sabor oxidado¹ | Días | Sabor oxidado¹ | Días | Sabor oxidado¹ |
| 122 | 4,0 | 11 | 8,1 | 7 | 2,3 |
| 145 | 8,1 | 20 | 13,9 | 14 | 7,5 |
| 164 | 6,6 | 30 | 19,0 | 21 | 15,3 |
| 183 | 16,9 | 39 | 23,0 | 28 | 24,4 |
| 201 | 19,3 | 48 | 31,6 | 35 | 36,9 |
| 224 | 21,2 | 56 | 33,2 | 42 | 43,4 |
| 245 | 28,2 | 61 | 37,2 | - | - |
¹ Medido con una escala sensorial en la que 0 = “nada oxidado” y 100 = “muy oxidado”.
4.1 Decisión del modelo que ajusta mejor a los datos
Se debe decidir para nuestros datos de estudio si el modelo cinético para la característica deseada es de orden cero o uno, o lo que es lo mismo, decidir cuál de los siguientes pares de datos poseen un mejor ajuste lineal:
4.1.1 Asumimos un modelo de orden cero
Si el comportamiento del Sabor oxidado (A) en función del tiempo posee cinética de orden cero, entonces el modelo \(A = A_o - kt\) posee un buen ajuste. Como k depende de la temperatura, se debe tener en cuenta cada una de las temperaturas de trabajo por separado (Figura 2).
#Gráfico de Sabor oxidado en función de los días de almacenamiento para las tres temperaturas.
ggplot(df, aes(x = Dias, y = Sabor_oxidado, color = as.factor(Temperatura))) +
geom_smooth(method = "lm") +
geom_point() +
geom_label(
x = 100,
y = 30,
label = "A == A[0] - k*t",
parse = TRUE,
inherit.aes = FALSE) +
labs(
title = "Sabor oxidado vs tiempo",
caption = "Figura 2: Comportamiento del sabor oxidado en función del tiempo para las muestras de mayonesa en estudio.\n La recta graficada corresponde al modelo cinético de orden cero.",
x = "Almacenamiento (días)",
y = "Sabor oxidado",
color = "Temperatura (ºC)",
linetype = NULL
) +
theme_minimal(base_size = 14) +
theme(legend.position = "top",
plot.caption = element_text(hjust = 0, size = 9)) 
# Ajustar modelos lineales de orden cero y extraer coeficientes + R²
r2_data_cero <- df %>%
group_by(Temperatura) %>%
nest() %>%
mutate(
model = map(data, ~ lm(Sabor_oxidado ~ Dias, data = .x)),
coef = map(model, ~ coef(.x)),
A_o = map_dbl(coef, 1),
k = map_dbl(coef, 2),
r2 = map_dbl(model, ~ glance(.x)$r.squared)
) %>%
select(Temperatura, A_o, k, r2)
# Mostrar tabla
kable(r2_data_cero, format = "markdown")| Temperatura | A_o | k | r2 |
|---|---|---|---|
| 20 | -21.275397 | 0.1972179 | 0.9335263 |
| 35 | 1.868343 | 0.5770626 | 0.9895385 |
| 45 | -8.646667 | 1.2359184 | 0.9862011 |
Para el modelo de orden cero, las tres condiciones poseen un muy buen ajuste, y en la práctica esto puede interpretarse de dos formas: O bien no es necesario complejizar la situación y plantear un modelo cinético de primer orden dado que el orden cero ajusta las condiciones del problema; o bien no tenemos datos suficientes para tomar esta decisión y el modelo puede confundirse con órdenes superiores.
En este punto es también necesario recalcar que el punto de corte del estudio de vida útil lo decide el investigador. Es decir, si un sabor oxidado mayor a 2 finaliza la vida útil, entonces para todas las temperaturas, el modelo de orden cero cubre el 100% del desplome de calidad sensorial.
A modo de ejercicio, estimaremos el modelo orden uno en la siguiente sección.
4.1.2 Asumimos un modelo de orden uno
Si el comportamiento del Sabor oxidado (A) en función del tiempo posee cinética de orden uno, entonces el modelo \(ln(A) = ln(A_o) - kt\) posee un buen ajuste. Como k depende de la temperatura, se debe tener en cuenta cada una de las temperaturas de trabajo por separado (Figura 3).
#Gráfico de Sabor oxidado en función de los días de almacenamiento para las tres temperaturas.
ggplot(df, aes(x = Dias, y = log(Sabor_oxidado), color = as.factor(Temperatura))) +
geom_smooth(method = "lm") +
geom_point() +
geom_label(
x = 100,
y = 3,
label = "ln(A) == ln(A[0]) - k*t",
parse = TRUE,
inherit.aes = FALSE) +
labs(
title = "ln(Sabor oxidado) vs tiempo",
caption = "Figura 3: Comportamiento del sabor oxidado en función del tiempo para las muestras de mayonesa en estudio. \n La recta graficada corresponde al modelo cinético de primer orden.",
x = "Almacenamiento (días)",
y = "ln(Sabor oxidado)",
color = "Temperatura (ºC)",
linetype = NULL
) +
theme_minimal(base_size = 14) +
theme(legend.position = "top",
plot.caption = element_text(hjust = 0, size = 9))
# Ajustar modelos lineales de orden uno y extraer coeficientes + R²
r2_data_uno <- df %>%
group_by(Temperatura) %>%
nest() %>%
mutate(
model = map(data, ~ lm(log(Sabor_oxidado) ~ Dias, data = .x)),
coef = map(model, ~ coef(.x)),
`log(A_o)` = map_dbl(coef, 1),
k = map_dbl(coef, 2),
r2 = map_dbl(model, ~ glance(.x)$r.squared)
) %>%
select(Temperatura, `log(A_o)`, k, r2)
# Mostrar tabla
kable(r2_data_uno, format = "markdown")| Temperatura | log(A_o) | k | r2 |
|---|---|---|---|
| 20 | -0.3658069 | 0.0156594 | 0.8960661 |
| 35 | 1.9653748 | 0.0287475 | 0.9511734 |
| 45 | 0.6980458 | 0.0813649 | 0.9199158 |
El ajuste del modelo cinético de primer orden presenta un peor ajuste para todas las temperaturas de trabajo, por ende, trabajaremos sobre el orden cero para este caso.
4.2 Estimación de la energía de activación.
Existen tres formas de estimar la energía de activación, a continuación se realizará de las tres formas posibles a partir de estos datos: Regresión lineal básica, regresión lineal con intervalos y regresión no lineal.
4.2.1 Regresión lineal básica
Este método sencillo retoma directamente desde los resultados de la Sección 4.1.1.
El mismo consiste en tomar para cada temperatura de trabajo la pendiente de la regresión lineal para el modelo cinético seleccionado. Recordando que \(A = A_o + kt\), y la ecuación de Arrhenius (Sección 3), se tiene que \(ln(k) = f(\frac{1}{T})\) con \(f\) una función lineal cuya pendiente es \(-\frac{E_A}{R}\). Por lo tanto, calculamos dichas transformaciones y obtenemos la pendiente de la recta que ajusta \(ln(k)\) vs \(\frac{1}{T}\), con T en ºK.
k_method1 <- r2_data_cero %>%
select(Temperatura, k) %>%
mutate(`1/T (1/ºK)` = (273.15+Temperatura)^(-1),
`ln(k)` = log(k))
kable(k_method1, format = "markdown")| Temperatura | k | 1/T (1/ºK) | ln(k) |
|---|---|---|---|
| 20 | 0.1972179 | 0.0034112 | -1.6234461 |
| 35 | 0.5770626 | 0.0032452 | -0.5498045 |
| 45 | 1.2359184 | 0.0031432 | 0.2118143 |
lm_Ea <- lm(`ln(k)` ~ `1/T (1/ºK)`, data = k_method1)
summary(lm_Ea)##
## Call:
## lm(formula = `ln(k)` ~ `1/T (1/ºK)`, data = k_method1)
##
## Residuals:
## 1 2 3
## 0.01575 -0.04138 0.02563
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21.5905 0.8739 24.71 0.0258 *
## `1/T (1/ºK)` -6809.7868 267.3731 -25.47 0.0250 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05116 on 1 degrees of freedom
## Multiple R-squared: 0.9985, Adjusted R-squared: 0.9969
## F-statistic: 648.7 on 1 and 1 DF, p-value: 0.02498Es el modelo lineal significativo (p = 0.025) con pendiente (\(-k\)) de -6809.79. Sin embargo, el intervalo de confianza (95%) de este parámetro es:
as.data.frame(confint(lm_Ea, level = 0.95)) %>%
mutate(`Mean value` = rowMeans(across(1:2))) #Mostrar también el promedio del intervalo de confianza. ## 2.5 % 97.5 % Mean value
## (Intercept) 10.48682 32.6942 21.59051
## `1/T (1/ºK)` -10207.08362 -3412.4901 -6809.78684Por lo tanto, tenemos un valor de \(E_AR^{-1}\) de \((6809.8 \pm 3397.3)ºK^{-1}\). Como se puede ver, el problema de este método de cálculo es su error. Al utilizar este camino estaremos cometiendo un error de casi el 50% el 95% de las veces.
4.2.2 Regresión lineal con intervalos
La forma trivial de mejorar el resultado anterior es agregar una mayor cantidad de valores de temperatura al estudio, sin embargo, esto suele ser costoso. Una estrategia matemática que mejora el error es tomar los extremos del intervalo de confianza para la variable k despejada en la Sección 4.1.1 anteriormente calculado como puntos de una misma temperatura y recalcular. Es decir, utilizar los siguientes datos para repetir la lógica de la Sección 4.2.1.
k_method2 <- df %>%
group_by(Temperatura) %>%
nest() %>%
mutate(
modelo = map(data, ~ lm(Sabor_oxidado ~ Dias, data = .x)),
coef = map(modelo, ~ coef(.x)["Dias"]),
conf = map(modelo, ~ confint(.x)["Dias", ])
) %>%
mutate(
k = map_dbl(coef, identity),
k_min = map_dbl(conf, 1),
k_max = map_dbl(conf, 2)
) %>%
select(Temperatura, k, k_min, k_max)
# Mostrar tabla
kable(k_method2, format = "markdown")| Temperatura | k | k_min | k_max |
|---|---|---|---|
| 20 | 0.1972179 | 0.1367181 | 0.2577177 |
| 35 | 0.5770626 | 0.5088524 | 0.6452728 |
| 45 | 1.2359184 | 1.0329687 | 1.4388680 |
Realizamos las transformaciones para graficar y determinar el modelo lineal de \(ln(k)\) vs \(\frac{1}{T}\).
k_method2long <- k_method2 %>% pivot_longer(-1,
values_to = "k",
names_to = "origen") %>%
mutate(`1/T (1/ºK)` = (273.15+Temperatura)^(-1),
`ln(k)` = log(k))
kable(k_method2long, format = "markdown")| Temperatura | origen | k | 1/T (1/ºK) | ln(k) |
|---|---|---|---|---|
| 20 | k | 0.1972179 | 0.0034112 | -1.6234461 |
| 20 | k_min | 0.1367181 | 0.0034112 | -1.9898345 |
| 20 | k_max | 0.2577177 | 0.0034112 | -1.3558903 |
| 35 | k | 0.5770626 | 0.0032452 | -0.5498045 |
| 35 | k_min | 0.5088524 | 0.0032452 | -0.6755972 |
| 35 | k_max | 0.6452728 | 0.0032452 | -0.4380821 |
| 45 | k | 1.2359184 | 0.0031432 | 0.2118143 |
| 45 | k_min | 1.0329687 | 0.0031432 | 0.0324369 |
| 45 | k_max | 1.4388680 | 0.0031432 | 0.3638567 |
#Gráfico de Sabor oxidado en función de los días de almacenamiento para las tres temperaturas.
ggplot(k_method2long, aes(x = `1/T (1/ºK)`,
y = `ln(k)`)) +
geom_smooth(method = "lm", width = 0.5,
alpha = 0.2) +
geom_point(aes(color = origen),
size = 2) +
geom_label(
x = 100,
y = 3,
label = "ln(A) == ln(A[0]) - k*t",
parse = TRUE,
inherit.aes = FALSE) +
labs(
title = "ln(k) vs 1/T",
caption = "Figura 4: Ajuste de ln(k) vs 1/T en el método de regresión lineal con intervalos.",
x = "ln(k)",
y = "1/T",
color = NULL,
linetype = NULL
) +
theme_minimal(base_size = 14) +
theme(legend.position = "none",
plot.caption = element_text(hjust = 0, size = 9))
lm_Ea2 <- lm(`ln(k)` ~ `1/T (1/ºK)`, data = k_method2long)
summary(lm_Ea2)##
## Call:
## lm(formula = `ln(k)` ~ `1/T (1/ºK)`, data = k_method2long)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.32106 -0.15010 0.02928 0.08386 0.31289
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21.89 2.01 10.89 1.21e-05 ***
## `1/T (1/ºK)` -6906.56 614.93 -11.23 9.90e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2038 on 7 degrees of freedom
## Multiple R-squared: 0.9474, Adjusted R-squared: 0.9399
## F-statistic: 126.1 on 1 and 7 DF, p-value: 9.902e-06Es el modelo lineal significativo (p = 0.025) con pendiente (\(-k\)) de -6809.79. Sin embargo, el intervalo de confianza (95%) de este parámetro es:
as.data.frame(confint(lm_Ea2, level = 0.95)) %>% mutate(`Mean value` = rowMeans(across(1:2))) #Mostrar el promedio también## 2.5 % 97.5 % Mean value
## (Intercept) 17.13855 26.64355 21.89105
## `1/T (1/ºK)` -8360.64459 -5452.48429 -6906.56444Por tanto, en este segundo método, tenemos un valor de \(E_AR^{-1}\) de \((6906.6 \pm 1454.0)ºK^{-1}\). Se tiene un error menor que en el caso anterior.
4.2.3 Regresión no lineal
Este método de cálculo suele ser el más exacto, pero requiere la utilización de un software que permita realizar la correlación no lineal.
El mismo utiliza los resultados de la sección 4.1.1 y la ecuación de Arrhenius explícita (Sección 3).
\[ \begin{array}{rl} \left\{ \begin{array}{l} \text{A}(t) = \text{A}_o + kt \quad \text{(Es orden cero)} \\ k = k_{\text{ref}} \, e^{-\frac{E_A}{R} \left( \frac{1}{T} - \frac{1}{T_{\text{ref}}} \right)} \quad \text{(Arrhenius)} \end{array} \right. & \Rightarrow \quad \text{A}(t) = \text{A}_o + k_{\text{ref}} \, e^{-\frac{E_A}{R} \left( \frac{1}{T} - \frac{1}{T_{\text{ref}}} \right)} t \end{array} \]
En programas estadísticos como R o Genstat, es posible utilizar como input la formula anterior con una serie de valores que lo cumplen. El programa itera sobre los valores y devuelve una salida como en el siguiente ejemplo para determinar las constantes \(E_AR^{-1}\), \(k_{ref}\) y \(A_o\).
#Convertir datos de temperatura a kelvin
df$Temp_K <- df$Temperatura + 273.15
#Fijar Tref en Kelvin
Tref <- 300
#Ajuste no lineal
nls_fit <- nls(
Sabor_oxidado ~ A_o + k_ref * exp(-EA_R * (1 / Temp_K - 1 / Tref)) * Dias,
data = df,
start = list(A_o = -10, k_ref = 0.56, EA_R = 6800) # valores iniciales razonables
)
#ver resultados:
summary(nls_fit)##
## Formula: Sabor_oxidado ~ A_o + k_ref * exp(-EA_R * (1/Temp_K - 1/Tref)) *
## Dias
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## A_o -1.19882 4.03246 -0.297 0.769844
## k_ref 0.21240 0.04321 4.915 0.000131 ***
## EA_R 8735.22606 658.10872 13.273 2.12e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.883 on 17 degrees of freedom
##
## Number of iterations to convergence: 7
## Achieved convergence tolerance: 8.689e-06Los intervalos de confianza de los parámetros determinados son:
as.data.frame(confint(nls_fit, level = 0.95)) %>% mutate(`Mean value` = rowMeans(across(1:2)))## 2.5% 97.5% Mean value
## A_o -9.6042956 7.167514e+00 -1.218391
## k_ref 0.1236815 3.006964e-01 0.212189
## EA_R 7586.4135818 1.041201e+04 8999.209295Por tanto, en este último método tenemos un valor de \(E_AR^{-1}\) de \((8999.2 \pm 1412.8)\). Si bien es un valor distinto a los anteriores, es el más riguroso desde un punto de vista estadístico/matemático y el que presenta un menor error de determinación.
4.2.4 Resumen de estimaciones
La siguiente tabla resume los cálculos obtenidos por los tres métodos propuestos para estimar la \(E_AR^{-1}\). En la misma se observa que el método con el que se obtiene la mejor estimación es el que parte de la regresión no lineal.
| Método | \(E_AR^{-1}\) (1/°K) | Intervalo de confianza (±) |
|---|---|---|
| Lineal (3 puntos) | 6809.8 | ±3397.3 |
| Lineal con intervalos de confianza (9 puntos) | 6906.6 | ±1454.0 |
| No lineal | 8999.2 | ±1412.8 |