Los sitios de metilación del ADN fluctuantes codifican el historial de crecimiento del tumor colorrectal.

Los tumores sólidos humanos se observan típicamente solo una vez, tras la resección. Por lo tanto, los parámetros de crecimiento tumoral, la edad y los modos de evolución deben inferirse a partir de una única instantánea de una estructura que ha estado creciendo durante meses o años. Aunque las líneas clonales pueden reconstruirse a partir de mutaciones genéticas, los errores en la metilación del ADN somático (marcas hereditarias que se acumulan durante la división celular [^[1], ^[2]]) ofrecen potencialmente una mayor resolución temporal a partir de una forma de datos ampliamente disponible y de bajo costo. Dado que los patrones de metilación se copian con pequeñas pero medibles tasas de error en cada replicación, actúan como relojes moleculares endógenos. Los primeros estudios explotaron esta propiedad para inferir el origen y la dinámica de crecimiento de los tumores en tumores colorrectales utilizando pequeños paneles de etiquetas CpG (citosina seguida de guanina) [^[3]–^[5]].

Más recientemente, Gabbutt et al. identificaron una clase de loci CpG fluctuantes (fCpG) cuyo estado de metilación no está constitutivamente metilado ni no metilado, sino que fluctúa en escalas de tiempo clínicamente relevantes debido a errores continuos de metilación y desmetilación [^[6]]. Estos loci fCpG sirven como trazadores de linaje de alta resolución cuyo divergencia codifica la historia reciente de la división celular del tejido muestreado [^[6], ^[7]]. El posterior desarrollo del método EVOFLUx demostró que la inferencia basada en fCpG de la dinámica evolutiva (incluida la tasa de crecimiento tumoral, la edad y las tasas de epimutación) es factible a escala clínica a partir de perfiles de metilación masiva en cánceres linfoides [^[8]].

El cáncer colorrectal (el tercer cáncer más común a nivel mundial y una de las principales causas de mortalidad relacionada con el cáncer [^[9]]) proporciona un sistema ideal para adaptar los métodos de inferencia fCpG a los tumores sólidos. Los tumores colorrectales crecen por fisión de glándulas (criptas), el contraparte neoplásico de la renovación normal de las criptas [^[10]]. Cada glándula tumoral se mantiene mediante una pequeña población de células madre cancerosas (CSCs) con capacidad proliferativa indefinida [^[11], ^[12]], análogas al compartimento de células madre en las criptas colónicas normales [^[13]]. Por lo tanto, las glándulas de tumores colorrectales conservan una organización jerárquica en la que unas pocas células madre generan la mayor parte de la población de la glándula, y el crecimiento tumoral a nivel tisular está impulsado por la fisión de estas unidades de glándula. Los análisis tanto de la diversidad genética como epigenética sugieren que la mayoría de la heterogeneidad genética intratumoral surge al principio de la expansión clonal final, y que la evolución interglandular posterior es efectivamente neutra [^[5], ^[14], ^[15]]. Según este modelo de "Big Bang", la historia de crecimiento tumoral puede inferirse a partir de datos espaciales de un solo punto en el tiempo, como los obtenidos de glándulas muestreadas de ubicaciones distantes en un tumor resecado. Un creciente cuerpo de investigación en biología computacional utiliza modelos basados en agentes e inferencia sin probabilidad en la evolución de tumores sólidos en general [^[16]–^[20]], y en el cáncer colorrectal en particular [^[3], ^[21]–^[24]]. Sin embargo, no existen métodos para explotar la rica información contenida en los datos fCpG estructurados espacialmente.

Aquí presentamos methdemon, un modelo basado en agentes del crecimiento de tumores colorrectales y la dinámica de la matriz de metilación fCpG a nivel de glándula, junto con methabc, un flujo de trabajo de inferencia ABC-SMC que estima los parámetros de crecimiento tumoral, principalmente las tasas de fisión de glándulas y las tasas de epimutación, a partir de matrices de metilación masiva de múltiples regiones. Aplicamos el marco para inferir los parámetros biológicos de 10 tumores colorrectales.

Resultados

Cohorte de pacientes y datos de metilación de CpG fluctuantes

Se perfilaron matrices de metilación de muestras de 10 tumores (7 carcinomas, 3 adenomas, que abarcan un rango de tamaños y estadios clínicos; Tabla 1) utilizando la plataforma Illumina Infinium Methylation EPIC BeadChip, lo que produjo mediciones de valores beta en aproximadamente 850 000 sitios CpG por muestra de glándula con muy alta pureza (Métodos). Se muestrearon ocho glándulas de dos regiones espacialmente separadas de cada tumor resecado (lados A y B), lo que proporcionó un conjunto de datos de 80 perfiles de metilación masiva a nivel de glándula con una estructura espacial conocida (Fig. 1).

Las ganancias y pérdidas estocásticas de metilación en sitios CpG individuales durante la división celular somática producen una variación informativa de linaje hereditaria en el estado de metilación del ADN [^[6]]. En los sitios CpG que no están constitutivamente metilados ni no metilados en toda la población, la fracción de metilación de cada glándula (la matriz fCpG) actúa como un código de barras de linaje intrínseco a la célula: las células genéticamente idénticas se separan en el estado de metilación a través de rondas repetidas de división, y la diversidad de la matriz resultante codifica el número de divisiones celulares desde el ancestro común más reciente. Estos loci CpG fluctuantes (fCpG) se han utilizado para reconstruir historias de linaje en tejidos normales y en cáncer [^[6]–^[8]].

Para identificar los loci fCpG, seguimos un procedimiento similar al de Gabbutt et al. [^[8]], utilizando datos de metilación de 138 muestras de tumores colorrectales en The Cancer Genome Atlas (TCGA) como una cohorte de referencia (Métodos). Este procedimiento identificó 1165 loci fCpG cuyos estados de metilación contienen información sobre la historia de las divisiones celulares en lugar de la identidad del tipo de célula. En toda la cohorte, las matrices fCpG por glándula mostraron la distribución característica en forma de W de los valores beta del sitio (Fig. 2A, 2B), lo que es consistente con las tasas de epimutación en el régimen intermedio donde los estados de metilación ancestrales se sobrescriben parcial pero no completamente.

La separación espacial predice la divergencia de la matriz fCpG en muestras de cáncer colorrectal de múltiples regiones

Para cuantificar la divergencia de metilación interglandular, calculamos la matriz de distancia interglandular D para cada tumor, definida como la diferencia cuadrática media en la fracción de metilación fCpG entre todos los pares de glándulas muestreadas (Ecuación 3; Métodos). Bajo la hipótesis de que las glándulas del mismo lado del tumor divergieron de un ancestro común más recientemente que las glándulas de lados opuestos, las distancias intra-lado deberían ser sistemáticamente menores que las distancias inter-lado.

Esta predicción se confirmó en toda la cohorte. La matriz de distancia interglandular para cada tumor mostró una clara estructura de bloques: los pares de glándulas del mismo lado son consistentemente más similares que los pares de lados opuestos (Fig. 2C, S1). En los 10 tumores, la distancia par por par intra-lado mediana fue significativamente menor que la distancia par por par inter-lado mediana (Fig. 3). Cada tumor exhibió el patrón esperado, con relaciones entre distancias inter-lado e intra-lado que oscilaron entre 1,12 (tumor D) y 4,71 (tumor I). Los tumores con relaciones más altas (I, U, M) tendieron a mostrar una estructura de bloques jerárquicos más clara en sus matrices de distancia, lo que es consistente con una mayor brecha entre los tiempos de coalescencia intra-lado y el tiempo de coalescencia de todo el tumor. La separación más débil observada en el tumor D (relación 1,12) puede reflejar su configuración de muestreo inusual (2 glándulas del lado A, 6 del lado B), lo que impide la estimación precisa de la distancia intra-lado en el lado minoritario.

La magnitud de las distancias interglandulares varió sustancialmente entre los pacientes (Fig. S1), lo que sugiere heterogeneidad en la historia de crecimiento o las tasas de epimutación de tumores individuales. Razonamos que esta variación codifica información sobre los parámetros evolutivos específicos del tumor que pueden recuperarse mediante un marco de inferencia calibrado, lo que motivó el desarrollo del modelo descrito en la siguiente sección.

Un modelo basado en agentes calibrado recapitula la dinámica fCpG multiglandular

Para inferir los parámetros de crecimiento específicos del tumor a partir de los datos fCpG observados, desarrollamos methdemon, un modelo basado en agentes del crecimiento de tumores colorrectales y la dinámica de la metilación fCpG a nivel de glándula. En este modelo, cada glándula se representa como un deme: una subpoblación local y bien mezclada que crece y se reproduce como una unidad. El modelo simula la expansión tumoral como un proceso de ramificación de fisión de glándulas, con cada glándula modelada como un deme de células madre cancerosas que experimentan metilación y desmetilación estocásticas en los loci fCpG (Fig. 4A). Antes de realizar la inferencia, caracterizamos cómo cada parámetro del modelo da forma a las estadísticas resumidas para determinar cuáles son identificables a partir de los datos de la matriz fCpG. Ejecutamos simulaciones de methdemon para una muestra de hipercubo latino de 200 conjuntos de parámetros extraídos de nuestra distribución previa de inferencia (Tabla 2), con todos los demás ajustes coincidentes con la configuración de inferencia (8 glándulas muestreadas, 2 lados, coincidentes con el tamaño promedio del tumor de la cohorte).

El barrido de parámetros reveló cómo las tasas de epimutación y (por división celular) controlan conjuntamente la forma y el sesgo de la distribución fCpG por glándula (Fig. 4B). A bajas tasas de epimutación combinadas, los sitios fCpG retienen su estado fundador y la distribución es marcadamente bimodal, de modo que todos los valores están cerca de 0 o 1. A altas tasas, la sobrescritura repetida del estado fundador impulsa la distribución hacia un pico unimodal cerca de la relación . A tasas intermedias, surge la característica forma de W, en la que la mayoría de los sitios permanecen cerca de 0 o 1, pero un subconjunto se ha desviado a valores intermedios.

La tasa de fisión controla la tasa por glándula a la que se producen nuevas glándulas durante la expansión y, por lo tanto, el número de divisiones celulares que separan las glándulas muestreadas. Las tasas más altas comprimen la expansión tumoral en menos divisiones celulares y reducen la divergencia fCpG interglandular, lo que da forma directamente a la magnitud y la estructura de bloques de la matriz de distancia interglandular (Fig. 4C). Debido a que la tasa de fisión por glándula entra en la dinámica solo a través del producto , la tasa de fisión y la capacidad de carga son estructuralmente no identificables: la reescalación de por un factor se compensa exactamente mediante la reescalación de por un factor . Confirmamos esta degeneración empíricamente con un análisis de sensibilidad. Las matrices de distancia a y difieren por una norma de Frobenius comparable en magnitud a la variación de otros parámetros dentro de la distribución previa. Surge un régimen de tamaño finito en , donde la capacidad del deme se acerca al número de linajes fundadores y aparece una deriva adicional. Por lo tanto, fijamos para la inferencia e interpretamos , en lugar de solo, como la cantidad biológicamente significativa.

Los análisis de sensibilidad establecieron que, si bien las tasas de epimutación y fisión son identificables a partir de la matriz de distancia interglandular (Métodos; Fig. S2), la tasa de mutación del gen impulsor y la ventaja selectiva dejan solo firmas débiles y confundidas. Debido a que las mutaciones del gen impulsor confieren una ventaja de fisión multiplicativa, el aumento de o produce perturbaciones similares y modestas que están dominadas por las señales más fuertes de fisión y epimutación. En consecuencia, aunque estos parámetros se mantuvieron en el modelo para la integridad estructural, sus distribuciones marginales posteriores rastrean de cerca sus distribuciones previas y no se informan como hallazgos biológicos.

