The Pennes bioheat equation with Caputo fractional derivative applied to the thermal treatment of ductal breast cancer

LATIN-AMERICAN JOURNAL OF COMPUTING (LAJC), Vol XIII, Issue 1, January - June 2026 La ecuación bio-calor de Pennes con derivada fraccionaria de Caputo aplicada al tratamiento térmico del cáncer ductal de mama The Pennes bioheat equation with Caputo fractional derivative applied to the thermal treatment of ductal breast cancer Eder Linares Vargas Universidad del Atlántico Facultad de Educación: Licenciatura en Matemática y Ciencias Naturales Barranquilla Colombia [email protected] https://orcid.org/0000-0002-4697-5773 Abstract—This article examines the Pennes bioheat equation in both its classical form and its extension using the Caputo fractional derivative to model tumor heating through magnetic hyperthermia with SPIONs. In the classical model (α = 1.0), simulations reach and maintain temperatures above 42 °C, consistent with the clinical and experimental results of Caizer et al., where nanoparticles raise and stabilize tissue within the therapeutic range. When incorporating the fractional derivative (α < 1.0), thermal memory effects emerge, allowing a more realistic description of tissue dynamics. Although the explicit L1 method exhibits numerical instability, the implicit L1 method provides stable and physically coherent solutions, showing slower and more localized heating for fractional orders, as expected in tissues with delayed diffusion. These fractional results computationally correspond to the three-dimensional simulations of Rahpeima & Lin, which report non-monotonic temperature patterns and diffusion dependent on SPION concentration. Overall, the implicit L1 method validates both the experimental behavior observed by Caizer and the numerical dynamics reported by Rahpeima & Lin, demonstrating that the fractional approach is promising for modeling tumor hyperthermia when stable numerical schemes are employed. Keywords— Penne’s equation, Magnetic hyperthermia, Caputo fractional derivative and Superparamagnetic nanoparticles (SPIONs) I. INTRODUCCIÓN El cáncer de mama constituye uno de los mayores desafíos para la salud pública, tanto a nivel mundial como en Colombia. Cada año se reportan más de un millón de casos nuevos en el planeta, lo que lo convierte en la neoplasia maligna más común entre las mujeres. En el país, su incidencia ha venido aumentando de manera constante, situación vinculada a elementos como el incremento de la esperanza de vida, la alta exposición a diversos factores de riesgo y las brechas en el acceso a servicios de detección temprana. La evidencia científica señala que esta enfermedad se presenta con mayor frecuencia en mujeres entre los 40 y 60 años, aunque puede manifestarse en cualquier momento de la vida. Particularmente, se ha establecido una correlación significativa entre ciertas características clínicas y la probabilidad de recurrencia del carcinoma ductal no metastásico. Tal como lo evidencian [1], este hallazgo refuerza la necesidad de implementar enfoques terapéuticos individualizados que consideren la variabilidad biológica y clínica de cada paciente. En la actualidad, el cáncer de mama es la primera causa de muerte por cáncer en las mujeres en Colombia, con una tasa de mortalidad que equivale a cerca del 15 % de los decesos por enfermedades oncológicas en mujeres mayores de 20 años. Aunque la incidencia nacional aún no alcanza los niveles reportados en países desarrollados —donde una de cada diez mujeres será diagnosticada con esta enfermedad a lo largo de su vida—, el crecimiento sostenido de casos, especialmente en zonas urbanas, evidencia una tendencia preocupante que pone de manifiesto el impacto clínico y social de esta patología [2]. Por otro lado, un gran número de pacientes que se han sometido a cirugías radicales, como la mastectomía, corren riesgos considerables de complicaciones postoperatorias. Según [3] dentro de estas se incluye la linfedema, algunas infecciones, dolor intenso o crónico y limitaciones funcionales de los miembros superiores. Estas consecuencias afectan la calidad de vida de las pacientes, lo que evidencia la importancia de mejorar los enfoques terapéuticos y de impulsar tratamientos menos invasivos y con efectos secundarios reducidos. En regiones con sistemas de salud más consolidados, como Europa Occidental, la implementación sistemática de programas de tamizaje mediante mamografía ha permitido diagnósticos más tempranos, lo cual se traduce en mejores tasas de sobrevida. En Colombia, pese a los esfuerzos por aumentar la disponibilidad de esta herramienta diagnóstica, continúan existiendo barreras significativas de acceso, especialmente en zonas rurales y en regiones con escasa oferta de servicios de salud, lo que reduce su efectividad en la detección temprana y el control de la enfermedad [4]. B. Last Name, G. Thomas, and D. Botto-Tobar, “Title of the article in English”, Latin-American Journal of Computing (LAJC), vol. 13, no. 1, 2026. Paralelamente, el uso terapéutico del calor en el tratamiento del cáncer mediante nanopartículas superparamagnéticas de óxido de hierro (SPIONs), como Fe₃O₄ o γ-Fe₂O₃ recubiertas con materiales biocompatibles (dextrano, PEG o sílice), ha cobrado gran relevancia en oncología. Estas nanopartículas, aprovechando el efecto EPR (Enhanced Permeability and Retention), son capaces de acumularse selectivamente en tejidos tumorales, permitiendo una liberación localizada de calor cuando son activadas por campos magnéticos alternos [5]. Este enfoque emergente fue impulsado por los estudios de [6], quienes evidenciaron los beneficios clínicos de la hipertermia en oncología, aunque también señalaron las limitaciones de métodos tradicionales como la diatermia en presencia de inflamación ganglionar, por el riesgo de hiperemia, obstrucción venosa y necrosis tisular. En respuesta a estas limitaciones, [7] propusieron el uso de nanopartículas magnéticas como una fuente de calor controlada, promoviendo así un abordaje más preciso y seguro para el tratamiento de tumores y metástasis [8]. II. MODELIZACIÓN MATEMÁTICA Y PERSPECTIVA FRACCIONARIA El creciente interés científico por la hipertermia magnética ha dado lugar a múltiples investigaciones in vitro e in vivo, acompañadas del desarrollo de modelos matemáticos que permitan predecir con precisión la distribución térmica en los tejidos afectados. Entre estos, el modelo más ampliamente utilizado está la ecuación bio-calor propuesta por [9] La cual, debido a su simple formulación —lo que la hace atractiva y correlaciona muy bien con datos experimentales—, ha sido adoptada como base para numerosas simulaciones térmicas en tejidos biológicos. Sin embargo, cuando se aplica este modelo al contexto clínico del cáncer ductal mamario, emergen importantes limitaciones. La ecuación de Pennes clásica asume un comportamiento térmico lineal, local y sin memoria, lo cual no refleja adecuadamente la naturaleza dinámica y estructuralmente heterogénea del tejido tumoral mamario, especialmente bajo exposiciones térmicas repetidas o moduladas. Esta falta de representación histórica en la conducción térmica es particularmente crítica en entornos biológicos donde la respuesta al calor puede presentar retardos o acoplamientos espacio-temporales no triviales. Ante esta problemática, en el presente trabajo se explora la formulación fraccionaria de la ecuación de Pennes, incorporando derivadas de Caputo. Esta extensión permite capturar efectos de memoria térmica y modelar de forma más realista la evolución de temperatura en tejidos tumorales. El enfoque fraccionario no solo mejora la fidelidad de las simulaciones térmicas, sino que también abre nuevas posibilidades para el diseño personalizado de tratamientos basados en hipertermia magnética, tal como lo sugieren estudios recientes como los de [10] III. MATERIALES Y MÉTODOS Se inicia con la solución numérica de la ecuación de biocalor propuesta por Pennes, se construyen los códigos en Python y se analizan las gráficas obtenidas. Seguidamente se resuelve numéricamente la ecuación de Pennes, pero esta vez desde la perspectiva del cálculo fraccionario propuesto por [11], se analizan las gráficas derivadas de los algoritmos en Python y se discuten los resultados contrastándolos con los obtenidos en el modelo de Pennes. IV. ECUACIÓN DE BIO-CALOR DE PENNES La ecuación de bio-calor de Pennes, formulada, describe cómo se distribuye y transfiere el calor en los tejidos biológicos humanos. Es ampliamente utilizada en estudios de termoterapia, hipertermia para tratamiento del cáncer, crioterapia y simulaciones médicas. V. FORMULACIÓN GENERAL DE LA ECUACIÓN DE PENNES ( ) ( ) b b b b m ext T c k T c w T T Q Q t  =   + − + +  , donde  es la densidad del tejido medido en (Kg/m3) c es el calor especifico del tejido J/ (Kg k) T es la temperatura del tejido (°C o K) k es la conductividad térmica del tejido (W/ mk) b  es la densidad de la sangre (Kg/m3) b c calor especifico de la sangre J/ (Kg k) b w es la tasa de perfusión sanguínea (1/seg) b T temperatura de la sangre arterial (°C o K) m Q es la tasa de generación metabólica de calor (W/m3) ext Q son las fuentes externas de calor (W/m3) Pennes propone una EDP, no lineal en ocasiones, parabólica y homogénea, con la que describe el efecto de la temperatura en tejidos biológicos. Es decir, mediante un enfoque matemático que integra diversos mecanismos de transferencia de calor. La ecuación, en su lado izquierdo representa la tasa de energía interna del tejido. El primer término del lado derecho modela la conducción térmica interna, reflejando la difusión del calor dentro del tejido debido a diferencias de temperatura. El segundo término introduce el efecto de la perfusión sanguínea, donde la sangre actúa como un vehículo de transporte térmico, intercambiando calor entre los vasos sanguíneos y el tejido circundante. El tercer término considera el calor generado por el metabolismo celular, una fuente constante en tejidos vivos. Finalmente, el cuarto término permite incorporar fuentes de calor externas, como las generadas por nanopartículas magnéticas en tratamientos de hipertermia, lo cual resulta especialmente relevante en el contexto del cáncer ductal de seno. La ecuación de bio-calor de Pennes ha convertido en una herramienta muy versátil en el campo de la bioingeniería, debido a sus múltiples aplicaciones en oncología, ingeniería de tejidos y tecnologías biomédicas. En el ámbito clínico, se utiliza ampliamente para simular tratamientos por hipertermia, permitiendo calcular con precisión cómo se eleva la temperatura en un tumor durante la exposición a LATIN-AMERICAN JOURNAL OF COMPUTING (LAJC), Vol XIII, Issue 1, January - June 2026 microondas, ultrasonido o campos magnéticos, mejorando así la efectividad terapéutica. Asimismo, esta ecuación es fundamental para modelar procedimientos de crioterapia localizada, donde se analiza la propagación del frío al aplicar agentes como nitrógeno líquido sobre tejidos afectados. En el diseño de implantes térmicamente activos o inteligentes, permite predecir cómo estos dispositivos alteran el entorno térmico del tejido que los rodea. También se aplica en la simulación de quemaduras, ayudando a comprender la propagación del calor en lesiones de primer y segundo grado. Finalmente, en medicina deportiva, la ecuación de Pennes facilita el estudio de la temperatura muscular durante el ejercicio, considerando el efecto de la perfusión sanguínea, lo cual es útil para crear planes de entrenamiento personalizado. VI. ECUACIÓN DE PENNES INTRODUCIENDO LA HIPERTERMIA MAGNÉTICA La hipertermia magnética es una técnica terapéutica empleada en oncología que se basa en elevar la temperatura de los tejidos tumorales mediante el uso de nanopartículas magnéticas (como los óxidos de hierro) previamente incorporadas en el tumor. Estas partículas se excitan mediante un campo magnético alterno, generando calor de forma localizada. Para este caso, la ecuación de Pennes se reescribe de esta forma: ( ) ( ) b b b b m mag T c k T c w T T Q Q t  =   + − + +  donde mag Q representa el calor generado por las nanopartículas. El calor generado por las NPs, se puede modelar a través de la ecuación ,, 2 0mag Q H f  = donde 0  es la permeabilidad magnética del vacío ,,  parte imaginaria de la susceptibilidad magnética H amplitud del campo magnético aplicado f frecuencia del campo magnético En general, el valor de mag Q depende del tipo de nanopartícula, la concentración y la frecuencia del campo. TABLE I. PARÁMETROS TÍPICOS EN TERAPIAS DE HIPERTERMIA MAGNÉTICA Parámetro Valor típico / Rango estimado Notas técnicas Tamaño de partícula 10 – 50 nm Influye en el super paramagnetismo [1213] Densidad del material ρ ~5200 kg/m³ (Fe₃O₄) Determina la cantidad de masa inyectada [14] Calor específico c 670 J/(kg·K) Valor específico para Fe₃O₄. Japan Atomic Energy Agency. (s. f.). Thermodynamic data Fe₃O₄ Conductividad térmica k ~6 W/(m·K) Superior a la del tejido humano (~0.5 W/m·K) [15] SAR (Specific Absorption Rate) 10 – 500 W/g Depende del campo magnético, tamaño, recubrimiento, etc [13;16] Campo magnético aplicado H 10 – 30 kA/m Típico en terapias clínicas [17] Frecuencia del campo ƒ 100 – 500 kHz Regulado por normativas de seguridad [17;18]) Concentración en tumor 1 – 10 mg/mL Varía según método de administración y acumulación [16] A. Calor generado: SAR y Qmag mag Q SAR C= donde C representa la concentración de SPIONs en el tejido (g/m3) mag Q es expresado en W/m3, lo cual es compatible con la ecuación de Pennes Por ejemplo, con un SAR = 200 W/g y una concentración C= 5mg/ml = 5x103 g/m3, como mag Q SAR C= entonces 6 3 200 5000 1 10 mag W Qm =  =  . Esto genera un calor altamente focalizado, que es capaz de aumentar la temperatura del tumor a 42-45°C, ideal para destruir las células cancerígenas sin destruir el tejido sano VII. SIMULACIÓN DEL CALENTAMIENTO TUMORAL MEDIANTE LA ECUACIÓN DE PENNES En el contexto clínico, el carcinoma ductal invasivo es un tipo común de cáncer localizado en la región mamaria. Para su tratamiento, se utiliza una técnica de hipertermia magnética basada en la inyección de nanopartículas súper paramagnéticas (SPIONs). Estas partículas se concentran en el tumor y, al aplicar un campo magnético variable, generan calor de forma localizada. El objetivo terapéutico es aumentar la temperatura del tejido tumoral entre 42 y 45 °C. Esta elevación térmica debe mantenerse durante un periodo de 30 a 60 minutos para inducir daño celular selectivo. El modelo que describe este proceso utiliza la ecuación bio-calor de Pennes. Esto permite simular y hasta optimizar la distribución térmica en el tejido afectado. Partiendo de la ecuación: ( ) ( ) b b b b m mag T c k T c w T T Q Q t  =   + − + +  , teniendo en cuenta la tabla 2. TABLE II. VALORES REFERENCIALES PARA TEJIDO MAMARIO Y TUMOR Parámetro Tumor ductal (valor estimado) Unidad ρ 1050 kg/m³. c 3600 J/(kg·K) k 0.5 W/(m·K) ρb 1060 kg/m³ cb 3770 J/(kg·K) ωb 0.002 1/s Tb 37 °C Qm 10,000 W/m³ Qmag 1 × 10⁶ (SAR × concentración) W/m³ B. Last Name, G. Thomas, and D. Botto-Tobar, “Title of the article in English”, Latin-American Journal of Computing (LAJC), vol. 13, no. 1, 2026. Suponemos un modelo esférico de tumor con las siguientes características: • Radio del tumor: R=1cm • Dominio en 3D o radial simétrico • Condiciones de contorno: temperatura fija en el exterior del tejido sano (T = 37 °C) A. Objetivo del modelo • Determinar el perfil térmico en el tumor. • Comprobar si se alcanza y mantiene una temperatura terapéutica. • Analizar el efecto de la perfusión sobre la disipación térmica. • Estudiar el efecto de usar la derivada fraccionaria de Caputo en la Ecuación de Pennes VIII. MÉTODO NUMÉRICO UTILIZADO IX. DIFERENCIAS FINITAS ESPACIALES Se utilizan aproximaciones por diferencias finitas centradas para calcular el laplaciano (es decir, la difusión del calor). Esto se hace con mallas uniformes (espaciado constante) sobre el dominio. 2 11 22 2 i i i T T T T xx +−  − +   , una vez discretizado el espacio, la ecuación de Pennes se transforma en un sistema de EDOs. Luego estas, se resuelven mediante Runge-Kutta. A continuación, se muestran las gráficas derivadas de los algoritmos hechos en Python. Fig. 1. La figura 1, muestra cómo varía la temperatura a lo largo del radio del tumor mamario después de 600 segundos (10 minutos) de tratamiento con hipertermia magnética. Se observa un gradiente térmico desde el centro hacia la periferia: la temperatura es más alta en la zona central del tumor y disminuye conforme se aproxima al tejido sano, donde se mantiene constante a 37 °C debido a la condición de frontera. Esta distribución es consecuencia del calentamiento interno generado por las nanopartículas magnéticas, así como del balance térmico con el entorno y la perfusión sanguínea. La línea roja discontinua representa el umbral terapéutico de 42 °C, temperatura mínima requerida para inducir daño térmico en células tumorales. Gran parte del tejido tumoral, específicamente en su núcleo, se supera dicho umbral, lo que es un claro indicio de que el modelo simulado es ideal para lograr el efecto terapéutico esperado. Sin embargo, en regiones cercanas al borde (más allá de 8–9 mm), la temperatura desciende por debajo del umbral, lo que sugiere que la hipertermia podría no ser completamente uniforme sin ajustes en el diseño del tratamiento. Fig. 2. La figura 2 muestra cómo varía la temperatura a lo largo del tiempo en el centro del tumor (r ≈ 0 mm) durante los 10 minutos de la simulación. Desde una temperatura inicial de 37 °C, se observa un incremento progresivo debido al efecto acumulativo del calor generado por el campo magnético sobre las nanopartículas. En los primeros 200– 300 segundos, el aumento es más pronunciado, seguido por una estabilización cerca de 44 °C, lo que indica un nuevo equilibrio térmico entre generación, conducción, y disipación de calor. El hecho de que la temperatura en el centro del tumor supere y mantenga un nivel superior al umbral terapéutico de 42 °C durante varios minutos es clínicamente relevante, debido a que este es el rango térmico esperado que permite la necrosis selectiva del tejido tumoral, logrando no afectar significativamente el tejido sano circundante. Este comportamiento refleja que el modelo es funcional para predecir la respuesta térmica en hipertermia, permitiendo simular distintos escenarios antes de la aplicación clínica real. Fig. 3. La figura 3 indica la biodistribución de temperatura en un tejido mediante el modelo bio-térmico de Pennes, que considera conducción, perfusión sanguínea y generación de calor; además de las temperaturas entre 37.00 °C y 37.10 °C, las cuales son clínicamente coherentes. Esto sugiere una fuente de calor suave y controlada. En cuanto a la distribución térmica, el máximo de temperatura se encuentra en el centro del dominio, como se esperaba debido a la fuente térmica gaussiana localizada ahí. Los bordes se mantienen en 37 °C, gracias a las condiciones de frontera fijas correctamente implementadas. La forma circular simétrica del gradiente térmico demuestra que la simulación es estable y físicamente realista. En cuanto al rango, este va de 0 a 10 mm en ambos ejes y representa un área de tejido de 1 cm², que es adecuado para simulaciones locales de hipertermia, láser terapéutico o ablación térmica. LATIN-AMERICAN JOURNAL OF COMPUTING (LAJC), Vol XIII, Issue 1, January - June 2026 Fig. 4. La figura 4 muestra con gran claridad una distribución térmica centrada y clínicamente significativa. En ella se observan contornos isotérmicos que delimitan regiones fisiológicamente relevantes: temperaturas superiores a 38 °C asociadas a activación metabólica, zonas de 39–40 °C vinculadas a vasodilatación y sensibilización celular, y una región interna que supera los 41–42 °C, indicativa de una hipertermia terapéutica efectiva. El gradiente de temperatura es simétrico y suave, confirmando que el calor se ha concentrado en el centro de la imagen con un patrón de disipación bien controlado. Esta imagen muestra una zona de tratamiento térmico altamente focalizado, donde se evidencia un núcleo sobrecalentado que está rodeado por anillos concéntricos donde se refleja una menor temperatura. La escala de color y los contornos permiten interpretar fácilmente la eficacia del calentamiento, mostrando que se ha alcanzado el rango térmico adecuado para inducción de necrosis celular selectiva. En conjunto, esta visualización constituye una representación gráfica precisa, robusta y clínicamente aplicable de un escenario de hipertermia oncológica localizada [19]. Fig. 5. La figura 5 muestra un Mapa de daño térmico acumulado (Ω ≥ 1= necrosis), lo que indica una simulación de la extensión del daño tisular causado por la exposición al calor. En esta representación bidimensional de 10 mm×10 mm, los colores varían desde el negro (menor daño) hasta el amarillo brillante (mayor daño). La barra de color a la derecha, titulada "log10(Ω) (Daño térmico acumulado)", muestra que los valores van desde −3.8 hasta −2.2. Un valor de Ω≥1 se asocia con necrosis, lo que implica que las zonas con los valores de log10(Ω) más altos (más cercanos a cero o positivos, aunque en este gráfico solo se muestran valores negativos) representan las áreas donde el daño es suficiente para causar la muerte celular. Se observa una concentración de daño severo en el centro del mapa (amarillo/naranja brillante), disminuyendo a medida que uno se aleja del punto central, lo que es consistente con una aplicación de calor localizada, y sugiere que el objetivo de la terapia de hipertermia, que es inducir daño térmico en una región específica, se está logrando. En menester aclarar que Ω, representa la integral de daño térmico. Esta integral se fundamenta en la cinética de Arrhenius, la cual explica cómo la temperatura influye en la velocidad de una reacción química, en este caso, la desnaturalización de proteínas y la muerte celular. En esencia, la integral de daño Ω representa la probabilidad acumulada de daño irreversible en el tejido. () 0 () E tRT t Ae d    −  =  donde Ω(t) representa el daño en el tumor en un tiempo t, A es el factor de frecuencia, el cual es constante y depende del tipo de tejido. Corresponde a la energía de activación, es decir, la cantidad mínima de energía requerida para que se produzca la desnaturalización de proteínas o el daño celular, y es propia de cada tipo de tejido. R representa la constante de los gases ideales. T(τ) es la temperatura absoluta del tejido (expresada en Kelvin) en el instante τ. Esta ecuación muestra el daño acumulado en un tejido en función exponencial de la temperatura y función lineal del tiempo. Es decir, para pequeños aumentos en la temperatura, se induce un aumento drástico en la tasa de daño, y una exposición prolongada a una temperatura específica también aumenta el daño acumulado. Fig. 6. En la figura 6, la imagen izquierda representa un corte transversal del tejido simulado (en el plano z = 5 mm) al finalizar 30 minutos de exposición térmica, y muestra la distribución de temperaturas alcanzadas. En ella, los colores más claros muestran zonas donde la temperatura es mayor, alcanzando hasta 42 °C en el centro, mientras que los bordes permanecen cercanos a los 37 °C (temperatura corporal). Se incluyen contornos que marcan los umbrales terapéuticos: 41.5 °C (azul claro), asociado a hipertermia útil en tratamientos oncológicos, y 42.0 °C (en verde lima), vinculado a daño térmico más severo. La figura permite identificar las regiones más calientes, esenciales para planificar terapias como la ablación térmica mediante SPIONs, con láser o radiofrecuencia. La imagen de la derecha refleja el daño térmico acumulado en el mismo corte transversal, mediante el modelo de Arrhenius, el cual se expresa como log₁₀(Ω), donde Ω es un indicador del daño celular. Una zona alcanza necrosis irreversible cuando Ω ≥ 1 (log₁₀(Ω) = 0), lo que permite diferenciar entre daño parcial y completo. Las zonas claras indican mayor daño térmico, mientras que las oscuras reflejan tejido no afectado. Esta imagen es crucial para validar si el tratamiento fue suficiente, estimar el margen de ablación y evaluar riesgos a tejidos circundantes. Ambas imágenes se complementan: la primera muestra las temperaturas alcanzadas y la segunda el efecto biológico acumulado, recordando que una alta temperatura durante poco tiempo puede no ser suficiente para inducir daño térmico permanente X. CÁLCULO FRACCIONARIO EN LA ECUACIÓN DE PENNES Integrar cálculo fraccionario en la ecuación de Pennes es una estrategia poderosa para modelar procesos térmicos anómalos, como los que pueden ocurrir en tejidos biológicos heterogéneos como un tumor ductal mamario. Para tratar de modelar el transporte térmico en un carcinoma ductal mamario, la derivada fraccionaria de Caputo es ideal. Esto se debe a que permite incorporar memoria térmica, que es relevante en tejidos con respuesta retardada al calor. Es compatible con condiciones iniciales clásicas como: B. Last Name, G. Thomas, and D. Botto-Tobar, “Title of the article in English”, Latin-American Journal of Computing (LAJC), vol. 13, no. 1, 2026. ( ) ( ) 0 01 ( ,0) ( ) b b b b m mag T c k T c w T T Q Q t T x T x      =   + − + +    = 1 1. Si Modelo Clásico de Pennes Si Modelo con efecto memoria   =→ → XI. MÉTODOS NUMÉRICOS EN PENNES CON CAPUTO XII. MÉTODO L1 EXPLÍCITO L1 en la ecuación de Caputo puede representar un coeficiente de diferencias finitas, usado en la aproximación numérica de la derivada fraccionaria. ( ) ( ) 111 1 0 11 (2 ) njj j TT Tn j n j tt    −−− + = −  − − −   −   Fig. 7. La figura 7 muestra cómo se distribuye la temperatura dentro del tumor al finalizar los 600 segundos de simulación para diferentes órdenes fraccionarios α. Se observa que cuanto mayor es el valor de α, más alta es la temperatura alcanzada en el centro del tumor, y más uniforme es la distribución térmica. En particular, para α=1.0, correspondiente al modelo clásico de Pennes, se logra alcanzar exactamente los 42 °C necesarios para un tratamiento térmico efectivo (umbral terapéutico). En cambio, para valores fraccionarios menores de α, el calor no se propaga tan eficientemente hacia los bordes del tumor, lo que refleja el comportamiento de memoria y retardo característico de los modelos fraccionarios. Fig. 8. En la figura 8 se representa la evolución temporal de la temperatura justo en el centro del tumor (r=0) para cada valor de α. Se aprecia que todos los casos parten desde la temperatura basal de 37 °C, y conforme pasa el tiempo, la temperatura se incrementa gradualmente. Cuanto mayor es el valor de α, más rápida y eficaz es la acumulación de calor. La evolución térmica es más lenta en los modelos fraccionarios (α<1), lo cual pone de manifiesto cómo la dinámica de transferencia térmica se ve afectada por el orden de la derivada de Caputo. Evidentemente el modelo clásico con α=1.0 alcanza los 42 °C necesarios para obtener la hipertermia terapéutica. A continuación, se presenta una tabla que recoge el análisis para cada valor de α. TABLE III. ANÁLISIS PARA CADA VALOR DE Α AL CABO DE 10 MINUTOS DE SIMULACIÓN TÉRMICA α Temp. Máx (°C) Temp. en el centro (°C) Temp. en el borde (°C) Temp. Mín (°C) 0.6 40.00 40.00 38.10 38.10 0.7 40.50 40.50 38.29 38.29 0.8 41.00 41.00 38.47 38.47 0.9 41.50 41.50 38.66 38.66 1.0 42.00 42.00 38.84 38.84 Sólo para α=1.0 se alcanzó el umbral terapéutico (42°C) en el centro del tumor. La introducción de la derivada fraccionaria de Caputo se debe a la limitación del modelo clásico de Pennes, que asume una respuesta térmica instantánea y sin memoria. La derivada de Caputo, al incorporar la historia térmica del tejido (memoria), busca ofrecer una representación más realista de la transferencia de calor en tejidos biológicos complejos y heterogéneos, como los tumores. Esto es clave cuando se pretende estudiar el efecto de usar las derivadas fraccionarias. El desafío inmediato al introducir la derivada fraccionaria con el método L1 explícito es la inestabilidad numérica severa observada en las gráficas de evolución térmica para α<1.0. Las soluciones se disparan a valores irrealistas después de un corto periodo de tiempo. Esto no permite completamente cumplir con los objetivos de determinar los perfiles térmicos que son estables, ni comprobar el mantener la temperatura terapéutica, o analizar cómo se disipa el calor en el contexto del cálculo fraccionario [20]. Se observa una similitud en el patrón de inestabilidad para todos los valores de α<1.0, donde la inestabilidad se retrasa ligeramente a medida que α se acerca a 1.0. Esto sugiere que el problema de estabilidad es inherente al esquema explícito para derivadas fraccionarias y no a un valor particular de α. En los primeros instantes de la simulación, antes de que aparezca la inestabilidad, las curvas de Caputo y la de Pennes clásica muestran una tendencia similar de aumento de temperatura desde el valor inicial. Ambas predicen que la temperatura terapéutica de 42∘C puede ser alcanzada en el centro del tumor. Las diferencias cruciales que la derivada de Caputo debe introducir, tales como una respuesta más lenta o amortiguada debido a la memoria, son afectadas por la inestabilidad. Sin embargo, si el método fuera estable, se esperaría que los perfiles y la evolución temporal de Caputo difieran de los clásicos, reflejando una dinámica de calor más compleja y potencialmente más precisa para tejidos con propiedades viscoelásticas. La inestabilidad de los métodos explícitos para las EDPs fraccionarias subraya la necesidad crítica de emplear métodos numéricos implícitos (como el Método L1 Implícito y predictor corrector para obtener soluciones estables y físicamente significativas. Sin soluciones estables, el potencial del modelo fraccionario para estudiar la biotransferencia de calor y optimizar los tratamientos de hipertermia no puede ser explotados cabalmente. Aunque la LATIN-AMERICAN JOURNAL OF COMPUTING (LAJC), Vol XIII, Issue 1, January - June 2026 formulación de Caputo es teóricamente superior para ciertos fenómenos, su aplicación práctica depende directamente de la robustez del algoritmo numérico. Los resultados obtenidos de las simulaciones de la evolución térmica en el centro del tumor, contrastando el modelo clásico de Pennes (α=1.0) con la formulación fraccionaria de Caputo (para α<1.0), revelan aspectos fundamentales sobre la dinámica de la transferencia de calor en tejidos biológicos. Un objetivo fundamental del modelo es determinar el perfil térmico en el tumor y comprobar si se alcanza y mantiene una temperatura terapéutica. En este sentido, la gráfica muestra que, independientemente del orden fraccionario, la temperatura en el centro del tumor inicialmente se eleva por encima del umbral terapéutico de 42∘C, lo que indica que la fuente de calor (generada por las SPIONs) es lo suficientemente potente para iniciar el proceso de calentamiento deseado. Sin embargo, la estabilidad y la evolución a largo plazo de esta temperatura difieren drásticamente entre los enfoques [21-23]. La discusión más evidente surge al analizar la estabilidad numérica de las soluciones. Mientras que la curva correspondiente al modelo clásico de Pennes (α=1.0) exhibe un comportamiento estable y físicamente plausible, con un ascenso gradual de la temperatura hacia un plateau, las simulaciones con la derivada fraccionaria de Caputo (para α=0.6,0.7,0.8,0.9) dan evidencia de una notable inestabilidad, caracterizada por picos y caídas abruptas y sin control de la temperatura después de un periodo de tiempo muy corto. La inestabilidad es un indicativo de que el método numérico explícito L1 explicito, no es el adecuado para resolver la ecuación de Pennes con la derivada de Caputo bajo las condiciones de discretización utilizadas. Esto genera una discusión crítica sobre la elección del esquema numérico, destacando la necesidad de pasos de tiempo significativamente más pequeños o, preferiblemente, la implementación de métodos implícitos que son inherentemente más estables para ecuaciones diferenciales fraccionarias y sistemas rígidos. El efecto de usar la derivada fraccionaria de Caputo, induce la inestabilidad observada cuando α<1.0, lo que impide una comparación directa y significativa del comportamiento de memoria térmica que el modelo pretende introducir. Si bien se esperaba que la derivada fraccionaria permitiera modelar con mayor realismo la evolución de temperatura en tejidos con propiedades viscoelásticas o heterogéneas, los resultados actuales no permiten evidenciar este efecto de manera concluyente debido a la divergencia de las soluciones. Sin embargo, se puede discutir que, en los breves periodos de estabilidad inicial, las curvas fraccionarias muestran una dinámica de calentamiento que podría diferir sutilmente de la clásica, lo que, de ser resuelto con un método estable, podría dar lugar a perfiles térmicos y evoluciones temporales más complejos y potencialmente más precisos para la predicción del daño tisular. La discusión se extiende a la capacidad de comprobar si se alcanza y mantiene una temperatura terapéutica y analizar el efecto de la perfusión sobre la disipación térmica. Aunque las temperaturas iniciales superan los 42∘C, la inestabilidad de las soluciones fraccionarias impide determinar si estas temperaturas se mantendrían durante el periodo terapéutico deseado (30-60 minutos) o cómo la perfusión sanguínea influiría en la disipación del calor en un modelo con memoria. El modelo clásico de Pennes, al ser estable, sí permite inferir que la temperatura terapéutica se mantiene en el centro del tumor. La falta de estabilidad en el enfoque de Caputo resalta que, para cumplir los objetivos del modelo en un contexto fraccionario, es imperativo desarrollar o emplear esquemas numéricos que garanticen la convergencia y la estabilidad de las soluciones a lo durante la simulación. En síntesis, mientras que el modelo clásico de Pennes demuestra la viabilidad de la hipertermia magnética y permite evaluar el perfil térmico y la disipación, la introducción de la derivada fraccionaria de Caputo, aunque teóricamente muy prometedora para modelar la memoria térmica, se ve limitada por los desafíos que se presentan en cuanto a la estabilidad numérica de los esquemas explícitos. Superar esta limitación con métodos implícitos es el siguiente paso crucial para validar y aprovechar el potencial predictivo de los modelos fraccionarios en la bio-transferencia de calor. XIII. MÉTODO L1 IMPLÍCITO Con el objetivo de mejorar la estabilidad numérica del modelo fraccionario de Pennes, se incorpora la formulación del método L1 en su versión implícita, en sustitución del esquema explícito utilizado inicialmente. A diferencia del método explícito —cuya estabilidad está fuertemente condicionada por el tamaño del paso temporal y conduce a divergencias para órdenes fraccionarios α < 1— el método L1 implícito permite resolver el término de memoria fraccionaria como un sistema acoplado, lo que extiende de manera significativa la región de estabilidad del algoritmo. Este ajuste metodológico resulta esencial para obtener perfiles térmicos físicamente plausibles y evitar oscilaciones numéricas no deseadas, especialmente en contextos donde la dinámica térmica está dominada por gradientes pronunciados y efectos de memoria térmica propios de la formulación fraccionaria de Caputo. Adicionalmente, se incorpora un esquema predictor– corrector adaptado a derivadas fraccionarias como mecanismo complementario al L1 implícito, con el fin de mejorar la precisión global del modelo y reducir errores acumulativos asociados al tratamiento del término histórico. Este método, basado en la estructura Adams–Bashforth–Moulton fraccionaria, permite realizar una predicción inicial de la solución seguida de un paso corrector que refina el valor calculado, proporcionando un equilibrio adecuado entre estabilidad numérica y precisión. La combinación del L1 implícito con el método predictor–corrector constituye así una estrategia robusta para el tratamiento de modelos termofraccionarios, permitiendo obtener simulaciones más consistentes y abriendo la posibilidad de estudiar escenarios clínicos con parámetros más realistas sin comprometer la estabilidad del cálculo. La ecuación bio-calor fraccionaria que se está resolviendo (Caputo en el tiempo, 𝟎<𝜶≤𝟏) es: 𝝆𝒄 𝑪𝑫𝒕𝜶𝑻(𝐱,𝒕)=𝒌𝛁𝟐𝑻(𝐱,𝒕)+𝝆𝒃𝒄𝒃𝝎𝒃(𝑻𝒃−𝑻(𝐱,𝒕)) +𝑸(𝐱,𝒕), con condiciones iniciales 𝑇(𝐱,0)=𝑇0(𝐱) y condiciones de frontera ( 𝑇=37∘C en contorno exterior). B. Last Name, G. Thomas, and D. Botto-Tobar, “Title of the article in English”, Latin-American Journal of Computing (LAJC), vol. 13, no. 1, 2026. A. Discretización temporal L1 Se define una malla temporal 𝒕𝒏=𝒏𝚫𝒕. Para 𝟎<𝜶<𝟏, los pesos L1 se definen como 𝒘𝒋=(𝒋+𝟏)𝟏−𝜶−𝒋𝟏−𝜶, 𝒋=𝟎,𝟏,… La aproximación L1 para la derivada de Caputo en el instante 𝒕𝒏+𝟏 es 𝑪𝑫𝒕𝜶𝑻(𝒕𝒏+𝟏)≈ 𝚫𝒕−𝜶 𝚪(𝟐−𝜶)∑𝒘𝒋 𝒏 𝒋=𝟎 (𝑻𝒏+𝟏−𝒋−𝑻𝒏−𝒋), donde 𝑻𝒎 denota la solución aproximada en 𝒕𝒎. Esta versión estándar L1 incluye la diferencia de cada intervalo. (Peso 𝒘𝟎=𝟏𝟏−𝜶−𝟎𝟏−𝜶=𝟏.) Para la derivación del sistema lineal por paso temporal se sustituye por aproximación L1 en la ecuación y se evalúa el término de difusión y perfusión en el nivel implícito 𝑛+1: 𝜌𝑐 Δ𝑡−𝛼 Γ(2−𝛼)∑𝑤𝑗 𝑛 𝑗=0 (𝑇𝑛+1−𝑗−𝑇𝑛−𝑗) =𝑘∇2𝑇𝑛+1+𝜌𝑏𝑐𝑏𝜔𝑏(𝑇𝑏−𝑇𝑛+1) +𝑄𝑛+1. Agrupando los términos dependientes de 𝑇𝑛+1 se obtiene una ecuación lineal del tipo 𝒜𝑇𝑛+1=𝑏𝑛+1, • Coeficiente masa (por punto espacial): 𝐶0=𝜌𝑐 Δ𝑡−𝛼 Γ(2−𝛼)𝑤0=𝜌𝑐 Δ𝑡−𝛼 Γ(2−𝛼), porque 𝑤0=1. • Lado derecho (términos conocidos con tiempos ≤ 𝑛): 𝑏𝑛+1=𝜌𝑐 Δ𝑡−𝛼 Γ(2−𝛼)(𝑤0𝑇𝑛+∑𝑤𝑗 𝑛 𝑗=1 (𝑇𝑛+1−𝑗 −𝑇𝑛−𝑗))+𝜌𝑏𝑐𝑏𝜔𝑏𝑇𝑏+𝑄𝑛+1. • Operador espacial (discretizado por diferencias finitas centradas): con la matriz Laplaciana 𝐿 tal que ∇2𝑇𝑛+1→𝐿𝑇𝑛+1. Por lo tanto, la matriz del sistema es: 𝒜=𝐶0𝐼+𝜌𝑏𝑐𝑏𝜔𝑏𝐼−𝑘𝐿, y resolviendo 𝒜𝑇𝑛+1=𝑏𝑛+1. Es válido aclarar que, el único término que involucra incógnitas en el sumatorio L1 es el término con 𝑤0(𝑇𝑛+1− 𝑇𝑛) (ya incorporado en 𝐶0). Los demás sumandos contienen 𝑇𝑚 con 𝑚≤𝑛 y por tanto son conocidos al resolver el paso 𝑛+1. Esto es lo que convierte a L1 en implícito respecto a los términos espaciales y a la contribución principal temporal, formando un sistema lineal (no no lineal si Q no depende de 𝑇). B. Discretización espacial Se usa diferencias finitas centradas para ∇2 en la malla radial. ∇2𝑇𝑖≈𝑇𝑖+1−2𝑇𝑖+𝑇𝑖−1 Δ𝑥2. 𝐿 es la matriz Laplaciana con condiciones de contorno aplicadas (Dirichlet 𝑇=37∘C en frontera). XIV. PREDICTOR–CORRECTOR ADAMS–BASHFORTH– MOULTON FRACCIONARIO En este caso se usa L1 implícito como esquema base de estabilidad (resolviendo el sistema lineal 𝒜{𝑛+1}=𝑏{𝑛+1} en cada paso) y emplear un predictor–corrector ABM fraccionario para obtener un predictor 𝑇pred 𝑛+1 más preciso del término fuente no lineal. En presencia de términos lineales (como en Pennes con Q conocido y constantes), el L1 implícito por sí solo es suficiente y el predictor–corrector aporta precisión adicional (menor error temporal) y puede acelerar convergencia en casos no lineales. A. Forma práctica del predictor–corrector (esquema Diethelm-type, adaptado): • Predictor (Adams–Bashforth fraccionario explícito): calcula una estimación 𝑇𝑃𝑛+1 usando datos previos 𝑇0,…,𝑇𝑛 y evaluaciones 𝑓𝑘 del lado derecho espacial 𝑓(𝑇𝑘):=𝑘∇2𝑇𝑘+𝜌𝑏𝑐𝑏𝜔𝑏(𝑇𝑏− 𝑇𝑘)+𝑄𝑘. Fórmula conceptual: 𝑇𝑃𝑛+1=𝑇0+Δ𝑡𝛼 Γ(𝛼+1)∑𝛽𝑘(𝑛+1) 𝑛 𝑘=0 𝑓(𝑡𝑘,𝑇𝑘), donde los coeficientes 𝛽𝑘(𝑛+1) son los pesos adecuados del integrador fraccionario (estos coeficientes pueden obtenerse por regla de Newton–Cotes fraccionaria; en implementaciones práctica se usan tablas o fórmulas estándar. • Corrector (Adams–Moulton fraccionario implícito): mejora la predicción resolviendo implícitamente: 𝑇𝑛+1=𝑇0+Δ𝑡𝛼 Γ(𝛼+1)(𝛽0(𝑛+1)𝑓(𝑡𝑛+1,𝑇𝑛+1) +∑𝛽𝑘+1 (𝑛+1) 𝑛 𝑘=0 𝑓(𝑡𝑘,𝑇𝑘)). En la práctica, el término 𝑓(𝑡𝑛+1,𝑇𝑛+1) contiene ∇2𝑇𝑛+1 y la perfusión en 𝑇𝑛+1, por lo que la corrección exige resolver de nuevo un sistema lineal muy parecido al del L1 implícito. Por eso una estrategia eficiente es: 1) Construir 𝑏𝑛+1 para L1 implícito (historia + fuentes). 2) Predictor: evaluar 𝑓(𝑡𝑘,𝑇𝑘) y construir 𝑇𝑃𝑛+1 por la fórmula explícita AB (rápida). LATIN-AMERICAN JOURNAL OF COMPUTING (LAJC), Vol XIII, Issue 1, January - June 2026 3) Corrector: usar 𝑇𝑃𝑛+1 para evaluar 𝑓(𝑡𝑛+1,𝑇𝑃𝑛+1) y resolver el sistema lineal para 𝑇𝑛+1 con RHS que incluya la porción correcta de ABM (esto es análogo a resolver 𝒜𝑇𝑛+1=𝑏corr 𝑛+1 donde 𝑏corr 𝑛+1 incorpora la evaluación en 𝑇𝑃𝑛+1). Si se desea, iterar predictor–corrector 1–2 veces (iteración de punto fijo) hasta tolerancia. Fig. 9. La gráfica temporal muestra un comportamiento térmico oscilatorio en el centro del tejido para α = 0.8. La temperatura sube rápidamente hasta unos 42 °C, luego desciende, vuelve a aumentar y finalmente se estabiliza. Este patrón no monótono es típico de los modelos fraccionarios, donde la memoria térmica del tejido genera retrasos y variaciones en la difusión del calor. En el contexto del tratamiento térmico del cáncer ductal de mama, esto indica que pueden presentarse fases alternadas de sobrecalentamiento y enfriamiento parcial, lo que debe considerarse al diseñar protocolos de hipertermia para evitar picos de temperatura que puedan dañar el tejido sano. Fig. 10. La gráfica muestra cómo evoluciona la temperatura en el centro del tejido para distintos órdenes fraccionarios α. Se observa que, cuando α es menor, el tejido se calienta más lentamente y alcanza temperaturas menores debido a un efecto de memoria térmica más fuerte. Por el contrario, valores de α cercanos a 1 reproducen el caso clásico de Pennes, con una difusión del calor más rápida y eficiente. Esto es importante para el tratamiento térmico del cáncer ductal de mama, ya que indica que la memoria térmica influye en la capacidad del tejido para acumular calor: con α pequeño, el calentamiento es más lento y podría requerir mayor tiempo o intensidad para lograr temperaturas terapéuticas seguras y efectivas en la zona tumoral. Fig. 11. La gráfica muestra el perfil radial de temperatura al final del calentamiento para distintos valores de α. Cuando α es menor, el gradiente térmico se reduce: la temperatura máxima en el centro es menor y el calor se distribuye menos hacia la periferia. En cambio, con α cercano a 1, la difusión es más eficiente y se alcanzan temperaturas más altas en todo el dominio, como en el modelo clásico sin memoria. Esto es importante en el tratamiento térmico del cáncer ductal de mama, porque la forma en que se distribuye el calor determina la eficacia y seguridad del procedimiento. Con α pequeño, el calor se concentra más en la zona central, favoreciendo una ablación focalizada, mientras que valores mayores de α generan un calentamiento más amplio y uniforme. XV. VALIDACIÓN EXPERIMENTAL O CLÍNICA Los resultados obtenidos con las derivadas fraccionarias de Caputo en el modelo de Bio-calor de Pennes y utilizando L1 implícito, son comparados con trabajos previos como el de Caizer et al [24], para validar la parte experimental clínica y con Rahpeima & Lin [25] para validar la parte numérica y computacional. XVI. ANÁLISIS Y VALIDACIÓN CON LOS RESULTADOS EXPERIMENTALES Correspondencia entre la dinámica térmica del modelo fraccionario y la hipertermia real Las curvas térmicas obtenidas con el modelo fraccionario muestran aumentos, retardos y estabilización de la temperatura según el valor de α. Esto coincide con lo observado por Caizer et al. [24], quienes demostraron que el calentamiento real con nanopartículas no es lineal, sino progresivo y regulado por la relajación magnética. La presencia de oscilaciones y memoria térmica para α < 1 es coherente con el comportamiento superparamagnético, donde la temperatura se estabiliza alrededor de 42–43 °C, igual que en los experimentos de Caizer. XVII. EFECTO DEL PARÁMETRO Α (MEMORIA TÉRMICA) Tus resultados indican que para α pequeño el calentamiento es más lento y alcanza temperaturas menores. Este comportamiento corresponde con los hallazgos experimentales de Caizer et al, donde el incremento térmico depende de la respuesta magnética y de la estructura de los nanobioconjugados. Esto demuestra que el tejido tiene un comportamiento visco-térmico con retardo, exactamente lo que el modelo fraccionario representa a través de α. XVIII. RITMO DE CALENTAMIENTO FRENTE A LOS TIEMPOS EXPERIMENTALES Caizer mostraron que llegar a ~43 °C depende de la concentración de nanopartículas: concentraciones altas calientan rápido y las bajas tardan más. En tu modelo, α cercano a 1 produce calentamiento rápido, mientras que α pequeño genera un ascenso lento. Esta relación es equivalente: α alto representa sistemas que calientan eficientemente, y α bajo imita sistemas con baja potencia térmica. Por tanto, el modelo fraccionario reproduce correctamente la dependencia entre concentración y ritmo de calentamiento observada en laboratorio. XIX. PERFIL ESPACIAL FINAL DE TEMPERATURA Los resultados muestran que con α bajo el calor se concentra en el centro del tumor y que con α≈1 se difunde más lejos. Esto coincide con Caizer et al., quienes describen