Manual de usuario para fraping
Gabriela Itzetl Medina-Ruíz 1,3 , Arturo Itzcóatl Medina-Ruiz 2 and Julio Morán 1 * 1. División de Neurociencias, Instituto de Fisiología Celular, Universidad Nacional Autónoma de México, Coyoacán, 04510, Ciudad de México 2. Facultad de Ciencias, Universidad Nacional Autónoma de México, Coyoacán, 04510, Ciudad de México, México 3. Posgrado en Ciencias Biológicas, Unidad de Posgrado, Edificio D, 1° Piso, Circuito de Posgrados, Ciudad Universitaria, Coyoacán, 04510, Ciudad de México, México
June 12, 2023
Introducción
La recuperación de fluorescencia después del fotoblanqueo (FRAP), una técnica ampliamente utilizada en biología y biomedicina, ofrece una oportunidad única para investigar la movilidad y la dinámica de las moléculas en células vivas. Proporciona información valiosa sobre la recuperación de moléculas marcadas con fluoróforos después de ser sometidas a fotoblanqueo.
Este manual presenta una metodología avanzada para el análisis de la recuperación de fluorescencia en la técnica de FRAP (Fluorescence Recovery After Photobleaching). Su objetivo principal es proporcionar a los investigadores una herramienta completa que les permita comprender y cuantificar la dinámica de las moléculas en células vivas mediante el uso de FRAP.
El contenido de este manual abarca los siguientes puntos clave:
Fundamentos de la técnica de FRAP: Se proporciona una descripción detallada de los principios fundamentales de la técnica de FRAP, incluyendo los conceptos básicos de fotoblanqueo, recuperación de fluorescencia y el uso de fluoróforos en experimentos de FRAP.
Modelos de recuperación de fluorescencia: Se introducen diferentes modelos matemáticos y estadísticos para describir la recuperación de la fluorescencia después del fotoblanqueo. Se explican los fundamentos teóricos de cada modelo y se muestran ejemplos de su aplicación en el análisis de datos de FRAP.
Ajuste de curvas y estimación de parámetros: Se presenta una metodología para ajustar curvas de recuperación de fluorescencia a los datos experimentales obtenidos mediante FRAP. Se discuten las técnicas de optimización y se proporcionan herramientas prácticas para estimar los parámetros del modelo de recuperación.
Análisis de resultados y evaluación de calidad: Se aborda la evaluación de la calidad del ajuste de las curvas de recuperación de fluorescencia mediante métodos estadísticos y visualización de resultados. Se discuten los criterios para determinar la bondad de ajuste y se presentan herramientas gráficas para analizar y comparar los resultados obtenidos.
Simulación de datos y validación del modelo: Se explica el proceso de simulación estocástica de curvas de recuperación de fluorescencia basado en los parámetros estimados. Se discute la importancia de la simulación en el análisis de FRAP y se proporcionan herramientas para validar el modelo de recuperación mediante comparación con datos reales.
En resumen, este manual ofrece una guía completa y práctica para el análisis de la recuperación de fluorescencia en la técnica de FRAP. Los investigadores encontrarán en este manual una referencia útil para comprender los fundamentos teóricos, aplicar modelos estadísticos, ajustar curvas, evaluar resultados y realizar simulaciones que permitan un análisis riguroso y cuantitativo de los datos obtenidos mediante FRAP.
Vamos a comenzar
En la paquetería fraping se implementan
las funciones necesarias para realizar el análisis de datos de FRAP. Sin
embargo, para que esta funcione correctamente, también es necesario
instalar la paquetería itz. Ambos son repositorios
de GitHub, y su instalación se lleva a
cabo mediante la función install_github de la paquetería devtools
de R.
devtools::install_github("artitzco/itz")
devtools::install_github("artitzco/fraping")De esta forma, las librerías quedarán instaladas en su equipo, por lo que cada que inicie un nuevo proyecto para el análisis de datos de FRAP, únicamente tendrá que importar las librerías correspondientes:
library(itz)
library(fraping)Comienza Ahora (opcional)
Si deseas comenzar a utilizar la librería fraping de inmediato,
puedes descargar un conjunto de datos de ejemplo utilizando la función
downloadData. Esta función descargará un archivo con datos
de FRAP en la carpeta “data” dentro de tu directorio de trabajo, que
puedes utilizar como punto de partida para practicar los métodos de
análisis que se describen en este manual. Ten en cuenta que esta
descarga es opcional y puedes utilizar tus propios datos de FRAP si los
tienes disponibles.
Si deseas utilizar los datos de ejemplo, ejecuta el siguiente código
en la consola de comandos de R:
fraping::downloadData(). Después de ejecutar esta función,
además de la carpeta “data”, encontrarás un script de prueba llamado
“test.fraping.R” en la raíz del directorio. Puedes utilizar este script
como punto de partida para explorar y analizar los datos de FRAP.
La función newDataFrap crea una base para un nuevo grupo
de datos de FRAP, la función sintetiza la información y devuelve una
lista que permite al usuario manipular sus datos con mayor facilidad.
Los parámetros que admite son el nombre del grupo y la intensidad del
fotoblanqueo en porcentaje. Para más detalles ir a newDataFrap.
Para evitar futuros problemas, el nombre del grupo debe cumplir con las siguientes recomendaciones:
- No debe incluir espacios,
- Debe inicial con una letra,
- Sólo se admiten los caracteres especiales “.” y “_“.
A <- newDataFrap("Neurites", 80)
B <- newDataFrap("Secondary.Neurites", 80)
C <- newDataFrap("Dendrites", 80)Recuerda seguir estas recomendaciones al crear nuevos grupos de datos de FRAP. Esto asegurará la consistencia y la correcta manipulación de los datos en futuros análisis.
Importar datos de FRAP
Una vez inicializados los grupos de datos, es hora de añadir la
información. Una forma de hacerlo es mediante la función
addDir del objeto newDataFrap. Esta función
requiere que el usuario organice las tablas de FRAP en tres ficheros en
su equipo: el primero con los datos de recuperación, el segundo con las
tablas de control y el tercero con la información del fondo.
La función addDir de la paquetería fraping acepta los
siguientes parámetros para agregar la información de los ficheros de
datos de FRAP:
seconds: Duración del experimento en segundos.photo: Índice del fotograma en el que se aplica el fotoblanqueo.areaColumn: Nombre de la columna que contiene los datos de área.inputColumn: Nombre de la columna que contiene los datos de intensidad.dirR: Ruta del fichero de datos de recuperación.dirC: Ruta del fichero de datos de control.dirB: Ruta del fichero de datos de fondo (opcional).
Es importante tener en cuenta que los datos del fondo son opcionales. Si decides omitir la información del fichero de datos de fondo, debes considerar que esto tendrá un impacto directo en la estandarización de los datos durante el análisis.
A continuación se muestra un ejemplo de cómo utilizar la función
addDir:
A$addDir(480, 4, "Area", "IntDen",
"data/GFP/GFP Secondary Neurite",
"data/GFP/GFP Secondary Neurite Control",
"data/GFP/GFP Secondary Neurite Back")
B$addDir(480, 4, "Area", "IntDen",
"data/LIFEACT GFP/LIFEACT GFP Secondary Neurite",
"data/LIFEACT GFP/LIFEACT GFP Secondary Neurite Control",
"data/LIFEACT GFP/LIFEACT GFP Secondary Neurite Back")
C$addDir(480, 4, "Area", "IntDen",
"data/LIFEACT GFP/LIFEACT GFP Dendrite",
"data/LIFEACT GFP/LIFEACT GFP Dendrite Control",
"data/LIFEACT GFP/LIFEACT GFP Dendrite Back")Al momento de crear un nuevo grupo de datos utilizando la función
newDataFrap, se asigna automáticamente un color aleatorio
al grupo. Este color se utiliza para diferenciarlo de los demás en las
gráficas. Si deseas reasignar un color específico al grupo, puedes
utilizar la función setColor del objeto
newDataFrap. Esta función acepta como parámetro el valor
del color en formato hex, o puedes elegir un color predefinido
del vector grDevices::colors() o utilizar la función color
de la paquetería itz.
A$setColor("#05753D")
B$setColor("#0000A6")
C$setColor("#E014F7")Antes de realizar un ajuste paramétrico, es posible que prefieras
visualizar los datos mediante gráficos para hacer una inspección visual.
Para lograr esto, puedes utilizar la función plotRecover,
que acepta como parámetros las bases creadas con
newDataFrap. Puedes personalizar los gráficos utilizando
diversos parámetros. Para obtener más detalles sobre esta función,
consulta plotRecover.
plotRecover(A, B, C, plot.lines=T)
plotRecover(A, B, C, plot.shadow=T, plot.mean=T)
plotRecover(A, B, C, plot.points=T)
Elegir un modelo de probabilidad
El modelo paramétrico propuesto para describir el comportamiento de la curva de curva de recuperación de la fluorescencia a partir del momento del fotoblanqueo, es el siguiente:
\[F(t)=f_{min}+\alpha\,(1-f_{min})\,p(t-t_0|\Theta),\] donde \(p\) es una función de distribución acumulada parametrizada por \(\Theta\), \(t_0\) es el tiempo del fotoblanqueo y \(f_{min}\) es el valor mínimo de la fluorescencia después del fotoblanqueo, también se define el valor teórico \(f_{max}\) como el valor máximo que puede alcanzar la fluorescencia después del fotoblanqueo, y está dado por: \[f_{max}= f_{min}+\alpha(1-f_{min}).\]
Puede observarse que la función \(F(t)\) está definida en el intervalo de tiempo \(t\in[t_0, \infty)\). De forma análoga, se define la función \(F^{^{AB}}(t)\) para representar a la curva de recuperación de la fluorescencia después del fotoblanqueo la cual está definida en el intervalo de tiempo \(t\in[0, \infty)\):
\[F^{^{AB}}(t)=f_{min}+\alpha\,(1-f_{min})\,p(t|\Theta).\]
En el apéndice de este documento se puede
encontrar más información acerca de la construcción del modelo. Para que
el usuario pueda realizar el ajuste paramétrico, para el posterior
análisis de datos sus datos, tendrá que ocuparse únicamente de elegir el
modelo de probabilidad, \(p\), más
adecuado. Mediante la función newFit es posible declarar un
nuevo modelo de probabilidad, la función recibe como parámetros el
nombre que se le quiera dar al modelo, la función de distribución
acumulada, el nombre de los parámetros de la función de probabilidad y
una lista con los rangos de valores para dichos parámetros. Para más
detalles ir a newFit.
Exp<-newFit("Exponential", pexp,"rate", list(c(0,1)))
Gam<-newFit("Gamma", pgamma, c("shape","rate"), list(c(0,5), c(0,5)))
Wei<-newFit("Weibull", pweibull, c("shape","scale"), list(c(0,5),c(0,2000)))Uno de los criterios más importantes para elegir entre distintos
modelos de probabilidad es mediante el error medio cuadrático, aquellos
modelos que tengan menor error deben ser mayormente preferidos. El error
medio cuadrático está dado por la siguiente expresión: \[E_i =\frac{1}{m_i}\sum_{j=1}^{m_i}
[f_i(t_j)-F(t_j)]^2,\] donde \(\{f_1,..,f_n\}\) es un grupo de datos
fluorescencia y \(\{t_1,t_2,\cdots,t_{m_i}\}\) son todos
aquellos valores mayores al tiempo de fotoblanqueo cuantificados para
\(f_i\). La función
compareFit comparar el error medio cuadrático de diferentes
modelos de probabilidad ajustados a un grupo de datos, por medio de una
gráfica Q-Q y el histograma del error de cada uno de los modelos. La
función recibe como parámetros al grupo de datos en cuestión, una lista
de los modelos que se requieran comparar, así como una serie de
parámetros que permiten personalizar el gráfico. Para más detalles ir a
compareFit.
compareFit(B, fit=l(Exp, Gam, Wei), col.lines=c("red","yellow"), lty.lines=c(2,1), lwd.lines=2, lwd.mean=2)
compareFit(C, fit=l(Exp, Gam, Wei), col.lines=c("red","yellow"), lty.lines=c(2,1), lwd.lines=2, lwd.mean=2)
Para más información acerca de la elección de un buen modelo de probabilidad ir al apéndice.
Ajustar un modelo paramétrico
Para este caso, en particular, se tiene que el modelo Weibull arroja
mejores resultados, por lo que se recurre a la función Wei,
para realizar el ajuste correspondiente de los datos, la función
Wei, aplicada en un grupo de datos específico, devuelve una
lista con la información obtenida por el ajuste paramétrico, así como la
estimación de la función de recuperación de la fluorescencia \(F(t)\), entre otros datos relevantes, para
más detalles ir a newFit.
Para realizar el gráfico del ajuste paramétrico de la recuperación de
la fluorescencia se utiliza la función plotFit la cual
recibe como parámetros, en primer lugar, a los grupos de datos por
ajustar y posteriormente el modelo de probabilidad que se haya elegido
para cada grupo de datos, entre parámetros opcionales que permiten
personalizar y dar formato a cada gráfica, ir a plotFit para saber más. Si el modelo de probabilidad
elegido es el mismo para todos los grupos de datos, basta con
especificarlo una única vez como se ve a continuación:
plotFit(B, C, fit=Wei, plot.lines=T)
plotFit(B, C, fit=Wei, plot.shadow=T, plot.mean=T, alp.shadow=0.2, ylim=c(0.2, 1), xdigits=0)
Con la función plotFit se pueden graficar estimaciones
futuras, de la curva de recuperación de la fluorescencia, con solo
modificar el rango del gráfico. Las gráficas obtenidas por las funciones
plotRecover y plotFit pueden ser apiladas
cambiando el parámetro opcional new.plot a
FALSE.
plotFit(B, fit=Wei, plot.shadow=T, alp.shadow=0.3, xlim=c(0,500), ylim=c(0.2, 1))
plotRecover(B, AB=T, plot.lines=T, new.plot=F)
Otro uso para la función plotFit es hacer una
comparación visual entre dos o más modelos de probabilidad sobre un
mismo grupo de datos.
plotRecover(B, plot.shadow =T, ylim=c(0.2, 1.2), xdigits=0,lwd.border=1, col="#808080")
plotFit(B, B, fit=l(Wei, Exp), plot.lines=T, new.plot=F, displacement=T, col=c("blue2","red2"))
Comparar conjuntos de datos
Análisis convencional
Por medio de la función compareParam es posible comparar
dos grupos de datos por medio de las variables asociadas al ajuste del
modelo paramétrico. Los parámetros que recibe esta función son los dos
grupos de datos, el modelo de probabilidad que ajusta a cada grupo de
datos, la variable que se desea comparar, entre otros valores que
permiten personalizar el gráfico. para saber más ir a compareParam. Los nombres de las variables
pueden obtenerse de la siguiente forma:
colnames(Wei(B)$table)## [1] "shape" "scale" "alpha" "fmax" "UF" "MF" "error"La función compareParam realiza un gráfico Q-Q de la
variable elegida, cambiando el parámetro new.plot=FALSE es
posible apilar los gráficos. Los detalles de la prueba \(t\mbox{-}student\) para la diferencia de
medias al \(95\%\) de confianza se
obtienen cambiando el parámetro return=TRUE (en
consecuencia, se omitirá el gráfico), el nivel de confianza y la
hipótesis alternativa pueden se modificados mediante los parámetros
alternative y conf.level respectivamente, para
más detalles ir a compareParam.
compareParam(B, C, fit=Wei, param="shape")
De esta forma es posible realizar un análisis convencional de la recuperación de la fluorescencia, es decir, comparar la Fracción Inmóvil (UF) y la Fracción Móvil (MF). Recuerde que \[UF=\frac{1-f_{max}}{1-f_{min}}\] \[MF=\frac{f_{max}-f_{min}}{1-f_{min}}\] \[UF+MF=1\]
compareParam(B, C, fit=Wei, param="UF", xlim=c(0,1), ylim=c(0,1),lty.lines=c(2,1), col.lines=c("white","red"), lwd.lines=2)
compareParam(B, C, fit=Wei, param="MF", new.plot=F, lwd.lines=2)
compareParam(B, C, fit=Wei, param="UF", return=T)#>>><<<##
## Welch Two Sample t-test
##
## data: UF_Weibull_Secondary.Neurites and UF_Weibull_Dendrites
## t = -0.51751, df = 4.2562, p-value = 0.6305
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.2846008 0.1933821
## sample estimates:
## mean of x mean of y
## 0.08361768 0.12922699compareParam(B, C, fit=Wei, param="MF", return=T)#>>><<<##
## Welch Two Sample t-test
##
## data: MF_Weibull_Secondary.Neurites and MF_Weibull_Dendrites
## t = 0.51751, df = 4.2562, p-value = 0.6305
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.1933821 0.2846008
## sample estimates:
## mean of x mean of y
## 0.9163823 0.8707730También, con la función compareParam es posible comparar
dos modelos de probabilidad, ajustados a un mismo grupo de datos, con
respecto al error que produce cada ajuste, puede notar que esta gráfica
ya se encuentra incluida dentro del mosaico que arroja la función
compareFit vista anteriormente.
compareParam(B, B, fit=l(Exp,Wei), param="error", col.lines=c("red", "yellow"), lwd.lines=2, lty.lines=c(2,1), ydigits=2, xdigits=2)
compareParam(B, B, fit=l(Exp,Wei), param="error", return=T)#>>><<<##
## Welch Two Sample t-test
##
## data: error_Exponential_Secondary.Neurites and error_Weibull_Secondary.Neurites
## t = 1.2307, df = 21.424, p-value = 0.2318
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.09746458 0.38089484
## sample estimates:
## mean of x mean of y
## 0.7222331 0.5805180Análisis del comportamiento medio
El análisis del comportamiento medio implica evaluar la recuperación promedio de la fluorescencia a lo largo del tiempo. La curva media se calcula de la siguiente manera:
\[\overline F(t)=\frac{1}{n}\sum_{i=1}^{n}F^{^{AB}}_i(t),\]
donde \(n\) es el tamaño del conjunto de datos y \(F^{AB}_i\) representa la curva de recuperación de la fluorescencia después del fotoblanqueo asociada a la cuantificación de la fluorescencia \(f_i\).
plotFit(B, fit=Wei, plot.lines=T, plot.mean=T, col.mean="purple",lwd.mean=2, xdigits=0, ydigits=2, alp.lines=0.5)
La función compareMean permite contrastar las curvas
medias de dos grupos de datos por medio de una prueba \(t\mbox{-}student\) para la diferencia de
medias aplicada de manera continua en un periodo de tiempo. La función
recibe como parámetros a los grupos de datos por comparar, el modelo
probabilidad que ajusta a los datos, entre otros valores que permiten
personalizar la gráfica.
El gráfico que devuelve representa al \(p\mbox{-}value\) de la prueba \(t\mbox{-}student\), la línea de referencia
indica el nivel de significancia. Cambiando el parámetro
return=TRUE se obtienen los detalles de la prueba
almacenados en una lista (omitiendo la realización de la gráfica), el
nivel de confianza y la hipótesis alternativa pueden ser modificados
mediante los parámetros alternative y
conf.level respectivamente, para más detalles ir a compareMean.
compareMean(B, C, fit=Wei, lty.lines=c(2,1), col.lines=c("red","purple"), ydigits=2)
Cambiando el parámetro p.value=FALSE se calculará la
diferencia de medias entre los dos grupos de datos, la línea de
referencia estará situada en cero. En este caso al cambiar el parámetro
return=TRUE se obtiene un vector con la diferencia de
medias entre ambas curvas.
compareMean(B, C, fit=Wei, p.value=F, lty.lines=c(2,1), col.lines=c("red","magenta"), ydigits=3)
Simulación
La implementación de métodos de simulación para el análisis de la recuperación de la fluorescencia es una de las principales aportaciones de este trabajo. La simulación permite incrementar el tamaño de una muestra sin alterar sus medidas de tendencia centra ni dispersión, esto es especialmente útil para los casos en los que no se cuenta con un tamaño de muestra suficientemente grande, recuerde que las pruebas de hipótesis que dependen de procesos límite son particularmente sensibles al tamaño de la muestra. Por otro lado, el proceso de simulación se retroalimenta con una elección adecuada del modelo de probabilidad, esto se explica más adelante.
La simulación del modelo de recuperación de la fluorescencia \(F(t)\) se logra gracias a la simulación de
los parámetros \(\Theta\) del modelo de
probabilidad \(p(t|\Theta)\), ir al apéndice para saber más. La función simFit
permite simular los datos ajustados a un modelo de probabilidad, recibe
como parámetro obligatorio a un grupo de datos ajustados a un modelo de
probabilidad y como parámetros opcionales el tamaño de la muestra
simulada y una semilla de simulación, para más detalles sobre la semilla
de simulación ir al apéndice.
Para el caso del ajuste Weibull la simulación de los datos se realiza
de la siguiente forma simFit(Wei(B)) y
simFit(Wei(C)), para más especificaciones sobre los datos
que arroja esta función ir simFit. Por otro lado,
las funciones plotFit, compareParam y
compareMean tienen implementado el método de simulación,
para acceder a esta función es necesario cambiar el parámetro
simulated=FALSE. Podemos replicar las gráficas hechas
anteriormente ahora con datos simulados:
plotFit(B, C, fit=Wei, simulated=T, plot.lines=T)
plotFit(B, C, fit=Wei, simulated=T, plot.shadow=T, plot.mean=T, alp.shadow=0.2, ylim=c(0.2, 1), xdigits=0)
plotFit(B, fit=Wei, simulated=T, plot.shadow=T, alp.shadow=0.3, xlim=c(0,500), ylim=c(0.2, 1))
plotRecover(B, AB=T, plot.lines=T, new.plot=F)
plotRecover(B, AB=T, plot.shadow=T,ylim=c(0.2, 1.2), xdigits=0,lwd.border=1, col="#808080")
plotFit(B, B, fit=l(Wei, Exp), simulated=T, plot.lines=T, new.plot=F, displacement=T, col=c("blue2","red2"))
compareParam(B, C, fit=Wei, simulated=T, param="shape")
compareParam(B, C, fit=Wei, param="UF", simulated=T, xlim=c(0,1), ylim=c(0,1),lty.lines=c(2,1), col.lines=c("white","red"), lwd.lines=2)
compareParam(B, C, fit=Wei, param="MF", simulated=T, new.plot=F, lwd.lines=2)
compareParam(B, C, fit=Wei, simulated=T, param="UF", return=T)##
## Welch Two Sample t-test
##
## data: UF_WeibullSim_Secondary.Neurites and UF_WeibullSim_Dendrites
## t = -1.4645, df = 97.842, p-value = 0.1463
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.0891179 0.0134351
## sample estimates:
## mean of x mean of y
## 0.07055279 0.10839419compareParam(B, C, fit=Wei, simulated=T, param="MF", return=T)##
## Welch Two Sample t-test
##
## data: MF_WeibullSim_Secondary.Neurites and MF_WeibullSim_Dendrites
## t = 1.4645, df = 97.842, p-value = 0.1463
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.0134351 0.0891179
## sample estimates:
## mean of x mean of y
## 0.9294472 0.8916058plotFit(B, fit=Wei, simulated=T, plot.lines=T, plot.mean=T, col.mean="purple",lwd.mean=2, xdigits=0, ydigits=2, alp.lines=0.5)
compareMean(B, C, fit=Wei, simulated=T, lty.lines=c(2,1), col.lines=c("red","purple"), ydigits=2)
compareMean(B, C, fit=Wei, simulated=T, p.value=F, lty.lines=c(2,1), col.lines=c("red","magenta"), ydigits=3)
Apéndice
El modelo de recuperación de fluorescencia
La cuantificación de la fluorescencia en el proceso de la técnica de FRAP sigue un patrón específico de recuperación. Inicialmente, la fluorescencia se mantiene constante, pero al momento del fotoblanqueo (también conocido como \(t_0\)), la intensidad de la fluorescencia disminuye hasta alcanzar un valor mínimo denominado \(f_{min}\). A partir de ese punto, se inicia un proceso gradual de recuperación de la fluorescencia que se caracteriza por dos condiciones fundamentales: en primer lugar, la recuperación no disminuye a medida que transcurre el tiempo, y en segundo lugar, la fluorescencia alcanza un límite máximo (\(f_{max}\)) que no puede superar el valor inicial antes del fotoblanqueo. Esta evolución se representa de manera gráfica para una mejor comprensión del proceso.
En ese sentido, el modelo propuesto para la recuperación de la fluorescencia debe cumplir con las características anteriores, el modelo general para la recuperación de la fluorescencia después del fotoblanqueo es el siguiente:
\[F(t)=f_{min}+\alpha\,(1-f_{min})\,p(t-t_0|\Theta)\]
En donde \(p\) es una función de distribución acumulada continua (CDF), parametrizada por \(\Theta\).
Observaciones-
La función \(F(t)\) está definida en el intervalo \([t_0,\infty)\)
-
El valor \(f_{max}\) está dado por \[f_{max} =\lim_{t\to\infty}F(t)\]
-
El valor \(\alpha\) está acotado por los valores \[\alpha\in[0,1]\]
-
La función de distribución acumulada \(p\) debe estar asociada a una variable aleatoria que con soporte en \((0,\infty)\), recuerde que la CDF se define de la siguiente manera:
\[p(t|\Theta)=\int_0^td(x|\Theta)\,\mbox{d}x\] donde \(d\) es una función de densidad paramétrica.
De forma análoga, se define la función \(F^{^{AB}}(t)\) para representar a la curva de recuperación de la fluorescencia después del fotoblanqueo la cual está definida en el intervalo de tiempo \(t\in[0, \infty)\):
\[F^{^{AB}}(t)=f_{min}+\alpha\,(1-f_{min})\,p(t|\Theta).\]
Probemos que \(F\) cumple con las características de la recuperación de la fluorescencia después del fotoblanqueo:
- \(F(t_0)=f_{min}\)
- \(F(t)\) es no decreciente
- \(f_{max}=f_{min}+\alpha\,(1-f_{min})\) and \(f_{max}\in[f_{min}, 1]\)
Se evalua
\[\begin{aligned} F(t_0) &=f_{min}+\alpha\,(1-f_{min})\,p(t_0-t_0|\Theta)\\ &=f_{min}+\alpha\,(1-f_{min})\,p(0|\Theta)\\ &=f_{min}+\alpha\,(1-f_{min})\int_0^0d(x|\Theta)\,\mbox{d}x\\ &=f_{min}+\alpha\,(1-f_{min})\,(0)\\ &=f_{min} \end{aligned}\] \[\therefore F(t_0)=f_{min}\]
Considere a \(t_1,t_2\in\mathbb R^+\) tales que \(t_1\leq t_2\), entonces
\[\begin{aligned} t_1-t_0\leq &\;t_2-t_0\\ \Rightarrow p(t_1-t_0|\Theta)\leq &\;p(t_2-t_0|\Theta)\\ f_{min}+\alpha\,(1-f_{min})\,p(t_1-t_0|\Theta)\leq &\;f_{min}+\alpha\,(1-f_{min})\,p(t_2-t_0|\Theta) \end{aligned}\] \[\therefore F(t_1)\leq F(t_2)\]
Lo anterior prueba que el modelo paramétrico \(F\) es no decreciente a lo largo del tiempo.
Se calcula el límite de \(F(t)\) cuando \(t\) tiende a infinito:
\[\begin{aligned} f_{max}=\lim_{t\to\infty}F(t) &=\lim_{t\to\infty}f_{min}+\alpha\,(1-f_{min})\,p(t-t_0|\Theta)\\ &=f_{min}+\alpha\,(1-f_{min})\lim_{t\to\infty}\,p(t-t_0|\Theta)\\ &=f_{min}+\alpha\,(1-f_{min})\lim_{t\to\infty}\int_0^{t-t_0}d(x|\Theta)\,\mbox{d}x\\ &=f_{min}+\alpha\,(1-f_{min})\int_0^{\infty}d(x|\Theta)\,\mbox{d}x\\ &= f_{min}+\alpha\,(1-f_{min})\,(1) \end{aligned}\]
\[\therefore f_{max}=f_{min}+\alpha\,(1-f_{min})\] Por otro lado se tiene que \(\alpha\in[0,1]\), entonces
\[\begin{aligned} &0\leq \alpha\leq 1\\ \Rightarrow &0\leq \alpha\,(1-f_{min})\leq 1-f_{min}\\ \Rightarrow &f_{min}\leq f_{min}+\alpha\,(1-f_{min})\leq f_{min}+1-f_{min}\\ \Rightarrow &f_{min}\leq f_{min}+\alpha\,(1-f_{min})\leq 1\\ \Rightarrow &f_{min}\leq f_{max}\leq 1\\ \end{aligned}\] \[\therefore f_{max}\in[f_{min},1]\]Dados un conjunto de datos de fluorescencia \(f\) y una función de distribución acumulada \(p(\bullet|\Theta)\), entonces los parámetros \((\alpha,\Theta)\) son estimados por medio del metodo de mínimos cuadrados, que consiste en minimizar el error medio cuadrático producido entre las curvas \(f\) y \(F\), la función del error medio cuadrático es la siguiente:
\[\begin{aligned}E(\alpha,\Theta)&=\sum_{j=1}^m\left[f(t_j)-F(t_j)\right]^2\\&=\sum_{j=1}^m\left[f(t_j)-f_{min}-\alpha\,(1-f_{min})\,p(t_j-t_0|\Theta)\right]^2\end{aligned}\] donde \(\{t_1,t_2,\cdots,t_m\}\) son todos aquellos tiempos mayores a \(t_0\) cuantificados para \(f\).
Elegir un modelo de probabilidad
La función de distribución acumulada \(p(t|\Theta)\) es la base para construir el
modelo paramétrico de recuperación de la fluorescencia \(F(t|\Theta)\). La elección de la función
\(p\) queda a libre criterio del
usuario, sin embargo, es indispensable que \(p\) tenga soporte en \((0,\infty)\), otro aspecto del que debe
encargarse el usuario es la asignación de los intervalos donde toman
valores los parámetros \(\Theta\), esto
se logra por medio de la función newFit.
Ejemplo:
Suponga que el usuario ha decidido utilizar una función de probabilidad Pareto, entonces \[p(t|\,\mbox{shape},\mbox{scale})=1-\left(\frac{\mbox{scale}}{t+\mbox{scale}}\right)^{\mbox{shape}}\] \[t>0,\;\;\mbox{scale}>0\;\;\&\;\;\mbox{shape}>0\] Se implementa la función de distribución acumulada Pareto y se declara su ajuste paramétrico asociado:
ppareto <- function(t, shape, scale){
ppar <- NULL
for (t in t)
ppar <- c(ppar, 1-(scale/(t+scale))^shape)
ppar
}
Par <- newFit("Pareto", ppareto, c("shape","scale"), list(c(0,a), c(0,b)))El código anterior establece que \(\mbox{shape}\in(0,a)\) y \(\mbox{scale}\in(0,b)\). Después de ajustar
el modelo Par a un grupo de datos el usuario debe verificar
que los valores shape y scale se encuentren
completamente contenidos en el interior de los intervalos \((0,a)\) y \((0,b)\) respectivamente, de lo contrario,
tendrá que modificar los valores de \(a\) y/o \(b\) y repetir el procedimiento. Lo anterior
se consigue mediante una inspección visual de las variables devueltas
por el ajuste: se deben ejecutar los códigos
Par(B)$table$shape y Par(B)$table$scale e
inspeccionar los resultados, donde B es un conjunto de
datos de FRAP obtenido por medio de la función
newDataFrap.
Las funciones creadas por newFit tienen el propósito de
encontrar los mejores valores \(\Theta\) tales que minimicen el error
cuadrático entre el modelo \(F(t|\Theta)\) y los datos reales de
recuperación de la fluorescencia. Por lo que el principal criterio para
elegir el modelo de probabilidad \(p\)
más adecuado debe ser: preferir el modelo que produzca el menor error
cuadrático medio, esto se logra, como se vio anterior mente, por medio
de la función compareFit.
compareFit(B, fit=l(Exp, Gam, Wei, Par), col.lines=c("red","yellow"), lty.lines=c(2,1), lwd.lines=2, lwd.mean=2)
Por otro lado, un segundo criterio para evaluar la viabilidad de un modelo de probabilidad \(p\) es viendo cómo se comportan las simulaciones de \(F\). Si existen discrepancias entre los datos reales de fluorescencia y las simulaciones de \(F\) entonces lo más razonable es descartar el modelo \(p\). Ya que esto implicaría que no hay una relación causal entre los parámetros \(\Theta\) y los datos de recuperación de la fluorescencia y, como se verá en el siguiente apéndice, las simulaciones de \(f_{max}\) son particularmente sensibles a las variaciones de \(\Theta\).
plotFit(B, B, fit=l(NULL,Exp), col=c("#916F91","#FF4D66"), plot.lines=c(F,T), plot.mean=c(T,F), lwd.mean=2, plot.shadow=c(T,F), lwd.border=1, alp.border=0.4, alp.shadow=0.3, ylim=c(0.2, 1), xdigits=0, simulated=T, seed=1992, main="Exponential simulation", cex.main=1.5)
plotFit(B, B, fit=l(NULL,Gam), col=c("#916F91","#4DB34D"), plot.lines=c(F,T), plot.mean=c(T,F), lwd.mean=2, plot.shadow=c(T,F), lwd.border=1, alp.border=0.4, alp.shadow=0.3, ylim=c(0.2, 1), xdigits=0, simulated=T, seed=1992, main="Gamma simulation", cex.main=1.5)
plotFit(B, B, fit=l(NULL,Wei), col=c("#916F91","#664DFF"), plot.lines=c(F,T), plot.mean=c(T,F), lwd.mean=2, plot.shadow=c(T,F), lwd.border=1, alp.border=0.4, alp.shadow=0.3, ylim=c(0.2, 1), xdigits=0, simulated=T, seed=1992, main="Weibull simulation", cex.main=1.5)
plotFit(B, B, fit=l(NULL,Par), col=c("#916F91","#B34DB3"), plot.lines=c(F,T), plot.mean=c(T,F), lwd.mean=2, plot.shadow=c(T,F), lwd.border=1, alp.border=0.4, alp.shadow=0.3, ylim=c(0.2, 1), xdigits=0, simulated=T, seed=1992, main="Pareto simulation", cex.main=1.5)
Teoría de la simulación estocástica
La simulación estocástica es un proceso mediante el cual se pretende reproducir el comportamiento de una variable aleatoria, la base de la simulación se sigue del siguiente teorema:
Sea \(X\) es una variable aleatoria y \(\mathbb F_X\) su función de distribución acumulada, sea \(U\) una variable aleatoria con distribución \(\mbox{Uniform}(0,1)\), entonces la variable aleatoria \(\mathbb F^{-1}_X(U)\) tiene la misma distribución que \(X\), es decir:
\[\mathbb F_X = \mathbb F_{\mathbb F^{-1}_X(U)}\]
Proof:
\[\begin{aligned}\mathbb F_{\mathbb F^{-1}_X(U)} &=\mathbb P\left[\mathbb F^{-1}_X(U)\leq t\right] =\mathbb P\left[\mathbb F_X\left(\mathbb F^{-1}_X(U)\right)\leq \mathbb F_X\left(t \right)\right]\\ &=\mathbb P\left[U\leq \mathbb F_X\left(t \right)\right] =\int_0^{\mathbb F_X\left(t \right)}f_U(u)\,\mbox{d}u =\int_0^{\mathbb F_X\left(t \right)}1\,\mbox{d}u\\ &=\left.u\right|_0^{\mathbb F_X\left(t \right)}={\mathbb F_X\left(t \right)}-0={\mathbb F_X\left(t \right)} \end{aligned}\] \[\therefore {\mathbb F_X\left(t \right)}=\mathbb F_{\mathbb F^{-1}_X(U)}\]
Nótese que como \(U\) es una variable uniforme, se tiene que \(f_U(u)=1,\;\; 0<u<1\). (\(f_U\) es la función de densidad de \(U\))
El teorema anterior permite simular cualquier variable aleatoria \(X\) siempre que se conozca su función de distribución acumulada \(F_X\). Sea \(\{u_1,u_2,\cdots,u_n\}\) una muestra aleatoria de valores observados de la variable aleatoria \(\mbox{Uniform}(0,1)\), entonces \(\{\mathbb F_X^{-1}(u_1),\mathbb F_X^{-1}(u_2),\cdots,\mathbb F_X^{-1}(u_n)\}\) es una muestra de valores simulados de la variable \(X\). En términos computacionales, para realizar el proceso de simulación de simulación, la generación de una muestra aleatoria uniforme se lleva a cabo por medio de algoritmos de generación de números pseudoaleatorios.
Habrá situaciones en las que se requiera simular el comportamiento de una muestra de datos \(x_1,x_2,\cdots,x_n\) de alguna distribución desconocida, por lo que no se dispone de función de distribución acumulada, en estos caso para realizar la simulación se utiliza a la función de distribución acumulada empírica (ecdf), la cual se define de la siguiente manera:
\[\widehat{\mathbb F}_{x_1,x_2,\cdots,x_n}(t)=\frac{1}{n}\sum_{i=1}^n\mathbb {I}_{x_i\leq t}\]
donde \(\mathbb{I}_{x_i\leq t}=\left\lbrace\begin{aligned} 1,\; x_i\leq t\\ 0,\; x_i> t\end{aligned}\right.\)