Las glándulas de tumores colorrectales crecen como un proceso de ramificación casi puro

Aplicamos el flujo de trabajo de inferencia ABC-SMC a cada tumor de la cohorte de forma independiente para obtener distribuciones posteriores de los cinco parámetros del modelo (Métodos; Tabla 2). El marco de inferencia tuvo éxito en reproducir las características cualitativas clave (compare la Fig. 5 con la Fig. 2) y los aspectos cuantitativos de los datos (Fig. S3; Métodos). Estos resultados de inferencia sugieren que los datos son ampliamente consistentes con la suposición subyacente a nuestro modelo basado en agentes, que los tumores colorrectales crecen a través de sucesivas fisiones de glándulas a partir de un solo progenitor. Sin embargo, los tumores reales pueden experimentar un período de renovación de estado estacionario posterior a la expansión en el que el nacimiento y la muerte de las células continúan dentro de las glándulas a capacidad de carga sin una fisión adicional. Una pregunta clave es si una fase de renovación tendría un impacto sustancial en los patrones de divergencia fCpG interglandular en los que se basa la inferencia.

Dos características de los datos indican que la fase de expansión domina la señal de divergencia interglandular. Primero, la clara estructura en bloques en las matrices de distancia interglandular (Fig. 2C; Fig. S1) refleja una historia de ramificación jerárquica: las glándulas del mismo lado del tumor son consistentemente más similares que las glándulas de lados opuestos, como se esperaría si la divergencia de fCpG se acumulara principalmente durante el proceso de ramificación que produjo la separación espacial. Dado que la renovación en estado estacionario actúa simétricamente sobre todas las glándulas, independientemente de sus relaciones espaciales, tendería a erosionar la estructura jerárquica que emerge debido a la ramificación. La pronunciada estructura en bloques en toda la cohorte sugiere, por lo tanto, que la expansión es el principal contribuyente a la señal observada.

En segundo lugar, las distribuciones de fCpG dentro de las glándulas individuales conservan una característica forma de W (Fig. 2A, 2B), lo que indica que las tasas de epimutación se encuentran en el régimen intermedio donde los estados de metilación ancestrales se sobrescriben parcialmente, pero no por completo. Una renovación prolongada a la capacidad de carga impulsaría progresivamente las distribuciones dentro de las glándulas hacia una distribución unimodal centrada en 0,5 a través de una deriva neutral continua y una epimutación.

Para cuantificar los efectos de la renovación post-expansión, realizamos un análisis de sensibilidad que comparó las matrices de distancia simuladas en diferentes fracciones de renovación (0%, 15%, 30%, 50%) en una cuadrícula de 200 conjuntos de parámetros muestreados de la distribución a priori (Métodos). La renovación aumentó modestamente la divergencia general (la mediana de la distancia media por pares aumentó un 30% con un 50% de renovación; Fig. S4A), pero armonizó progresivamente la estructura en bloques entre lados y dentro de los lados (relación mediana: 1,14 con 0% a 1,04 con 30%; Fig. S4B, S4C). Las relaciones entre lados y dentro de los lados observadas en la cohorte (mediana 2,24) superan con creces la mediana en las simulaciones, incluso para una renovación del 0%, lo que sugiere que la expansión es el principal impulsor de la divergencia observada.

Las tasas de fisión glandular inferidas varían ampliamente entre los pacientes

La aplicación del marco de inferencia ABC-SMC a los 10 tumores produjo distribuciones posteriores para la tasa de fisión, la tasa de metilación y la tasa de desmetilación de cada tumor (Tabla 3). Las distribuciones posteriores de los parámetros restantes del modelo (es decir, la tasa de mutación del gen impulsor y la ventaja selectiva) se mantuvieron amplias en todos los tumores, lo que es consistente con una dinámica interglandular efectivamente neutral durante la expansión, y estos parámetros se excluyen, por lo tanto, de una mayor interpretación.

Las tasas de fisión por célula inferidas medianas abarcaron un factor de aproximadamente ocho, o casi una magnitud completa, en toda la cohorte, desde 2,2 × 10−3 célula−1 célula div−1 (tumor M) hasta 1,7 × 10−2 célula−1 célula div−1 (tumor U) (Fig. 6A). La tasa de metilación mediana abarcó 6,5 veces (rango: 7,5 × 10−4 a 4,8 × 10−3 célula div−1; Fig. 6B), comparable a la dispersión de la tasa de fisión, mientras que las tasas de desmetilación fueron más consistentes (4,1 veces, rango: 7,9 × 10−4 a 3,3 × 10−3 célula div−1; Fig. 6B, 7). La tasa de desmetilación relativamente conservada es consistente con la fidelidad del mantenimiento de CpG que actúa como una propiedad celular estable, mientras que las tasas de fisión varían más ampliamente y codifican las historias de crecimiento específicas del tumor. Tenga en cuenta, sin embargo, que los intervalos de credibilidad posteriores individuales del 90% para la tasa de fisión abarcan de 1,3 a 1,7 órdenes de magnitud dentro de cada tumor (Tabla 3), lo que es más amplio que la dispersión entre tumores de las medianas posteriores; por lo tanto, las comparaciones entre cohortes deben interpretarse como ordinales en lugar de como separaciones puntuales precisas.

Estimación de las edades mitóticas del tumor

Dado que methdemon utiliza el algoritmo de Gillespie [^[25]] para programar los eventos, todas las tasas inferidas están en unidades de (división celular)−1, y el marco, por lo tanto, estima directamente la edad mitótica de cada tumor, es decir, el número de divisiones de células madre cancerosas que ocurrieron durante la expansión. Multiplicar la tasa de fisión por célula por la capacidad de carga da la tasa de fisión por glándula. Bajo un crecimiento exponencial desde una sola glándula fundadora hasta n glándulas, la edad mitótica es (Métodos):

se estimó asumiendo una geometría esférica y un área de sección transversal de la glándula derivada de las dimensiones típicas de las glándulas de tumores colorrectales. Se obtuvieron distribuciones posteriores completas de n evaluando la ecuación 1 en cada muestra posterior de .

Las edades mitóticas medianas oscilaron entre 10 divisiones celulares (tumor U) y 73 divisiones celulares (tumor M) en toda la cohorte (Fig. 8). La conversión a tiempo calendario requiere una estimación de la tasa de división de las células madre cancerosas, que no se conoce con precisión en los tumores colorrectales. Utilizando las estimaciones de la literatura que abarcan de una vez por semana a una vez por mes [^[26], ^[27]], las edades mitóticas inferidas corresponden a tiempos de expansión tumoral del orden de varios meses a varios años. El amplio rango entre tumores se debe principalmente a la variación en las tasas de fisión, mientras que la incertidumbre dentro de cada tumor refleja tanto la incertidumbre posterior en n como la tasa de división de las células madre cancerosas desconocida. Tenga en cuenta que la tasa de división discutida aquí se refiere únicamente a las divisiones simétricas de las células madre cancerosas.

Asociaciones exploratorias con características clínicas

Los tres adenomas exhibieron tasas de metilación más bajas y tasas de desmetilación más altas que los siete carcinomas (prueba de Mann-Whitney en ambos casos; Fig. S5A), con diferencias de aproximadamente dos veces en promedio. No se observó una correlación significativa entre la tasa de fisión inferida y el tamaño del tumor (Spearman ), la edad del paciente o el estadio clínico (, solo carcinomas) (Fig. S5B). Se observó una tendencia entre la edad del paciente y la tasa de desmetilación (Spearman ), pero esto no alcanzó una significación nominal. Ni la edad del paciente ni el tamaño del tumor difirieron significativamente entre los adenomas y los carcinomas en esta cohorte (prueba de Mann-Whitney y , respectivamente), lo que sugiere que es poco probable que exista confusión por estas variables, aunque el tamaño de la muestra impide una evaluación definitiva. En particular, el amplio rango de variación en las tasas de fisión se conserva al restringirse a los carcinomas solos (; rango 0,0022–0,017; CV = 0,78), lo que demuestra que este hallazgo no depende de la agrupación de adenomas y carcinomas. Estas asociaciones exploratorias se basan en grupos pequeños y desequilibrados ( vs ) sin corrección para comparaciones múltiples, y deben considerarse observaciones generadoras de hipótesis que requieren validación en cohortes más grandes.

Cohorte de pacientes y datos de metilación de CpG fluctuantes

Se perfilaron matrices de metilación de muestras de 10 tumores (7 carcinomas, 3 adenomas, que abarcan un rango de tamaños y estadios clínicos; Tabla 1) utilizando la plataforma Illumina Infinium Methylation EPIC BeadChip, lo que produjo mediciones de valores beta en aproximadamente 850 000 sitios CpG por muestra de glándula con muy alta pureza (Métodos). Se muestrearon ocho glándulas de dos regiones separadas espacialmente de cada tumor resecado (lados A y B), lo que proporciona un conjunto de datos de 80 perfiles de metilación a nivel de glándula con una estructura espacial conocida (Fig. 1).

Las ganancias y pérdidas estocásticas de metilación en sitios CpG individuales durante la división celular somática producen una variación hereditaria e informativa para la línea de células en el estado de metilación del ADN [^[6]]. En los sitios CpG que no están constitutivamente metilados ni no metilados en toda la población, la fracción de metilación de cada glándula (el arreglo de fCpG) actúa como una huella digital de línea de células intrínseca: las células genéticamente idénticas se separan en el estado de metilación a través de rondas repetidas de división, y la diversidad del arreglo resultante codifica el número de divisiones celulares desde el ancestro común más reciente. Estos loci de CpG fluctuantes (fCpG) se han utilizado para reconstruir historias de linaje en tejidos normales y en el cáncer [^[6]–^[8]].

Para identificar los loci de fCpG, seguimos un procedimiento similar al de Gabbutt et al. [^[8]], utilizando datos de metilación de 138 muestras de tumores colorrectales en The Cancer Genome Atlas (TCGA) como una cohorte de referencia (Métodos). Este procedimiento identificó 1165 loci de fCpG cuyos estados de metilación llevan información sobre la historia de las divisiones celulares en lugar de la identidad del tipo de célula. En toda la cohorte, los arreglos de fCpG por glándula mostraron la característica distribución en forma de W de los valores beta del sitio (Fig. 2A, 2B), lo que es consistente con las tasas de epimutación en el régimen intermedio donde los estados de metilación ancestrales se sobrescriben parcialmente, pero no por completo.

La separación espacial predice la divergencia del arreglo de fCpG en muestras de cáncer colorrectal multirregión

Para cuantificar la divergencia de metilación interglandular, calculamos la matriz de distancia interglandular D para cada tumor, definida como la diferencia cuadrática media en la fracción de metilación de fCpG entre todos los pares de glándulas muestreadas (Ecuación 3; Métodos). Bajo la hipótesis de que las glándulas del mismo lado del tumor divergieron de un ancestro común más recientemente que las glándulas de lados opuestos, las distancias entre lados deben ser sistemáticamente menores que las distancias entre lados.

Esta predicción se confirmó en toda la cohorte. La matriz de distancia interglandular para cada tumor mostró una clara estructura en bloques: los pares de glándulas del mismo lado son consistentemente más similares que los pares de lados opuestos (Fig. 2C, S1). En los 10 tumores, la distancia por pares entre lados mediana fue significativamente menor que la distancia entre lados mediana (Fig. 3). Cada tumor exhibió el patrón esperado, con relaciones de distancia entre lados y dentro de los lados que oscilaron entre 1,12 (tumor D) y 4,71 (tumor I). Los tumores con relaciones más altas (I, U, M) tendieron a mostrar una estructura en bloques jerárquicos más clara en sus matrices de distancia, lo que es consistente con una mayor brecha entre los tiempos de coalescencia dentro de los lados y el tiempo de coalescencia de todo el tumor. La separación más débil observada en el tumor D (relación 1,12) puede reflejar su configuración de muestreo inusual (2 glándulas del lado A, 6 del lado B), lo que impide una estimación precisa de la distancia dentro del lado en el lado minoritario.

La magnitud de las distancias interglandulares varió sustancialmente entre los pacientes (Fig. S1), lo que sugiere heterogeneidad en la historia de crecimiento o las tasas de epimutación de los tumores individuales. Razonamos que esta variación codifica información sobre los parámetros evolutivos específicos del tumor que se pueden recuperar mediante un marco de inferencia calibrado, lo que motivó el desarrollo del modelo descrito en la siguiente sección.

Un modelo basado en agentes calibrado recapitula la dinámica de fCpG multiglándula

Para inferir los parámetros de crecimiento específicos del tumor a partir de los datos de fCpG observados, desarrollamos methdemon, un modelo basado en agentes del crecimiento del tumor colorrectal y la dinámica de metilación de fCpG. En este modelo, cada glándula se representa como un deme: una subpoblación local bien mezclada que crece y se reproduce como una unidad. El modelo simula la expansión del tumor como un proceso de ramificación de fisiones de glándulas, con cada glándula modelada como un deme de n células madre cancerosas que experimentan metilación y desmetilación estocásticas en los loci de fCpG (Fig. 4A). Antes de realizar la inferencia, caracterizamos cómo cada parámetro del modelo da forma a las estadísticas resumidas para determinar cuáles son identificables a partir de los datos del arreglo de fCpG. Ejecutamos simulaciones de methdemon para una muestra de hipercubo latino de 200 conjuntos de parámetros extraídos de nuestra distribución a priori de inferencia (Tabla 2), con todos los demás ajustes coincidentes con la configuración de inferencia (8 glándulas muestreadas, 2 lados, n coincidente con el tamaño promedio del tumor de la cohorte).

El barrido de parámetros reveló cómo las tasas de epimutación μ y ν (por división celular) controlan conjuntamente la forma y el sesgo de la distribución de fCpG por glándula (Fig. 4B). A bajas tasas de epimutación combinadas, los loci de fCpG retienen su estado fundador y la distribución es marcadamente bimodal, de modo que todos los valores β están cerca de 0 o 1. A altas tasas, la sobrescritura repetida del estado fundador impulsa la distribución hacia un pico unimodal cerca de la relación 0,5. A tasas intermedias, surge la característica forma de W, en la que la mayoría de los sitios permanecen cerca de 0 o 1, pero un subconjunto se ha desviado hacia valores intermedios.

La tasa de fisión controla la tasa por glándula a la que se producen nuevas glándulas durante la expansión y, por lo tanto, el número de divisiones celulares que separan las glándulas muestreadas. Una tasa más alta comprime la expansión del tumor en menos divisiones celulares y reduce la divergencia de fCpG interglandular, lo que influye directamente en la magnitud y la estructura en bloques de la matriz de distancia interglandular (Fig. 4C). Dado que la tasa de fisión por glándula solo entra en la dinámica a través del producto , la tasa de fisión y la capacidad de carga son estructuralmente no identificables: escalar por un factor se compensa exactamente escalando por . Confirmamos esta degeneración empíricamente con un análisis de sensibilidad. Las matrices de distancia en y difirieron por una norma de Frobenius comparable en magnitud a la variación de otros parámetros dentro del rango previo. Surge un régimen de tamaño finito en , donde la capacidad del grupo se acerca al número de linajes fundadores y aparece una deriva adicional. Por lo tanto, fijamos para la inferencia e interpretamos , en lugar de solo , como la cantidad biológicamente significativa.

Los análisis de sensibilidad establecieron que, si bien la epimutación y las tasas de fisión son identificables a partir de la matriz de distancia interglandular (Métodos; Fig. S2), la tasa de mutación del gen impulsor y la ventaja selectiva dejan solo firmas débiles y mutuamente confundidas. Dado que las mutaciones del gen impulsor confieren una ventaja de fisión multiplicativa, aumentar o produce perturbaciones similares y modestas que están dominadas por las señales más fuertes de fisión y epimutación. En consecuencia, aunque estos parámetros se mantuvieron en el modelo para la integridad estructural, sus distribuciones posteriores marginales siguen de cerca sus distribuciones previas y no se informan como hallazgos biológicos.

Las glándulas tumorales de cáncer colorrectal crecen como un proceso de ramificación casi puro

Aplicamos el flujo de trabajo de inferencia ABC-SMC a cada tumor de la cohorte de forma independiente para obtener distribuciones posteriores de los cinco parámetros del modelo (Métodos; Tabla 2). El marco de inferencia tuvo éxito en reproducir características cualitativas clave (compare la Fig. 5 con la Fig. 2) y aspectos cuantitativos de los datos (Fig. S3; Métodos). Estos resultados de inferencia sugieren que los datos son ampliamente consistentes con la suposición subyacente a nuestro modelo basado en agentes, que los tumores colorrectales crecen a través de sucesivas fisiones glandulares a partir de una sola célula progenitora. Sin embargo, los tumores reales pueden experimentar un período de renovación en estado estacionario post-expansión en el que la proliferación y la muerte celular continúan dentro de las glándulas a capacidad de carga sin una fisión adicional. Una pregunta clave es si una fase de renovación de este tipo alteraría sustancialmente los patrones de divergencia de fCpG interglandular en los que se basa la inferencia.

Dos características de los datos indican que la fase de expansión domina la señal de divergencia interglandular. Primero, la clara estructura en bloques en las matrices de distancia interglandular (Fig. 2C; Fig. S1) refleja una historia de ramificación jerárquica: las glándulas del mismo lado del tumor son consistentemente más similares que las glándulas de lados opuestos, como se espera si la divergencia de fCpG se acumula principalmente durante el proceso de ramificación que produjo la separación espacial. Dado que la renovación en estado estacionario actúa simétricamente sobre todas las glándulas independientemente de sus relaciones espaciales, tendería a erosionar la estructura jerárquica que surge debido a la ramificación. La pronunciada estructura en bloques en toda la cohorte sugiere, por lo tanto, que la expansión es el principal contribuyente a la señal observada.

En segundo lugar, las distribuciones de fCpG dentro de las glándulas individuales conservan una característica forma en W (Fig. 2A, 2B), lo que indica que las tasas de epimutación se encuentran en el régimen intermedio donde los estados de metilación ancestrales se sobrescriben parcialmente pero no por completo. Una renovación prolongada a capacidad de carga impulsaría progresivamente las distribuciones dentro de las glándulas hacia una distribución unimodal centrada en 0,5 a través de una deriva y una epimutación continuas y neutras.

Para cuantificar los efectos de la renovación post-expansión, realizamos un análisis de sensibilidad que comparó las matrices de distancia simuladas en diferentes fracciones de renovación (0 %, 15 %, 30 %, 50 %) en una cuadrícula de 200 conjuntos de parámetros muestreados del rango previo (Métodos). La renovación aumentó modestamente la divergencia general (la mediana de la distancia media por pares +30 % a una renovación del 50 %; Fig. S4A) pero armonizó progresivamente la estructura en bloques entre lados y dentro de los lados (relación mediana: 1,14 al 0 % a 1,04 al 30 %; Fig. S4B, S4C). Las relaciones entre lados/dentro de los lados observadas en la cohorte (mediana 2,24) superan con creces la mediana en las simulaciones incluso para una renovación del 0 %, lo que sugiere que la expansión es el principal impulsor de la divergencia observada.

Las tasas de fisión glandular inferidas varían ampliamente entre los pacientes

Aplicar el marco de inferencia ABC-SMC a los 10 tumores produjo distribuciones posteriores para la tasa de fisión, la tasa de metilación y la tasa de desmetilación de cada tumor (Tabla 3). Las distribuciones posteriores de los parámetros restantes del modelo (es decir, la tasa de mutación del gen impulsor y la ventaja selectiva) se mantuvieron amplias en todos los tumores, lo que es consistente con una dinámica interglandular efectivamente neutra durante la expansión, y estos parámetros se excluyen, por lo tanto, de una mayor interpretación.

Las tasas de fisión por célula inferidas medianas abarcaron un factor de aproximadamente ocho, o casi una magnitud completa, en toda la cohorte, desde 2,2 × 10−3 célula−1 división celular−1 (tumor M) hasta 1,7 × 10−2 célula−1 división celular−1 (tumor U) (Fig. 6A). La tasa de metilación mediana abarcó 6,5 veces (rango: 7,5 × 10−4 a 4,8 × 10−3 célula div−1; Fig. 6B), comparable a la dispersión de la tasa de fisión, mientras que las tasas de desmetilación fueron más consistentes (4,1 veces, rango: 7,9 × 10−4 a 3,3 × 10−3 célula div−1; Fig. 6B, 7). La tasa de desmetilación relativamente conservada es consistente con la fidelidad del mantenimiento de CpG que actúa como una propiedad celular biológica estable, mientras que las tasas de fisión varían más ampliamente y codifican las historias de crecimiento específicas del tumor. Tenga en cuenta, sin embargo, que los intervalos de credibilidad posteriores individuales del 90 % para la tasa de fisión abarcan de 1,3 a 1,7 órdenes de magnitud dentro de cada tumor (Tabla 3), lo que es más amplio que la dispersión entre tumores de las medianas posteriores; por lo tanto, las comparaciones entre cohortes deben interpretarse como ordinales en lugar de como separaciones puntuales precisas.

Estimación de las edades mitóticas de los tumores

Dado que methdemon utiliza el algoritmo de Gillespie [^[25]] para programar los eventos, todas las tasas inferidas están en unidades de (división celular)−1, y el marco, por lo tanto, estima directamente la edad mitótica de cada tumor, es decir, el número de divisiones celulares de las células madre cancerosas que ocurrieron durante la expansión. Multiplicar la tasa de fisión por célula por la capacidad de carga da la tasa de fisión por glándula. Bajo un crecimiento exponencial desde una sola glándula fundadora hasta glándulas, la edad mitótica es (Métodos):

se estimó asumiendo una geometría esférica y un área de sección transversal de la glándula derivada de las dimensiones típicas de las glándulas de cáncer colorrectal. Se obtuvieron distribuciones posteriores completas sobre evaluando la ecuación 1 en cada muestra posterior de .

Las edades mitóticas medianas variaron de 10 divisiones celulares (tumor U) a 73 divisiones celulares (tumor M) en toda la cohorte (Fig. 8). La conversión a tiempo calendario requiere una estimación de la tasa de división celular de las células madre cancerosas, que no se conoce con precisión en los tumores colorrectales. Utilizando estimaciones de la literatura que abarcan de una vez por semana a una vez por mes [^[26], ^[27]], las edades mitóticas inferidas corresponden a tiempos de expansión tumoral del orden de varios meses a varios años. El amplio rango en los tumores está impulsado principalmente por la variación en las tasas de fisión, mientras que la incertidumbre dentro de cada tumor refleja tanto la incertidumbre posterior en como la tasa de división de las células madre cancerosas desconocida. Tenga en cuenta que la tasa de división discutida aquí se refiere únicamente a las divisiones simétricas de las células madre cancerosas.

Asociaciones exploratorias con características clínicas

Los tres adenomas exhibieron tasas de metilación más bajas y tasas de desmetilación más altas que los siete carcinomas (Mann-Whitney en ambos casos; Fig. S5A), con aproximadamente dos veces de diferencia en promedio. No se observó una correlación significativa entre la tasa de fisión inferida y el tamaño del tumor (Spearman ), la edad del paciente o el estadio clínico (, carcinomas solamente) (Fig. S5B). Se observó una tendencia entre la edad del paciente y la tasa de desmetilación (Spearman ), pero esto no alcanzó una significación nominal. Ni la edad del paciente ni el tamaño del tumor difirieron significativamente entre los adenomas y los carcinomas en esta cohorte (Mann-Whitney y , respectivamente), lo que reduce la preocupación por la confusión por estas variables, aunque el tamaño de la muestra impide una evaluación definitiva. En particular, el amplio rango de variación en las tasas de fisión se conserva cuando el análisis se restringe a los carcinomas solos (; rango 0,0022–0,017; CV = 0,78), lo que demuestra que este hallazgo no depende de la agrupación de adenomas y carcinomas.

No se observó una correlación significativa entre la tasa de fisión inferida y el tamaño del tumor en la resección. Si bien esto puede reflejar parcialmente una potencia estadística limitada, el resultado también es consistente con la expectativa de que el tamaño del tumor en un momento dado es un indicador deficiente de la tasa de crecimiento: los tumores del mismo tamaño pueden haber estado creciendo durante diferentes duraciones, y el tamaño en la resección se ve afectado por el momento de la detección clínica. Además, la tasa de fisión, tal como se infiere aquí, refleja la dinámica de la fase de expansión, mientras que el tamaño en la resección integra el crecimiento, la renovación y cualquier período de latencia o regresión.

SPANISH TRANSLATION:
La tasa de fisión controla la tasa por glándula a la que se producen nuevas glándulas durante la expansión y, por lo tanto, el número de divisiones celulares que separan las glándulas muestreadas. Una tasa más alta comprime la expansión del tumor en menos divisiones celulares y reduce la divergencia de fCpG interglandular, lo que influye directamente en la magnitud y la estructura en bloques de la matriz de distancia interglandular (Fig. 4C). Dado que la tasa de fisión por glándula solo entra en la dinámica a través del producto , la tasa de fisión y la capacidad de carga son estructuralmente no identificables: escalar por un factor se compensa exactamente escalando por . Confirmamos esta degeneración empíricamente con un análisis de sensibilidad. Las matrices de distancia en y difirieron por una norma de Frobenius comparable en magnitud a la variación de otros parámetros dentro del rango previo. Surge un régimen de tamaño finito en , donde la capacidad del grupo se acerca al número de linajes fundadores y aparece una deriva adicional. Por lo tanto, fijamos para la inferencia e interpretamos , en lugar de solo , como la cantidad biológicamente significativa.

Los análisis de sensibilidad establecieron que, si bien la epimutación y las tasas de fisión son identificables a partir de la matriz de distancia interglandular (Métodos; Fig. S2), la tasa de mutación del gen impulsor y la ventaja selectiva dejan solo firmas débiles y mutuamente confundidas. Dado que las mutaciones del gen impulsor confieren una ventaja de fisión multiplicativa, aumentar o produce perturbaciones similares y modestas que están dominadas por las señales más fuertes de fisión y epimutación. En consecuencia, aunque estos parámetros se mantuvieron en el modelo para la integridad estructural, sus distribuciones posteriores marginales siguen de cerca sus distribuciones previas y no se informan como hallazgos biológicos.

Las glándulas tumorales de cáncer colorrectal crecen como un proceso de ramificación casi puro

Aplicamos el flujo de trabajo de inferencia ABC-SMC a cada tumor de la cohorte de forma independiente para obtener distribuciones posteriores de los cinco parámetros del modelo (Métodos; Tabla 2). El marco de inferencia tuvo éxito en reproducir características cualitativas clave (compare la Fig. 5 con la Fig. 2) y aspectos cuantitativos de los datos (Fig. S3; Métodos). Estos resultados de inferencia sugieren que los datos son ampliamente consistentes con la suposición subyacente a nuestro modelo basado en agentes, que los tumores colorrectales crecen a través de sucesivas fisiones glandulares a partir de una sola célula progenitora. Sin embargo, los tumores reales pueden experimentar un período de renovación en estado estacionario post-expansión en el que la proliferación y la muerte celular continúan dentro de las glándulas a capacidad de carga sin una fisión adicional. Una pregunta clave es si una fase de renovación de este tipo alteraría sustancialmente los patrones de divergencia de fCpG interglandular en los que se basa la inferencia.

Dos características de los datos indican que la fase de expansión domina la señal de divergencia interglandular. Primero, la clara estructura en bloques en las matrices de distancia interglandular (Fig. 2C; Fig. S1) refleja una historia de ramificación jerárquica: las glándulas del mismo lado del tumor son consistentemente más similares que las glándulas de lados opuestos, como se espera si la divergencia de fCpG se acumula principalmente durante el proceso de ramificación que produjo la separación espacial. Dado que la renovación en estado estacionario actúa simétricamente sobre todas las glándulas independientemente de sus relaciones espaciales, tendería a erosionar la estructura jerárquica que surge debido a la ramificación. La pronunciada estructura en bloques en toda la cohorte sugiere, por lo tanto, que la expansión es el principal contribuyente a la señal observada.

En segundo lugar, las distribuciones de fCpG dentro de las glándulas individuales conservan una característica forma en W (Fig. 2A, 2B), lo que indica que las tasas de epimutación se encuentran en el régimen intermedio donde los estados de metilación ancestrales se sobrescriben parcialmente pero no por completo. Una renovación prolongada a capacidad de carga impulsaría progresivamente las distribuciones dentro de las glándulas hacia una distribución unimodal centrada en 0,5 a través de una deriva y una epimutación continuas y neutras.

Para cuantificar los efectos de la renovación post-expansión, realizamos un análisis de sensibilidad que comparó las matrices de distancia simuladas en diferentes fracciones de renovación (0 %, 15 %, 30 %, 50 %) en una cuadrícula de 200 conjuntos de parámetros muestreados del rango previo (Métodos). La renovación aumentó modestamente la divergencia general (la mediana de la distancia media por pares +30 % a una renovación del 50 %; Fig. S4A) pero armonizó progresivamente la estructura en bloques entre lados y dentro de los lados (relación mediana: 1,14 al 0 % a 1,04 al 30 %; Fig. S4B, S4C). Las relaciones entre lados/dentro de los lados observadas en la cohorte (mediana 2,24) superan con creces la mediana en las simulaciones incluso para una renovación del 0 %, lo que sugiere que la expansión es el principal impulsor de la divergencia observada.

Las tasas de fisión glandular inferidas varían ampliamente entre los pacientes

Aplicar el marco de inferencia ABC-SMC a los 10 tumores produjo distribuciones posteriores para la tasa de fisión, la tasa de metilación y la tasa de desmetilación de cada tumor (Tabla 3). Las distribuciones posteriores de los parámetros restantes del modelo (es decir, la tasa de mutación del gen impulsor y la ventaja selectiva) se mantuvieron amplias en todos los tumores, lo que es consistente con una dinámica interglandular efectivamente neutra durante la expansión, y estos parámetros se excluyen, por lo tanto, de una mayor interpretación.

Las tasas de fisión por célula inferidas medianas abarcaron un factor de aproximadamente ocho, o casi una magnitud completa, en toda la cohorte, desde 2,2 × 10−3 célula−1 división celular−1 (tumor M) hasta 1,7 × 10−2 célula−1 división celular−1 (tumor U) (Fig. 6A). La tasa de metilación mediana abarcó 6,5 veces (rango: 7,5 × 10−4 a 4,8 × 10−3 célula div−1; Fig. 6B), comparable a la dispersión de la tasa de fisión, mientras que las tasas de desmetilación fueron más consistentes (4,1 veces, rango: 7,9 × 10−4 a 3,3 × 10−3 célula div−1; Fig. 6B, 7). La tasa de desmetilación relativamente conservada es consistente con la fidelidad del mantenimiento de CpG que actúa como una propiedad celular biológica estable, mientras que las tasas de fisión varían más ampliamente y codifican las historias de crecimiento específicas del tumor. Tenga en cuenta, sin embargo, que los intervalos de credibilidad posteriores individuales del 90 % para la tasa de fisión abarcan de 1,3 a 1,7 órdenes de magnitud dentro de cada tumor (Tabla 3), lo que es más amplio que la dispersión entre tumores de las medianas posteriores; por lo tanto, las comparaciones entre cohortes deben interpretarse como ordinales en lugar de como separaciones puntuales precisas.

Estimación de las edades mitóticas de los tumores

Dado que methdemon utiliza el algoritmo de Gillespie [^[25]] para programar los eventos, todas las tasas inferidas están en unidades de (división celular)−1, y el marco, por lo tanto, estima directamente la edad mitótica de cada tumor, es decir, el número de divisiones celulares de las células madre cancerosas que ocurrieron durante la expansión. Multiplicar la tasa de fisión por célula por la capacidad de carga da la tasa de fisión por glándula. Bajo un crecimiento exponencial desde una sola glándula fundadora hasta glándulas, la edad mitótica es (Métodos):

se estimó asumiendo una geometría esférica y un área de sección transversal de la glándula derivada de las dimensiones típicas de las glándulas de cáncer colorrectal. Se obtuvieron distribuciones posteriores completas sobre evaluando la ecuación 1 en cada muestra posterior de .

Las edades mitóticas medianas variaron de 10 divisiones celulares (tumor U) a 73 divisiones celulares (tumor M) en toda la cohorte (Fig. 8). La conversión a tiempo calendario requiere una estimación de la tasa de división celular de las células madre cancerosas, que no se conoce con precisión en los tumores colorrectales. Utilizando estimaciones de la literatura que abarcan de una vez por semana a una vez por mes [^[26], ^[27]], las edades mitóticas inferidas corresponden a tiempos de expansión tumoral del orden de varios meses a varios años. El amplio rango en los tumores está impulsado principalmente por la variación en las tasas de fisión, mientras que la incertidumbre dentro de cada tumor refleja tanto la incertidumbre posterior en como la tasa de división de las células madre cancerosas desconocida. Tenga en cuenta que la tasa de división discutida aquí se refiere únicamente a las divisiones simétricas de las células madre cancerosas.

Asociaciones exploratorias con características clínicas

Los tres adenomas exhibieron tasas de metilación más bajas y tasas de desmetilación más altas que los siete carcinomas (Mann-Whitney en ambos casos; Fig. S5A), con aproximadamente dos veces de diferencia en promedio. No se observó una correlación significativa entre la tasa de fisión inferida y el tamaño del tumor (Spearman ), la edad del paciente o el estadio clínico (, carcinomas solamente) (Fig. S5B). Se observó una tendencia entre la edad del paciente y la tasa de desmetilación (Spearman ), pero esto no alcanzó una significación nominal. Ni la edad del paciente ni el tamaño del tumor difirieron significativamente entre los adenomas y los carcinomas en

Es importante señalar varias limitaciones. Quizás lo más relevante es que el tamaño de la cohorte, compuesto por diez tumores, limita la potencia estadística para detectar asociaciones con características clínicas. El estudio también carece de un control basal de los criptas cólicas normales de los mismos pacientes. Perfiles de metilación de criptas normales emparejadas obtenidos de los márgenes quirúrgicos proporcionarían una línea de base específica del paciente y confirmarían que los parámetros inferidos son específicos del crecimiento tumoral. Nuestro hallazgo de que las tasas de epimutación tumoral son considerablemente más altas que las estimadas a partir de estudios de envejecimiento normal [^[1], ^[2]] es, no obstante, consistente con la dinámica acelerada durante la expansión neoplásica y también refleja nuestra selección de sitios fCpG que fluctúan en escalas de tiempo relevantes para el tumor.

La tasa de mutación impulsora y la ventaja selectiva resultaron ser no identificables a partir de los datos disponibles de la matriz fCpG. Por lo tanto, si bien nuestros resultados son consistentes con una dinámica interglandular efectivamente neutra [^[14], ^[15]], no excluyen una selección intraglandular débil a moderada. La capacidad de EVOFLUx para detectar cambios selectivos subclonales en cánceres linfoides [^[8]] implica el potencial de un muestreo más denso o dirigido dentro de las glándulas para probar la selección en tumores sólidos.

Nuestro modelo asume una población de células madre bien mezclada dentro de cada glándula y no captura la estructura espacial dentro de las glándulas individuales. Si bien las matrices obtenidas a partir de muestras masivas hacen que esta suposición sea razonable, uno podría preguntarse sobre la influencia de la estructura espacial de una glándula en su dinámica poblacional interna. Restringir el nicho de las células madre a una estructura poblacional cuasi-unidimensional en la base de una glándula podría cambiar el tamaño efectivo de la población en relación con un modelo bien mezclado del mismo tamaño censal, al influir en la tasa de deriva y fijación dentro de las glándulas [^[28]]. Se necesitaría un modelo más complejo y computacionalmente intensivo para explorar estos efectos.

Dos salvedades limitan la interpretación de nuestras edades mitóticas inferidas. En primer lugar, nuestro modelo cuenta solo las divisiones simétricas de las células madre a lo largo del árbol de expansión. Las células madre de cáncer colorrectal también se dividen asimétricamente y se replican a su máxima capacidad, y dado que la epimutación en el modelo está acoplada a la división, estos eventos también generan cambios en los sitios fCpG. Por lo tanto, las tasas inferidas y deben interpretarse como tasas de epimutación efectivas por división simétrica, análogas a la tasa de mutación efectiva por división en modelos neutrales de la evolución tumoral [^[15]]. La recuperación de una estimación del número total de divisiones requeriría conocer la fracción de divisiones de células madre que son simétricas, lo cual no es identificable a partir de los datos fCpG masivos. La segunda salvedad es que la conversión de las divisiones celulares a tiempo calendario introduce una incertidumbre sustancial porque las tasas de división de las células madre cancerosas no se conocen con precisión y pueden variar entre los tumores [^[26], ^[27]]. Una reescalada basada en estimaciones de las tasas de división de las CSC sitúa, no obstante, los tiempos de expansión tumoral dentro de un rango clínicamente plausible de varios meses a varios años [^[14]].

El diseño modular de methdemon y methabc facilita la adaptación a otros tipos de cáncer para los cuales se pueden obtener matrices de metilación resueltas espacialmente y se pueden construir modelos mecanicistas apropiados. Dado que los datos de entrada son relativamente económicos, en principio, el enfoque podría aplicarse a escala de cohorte para determinar si la historia de crecimiento proporciona información pronóstica. El marco también podría integrarse con datos genómicos emparejados para proporcionar una validación ortogonal de los parámetros de crecimiento inferidos mediante la comparación con filogenias basadas en mutaciones somáticas.

Métodos

Cohorte de pacientes y recolección de muestras

Las muestras de tumor (Tabla 1) fueron tejidos sobrantes obtenidos durante la atención clínica de rutina en la Facultad de Medicina de la Universidad del Sur de California, con la aprobación de la Junta de Revisión Institucional local. Las muestras de tumor se examinaron frescas y se tomaron fragmentos de aproximadamente 0,5 cm3 de lados opuestos del tumor. De cada fragmento, se aislaron glándulas tumorales individuales (~10 000 células) mediante un método de lavado con EDTA [^[29]].

Procesamiento de la matriz de metilación

La metilación del ADN de las glándulas (8 glándulas por tumor) se midió con matrices de cuentas EPIC (Illumina) utilizando el protocolo Restore y los protocolos del fabricante [^[30]]. Los archivos IDAT se procesaron utilizando la función de normalización noob en el paquete R minfi [^[31]].

Identificación de loci CpG fluctuantes

Los loci CpG fluctuantes (fCpG) se identificaron como el subconjunto de sitios CpG cuyo estado de metilación no está constitutivamente metilado ni no metilado en toda la población, y que muestran la menor coherencia entre las muestras, lo que los hace informativos sobre la dinámica de la población celular en lugar de la identidad del tipo de célula.

Identificamos los loci fCpG utilizando datos de la matriz de metilación del conjunto de datos de cáncer colorrectal de The Cancer Genome Atlas (TCGA), siguiendo una adaptación del procedimiento descrito en Gabbutt et al. [^[8]]. Brevemente, empleamos los 138 tumores colorrectales con una pureza de al menos 0,6 (evaluado mediante ABSOLUTE [^[32]]) y las 38 muestras de tejido de colon normal como nuestro conjunto de descubrimiento. Retuvimos los CpGs con una alta desviación estándar en las muestras de tumor masivas, con una metilación media de aproximadamente 0,5 y con una alta puntuación de Laplace (es decir, CpGs que no preservaron el gráfico de vecinos más cercanos). A diferencia de Gabbutt et al. [^[8]], también excluimos los CpGs que tenían una alta distancia absoluta de 0,5 en las muestras de colon normal. Este procedimiento identificó 1258 fCpGs específicos de los tumores colorrectales, de los cuales 1164 estaban presentes en todas nuestras muestras de pacientes y, por lo tanto, se utilizaron para todos los análisis posteriores. El uso de un gran conjunto de referencia para la identificación de loci garantiza la solidez frente a las limitaciones del tamaño de la muestra de los conjuntos de pacientes individuales.

methdemon: un modelo basado en demes de matrices de metilación fluctuantes en el cáncer colorrectal

Estructura del modelo

Desarrollamos methdemon, un modelo basado en agentes escrito en C++ para simular el crecimiento de un tumor colorrectal y la dinámica correspondiente de las matrices de metilación fCpG en las glándulas tumorales. El modelo representa un tumor como una colección de glándulas (demes), cada una de las cuales contiene una población de células madre cancerosas con potencial proliferativo hasta una capacidad máxima . El supuesto de mezcla completa dentro de cada glándula es consistente con el secuenciamiento masivo de muestras de glándulas y refleja la organización jerárquica de los tumores colorrectales, en la que un pequeño compartimento de células madre mantiene la población de glándulas [^[12], ^[13]]. Esta estructura del modelo está inspirada en demon, un modelo de oncología basado en demes [^[33], ^[34]] utilizado en estudios previos de la evolución tumoral [^[35]–^[37]].

La programación de eventos sigue el algoritmo de Gillespie [^[25]], en el cual primero se selecciona un deme con una probabilidad proporcional a su tasa de eventos total, seguido de la selección de una célula dentro de ese deme. A la capacidad máxima, la dinámica celular dentro de una glándula sigue un proceso de Moran, con eventos que consisten en el nacimiento celular, la muerte celular y la fisión de la glándula. Tras la división celular, cada sitio fCpG metilado en la célula hija sufre desmetilación de forma independiente con una probabilidad por sitio por división, y cada sitio no metilado sufre metilación con una probabilidad por sitio por división. Los eventos de epimutación están, por lo tanto, acoplados a la división celular, lo que es consistente con el origen predominantemente acoplado a la replicación de los errores de metilación [^[1], ^[2]].

Las mutaciones impulsoras surgen con una probabilidad por célula por división y confieren una ventaja proliferativa multiplicativa a su portador. Cada célula rastrea su matriz fCpG como un vector binario de longitud , y cada deme rastrea la matriz fCpG promedio en su población celular como un vector de fracciones de metilación.

Fisión de la glándula

La fisión de la glándula se modela como un proceso de ramificación espacial neutral, consistente con el mecanismo predominante del crecimiento del tumor colorrectal [^[10]] y con los hallazgos que respaldan la dinámica interglandular neutral [^[14]]. Un evento de fisión distribuye la población de la glándula aleatoriamente en dos mitades iguales, tras lo cual cada glándula hija repobla de forma independiente hasta la capacidad máxima. La probabilidad por célula de fisión determina la velocidad a la que se producen nuevas glándulas.

Para centrar los recursos computacionales en el subconjunto de glándulas que se muestrean finalmente, el modelo distingue entre fisiones rastreadas y no rastreadas. Las fisiones no rastreadas producen glándulas que no se siguen hasta el final de la simulación; contribuyen a la dinámica de la población del tumor, pero no a la salida final. La probabilidad de que un evento de fisión se rastree es

Esto asegura que el número esperado de eventos de fisión antes de una fisión rastreada sea de aproximadamente la mitad del número total medio de fisiones. Además, la distribución de los intervalos de fisión en el tiempo sigue la distribución de Poisson con una media , como se deriva en [^[38]].

Parámetros del modelo y condiciones de parada

La lista completa de parámetros del modelo se proporciona en la Tabla 4. La simulación termina cuando el número medio de fisiones por deme rastreado alcanza un valor objetivo , correspondiente a un tumor de glándulas en expectativa. Se puede activar una fase opcional de equilibrio de rotación post-crecimiento, en la que el nacimiento y la muerte celular continúan a la capacidad máxima sin una fisión adicional, lo que aproxima un régimen de crecimiento saturado.

methabc: flujo de inferencia de computación bayesiana aproximada

Marco ABC-SMC

La inferencia de parámetros se realizó utilizando la computación bayesiana aproximada con Monte Carlo secuencial (ABC-SMC), implementada en un clúster de alto rendimiento con el paquete pyabc [^[39], ^[40]]. ABC-SMC refina iterativamente una población de partículas de parámetros a través de generaciones sucesivas, aceptando partículas cuyos resultados simulados se encuentran dentro de una tolerancia de los datos observados. El umbral de tolerancia se determinó de forma adaptativa en cada generación utilizando el método SilkOptimalEpsilon [^[41]]. Cada generación utilizó 1000 partículas con una tasa de aceptación de aproximadamente el 2%. La inferencia se ejecutó durante 9 a 13 generaciones hasta que el umbral de tolerancia y las distribuciones marginales posteriores se estabilizaron entre generaciones sucesivas.

Todos los parámetros, excepto la ventaja selectiva, se extrajeron de una distribución -transformada para garantizar una exploración uniforme del espacio de parámetros en múltiples órdenes de magnitud. Las distribuciones previas se proporcionan en la Tabla 2.

La capacidad del deme se fijó en para todas las ejecuciones de inferencia, lo que es consistente con las estimaciones de la fracción de células madre de tumores colorrectales de ~1% [^[11], ^[12]] dadas las dimensiones típicas de las glándulas de ~104 células. Dado que la tasa de fisión por glándula entra en el modelo como el producto , el valor específico de es una convención de modelado en lugar de un parámetro libre, y es la cantidad biológicamente significativa inferida (§).

Estadísticos de resumen

Se utilizaron dos estadísticos de resumen complementarios para comparar los datos de la matriz fCpG simulados y observados.

Matriz de distancia interglandular.

Para un tumor que consta de glándulas muestreadas con matrices fCpG promedio , la matriz de distancia interglandular D se define como

donde es el número de sitios fCpG. La diferencia al cuadrado enfatiza las divergencias interglandulares grandes y reduce la influencia de las pequeñas fluctuaciones estocásticas. Se comparan dos matrices de distancia interglandular y utilizando la norma de Frobenius normalizada de su diferencia:

donde el factor de 1/2 tiene en cuenta la simetría de la matriz de distancia.

Distancia de la distribución fCpG por glándula.

Para capturar información sobre la forma y el sesgo de las distribuciones fCpG individuales de las glándulas, lo que codifica las tasas relativas de metilación y desmetilación, también calculamos la distancia media de Wasserstein entre las distribuciones fCpG correspondientes en los datos simulados y observados:

donde es el conjunto de todas las distribuciones conjuntas con marginales y .

La distancia total utilizada para la aceptación de ABC fue una suma de los dos estadísticos:

donde es la distancia de Wasserstein media por glándula. La contribución relativa de los dos componentes se gestionó mediante pyabc.AdaptiveAggregatedDistance, que reescala cada componente en cada generación ABC-SMC para que ambos contribuyan aproximadamente por igual al umbral de aceptación. Este esquema adaptativo elimina la necesidad de un parámetro de ponderación fijo y se adapta a los cambios de escala a medida que la distribución posterior se concentra durante las generaciones sucesivas.

Comprobaciones predictivas posteriores

Realizamos comprobaciones predictivas posteriores para evaluar si el modelo inferido reproduce adecuadamente los datos observados. Para cada tumor, se muestrearon 100 conjuntos de parámetros de la distribución posterior ABC de la generación final. Para cada muestra, se realizó una simulación de tumor utilizando methdemon con la misma configuración que la ejecución de inferencia correspondiente (8 glándulas muestreadas, 2 lados). Se calculó la matriz de distancia interglandular simulada para cada réplica y se comparó la distribución de las distancias por pares simuladas con los valores observados. Se consideró que un ajuste del modelo era adecuado si la distancia media interglandular observada se encontraba dentro del 90% central de la distribución predictiva posterior. La cobertura por pares se definió como la fracción de distancias interglandulares observadas que se encontraban dentro del intervalo predictivo posterior del 90% elemento por elemento.

Para los 10 tumores del estudio, la distancia interglandular observada se encontró dentro del 90% central de la distribución predictiva posterior (Fig. S3). La cobertura por pares (la fracción de distancias interglandulares que se encuentran dentro del intervalo predictivo posterior del 90%) superó el 96% para todos los tumores (rango: 96-100%). Las distancias por pares dentro del mismo lado mostraron una cobertura ligeramente menor (75-100%) que las distancias entre lados (100% para todos los tumores), lo que sugiere que el modelo captura la magnitud general y la estructura jerárquica de la divergencia de fCpG, aunque la variación intra-lado puede estar subestimada en una minoría de los tumores.

Recuperación de parámetros sintéticos

Para validar que el flujo de trabajo de inferencia recupera los parámetros identificables, se realizaron experimentos de recuperación de parámetros sintéticos. Se extrajo un conjunto de vectores de parámetros verdaderos de la distribución previa, y para cada uno, se simuló un tumor sintético utilizando methdemon. Luego, se aplicó el flujo de trabajo completo de inferencia ABC-SMC a cada conjunto de datos sintéticos para obtener una distribución posterior. La recuperación se evaluó comparando la mediana posterior con el valor verdadero para cada parámetro y reportando la fracción de experimentos en los que el valor verdadero se encontraba dentro del intervalo creíble del 90% (Fig. S2). Las tasas de epimutación y la tasa de fisión por glándula se recuperaron de manera confiable; y no se recuperaron, lo que es consistente con el análisis de sensibilidad descrito anteriormente.

Estimación de la edad mitótica

La edad mitótica de cada tumor se estimó a partir de la distribución posterior de la tasa de fisión inferida de la siguiente manera.

Bajo el algoritmo de Gillespie, todas las tasas en methdemon se expresan por división celular. La tasa de fisión por glándula es , donde es la tasa de fisión por célula inferida y es la capacidad de carga del deme. Bajo un crecimiento exponencial a partir de una sola glándula fundadora, el número esperado de glándulas después de divisiones celulares es , lo que da

donde es el número total de glándulas en el tumor en el momento de la resección.

Dado que los recuentos de glándulas no están disponibles directamente, se estimó a partir del diámetro tumoral aproximado mediante el supuesto de una geometría esférica:

donde se derivó de las dimensiones típicas de las glándulas de tumores colorrectales: aproximadamente 70 μm de diámetro y 450 μm de profundidad [^[13], ^[26]].

Para cada muestra posterior de , calculamos , lo que da como resultado una distribución posterior completa de la edad mitótica. La conversión a tiempo calendario se realizó dividiendo por las tasas de división de CSC asumidas que abarcan desde una división por semana hasta una división por mes, siguiendo las estimaciones publicadas [^[26], ^[27]].

El estimador debe interpretarse como la edad mitótica simétrica de la fase de expansión: methdemon contiene solo divisiones simétricas, por lo que las contribuciones asimétricas y post-expansión se absorben en las tasas inferidas en lugar de resolverse como divisiones adicionales.

Como comprobación cruzada independiente, también estimamos las edades mitóticas a partir de las distribuciones de fCpG intra-glandulares utilizando una aproximación de dos estados a la dinámica de metilación de la cadena de Markov descrita en [^[2]]. Para cada glándula, se calculó la polarización media, definida como

en todos los loci de fCpG (Fig. S6A). Esta cantidad mide qué tan lejos se ha desviado la matriz de metilación de los estados completamente metilados o no metilados. Tratando cada locus como un sistema de dos estados en el que los eventos de metilación y desmetilación ocurren independientemente en cada división celular, la polarización disminuye geométricamente, y la edad mitótica se puede estimar como

donde es la polarización máxima observada en todo el estudio y son las tasas de epimutación medianas posteriores para ese tumor (Fig. S6B).

Análisis de sensibilidad del recambio

Para evaluar si la inclusión de un recambio en estado estacionario alteraría sustancialmente los patrones de divergencia interglandular en los que se basa la inferencia, realizamos una búsqueda en la cuadrícula en el espacio de parámetros del modelo en cuatro longitudes diferentes de la fase de recambio post-crecimiento (0%, 15%, 30%, 50%).

Se extrajeron 200 conjuntos de parámetros de la distribución previa de inferencia mediante el muestreo de hipercubo latino. Para cada conjunto de parámetros, se ejecutó methdemon en las cuatro fracciones de recambio utilizando la misma semilla aleatoria, de modo que la fase de expansión fuera idéntica y la única diferencia fuera la duración del período de recambio posterior. Para cada ejecución, se calculó la matriz de distancia interglandular de 8 × 8 y se resumió mediante (i) la distancia media por pares, (ii) la relación de distancias entre lados y dentro del mismo lado, y (iii) la distancia de Frobenius entre la matriz de distancia en cada fracción de recambio y la línea de base del 30%.

Análisis estadístico

Todas las pruebas estadísticas fueron bicaudales con un umbral de significación de . Dado el pequeño tamaño del estudio, los tamaños del efecto se informan junto con los valores de p y los resultados se interpretan con la precaución adecuada con respecto al poder estadístico. La comparación entre las distancias dentro del mismo lado y entre lados se evaluó utilizando una prueba de rango con signo de Wilcoxon en los 10 tumores. La correlación de rango de Spearman se utilizó para evaluar la relación entre los valores de parámetros medianos posteriores inferidos y las variables clínicas continuas (tamaño del tumor, edad del paciente). Se utilizaron pruebas de Mann-Whitney para comparar los valores de parámetros inferidos entre los grupos categóricos (sexo, adenoma versus carcinoma). La consistencia de las tasas de epimutación inferidas entre los pacientes se cuantificó utilizando el coeficiente de variación (CV) de las medianas posteriores. Todos los análisis estadísticos se realizaron en Python 3.12 utilizando scipy v1.17.1 (módulo scipy.stats).

Cohorte de pacientes y recolección de muestras

Las muestras de tumor (Tabla 1) fueron tejidos sobrantes obtenidos durante la atención clínica de rutina en la Facultad de Medicina de la Universidad del Sur de California con la aprobación del Comité de Revisión Institucional local. Las muestras de tumor se examinaron frescas y se tomaron fragmentos de aproximadamente 0,5 cm3 de los lados opuestos del tumor. De cada fragmento, se aislaron glándulas tumorales individuales (~10 000 células) con un método de lavado con EDTA [^[29]].

Procesamiento de la matriz de metilación

El ADN de metilación de las glándulas (8 glándulas por tumor) se midió con matrices EPIC (Illumina) utilizando el protocolo Restore y los protocolos del fabricante [^[30]]. Los archivos IDAT se procesaron utilizando la función de normalización noob en el paquete R minfi [^[31]].

Identificación de loci de CpG fluctuantes

Los loci de CpG fluctuantes (fCpG) se identificaron como el subconjunto de sitios de CpG cuyo estado de metilación no es constitutivamente metilado ni no metilado en toda la población, y que muestran la menor coherencia entre las muestras, lo que los hace informativos sobre la dinámica de la población celular en lugar de la identidad del tipo de célula.

Identificamos los loci de fCpG utilizando datos de la matriz de metilación del estudio del cáncer del genoma (TCGA) de la cohorte de cáncer colorrectal, siguiendo una adaptación del procedimiento descrito en Gabbutt et al. [^[8]]. Brevemente, empleamos los 138 tumores colorrectales con una pureza de al menos 0,6 (evaluado mediante ABSOLUTE [^[32]]) y las 38 muestras de tejido de colon normal como nuestro conjunto de descubrimiento. Retuvimos los CpGs con una alta desviación estándar en todas las muestras de tumor, con una metilación media de aproximadamente 0,5 y con una alta puntuación de Laplace (es decir, CpGs que no preservan el gráfico de vecinos más cercanos). A diferencia de Gabbutt et al. [^[8]], también excluimos los CpGs que tenían una alta distancia absoluta de 0,5 en las muestras de colon normal. Este procedimiento identificó 1258 fCpGs específicos de los tumores colorrectales, de los cuales 1164 estaban presentes en todas nuestras muestras de pacientes y, por lo tanto, se utilizaron para todos los análisis posteriores. El uso de una gran cohorte de referencia para la identificación de loci garantiza la solidez frente a las limitaciones del tamaño de la muestra de las cohortes de pacientes individuales.

methdemon: un modelo basado en demes de matrices de metilación fluctuantes en el cáncer colorrectal

Estructura del modelo

Desarrollamos methdemon, un modelo basado en agentes escrito en C++ para simular el crecimiento de un tumor colorrectal y la dinámica correspondiente de las matrices de metilación de fCpG en las glándulas tumorales. El modelo representa un tumor como una colección de glándulas (demes), cada una de las cuales contiene una población de células madre cancerosas con potencial proliferativo hasta una capacidad de carga máxima. El supuesto de mezcla homogénea dentro de cada glándula es consistente con el secuenciamiento masivo de muestras de glándulas y refleja la organización jerárquica de los tumores colorrectales, en la que un pequeño compartimento de células madre mantiene la población de glándulas [^[12], ^[13]]. Esta estructura del modelo está inspirada en demon, un modelo de oncología basado en demes [^[33], ^[34]] utilizado en estudios anteriores de la evolución tumoral [^[35]–^[37]].

La programación de eventos sigue el algoritmo de Gillespie [^[25]], en el que primero se selecciona un deme con una probabilidad proporcional a su tasa de eventos total, seguido de la selección de una célula dentro de ese deme. En la capacidad de carga, la dinámica celular dentro de una glándula sigue un proceso de Moran, con eventos que consisten en el nacimiento de células, la muerte de células y la fisión de glándulas. Tras la división celular, cada sitio de fCpG metilado en la célula hija sufre desmetilación de forma independiente con una probabilidad de por sitio por división, y cada sitio no metilado sufre metilación con una probabilidad de por sitio por división. Los eventos de epimutación están, por lo tanto, acoplados a la división celular, lo que es consistente con el origen predominantemente acoplado a la replicación de los errores de metilación [^[1], ^[2]].

Las mutaciones impulsoras surgen con una probabilidad de por célula por división y confieren una ventaja proliferativa multiplicativa a su portador. Cada célula rastrea su matriz de fCpG como un vector binario de longitud , y cada deme rastrea la matriz de fCpG promedio en toda su población celular como un vector de fracciones de metilación.

Fisión de glándulas

La fisión de glándulas se modela como un proceso de ramificación espacial neutral, lo que es consistente con el mecanismo predominante del crecimiento del tumor colorrectal [^[10]] y con los hallazgos que respaldan la dinámica interglandular neutral [^[14]]. Un evento de fisión distribuye la población de glándulas aleatoriamente en dos mitades iguales, después de lo cual cada glándula hija repobla independientemente hasta la capacidad de carga. La probabilidad de fisión por célula determina la velocidad a la que se producen nuevas glándulas.

Para centrar los recursos computacionales en el subconjunto de glándulas que se muestrean en última instancia, el modelo distingue entre fisiones rastreadas y no rastreadas. Las fisiones no rastreadas producen glándulas que no se siguen hasta el final de la simulación; contribuyen a la dinámica de la población del tumor, pero no a la salida final. La probabilidad de que un evento de fisión se rastree es

Esto asegura que el número esperado de eventos de fisión antes de una fisión rastreada sea aproximadamente la mitad del número total medio de fisiones. Además, la distribución de los intervalos de fisión en el tiempo sigue la distribución de Poisson con una media , como se deriva en [^[38]].

Parámetros del modelo y condiciones de parada

La lista completa de parámetros del modelo se presenta en la Tabla 4. La simulación se detiene cuando el número medio de fisiones por deme rastreado alcanza un valor objetivo, correspondiente a un tumor de glándulas en promedio. Se puede activar una fase opcional de recambio en estado estacionario post-crecimiento, en la que el nacimiento y la muerte celular continúan a la capacidad de carga sin nuevas fisiones, lo que aproxima un régimen de crecimiento saturado.

Estructura del modelo

Desarrollamos methdemon, un modelo basado en agentes escrito en C++ para simular el crecimiento de un tumor colorrectal y la dinámica correspondiente de las matrices de metilación de fCpG en las glándulas tumorales. El modelo representa un tumor como una colección de glándulas (demes), cada una de las cuales contiene una población de células madre cancerosas con potencial proliferativo hasta una capacidad de carga máxima. El supuesto de mezcla homogénea dentro de cada glándula es consistente con el secuenciamiento masivo de muestras de glándulas y refleja la organización jerárquica de los tumores colorrectales, en la que un pequeño compartimento de células madre mantiene la población de glándulas [^[12], ^[13]]. Esta estructura del modelo está inspirada en demon, un modelo de oncología basado en demes [^[33], ^[34]] utilizado en estudios previos de la evolución tumoral [^[35]–^[37]].

La programación de eventos sigue el algoritmo de Gillespie [^[25]], en el que primero se selecciona un deme con una probabilidad proporcional a su tasa de eventos total, seguido de la selección de una célula dentro de ese deme. A la capacidad de carga, la dinámica celular dentro de una glándula sigue un proceso de Moran, con eventos que consisten en el nacimiento celular, la muerte celular y la fisión de la glándula. Tras la división celular, cada sitio de fCpG metilado en la célula hija se desmetila de forma independiente con una probabilidad por sitio por división, y cada sitio no metilado se metila con una probabilidad por sitio por división. Los eventos de epimutación están, por lo tanto, acoplados a la división celular, lo que es consistente con el origen predominantemente acoplado a la replicación de los errores de metilación [^[1], ^[2]].

Las mutaciones de genes impulsores surgen con una probabilidad por célula por división y confieren una ventaja proliferativa multiplicativa al portador. Cada célula rastrea su matriz de fCpG como un vector binario de longitud , y cada deme rastrea la matriz de fCpG promedio en su población celular como un vector de fracciones de metilación.

Fisión de glándulas

La fisión de glándulas se modela como un proceso de ramificación espacial neutral, lo que es consistente con el mecanismo predominante del crecimiento del tumor colorrectal [^[10]] y con los hallazgos que respaldan la dinámica interglandular neutral [^[14]]. Un evento de fisión distribuye la población de la glándula aleatoriamente en dos mitades iguales, tras lo cual cada glándula hija repobla independientemente hasta la capacidad de carga. La probabilidad de fisión por célula determina la velocidad a la que se producen nuevas glándulas.

Para centrar los recursos computacionales en el subconjunto de glándulas que se muestrean finalmente, el modelo distingue entre fisiones rastreadas y no rastreadas. Las fisiones no rastreadas producen glándulas que no se siguen hasta el final de la simulación; contribuyen a la dinámica de la población del tumor, pero no a la salida final. La probabilidad de que un evento de fisión se rastree es

Esto asegura que el número esperado de eventos de fisión antes de una fisión rastreada sea aproximadamente la mitad del número total medio de fisiones. Además, la distribución de los intervalos de fisión en el tiempo sigue la distribución de Poisson con una media , como se deriva en [^[38]].

Parámetros del modelo y condiciones de parada

La lista completa de parámetros del modelo se presenta en la Tabla 4. La simulación se detiene cuando el número medio de fisiones por deme rastreado alcanza un valor objetivo, correspondiente a un tumor de glándulas en promedio. Se puede activar una fase opcional de recambio en estado estacionario post-crecimiento, en la que el nacimiento y la muerte celular continúan a la capacidad de carga sin nuevas fisiones, lo que aproxima un régimen de crecimiento saturado.

methabc: flujo de trabajo de inferencia de computación bayesiana aproximada

Marco ABC-SMC

La inferencia de parámetros se realizó utilizando la computación bayesiana aproximada con Monte Carlo secuencial (ABC-SMC), implementada en un clúster de alto rendimiento con el paquete pyabc [^[39], ^[40]]. ABC-SMC refina iterativamente una población de partículas de parámetros a través de generaciones sucesivas, aceptando partículas cuyos resultados simulados se encuentran dentro de una tolerancia de los datos observados. El umbral de tolerancia se determinó de forma adaptativa en cada generación utilizando el método SilkOptimalEpsilon [^[41]]. Cada generación utilizó 1000 partículas con una tasa de aceptación de aproximadamente el 2%. La inferencia se ejecutó durante 9-13 generaciones hasta que el umbral de tolerancia y las distribuciones posteriores marginales se estabilizaron entre generaciones sucesivas.

Todos los parámetros, excepto la ventaja selectiva, se extrajeron de una distribución -transformada para garantizar una exploración uniforme del espacio de parámetros en múltiples órdenes de magnitud. Las distribuciones previas se dan en la Tabla 2.

La capacidad de carga del deme se fijó en para todas las ejecuciones de inferencia, lo que es consistente con las estimaciones de la fracción de células madre del tumor colorrectal de aproximadamente el 1% [^[11], ^[12]] dadas las dimensiones típicas de las glándulas de aproximadamente 104 células. Dado que la tasa de fisión por glándula entra en el modelo como el producto , el valor específico de es una convención de modelado en lugar de un parámetro libre, y es la cantidad inferida biológicamente significativa (§).

Estadísticos de resumen

Se utilizaron dos estadísticos de resumen complementarios para comparar los datos de la matriz de fCpG simulados y observados.

Matriz de distancia interglandular.

Para un tumor que consta de glándulas muestreadas con matrices de fCpG promedio , la matriz de distancia interglandular D se define como

donde es el número de sitios de fCpG. La diferencia al cuadrado enfatiza las grandes divergencias interglandulares y reduce la influencia de las pequeñas fluctuaciones estocásticas. Se comparan dos matrices de distancia interglandular y utilizando la norma de Frobenius normalizada de su diferencia:

donde el factor de 1/2 tiene en cuenta la simetría de la matriz de distancia.

Distancia de la distribución de fCpG por glándula.

Para capturar información sobre la forma y el sesgo de las distribuciones de fCpG individuales de las glándulas, lo que codifica las tasas relativas de metilación y desmetilación, también calculamos la distancia de Wasserstein media entre las distribuciones de fCpG correspondientes de las glándulas en los datos simulados y observados:

donde es el conjunto de todas las distribuciones conjuntas con marginales y .

La distancia total utilizada para la aceptación de ABC fue una suma de los dos estadísticos:

donde es la distancia de Wasserstein media por glándula. La contribución relativa de los dos componentes se gestionó mediante pyabc.AdaptiveAggregatedDistance, que reescala cada componente en cada generación de ABC-SMC para que ambos contribuyan aproximadamente por igual al umbral de aceptación. Este esquema adaptativo elimina la necesidad de un parámetro de ponderación fijo y se adapta a los cambios de escala a medida que la distribución posterior se concentra durante las generaciones sucesivas.

Comprobaciones predictivas posteriores

Realizamos comprobaciones predictivas posteriores para evaluar si el modelo inferido reproduce adecuadamente los datos observados. Para cada tumor, se muestrearon 100 conjuntos de parámetros del conjunto de distribución posterior de la última generación de ABC. Para cada , se ejecutó una simulación de tumor utilizando methdemon con la misma configuración que la ejecución de inferencia correspondiente (8 glándulas muestreadas, 2 lados). Se calculó la matriz de distancia interglandular para cada réplica, y se comparó la distribución de las distancias por pares simuladas con los valores observados. Se consideró que un ajuste del modelo era adecuado si la distancia interglandular media observada se encontraba dentro del 90% central de la distribución predictiva posterior. La cobertura por pares se definió como la fracción de distancias interglandulares observadas que caen dentro del intervalo predictivo posterior del 90% por elemento.

Para los 10 tumores del grupo, la distancia interglandular observada se encontró dentro del 90% central de la distribución predictiva posterior (Fig. S3). La cobertura por pares (la fracción de distancias interglandulares que caen dentro del intervalo predictivo posterior del 90%) superó el 96% para todos los tumores (rango: 96-100%). Las distancias por pares dentro del mismo lado mostraron una cobertura ligeramente menor (75-100%) que las distancias entre lados (100% para todos los tumores), lo que sugiere que el modelo captura la magnitud general y la estructura jerárquica de la divergencia de fCpG, aunque la variación intra-lado puede estar subestimada en una minoría de los tumores.

Recuperación de parámetros sintéticos

Para validar que el flujo de trabajo de inferencia recupera los parámetros identificables, se realizaron experimentos de recuperación de parámetros sintéticos. Se extrajo un conjunto de vectores de parámetros verdaderos del conjunto de distribución previa, y para cada se simuló un tumor sintético utilizando methdemon. A continuación, se aplicó todo el flujo de trabajo de inferencia de ABC-SMC a cada conjunto de datos sintéticos para obtener una distribución posterior . La recuperación se evaluó comparando la mediana posterior con el valor verdadero para cada parámetro, y informando de la fracción de experimentos en los que el valor verdadero se encontraba dentro del intervalo creíble del 90% (Fig. S2). Las tasas de epimutación y la tasa de fisión por glándula se recuperaron de forma fiable; y no se recuperaron, lo que es consistente con el análisis de sensibilidad descrito anteriormente.

Marco ABC-SMC

La inferencia de parámetros se realizó utilizando la computación bayesiana aproximada con Monte Carlo secuencial (ABC-SMC), implementada en un clúster de alto rendimiento con el paquete pyabc [^[39], ^[40]]. ABC-SMC refina iterativamente una población de partículas de parámetros a través de generaciones sucesivas, aceptando partículas cuyos resultados simulados se encuentran dentro de una tolerancia de los datos observados. El umbral de tolerancia se determinó de forma adaptativa en cada generación utilizando el método SilkOptimalEpsilon [^[41]]. Cada generación utilizó 1000 partículas con una tasa de aceptación de aproximadamente el 2%. La inferencia se ejecutó durante 9-13 generaciones hasta que el umbral de tolerancia y las distribuciones posteriores marginales se estabilizaron entre generaciones sucesivas.

Todos los parámetros, excepto la ventaja selectiva, se extrajeron de una distribución -transformada para garantizar una exploración uniforme del espacio de parámetros en múltiples órdenes de magnitud. Las distribuciones previas se dan en la Tabla 2.

La capacidad de carga del deme se fijó en para todas las ejecuciones de inferencia, lo que es consistente con las estimaciones de la fracción de células madre del tumor colorrectal de aproximadamente el 1% [^[11], ^[12]] dadas las dimensiones típicas de las glándulas de aproximadamente 104 células. Dado que la tasa de fisión por glándula entra en el modelo como el producto , el valor específico de es una convención de modelado en lugar de un parámetro libre, y es la cantidad inferida biológicamente significativa (§).

Estadísticos de resumen

Se utilizaron dos estadísticos de resumen complementarios para comparar los datos de la matriz de fCpG simulados y observados.

Matriz de distancia interglandular.

Para un tumor que consta de glándulas muestreadas con matrices de fCpG promedio , la matriz de distancia interglandular D se define como

donde es el número de sitios de fCpG. La diferencia al cuadrado enfatiza las grandes divergencias interglandulares y reduce la influencia de las pequeñas fluctuaciones estocásticas. Se comparan dos matrices de distancia interglandular y utilizando la norma de Frobenius normalizada de su diferencia:

donde el factor de 1/2 tiene en cuenta la simetría de la matriz de distancia.

Distancia de la distribución de fCpG por glándula.

Para capturar información sobre la forma y el sesgo de las distribuciones de fCpG individuales de las glándulas, lo que codifica las tasas relativas de metilación y desmetilación, también calculamos la distancia de Wasserstein media entre las distribuciones de fCpG correspondientes de las glándulas en los datos simulados y observados:

donde es el conjunto de todas las distribuciones conjuntas con marginales y .

La distancia total utilizada para la aceptación de ABC fue una suma de los dos estadísticos:

donde es la distancia de Wasserstein media por glándula. La contribución relativa de los dos componentes se gestionó mediante pyabc.AdaptiveAggregatedDistance, que reescala cada componente en cada generación de ABC-SMC para que ambos contribuyan aproximadamente por igual al umbral de aceptación. Este esquema adaptativo elimina la necesidad de un parámetro de ponderación fijo y se adapta a los cambios de escala a medida que la distribución posterior se concentra durante las generaciones sucesivas.

Matriz de distancia interglandular.

Para un tumor que consiste en glándulas muestreadas con matrices fCpG promedio, la matriz de distancia interglandular D se define como:

donde es el número de sitios fCpG. La diferencia al cuadrado enfatiza las grandes divergencias interglandulares y reduce la influencia de pequeñas fluctuaciones estocásticas. Dos matrices de distancia interglandular, y , se comparan utilizando la norma de Frobenius normalizada de su diferencia:

donde el factor de 1/2 tiene en cuenta la simetría de la matriz de distancia.

Distancia de distribución fCpG por glándula.

Para capturar información sobre la forma y el sesgo de las distribuciones fCpG individuales de cada glándula, lo que codifica las tasas relativas de metilación y desmetilación, también calculamos la distancia de Wasserstein promedio entre las distribuciones fCpG correspondientes de las glándulas en los datos simulados y observados:

donde es el conjunto de todas las distribuciones conjuntas con marginales y , y .

La distancia total utilizada para la aceptación de ABC fue una suma de las dos estadísticas:

donde es la distancia de Wasserstein promedio por glándula. La contribución relativa de los dos componentes se gestionó mediante pyabc.AdaptiveAggregatedDistance, que reescala cada componente en cada generación de ABC-SMC para que ambos contribuyan aproximadamente por igual al umbral de aceptación. Este esquema adaptativo elimina la necesidad de un parámetro de ponderación fijo y se adapta a los cambios de escala a medida que la distribución posterior se concentra durante las generaciones sucesivas.

Comprobaciones predictivas posteriores.

Realizamos comprobaciones predictivas posteriores para evaluar si el modelo inferido reproduce adecuadamente los datos observados. Para cada tumor, se muestrearon 100 conjuntos de parámetros de la distribución posterior de ABC de la generación final. Para cada muestreado, se realizó una simulación de tumor utilizando methdemon con la misma configuración que la ejecución de inferencia correspondiente (8 glándulas muestreadas, 2 lados). Se calculó la matriz de distancia interglandular simulada para cada réplica, y se comparó la distribución de las distancias simuladas por pares con los valores observados. Se consideró que un ajuste del modelo era adecuado si la distancia interglandular media observada se encontraba dentro del 90% central de la distribución predictiva posterior. La cobertura por pares se definió como la fracción de distancias interglandulares observadas que se encuentran dentro del intervalo predictivo posterior del 90% elemento por elemento.

Para los 10 tumores del estudio, la distancia interglandular observada se encontró dentro del 90% central de la distribución predictiva posterior (Fig. S3). La cobertura por pares (la fracción de distancias interglandulares que se encuentran dentro del intervalo predictivo posterior del 90%) superó el 96% para todos los tumores (rango: 96-100%). Las distancias por pares dentro del mismo lado mostraron una cobertura ligeramente menor (75-100%) que las distancias entre lados (100% para todos los tumores), lo que sugiere que el modelo captura la magnitud general y la estructura jerárquica de la divergencia fCpG, aunque la variación intra-lado puede estar subestimada en una minoría de los tumores.

Recuperación de parámetros sintéticos.

Para validar que el flujo de trabajo de inferencia recupera los parámetros identificables, se realizaron experimentos de recuperación de parámetros sintéticos. Se extrajo un conjunto de vectores de parámetros verdaderos de la distribución previa, y para cada se simuló un tumor sintético utilizando methdemon. Luego, se aplicó el flujo de trabajo completo de inferencia ABC-SMC a cada conjunto de datos sintéticos para obtener una distribución posterior . La recuperación se evaluó comparando la mediana posterior con el valor verdadero para cada parámetro, y reportando la fracción de experimentos en los que el valor verdadero se encontraba dentro del intervalo creíble del 90% (Fig. S2). Las tasas de epimutación y la tasa de fisión por glándula se recuperaron de manera confiable; y no, lo que es consistente con el análisis de sensibilidad descrito anteriormente.

Estimación de la edad mitótica.

La edad mitótica de cada tumor se estimó a partir de la distribución posterior inferida de la tasa de fisión de la siguiente manera.

Bajo el algoritmo de Gillespie, todas las tasas en methdemon se expresan por división celular. La tasa de fisión por glándula es , donde es la tasa de fisión por célula inferida y es la capacidad de carga del deme. Bajo un crecimiento exponencial a partir de una sola glándula fundadora, el número esperado de glándulas después de divisiones celulares es , lo que da

donde es el número total de glándulas en el tumor en el momento de la resección.

Dado que los recuentos de glándulas no están disponibles directamente, se estimó a partir del diámetro tumoral aproximado asumiendo una geometría esférica:

donde se derivó de las dimensiones típicas de las glándulas de los tumores colorrectales: aproximadamente 70 μm de diámetro y 450 μm de profundidad [^[13], ^[26]].

Para cada muestra posterior de , calculamos , lo que da como resultado una distribución posterior completa de la edad mitótica. La conversión a tiempo calendario se realizó dividiendo por las tasas de división de las células madre cancerosas (CSC) asumidas, que abarcan desde una división por semana hasta una división por mes, siguiendo las estimaciones publicadas [^[26], ^[27]].

El estimador debe interpretarse como la edad mitótica simétrica de la fase de expansión: methdemon contiene solo divisiones simétricas, por lo que las contribuciones asimétricas y post-expansión se absorben en las tasas inferidas en lugar de resolverse como divisiones adicionales.

Como comprobación cruzada independiente, también estimamos las edades mitóticas a partir de las distribuciones fCpG intra-glandulares utilizando una aproximación de dos estados a la dinámica de metilación de la cadena de Markov descrita en [^[2]]. Para cada glándula, se calculó la polarización media, definida como

en todos los loci fCpG (Fig. S6A). Esta cantidad mide qué tan lejos se ha desviado la matriz de metilación de los estados completamente metilados o no metilados. Tratando cada locus como un sistema de dos estados en el que los eventos de metilación y desmetilación ocurren independientemente en cada división celular, la polarización disminuye geométricamente, y la edad mitótica se puede estimar como

donde es la polarización máxima observada en todo el estudio y son las tasas de epimutación medianas posteriores para ese tumor (Fig. S6B).

Análisis de sensibilidad del recambio.

Para evaluar si la inclusión de un recambio en estado estacionario alteraría sustancialmente los patrones de divergencia interglandular en los que se basa la inferencia, realizamos una búsqueda en la cuadrícula en el espacio de parámetros del modelo en cuatro longitudes diferentes de la fase de recambio post-crecimiento (0%, 15%, 30%, 50%).

Se extrajeron 200 conjuntos de parámetros de la distribución previa de la inferencia mediante el muestreo de hipercubo latino. Para cada conjunto de parámetros, se ejecutó methdemon en las cuatro fracciones de recambio utilizando la misma semilla aleatoria, de modo que la fase de expansión fuera idéntica y la única diferencia fuera la duración del período de recambio posterior. Para cada ejecución, se calculó la matriz de distancia interglandular de 8 × 8 y se resumió mediante (i) la distancia media por pares, (ii) la relación de distancias entre lados y dentro del mismo lado, y (iii) la distancia de Frobenius entre la matriz de distancia en cada fracción de recambio y la línea de base del 30%.

Análisis estadístico.

Todas las pruebas estadísticas fueron bicaudales con un umbral de significación de . Dado el pequeño tamaño de la cohorte, se informan los tamaños del efecto junto con los valores p y los resultados se interpretan con la precaución adecuada con respecto al poder estadístico. La comparación entre las distancias dentro del mismo lado y entre lados se evaluó utilizando una prueba de rango con signo de Wilcoxon en los 10 tumores. La correlación de rango de Spearman se utilizó para evaluar la relación entre los valores de parámetros medianos posteriores inferidos y las variables clínicas continuas (tamaño del tumor, edad del paciente). Se utilizaron pruebas de Mann-Whitney para comparar los valores de parámetros inferidos entre los grupos categóricos (sexo, adenoma frente a carcinoma). La consistencia de las tasas de epimutación inferidas entre los pacientes se cuantificó utilizando el coeficiente de variación (CV) de las medianas posteriores. Todos los análisis estadísticos se realizaron en Python 3.12 utilizando scipy v1.17.1 (módulo scipy.stats).