Implementación en elementos finitos de modelos de remodelación ósea locales y no locales con activación discreta de BMUs

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo Fin de Grado Grado en Ingeniería de Tecnologías Industriales Implementación en elementos finitos de modelos de remodelación ósea locales y no locales con activación discreta de BMUs Autor: Pablo Manchado Morales Tutor: Francisco Javier Martínez Reina Tutor: José Luis Calvo Gallego Dpto. Ingeniería Mecánica y Fabricación Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2021 Trabajo Fin de Grado Grado en Ingeniería de Tecnologías Industriales Implementación en elementos finitos de modelos de remodelación ósea locales y no locales con activación discreta de BMUs Autor: Pablo Manchado Morales Tutores: Francisco Javier Martínez Reina Profesor Titular José Luis Calvo Gallego Profesor Ayudante Doctor Dpto. Ingeniería Mecánica y Fabricación Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2021 Trabajo Fin de Grado: Implementación en elementos finitos de modelos de remodelación ósea locales y no locales con activación discreta de BMUs Autor: Pablo Manchado Morales Tutor: Francisco Javier Martínez Reina Tutor: José Luis Calvo Gallego El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: Acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: A mis padres y hermano, por el apoyo constante. A Jaime y Ana, por los años compartidos. Agradecimientos Este trabajo simboliza el fin de una etapa que he compartido con personas a la que me gustaría agradecer su presencia en ella. En primer lugar, a mi familia por brindarme la oportunidad de poder cumplir mis metas y motivarme para lograrlas. A mis padres y a mi hermano, porque a pesar de estar lejos siempre han sido apoyo para continuar. A mis amigos, por el ánimo y apoyo incondicional. En especial a Miguel y Antonio por las visitas a Sevilla y a Fran, porque aunque seas un desastre me despertaste a tiempo para asistir al último examen de la carrera. A Jaime y Ana, con los que más tiempo he estado durante estos años y los que más me han tenido que soportar. Espero que sigamos compartiendo muchos años más. Por último, a todos los profesores que han contribuido a mi formación a lo largo de mi vida. Sobre todo, a Javier y José Luis por las horas dedicadas y el esfuerzo realizado para que este trabajo salga adelante. Pablo Manchado Morales Sevilla, 2021 III XÍndice 3.2 Limitaciones del modelo de poblaciones celulares 26 3.3 Implementación en Abaqus 27 3.3.1 Estructura y comentarios sobre el código principal 27 3.3.2 Modelo simulado 29 4 Modelo local 31 4.1 Modificaciones 31 4.1.1 Umbrales 31 4.1.2 Concentración de TGF-β32 4.1.3 Suavizado 33 4.2 Análisis de las simulaciones del Modelo local 37 4.2.1 Simulaciones iniciales 37 4.2.2 Variación de la velocidad de degradación TGF-β41 4.2.3 Variación velocidad Apoptosis OCa48 4.2.4 Variación del factor regulador Πdam act 53 4.2.5 Variación conjunta de los parámetros estudiados 54 4.2.6 Suavizado 57 5 Modelo no local 63 5.1 Implementación de las restricciones 64 5.2 Resultados y análisis de las simulaciones 65 6 Conclusiones y líneas futuras 71 6.1 Conclusiones 71 6.2 Líneas futuras 72 Apéndice A Parámetros empleados 73 Bibliografía 75 Notación BMU Basic Multicellular Unit MEF Método de los elementos finitos ROI Remodelación Ósea Interna ROE Remodelación Ósea Externa PMO Osteoporosis Postmenopaúsica PI Punto de integración OCpOsteoclastos precursores OCaOsteoclastos Activos OBaOsteoblastos Activos OBuOsteoblastos sin comprometer OBpOsteoblastos precursores OtOsteocitos TGF −βFactor de crecimiento transformante beta OPG Proteína osteoprotegerina RANK Receptor activador del factor nuclear κβ RANKL Ligando del receptor activador para factor nuclear κβ PTH Hormona paratiroidea ΠA act Función activadora de A ΠA rep Función represora de A fbm Fracción volumétrica de matriz ósea SED Stain Energy Density MES Minimally Effective Strain RBV Resorbed bone volume FBV Formed bone volume RVE Representative volume element PIoPunto de integración origen PIvPunto de integración vecino XI 1 Introducción La función principal que desempeñan los huesos es estructural, deben ser capaces de soportar las distintas cargas a las que se someten diariamente, para ello presentan un mecanismo que les permite "identificar" los estímulos provocados por esas cargas y modificar su estructura para obtener la distribución óptima. Este proceso es conocido como remodelación ósea y es llevado a cabo por una asociación de células llamada basic multicellular unit (BMU). Una BMU está constituida principalmente por osteoclastos y osteoblastos que actúan de manera coordinada reabsorbiendo y generando masa ósea, respectivamente. Esta coordinación se consigue gracias a la participación de los osteocitos, otro tipo de células que se encargan de "sentir" los estímulos mecánicos y activar la actuación de los anteriores. El proceso tiene como resultado una posible variación de la porosidad, de la anisotropía del hueso y una reducción de la densidad de microgrietas y del contenido mineral. Esto se debe a que durante este proceso el tejido antiguo, más dañado y con alto contenido mineral, será sustituido por un tejido intacto, llamado osteoide, que se encuentra formado únicamente por agua y colágeno [1]. 1.1 Motivación La motivación de este trabajo reside en el interés por el método de los elementos finitos (MEF) y por la bioingeniería. El MEF permite obtener soluciones aproximadas de ecuaciones diferenciales parciales asociadas a problemas físicos y presenta gran aplicabilidad en distintos campos, entre ellos se encuentra la bioingeniería. Esta disciplina aplica conceptos y métodos físico-matemáticos para abordar el estudio de los organismos vivos. A estos intereses personales se les suma la asociación directa entre gran cantidad de trastornos óseos y desequilibrios durante el proceso de remodelación. Algunas de estas enfermedades son la enfermedad de Paget o la osteoporosis, sobre la que se enfoca el trabajo. En periodo de vida de una BMU se dan distintas etapas: activación, reabsorción y formación, seguido de un periodo de reposo conocido como periodo de quiescencia [ 2 ]. La descoordinación de estas etapas ocasiona graves consecuencias, por ello el estudio del comportamiento de las distintas poblaciones celulares que intervienen en la remodelación cobra una gran importancia para el entendimiento de estas enfermedades y su tratamiento. 1.1.1 Osteoporosis La osteoporosis provoca una disminución de la densidad de masa ósea. Esto genera que los huesos se vuelvan más porosos, aumentando el número y tamaño de las cavidades internas, y a su vez la fragilidad y riesgo de fractura. Esta enfermedad es silenciosa porque no manifiesta síntomas hasta una pérdida significativa de densidad ósea como para que aparezcan fracturas. A partir de ese instante las consecuencias de la 1 2Capítulo 1. Introducción osteopososis serán evidentes para el paciente, el equipo médico y los responsables de financiar el sistema de salud. En 2002, se calculó que el coste combinado de las fracturas osteoporóticas en los Estados Unidos ascendía a 20.000 millones de dólares por año. Las fracturas osteoporóticas, también conocidas como fracturas traumáticas mínimas o leves, más frecuentes son vertebrales, de las muñecas y de cadera. Esta última tiene una importancia significativa ya que se considera un acontecimiento grave porque requiere intervención quirúrgica e ingreso hospitalario. Principalmente afecta a mujeres después de la menopausia pero también puede aparecer antes o afectar a hombres, adolescentes e incluso niños. En concreto, en España aproximadamente 2 millones de mujeres padecen osteoporosis con una prevalencia en la población postmenopáusica del 25 % (1 de cada 4). Se estima que esta enfermedad es la causante de unas 25.000 fracturas cada año y aproximadamente, 1 de cada 3 mujeres y 1 de cada 5 hombres mayores de 50 años sufrirá una fractura osteoporótica en su vida [ 3 ]. En la Figura 1.1 se muestran unas gráficas sobre las fracturas producidas por fragilidad: Figura 1.1 Epidemiología de las fracturas vertebral, de cadera y de muñeca con la edad. Adaptado de [4]. En la actualidad, a pesar de ser una enfermedad silenciosa se dispone de numerosos recursos para una detección precoz que permita una buena adaptación del tratamiento e incluso existen hábitos que pueden mejorar la calidad ósea. Entre estos se encuentran: la alimentación (dieta rica en calcio y vitamita D), abandono del consumo de tabaco y alcohol, o la práctica de ejercicio físico evitando el riesgo de caídas. Uno de los recursos disponibles es el tratamiento con fármacos, estos buscan tanto el aumento de la masa como la disminución de las fracturas. Los fármacos utilizados pueden dividirse en tres tipos según [5]: - Antirresortivos: actúan sobre los osteoclastos inhibiendo o retrasando el proceso, lo que provoca una estabilización e incluso un aumento discreto de la masa. En este grupo se encuentran los estrógenos, bifosfonatos, calcitonina y denosumab. - Osteoformadores: actúan sobre los osteoblastos produciendo un aumento mayor de la densidad ósea. Los más característicos son los fluoruros, paratohormona (PTH), factores de crecimiento y fosfanatos. 1.2 Objetivos 3 - El grupo restante se asocia a fármacos que mezclan los dos mecanismos anteriores, destacan: anabolizantes, vitamina D o calcitriol y progestágenos. Como se ha reflejado, la mayoría de los fármacos usados frente a estas enfermedades se basan en la interacción directa con las células que desarrollan la remodelación, o de manera indirecta, interactuando con alguno de sus factores de crecimiento o regulación. Por tanto, el estudio de esas poblaciones celulares juega un papel fundamental para la adaptación y mejora de los tratamientos. Llegados a este punto, los modelos matemáticos generados facilitan el estudio y permiten probar hipótesis experimentales y teóricas, entre ellos destacan los propuestos por Pivonka et al. [ 2 ], [ 6 ], Scheiner et al. [ 7 ] y el desarrollado por Martínez-Reina et al. [ 8 ], siendo este último en el que se basa el trabajo. 1.2 Objetivos El trabajo desarrollado presenta como fin la obtención de simulaciones realistas del proceso de remodelación ósea interna. Para ello se ha estudiado el proceso biológico, se ha revisado la literatura actual y se ha profundizado en los conocimientos y manejo del programa Abaqus, en concreto, en el uso de subrutinas. Se parte de un código que implementa las ecuaciones desarrolladas en el modelo matemático mencionado anteriormente. A continuación se implementa dicho modelo en un código de Elementos Finitos, en el que se discretiza el dominio de estudio en elementos finitos. Estos presentan puntos de integración (PI), que son justamente en los que se plantean las ecuaciones del modelo matemático. De esta forma en cada PI puede activarse una BMU y los resultados más relevantes que se obtendrán serán las cantidades de tejido óseo reabsorbido y formado. Tanto código de partida como el modelo matemático desarrollan un estudio a nivel local atendiendo únicamente a cada PI de manera independiente. En ambos se tiene en cuenta la reabsorción y formación de tejido óseo, pero no la secuenciación entre ambas fases, ni los sucesos que tienen lugar en el entorno cercano a cada PI. En una primera aproximación se intenta obtener simulaciones que presenten las distintas etapas de vida de una BMU indicadas anteriormente. Para ello se implementarán las condiciones necesarias para que tengan lugar las fases de reabsorción y formación. Este análisis se realiza a nivel local, sin tener en cuenta las distintas BMUs presentes en el dominio. También se realiza el estudio de los efectos producidos por la modificación de distintos parámetros o la relación entre ellos. A continuación se realiza un estudio no local. En este se implementan las distintas condiciones que deben darse en el entorno cercano a una BMU para que se produzca su activación, las cuales son: presentar una concentración de las células que desarrollan la reabsorción superior a un umbral, que esta concentración sea la máxima de su entorno, y que no haya ninguna otra BMU activa en el entorno. Estas condiciones se trasladarán a elementos finitos estudiando el dominio cercano a cada punto de integración e implementando las restricciones mencionadas. Las modificaciones generadas en el modelo local se trasladarán a este, obteniéndose mejores aproximaciones del proceso tanto a nivel local como global. El código de partida implementa el modelo matemático en Abaqus. Este programa aplica el MEF permitiendo la simulación y análisis del proceso de remodelación. 1.3 Estructura del documento En primer lugar, en el Capítulo 2 se desarrollan los fundamentos biológicos necesarios para la compresión del trabajo realizado, explicando nociones básicas sobre el hueso y profundizando en el proceso de remodelación ósea y la enfermedad de la osteoporosis. 4Capítulo 1. Introducción En el Capítulo 3 se explicará el modelo matemático desarrollado por Martínez-Reina et al. [ 8 ] en el que se basa el código y algunos aspectos sobre el código partida y modelo simulado. Los siguientes capítulos exponen el modelo local, Capítulo 4, y el modelo no local, Capítulo 5, que se han generado durante el desarrollo del trabajo. En ellos se indican las distintas modificaciones llevadas a cabo con la finalidad pretendida y se comentan los resultados obtenidos. Finalmente, en el Capítulo 6 se expondrán las conclusiones obtenidas de los dos modelos, y se indicarán las futuras líneas de estudio. 2 Biología Ósea El esqueleto se encuentra formado por el conjunto de los huesos. Estos desempeñan diversas funciones esenciales como la estructural : dan soporte y conforman el cuerpo. Simultáneamente permiten, junto con los músculos, la transmisión de cargas para realizar movimientos (función locomotora) y protegen órganos internos como el corazón. Desarrollan otras funciones menos intuitivas entre las que destaca su actuación como reserva de minerales como el calcio o el fósforo y el almacenamiento de médula ósea , crucial para el desarrollo de algunas funciones fisiológicas. Los huesos se encuentran constituidos por distintos tejidos y células, también presentan un mecanismo que les permite "renovarse". En este capítulo se profundizará en estos aspectos fundamentales para comprender el trabajo realizado. 2.1 Clasificación huesos Atendiendo a tamaño y forma los huesos se pueden clasificar según [10] en: 2.1.1 Huesos largos En este tipo de huesos una dimensión predomina sobre las otras. Son duros y densos, brindan resistencia estructura y movilidad, algunos ejemplos son el fémur o el radio. Presentan una diáfisis, parte intermedia del hueso, y dos extremos llamados epífisis, esta parte es la encargada de interactuar con otros huesos creando articulaciones [ 11 ]. La parte que une la diáfisis con la epífisis se conoce como metáfisis, estas partes se indican en la Figura 2.1: Figura 2.1 Partes de un fémur [12]. 5 6Capítulo 2. Biología Ósea 2.1.2 Huesos cortos No presentan una dimensión dominante sino que son todas similares. Suelen encontrarse agrupados dando estabilidad y soporte y permiten en movimiento de las articulaciones. Algunos ejemplos son: el carpo y el tarso (huesos que componen mano y pie, respectivamente). 2.1.3 Huesos planos Son aquellos de pequeño espesor, es decir, presentan dos dimensiones considerables frente a la tercera. En algunos casos actúan protegiendo a tejidos blandos como el cráneo, y en otros permiten anclar músculos en su amplia superficie como la escápula. 2.1.4 Huesos irregulares Son aquellos difíciles de caracterizar por su forma y no entran en ninguna de las otras categorías. Este es el caso de los huesos que conforman el oído. 2.1.5 Huesos sesamoideos Son huesos pequeños y redondos. Estos se encuentran en zonas de alta fricción, tensión o esfuerzo físico, ya que, su misión es incrementar la función palanca de los músculos. En concreto están presentes en las articulaciones. En la Figura 2.2 se muestran algunos ejemplos de los tipos indicados: Figura 2.2 Ejemplos de cada uno de los tipos de huesos mencionados [10]. 2.2 Tejido Óseo 7 2.2 Tejido Óseo Las propiedades mecánicas del tejido óseo dependen tanto de los elementos que lo componen como la disposición de los mismos, por ello resulta interesante profundizar en su composición. Al igual que el resto de tejidos conjuntivos, está formado por una matriz y un conjunto de células. 2.2.1 Matriz Ósea Está compuesta por varias fases: una orgánica (principalmente formada por colágeno, 20 % de la masa ósea) responsable de su rigidez, viscoelasticidad y tenacidad; una fase mineral (70 % de la masa ósea) formada por hidroxiapatita ( Ca10(PO4)6(OH)2 ) para el refuerzo estructural, rigidez y homeostasis mineral [13]. El contenido restante es agua (10 % de la masa ósea). Las células encargadas de la formación depositan osteoide, que es un material formado por la fase orgánica y agua. Posteriormente, durante el proceso de mineralización, parte de este agua es sustituida por mineral. 2.2.2 Células óseas y su función Existen cuatro tipos de células óseas maduras con sus funciones definidas. Los osteoclastos, procedentes de células madre hematopoyéticas y los osteoblastos, procentes del linaje de las células madre mesenquimales. Los dos tipos de células restantes son los osteocitos y las células de borde, ambos proceden de los osteoblastos de los que se diferencian para realizar otra función concreta. Osteoclastos Tienen gran tamaño, ya que, se tratan de células multinucleadas (presentan entre 3 y 20 núcleos, véase Figura 2.3) y desempeñan la función de reabsorción ósea. Como se ha mencionado, pertenecen al linaje de las células madre hematopoyéticas y surgen de la unión de varios osteoclastos precursores. Posterior a esta fusión y, estando varios osteoclastos maduros activos simultáneamente, comienzan a segregar ácidos que disuelven las fases orgánica e inorgánica de la matriz ósea. La velocidad de avance es de unas decenas de micras al día. Una vez cumplida su función los osteoclastos sufren apoptosis (muerte celular programada). Como consecuencia de dicha actuación aparece un hueco en la matriz ósea, denominado cavidad de reabsorción. Figura 2.3 Osteoclastos [14]. Osteoblastos Son las células encargadas de generar la matriz ósea, aunque segregan únicamente la fase orgánica. A diferencia de las anteriores son células mononucleadas, pequeñas, de forma cuboidal (Figura 2.4) y 14 Capítulo 2. Biología Ósea Figura 2.11 Comparativa de un hueso normal y otro con osteoporosis. Adaptada de [31]. Atendiendo a la causa productora de la osteoporosis, se pueden considerar dos grandes grupos: Osteoporosis primaria o involutiva Es el grupo más frecuente y su causa se asocia al transcurso de los años. Dentro de este grupo podemos distinguir [32]: - Tipo I u Osteoporosis postmenopaúsica (PMO): se da en las mujeres entre 50 y 75 años, como consecuencia del cese de la actividad ovárica. Afecta especialmente al hueso trabecular, siendo características las fracturas vertebrales por aplastamiento, las de la extremidad distal del antebrazo. - Tipo II o senil: se produce en ambos sexos y a edades más avanzadas (mayores de 70 años). Esta pérdida no es tan acelerada como la anterior y afecta a ambos tipos de hueso por igual. La fractura característica es la de cadera, además de las de pelvis, húmero proximal y vértebras. Es necesario recalcar que esta división en dos tipos no indica que sean dos procesos diferentes. Por el contrario, muestra momentos distintos de un mismo proceso que evoluciona con el paso de los años. Osteoporosis secundarias Se usa esta denominación cuando existe una causa capaz de producir el trastorno, con independencia de la menopausia y la edad. Hay diversos causantes entre los que destacan: algunas enfermedades endocrinas, gastrointestinales, hematológicas o conectivopatías, así como la inmovilización prolongada o el uso de distintos fármacos. 2.5.1 Osteoporosis Postmenopaúsica Dentro de la clasificación anterior esta es la más frecuente de todas. Se debe a una combinación del envejecimiento y de la aparición de la menopausia. El cese de la función ovárica provoca una reducción de estrógenos se acompaña de una pérdida ósea acelerada. El tratamiento con estrógenos permite revertir esta situación [33]. Existen otros tratamientos farmacológicos para este tipo de osteoporosis, entre los que se encuentran el uso de bifosfonatos o el denosumab. Los bifosfonatos inhiben la actividad osteoclástica 2.5 Osteoporosis 15 mediante la alteración de las interacciones entre proteínas, bloqueo de señales celulares y aumentando la apoptosis de los osteoclastos. El denosumab es un anticuerpo monoclonal (tienen como objetivo un único anticuerpo específico) que interactúa con el ligando RANKL impidiendo la unión RANK-RANKL, dificultado la diferenciación de los osteoclastos activos y por tanto la reabsorción. La administración de este último fármaco es estudiado en el modelo base del trabajo. El tratamiento habitual aprobado por la Organización Mundial de la Salud para combatir la osteoporosis postmenopáusica consiste en dosis de 60 mg inyectadas cada 6 meses. En mujeres que presentan osteoporosis postmenopaúsica se ha demostrado que provoca un aumento de la masa ósea y una disminución del riesgo de sufrir fracturas [8]. En este capítulo se han comentado los conocimientos biológicos básicos y necesarios para el entendimiento del trabajo. A continuación, se expone el modelo matemático que modela el proceso de remodelación ósea interna y algunos mecanismos que influyen en el mismo. 3 Modelo Matemático Dentro del conjunto de los modelos matemáticos que se han propuestos para simular la remodelación ósea interna se distinguen dos grupos: fenomenológicos y mecanobiológicos. Los primeros también son conocidos como modelos de "caja negra", esto se debe a que establecen una relación matemática entre variables de entrada y unas variables de salida, o una variación de las mismas. En concreto, las variables de entrada serían las cargas mecánicas soportadas por el hueso y las de salida las propiedades mecánicas. Los modelos fenomenológicos no describen los procesos biológicos que tienen lugar durante la remodelación, simplemente describen la relación causa-efecto mediante un conjunto de ecuaciones matemáticas. Por el contrario, los modelos mecanobiológicos intentan incluir los procesos biológicos que tienen lugar dentro del hueso, obteniendo más información de lo que ocurre y no solo los cambios producidos en las propiedades mecánicas. Durante el desarrollo de este capítulo se expondrá el modelo de Martínez-Reina et al. [ 8 ], se explicarán algunos aspectos del código que lo implementa y del modelo de Abaqus simulado. Se indican las funciones que implementan cada uno de los mecanismos expuestos, pero el código no se adjunta por motivos de propiedad intelectual. 3.1 Desarrollo del modelo Este modelo fue presentando inicialmente por Pivonka et al. [ 2 ] y ha ido sufriendo modificaciones. Algunas de ellas son: la inclusión del proceso de mineralización [ 1 ], la influencia de la superficie ósea [ 34 ] o la inclusión del daño [ 35 ]. El modelo que se presenta a continuación se ha obtenido de una de las publicaciones recientes [8] e incluye las distintas modificaciones indicadas. En primer lugar se explican los mecanismos bioquímicos que ocurren. Estos permiten la diferenciación de las células mencionadas en el Apartado 2.2.2. Los osteoclastos derivan de las células madre hematopoyéticas, en el modelo se denominan osteoclastos precursores (OCp), los cuales se diferenciarán a osteoclastos activos ( OCa ) gracias a la unión receptor-ligando (RANK-RANKL). Al finalizar el proceso de reabsorción sufrirán apoptosis. Por otro lado, los osteoblastos activos ( OBa ) provienen de las células madre mesenquimales capaces de diferenciarse en varios tipos de células. En el modelo esas células mesenquimales son denominadas osteoblastos sin comprometer ( OBu ) y, una vez se diferencian a osteoblastos, se conocen como osteoblastos precursores ( OCp ), los cuales podrán diferenciarse a activos. Esta diferenciación tiene lugar gracias a factores de crecimiento (como el TGF- β ) o de regulación (como el PTH). Al terminar la reabsorción sufrirán apoptosis o se diferenciarán a osteocitos (Ot). El proceso de remodelación se puede describir como un balance celular utilizando las concentraciones celulares indicadas anteriormente: osteoblastos precursores ( OBp ), osteoblastos activos ( OBa ), osteocitos ( Ot ), osteoclastos activos ( OCa ), osteoclastos sin comprometer ( OCu ), osteoclastos precursores (OCp). Las concentraciones de los dos últimos se consideran constantes. 17 18 Capítulo 3. Modelo Matemático dOBp dt =DOBu·ΠT GF−β act ·OBu+POBp·Πψbm act ·OBp−DOBp·ΠT GF-β rep ·OBp(3.1) dOBa dt =DOBp·ΠT GF−β rep ·OBp−AOBa·OBa(3.2) dOCa dt =DOCp·ΠRANKL act ·OCp−AOCa·ΠTGF-β act ·OCa(3.3) dOt dt =η·d fbm dt (3.4) siendo: -DOBu , DOBp y DOCp las tasas o velocidades de diferenciación de OBu en OBp , OBp en OBa y OCpen OCa, respectivamente. -POBpes la velocidad de proliferación de OBp. -AOBayAOCason las velocidades de apoptosis de OBayOCa, respectivamente. - Las variables ΠTGF−β act y ΠTGF−β rep representan las funciones activadoras y represoras relacionadas con la unión entre TGF-βy su receptor. - La variable ΠRANKL act representa la función activadora de la unión RANK-RANKL. - La función Πψbm act cuantifica la influencia del estímulo mecánico en la proliferación de OBp. -fbm es la fracción volumétrica de matriz ósea, que se desarrollará a continuación, y η es la concentración de osteocitos en el volumen de la matriz ósea. Figura 3.1 Representación esquemática del modelo con las distintas fases de diferenciación de las células, su interacción con los factores de crecimiento y regulación y la inyección externa del fármaco denosumab [8]. 3.1 Desarrollo del modelo 19 En modelo se suponen proporcionales el volumen reabsorbido y la concentración de osteoclastos activos así como el volumen formado y la concentración de osteoblastos activos, quedando la variación de la fracción volumétrica de matriz ósea: d fbm dt =−kres ·OCa+kf orm ·OBa(3.5) siendo kres y kf orm la velocidad de reabsorción de matriz ósea y la velocidad de formación de osteoide, respectivamente. En la Figura 3.1 se muestra un esquema del proceso y en el Apéndice A se incluye el valor de los parámetros del modelo. A continuación, se desarrollan algunos de los términos presentes en las ecuaciones (3.1-3.4). 3.1.1 Mecanismo de regulación bioquímica Esta regulación se tiene en cuenta mediante las funciones activadoras y represoras comentadas anteriormente. En primer lugar, la concentración de TGF- β almacenada en la matriz y liberada durante la reabsorción por los osteoclastos se calcula [2]: TGF-β=αT GF-β·kres ·OCa ˜ DTGF-β (3.6) siendo αTGF-β la concentración de TGF- β en la matriz ósea y ˜ DTGF-β la velocidad de degradación del TGF- β . Esta concentración es usada para definir las funciones activadora y represora de las ecuaciones (3.1-3.3): ΠTGF-β act =T GF-β KTGF−β act +T GF-β(3.7) ΠTGF-β rep =TGF-β 1−TGF-β/KT GF-β rep (3.8) donde KTGF-β act y KTGF-β rep las constantes de activación y represión, respectivamente. El cálculo de estas variables se realiza dentro de las rutinas de Matlab yAbaqus con la función: FU_TGFBE. Para la obtención de la función activadora del RANKL se necesitan algunas concentraciones, estas se obtienen de las ecuaciones propuestas en [34]: RANK =RRANK ·OCp(3.9) PTH =PPT H +βPT H ˜ DPTH (3.10) OPG =βOPG ·OBa·ΠPTH rep ˜ DOPG +βOPG·OBa·ΠPT H rep OPGmax (3.11) La capacidad portante de RANKL en la superficide de los OBpyOt(RANKLe f f ) se calcula: RANKLe f f =RRANKL ·OBp·ΠPT H act +RRANKL ·Ot·Πdam act (3.12) Como se ha indicado, el denosumab actúa compitiedo con el RANK y OPG por unirse al ligando RANKL. Una de las últimas modificaciones del modelo es la inclusión del término asociado al fármaco en la ecuación (3.13) realizado por Martínez-Reina et al. [9]. 20 Capítulo 3. Modelo Matemático RANKL =RANKLe f f ·βRANKL +PRANKL βRANKL +˜ DRANKL ·RANKLe f f · "1+OPG Kd,[RANKL−OPG] +RANK Kd,[RANKL−RANK] +ζ·Cp,den Kd,[RANKL−den]#−1 (3.13) A continuación, se indican los significados de los términos presentes en las ecuaciones: -RRANK es la capacidad de producción de RANK por parte de los OCp. -βOPG , βRANKL y βPTH son las velocidades de producción endógenas de OPG , RANKL y PTH , respectivamente. -ΠPTH rep y ΠPTH act las funciones activadora y represora de la hormona paratiroidea PTH . Estas se calculan de manera similar a las ecuaciones (3.7) y (3.8): ΠPTH act =PT H KPTH act +PT H (3.14) ΠPTH rep =PTH 1−PTH/KPT H rep (3.15) siendo KPTH act yKPT H rep las constantes de activación y represión, respectivamente. -˜ DOPG , ˜ DRANKL y ˜ DPTH son las velocidades de degradación del OPG , RANKL y PTH , respectivamente. -RANKe f f es la capacidad portante de RANKL , tal como se ha definido anteriormente en la ecuación (3.12). -Πdam act factor que regula la producción de RANKL en función del daño. Este se calcula: Πdam act =dξ dξ+δξ 50 (3.16) donde d es el daño, ξ el factor de forma y δ50 el valor del daño para el que se produce un 50 % de la respuesta máxima. El valor de estos parámetros se indica en el Apéndice A adjunto al final del trabajo y en el código se modela con la función: FU_DAM_RANKL. -PRANKL producción extra de RANKL debida al desuso o a la enfermedad, cada uno de los cuales contribuye con un término a su cálculo: PRANKL =Pmech RANKL +PPMO RANKL (3.17) El primer término se explicará en el Apartado 3.1.2 y está asociado al estímulo mecánico, mientras que el segundo término es consecuencia de la menopausia. Como se ha comentado antes esta favorece la producción de RANKL, que permite una mayor diferenciación de OCa . Además se debe tener en cuenta que este efecto se incrementa con el tiempo. Para modelar este efecto Martínez-Reina et al. [35] usan una función sigmoidal: PPMO RANKL(t) = PPMO,max RANKL ·(t−tonset)2 (t−tonset )2+δ2 PMO para t>tonset (3.18) 3.1 Desarrollo del modelo 21 siendo PPMO,max RANKL la velocidad de producción de RANKL máximo asociado a la osteoporosis postmenopaúsica, tonset indica el comienzo de la enfermedad y δPMO es una constante que indica el instante temporal en el que se alcanza el 50 % de la producción máxima ( PPMO,max RANKL ). Se incluye en el código con la función: FU _ LOADCASE esta función realiza distintas acciones dependiendo del valor del parámetro LoadCase. -Kd,[RANKL−OPG] , Kd,[RANKL−RANK] y Kd,[RANKL−den] son las constantes de disociación entre el RANKL yOPG,RANK y el fármaco denosumab. -Cp,den es la concentración de denosumab en el compartimento central (suero sanguíneo) y ζ es el factor de accesibilidad, que en el modelo farmacocinético controla la cantidad de fármaco del compartimento central que está disponible en el compartimento hueso. Finalmente, a partir de las ecuaciones anteriores, se puede obtener la función activadora de la ecuación (3.3): ΠRANKL act =RANKL ·RANK Kd,[RANKL−RANK]+RANKL ·RANK (3.19) Por último, hay que indicar que para el cálculo de la concentración y funciones asociadas al PTH se usa: FU_PTH y para RANK, OPG y RANKL se usa: FU_RANKL. 3.1.2 Mecanismo de regulación biomecánica El modelo incluye la respuesta biológica al estímulo mecánico propuesta por Frost [ 22 ]. Esta teoría postula la existencia de cuatro zonas o windows en función del valor de la tensión minima efectiva (Minimally Effective Strain, MES). En la Figura 3.2 se representan las cuatro zonas mencionadas. Figura 3.2 Teoría del Mecanostato [22]. - 1) Disuse window (DW) Estado de desuso, se produce una reabsorción neta. - 2) Adapted window (AW) Zona de adaptación, no se produce modificación. - 3) Mild overload window (MOW) Zona de sobrecarga leve, se produce una formación neta. - 4) Pathologic overload window (POW) Zona de sobrecarga alta, fractura. 22 Capítulo 3. Modelo Matemático En el modelo para un estado estímulo mecánico alto se producirá un incremento de OBa , esto se ve reflejado en la ecuación (3.1) por el término Πψbm act , que provocará un aumento directo de la proliferación de OBp e indirecto en la concentración de OBa . En caso de un estímulo bajo, se producirá un aumento en la producción de RANKL (a partir del término PRANKL ) que producirá una mayor cantidad de OCa. Ambos parámetros mencionados dependen de la densidad de energía de deformación (SED, ψbm ) presente en la zona del hueso. Esta, para estados uniaxiales de deformación se calcula: ψbm =1 2·E·MES2(3.20) siendo E el módulo de Young, dependiente de la densidad y del daño acumulado. En primer lugar se indica el cálculo de la densidad. Esta cambia como consecuencia de la variación de la porosidad y por la mineralización: ρt=m Vbm =ρm·vm+ρo·vo+ρw·vm(3.21) donde m y Vbm son la masa y el volumen ocupado por la matriz ósea, y vm , vo , vw , ρm , ρo y ρw son los volúmenes específicos y densidades de las fases presentes: mineral, orgánica y agua. El volumen de la fase orgánica se asume constante ( vo =3/7) [ 36 ]. Por el contrario, los volúmenes de las fases mineral y agua ( vm y vw ) no son constantes, el cálculo de la variación se detallará posteriormente (Apartado 3.1.4). Finalmente, la densidad del tejido óseo se calcula: ρ=m Vbm ·Vbm VRV E =ρt·fbm (3.22) A continuación se obtiene el daño acumulado en la matriz ósea. La aparición del daño se debe a la fatiga y este daño se reduce con la reabsorción del tejido y la generación del nuevo, por tanto: ˙ d=˙ dA−˙ dR(3.23) siendo ˙ dA la velocidad de acumulación de daño por fatiga y ˙ dR la velocidad de reparación de daño por remodelación. Se asume que el daño está uniformemente distribuido en un volumen representativo del elemento, así la cantidad de daño reparado es: ˙ dR=d·kres ·OCa fbm (3.24) El daño producido es evaluado a partir del trabajo desarrollado por García et al. [ 37 ] y Pattin et al. [ 38 ]. Los ensayos experimentales de fatiga proporcionan la evolución del daño con la tensión o deformación y el número de ciclos: Nf=Ki εδii=c(compresión),t(tensión) (3.25) donde Nf es la vida a fatiga en número de ciclos, Ki una constante y ε la deformación. El valor de estos parámetros ha sido determinado en varios estudios [ 37 ] y se incluyen en el Apéndice A localizado a final del trabajo. Con las ecuaciones anteriores se puede obtener el valor del daño: dc=−1 γcln(1−CcεδcN)(3.26) dt=1−1 Ct2 ln(eCt2−Ct1εδtN)1 γt(3.27) 3.1 Desarrollo del modelo 23 siendo N el número de ciclos y: γc=−5.238 ·(ε−6100)·10−3+7; Cc=1−e−γc Kc para compresión γt=−0.018 ·(ε−4100)+12; Ct1=eCt2−1 Kt para tracción Finalmente, a partir de la correlación de Beaupré se puede obtener el valor del módulo de Young: E(ρ,d) = 2014 ·ρ2.5·(1−d)para ρ<1.2g/cm3 1763 ·ρ3.2·(1−d)para ρ≥1.2g/cm3(3.28) donde se ha multiplicado en la correlación de Beaupré por (1-d) para tener en cuenta la degradación de las propiedades por el daño. Con los desarrollos indicados queda totalmente definida la SED y se continúa indicando los términos que implementan la Teoría del Mecanostato de Frost [ 22 ] en el modelo. Para estados de desuso: Pmech RANKL =   Pmech,max RANKL ·1−ψbm ψrpara ψbm <ψr 0para ψbm ≥ψr (3.29) siendo Pmech,max RANKL el valor máximo de producción de RANKL por desuso y ψr es el valor de la energía de deformación límite que provoca una producción de RANKL por desuso. Para los estados de sobrecarga: Πψbm act =˜ Πmech act,OBp·1+λ·ψbm ˜ ψbm −1λ=0para ψbm ≤˜ ψbm (Carga normal) λ=1.25 para ψbm >˜ ψbm (Sobrecarga) (3.30) donde ˜ ψbm es el valor límite de SED entre una carga normal y una sobrecarga. En el código empleado el daño se calcula a partir de la función: FU _ DANO, las propiedades mecánicas con MECH _ PROP y las variables que modelan la Teoría del Mecanostato con FU _ MECH. 3.1.3 Degradación de las propiedades a fatiga por el contenido mineral El aumento de la fase mineral provoca un aumento de la rigidez, pero también de la fragilidad, afectando a las propiedades a fatiga del hueso. La variación de estas propiedades de fatiga con el contenido de mineral se modela con las propuestas realizadas por Matínez-Reina et al. [ 39 ] bajo las cuales las propiedades a fatiga se degradan con el aumento del contenido mineral: - La forma de las curvas d-N, modeladas con las ecuaciones (3.26) y (3.27) se mantienen a pesar de los cambios en el contenido mineral. - Solo la vida a fatiga se ve afectada por la variación del contenido mineral, variando Kt y manteniendo constante el exponente δt en la ecuación (3.25). Esto modifica Ct1 (ecuación (3.1.2)) para aumentos de Ktse alarga la vida a fatiga y se reduce el daño acumulado. - El límite a fatiga se establece en 107 ciclos. Normalmente, este límite se alcanza para un valor último de la fracción de la deformación a tracción εu/β , donde β depende del tipo de 30 Capítulo 3. Modelo Matemático En este capítulo se ha expuesto el modelo matemático y se han desarrollado los mecanismos que forman parte del mismo. También se ha comentado alguna de las limitaciones que presenta el modelo que pueden servir para estudios futuros. Además, se indican algunos aspectos de su implementación en el programa de Abaqus, así como algunos comentarios sobre el modelo simulado. A continuación, se exponen los dos modelos generados y los resultados obtenidos. 4 Modelo local En primer lugar, destacar que el modelo de partida asume el tejido reabsorbido y formado directamente proporcional con la concentración de OCa y OBa , respectivamente. Pero para cualquier concentración de las células anteriores no se producen estos procesos, es decir, hay periodos temporales en la que la concentración de las poblaciones celulares no es suficiente para reabsorber o formar tejido óseo. Durante esos periodos los procesos de reabsorción y formación están desactivados, dando lugar a unas fases de reposo o quiescencia. El objetivo principal de este capítulo es la obtención de las fases que tienen lugar en el periodo de vida de una BMU. Para ello se realiza un estudio independiente de cada BMU sin tener en cuenta las demás BMUs presentes en el dominio cercano. Y se implementan algunas modificaciones que asemejan en mayor medida las simulaciones con el desarrollo real del proceso. Finalmente, se realiza el estudio del efecto de la variación de algunos parámetros que intervienen para un correcto ajuste. 4.1 Modificaciones 4.1.1 Umbrales Como se ha comentado en el Capítulo 2, para el comienzo de la reabsorción o formación es necesario un reclutamiento de OCa u OBa , respectivamente. Esto tiene especial importancia ya que, hasta que no se alcanza una cierta concentración de dichas poblaciones celulares no tienen lugar dichos procesos. Además, durante el periodo de vida de las BMUs se distinguen tres fases: reabsorción, formación y quiescencia o periodo de reposo, esta última fase en el código de partida solo tendría lugar si las concentraciones de OCayOBafuesen nulas, lo cual no es realista. El cambio que se plantea en el presente trabajo consiste en introducir unos umbrales para la concentración de las poblaciones celulares que llevan a cabo la reabsorción y formación, en concreto, uno para los OCa y otro para los OBa . Para valores inferiores a estos umbrales se considera que no se reabsorbe o forma tejido óseo mientras que, si se supera alguno de los umbrales, se activa el proceso correspondiente. Esto se implementa a partir de la definición de dos parámetros ( OBaTH y OCaTH ), que serán los umbrales, y dos variables ( FOBa y FOCa ), que actuarán como bandera, tomando valor unidad o nulo en función de si las concentraciones superan los umbrales establecidos o no. Al introducir este cambio se deben modificar las variables dependientes del volumen reabsorbido o formado. En este caso afecta a: la fracción volumétrica de masa ósea ( fbm ), la concentración del TGF-βy el daño reparado ( ˙ dR). En primer lugar se analiza la fracción volumétrica de masa ósea ( fbm ). El cálculo de esta se obtiene con la ecuación (3.5) indicada en el capítulo anterior. En el segundo miembro de esta ecuación aparecen dos sumandos. El primero corresponde al volumen óseo reabsorbido y el segundo al 31 32 Capítulo 4. Modelo local volumen óseo formado. Con la modificación realizada es necesario la introducción de las banderas que indican si los procesos de reabsorción o formación están activos. La bandera de reabsorción multiplicará al primer sumando y la de formación al segundo. La ecuación queda: d fbm dt =−kres ·OCa·FOCa+kf orm ·OBa·FOBa(4.1) La concentración de TGF- β se obtiene de la ecuación (3.6). En el segundo miembro de la ecuación aparece el volumen reabsorbido multiplicado por la concentración de TGF- β presente en la matriz ósea ( α ). La aparición del volumen reabsorbido se debe a la aceptación de que diversos tipos de factores de crecimiento, entre los que está el TGF- β , se almacenan en la matriz ósea y son liberados cuando el tejido es reabsorbido por los OCa [ 2 ]. Para una correcta implementación de la modificación es necesario multiplicar este término por la bandera de reabsorción, quedando la ecuación: TGF-β=αT GF-β·kres ·OCa·FOCa ˜ DTGF-β (4.2) Finalmente, se analiza la ecuación (3.24) que modela el daño reparado ( ˙ dR ). Debido a que el tejido que se reabsorbe es el más antiguo, con mayor contenido mineral y más dañado, el daño reparado será proporcional al volumen reabsorbido. De nuevo, es necesario multiplicar el término asociado al volumen reabsorbido por la bandera de reabsorción. La ecuación queda entonces: ˙ dR=d·kres ·OCa·FOCa fbm (4.3) De esta forma, si el proceso está activo, los términos asociados al volumen reabsorbido o formado se multiplican por la unidad quedando los términos igual que en las ecuaciones iniciales. Por el contrario, si el proceso está inactivo, toman valor nulo y anulan el término asociado a la reabsorción o formación de tejido. Hay que destacar que este cambio no supone una modificación del modelo matemático, pero sí de las líneas de código que implementan dichas ecuaciones. 4.1.2 Concentración de TGF-β La ecuación diferencial que gobierna la concentración de TGF-βes: dT GF-β dt =α·kres ·OCa·FOCa−˜ DTGF-β·TGF-β+PTGF-β(4.4) El TGF- β se libera de la matriz ósea por el proceso de reabsorción y se degrada a una velocidad ˜ DTGF−β . El término PT GF-β representa un posible aporte externo de TGF- β . Así, la ecuación (4.2) corresponde a una situación en la que no hay aporte externo de TGF- β y además se consigue el equilibrio en la concentración de TGF- β de forma instantánea. En otras palabras, si en la ecuación diferencial (4.4) hacemos dTGF- β /dt = 0 y resolvemos la ecuación resultante en TGF-beta, se obtiene la ecuación (4.2). De esta forma la liberación de TGF- β y su degradación ocurren de forma continua y siempre en los valores de equilibrio. Sin embargo, si la reabsorción ocurre de forma episódica, hay que evaluar el balance de ambos procesos de forma más exhaustiva, volviendo a utilizar para ello la ecuación diferencial que gobierna la concentración de TGF- β . Esta se integra con un algoritmo de Euler explícito, de la siguiente forma: TGF-β(t) = T GF-β(t−1)+αTGF-β·kres ·OCa·FOCa−TGF-β(t−1)·˜ DTGF-β(4.5) 4.1 Modificaciones 33 En el segundo miembro de dicha ecuación, de izquierda a derecha, el primer sumando es la concentración de TGF- β almacenada en el anterior incremento, el segundo representa el TGF- β liberado de la matriz ósea por la reabsorción de los osteoclastos y el tercero representa la degradación del TGF-β. Durante las primeras simulaciones que se realizaron para comprobar el funcionamiento de esta modificación aparecieron problemas en la convergencia de los resultados. El error se encontraba en la concentración de TGF- β . Para ciertos días de la simulación presentaba valores negativos lo cual es imposible, ya que, es la concentración de una proteína. Para solventarlo se añadió una condición simple que impone valor nulo a la concentración de TGF-βsi esta toma valor inferior a cero. 4.1.3 Suavizado La siguiente modificación consiste en realizar un "suavizado" de las concentraciones celulares y algunos factores biológicos partícipes del proceso. En concreto, se define una vecindad en torno a cada punto de integración y se realiza un promediado dependiente de la distancia entre los puntos de integración pertenecientes a la vecindad y el punto de integración origen ( PIo ), dando una mayor importancia a las contribuciones de los puntos de integración más cercanos. La justificación a nivel biológico de esta modificación reside en la incoherencia que se produce si en una zona del hueso hay una transición brusca en la concentración de algún tipo de célula o factor biológico. Es decir, si en un punto del hueso la concentración de OCa es máxima, en sus proximidades no puede presentarse zonas con concentraciones nulas ni reducciones bruscas. Debe producirse una transición suave y continua. Este cambio permite obtener dicha transición deseada entre los elementos que conforman el dominio, y es una primera aproximación para la obtención del Non-local Model. En el proceso se definen distintos subdominios abordando el estudio de manera no local pero no se atiende a los sucesos que tienen lugar en los mismos (zonas ya activas de reabsorción). Para la implementación del cambio se parte del algoritmo NAM (Neighbourhood averaging method) propuesto por Calvo-Gallego et al. [ 42 ]. Este algoritmo promedia las tensiones de cada vecindad y con estos valores se aplican los mecanismos desarrollados en el Capítulo 3 para obtener las variables de estado (STATEV). Con este algoritmo se evitan ciertos fenómenos de inestabilidad como el problema conocido como checkerboard y problemas de convergencia en los resultados. La modificación realizada cambia el proceso indicado, ya que, el "suavizado" se realiza a las variables de estado ya calculadas. En primer lugar se define la vecindad de cada punto. Esta será una esfera con un radio especificado. En la Figura 4.1 se representa el subdominio definido para un punto de integración. Figura 4.1 Representación de la vecindad de un punto dentro del dominio completo [43]. 34 Capítulo 4. Modelo local Durante el proceso de "suavizado" se realiza la llamada a varias funciones que trabajan con las coordenadas de los puntos de integración y las vecindades definidas. Las coordenadas se obtienen de la variable COORDS generada en SDVINI, subrutina de Abaqus comentada anteriormente. A continuación, se explica la finalidad de esas funciones que facilitará el entendimiento del proceso global. Función: NEIGHBOURHOOD1 Esta función calcula el número de puntos de integración que hay en la vecindad de cada punto y devuelve el valor máximo en la variable (NGP_NEIGHB2). Para determinar si un punto pertenece a una vecindad, con las coordenadas de los puntos conocidas, se calcula la distancia entre el punto origen (centro de la esfera) y el resto de puntos. Si el valor obtenido es inferior al radio definido, el punto de integración pertenecerá a dicho subdominio. La función realiza dos bucles anidados recorriendo los puntos de integración, el bucle exterior indicará el punto origen y el interior recorrerá el resto de puntos. En el bucle se calcula la distancia entre los dos puntos y, si se cumple la condición indicada anteriormente, se incrementa una variable que actúa como contador (SUMAT). Esta variable es un vector de dimensión igual al número total de puntos de integración presentes en el dominio completo y cada componente se asocia a un punto de integración. Al final las componentes de dicho vector serán el número de puntos de integración presentes en la vecindad asociada a cada punto. En la Figura 4.2 se refleja dicho vector. Finalmente, para obtener el valor máximo del vector SUMAT se utiliza la función MAXVAL. Figura 4.2 Vector que contiene el número de puntos de integración vecinos ( NPIv ) asociados a cada punto de integración (PIo). Siendo n el número total de PI en el dominio. El valor obtenido permite el correcto dimensionamiento de la variable NEIGHB, cuya definición y función se desarrollará a continuación. La llamada a NEIGHBOURHOOD1 tiene lugar en UEXTERNALDB justo después de SDVINI (Véase Figura 3.5). Función: NEIGHBOURHOOD2 Realiza un proceso similar al anterior recorriendo los puntos de integración con un doble bucle anidado calculando la distancia, pero en este caso si se cumple la condición de la distancia se guarda en una matriz (NEIGHB) el índice del punto de integración vecino ( PIv ) y la distancia entre el PIo y el PIv ( DPIo−PIv ). El resultado obtenido es una matriz que en cada fila contiene los puntos de integración vecinos y sus distancias. El número de los puntos de integración se encuentra en las columnas impares y la distancia en las pares. Si una vecindad no contiene más puntos la fila se completa con ceros. En la Figura 4.3 se representa una fila de dicha matriz de ejemplo. La matriz NEIGHB se define con las siguientes dimensiones: filas igual al número total de puntos de integración presentes y columnas de dimensión el doble del número máximo de puntos de integración presentes en una vecindad. Esto se debe a que guarda el índice asociado al punto y la distancia. También se realiza su llamada en UEXTERNALDB justo después de SDVINI (Véase Figura 3.5), pero se incluye una condición para que solo se use al comienzo del primer STEP de la simulación. 4.1 Modificaciones 35 Figura 4.3 Fila de la matriz NEIGHB asociada al PI "i". Siendo nm el número máximo de PI en una vecindad. Función: AVERAGE Con esta se obtiene el promedio de las variables en función de la distancia. Para ello, dado un punto de integración, p, el conjunto de puntos que forman parte de su vecindad, Np , y la variable que se "suaviza", x, el calculo se realiza: ¯xp=∑Np i=ixi·wi(dip) ∑Np i=iwi(dip) siendo wi(dip) una función peso dependiente de la distancia entre el punto pi y el punto origen p . Esta se define a partir de una función sigmoidal: wi(dip) = 2·Rγ Rγ+dγ ip −1 donde R es el radio que define la vecindad y γ es un factor de forma. El valor de los parámetros se adjunta en el Apéndice A. Desarrollo del proceso Con las funciones anteriores definidas se explica el proceso global seguido. En primer lugar, se realiza un análisis que permita obtener el número máximo de puntos de integración presentes en una vecindad para dimensionar correctamente la matriz NEIGHB. Esto permite reducir el tiempo computacional y ahorrar memoria. Después de este análisis no será necesaria la llamada a la función NEIGHBOURHOOD1, a no ser que se cambie la malla de elementos finitos. El "suavizado" tiene lugar en la UMAT, estando ya definida la matriz NEIGHB. Para este proceso se realiza un bucle que recorra el número de variables que se desea promediar. Antes de producirse la llamada a la función AVERAGE, se asigna el radio y factor de forma asociado a cada variable. Finalmente, se realiza la llamada mencionada y se guarda el valor devuelto en la variable de estado correspondiente. Al final de la UMAT se guardan las STATEV en una variable temporal, SDV1. Esta variable tendrá una fila asociada a cada PI presente en el dominio, y una columna por cada variable de estado, o STATEV, que se promedia. Su definición se realiza para poder disponer del valor de las STATEV de los PIven cada incremento. Es interesante destacar, que el programa Abaqus almacena el valor de las STATEV al final del incremento solo si los resultados han convergido. Si esto no ocurre, repite el incremento con los valores del equilibrio calculados en el paso anterior. Por ello, al terminar un incremento con los resultados convergidos, en UEXTERNALDB (LOP=2) (ver Figura 3.5), se almacenará el valor de SDV1 en otra variable, SDV2, con la que trabaja la función AVERAGE. Esto se realiza debido a que si en un incremento los resultados no convergiesen, y trabajase solo con SDV1, el incremento comenzaría de nuevo pero los valores almacenados en la variable temporal SVD1 no corresponderían a los de equilibrio calculados en el incremento anterior. Al trabajar con la variable SDV2 si los 36 Capítulo 4. Modelo local resultados nos convergen, se sobreescribe la variable SDV1 pero se mantienen los valores del último estado de equilibrio guardados en la primera. En la Figura 4.4 se muestra un esquema del proceso de "suavizado". Figura 4.4 Diagrama de flujo del proceso de "suavizado". 4.2 Análisis de las simulaciones del Modelo local 37 4.2 Análisis de las simulaciones del Modelo local Para el análisis del modelo se han estudiado los efectos de la variación de distintos parámetros, en concreto: la velocidad de degradación del TGF- β ( ˜ DTGF-β ); la velocidad de apoptosis de los OCa (AOCa) y el factor que regula la producción de RANKL debido al daño (Πdam act ). 4.2.1 Simulaciones iniciales En primer lugar, se realizan simulaciones con las modificaciones comentadas variando las condiciones iniciales de las concentraciones celulares implicadas en la reabsorción y formación de tejido óseo y el valor de los umbrales. En este apartado se incluirán varias figuras que representan la evolución de las variables de estudio para cinco elementos situados en una dirección radial del cilindro. Estos se encuentran situados por orden de numeración, siendo el elemento 4025 interior y el 4153 el exterior. El hecho de no representar la evolución de las variables en todos los elementos se debe a que en los resultados obtenidos se distinguen cinco comportamientos distintos. Sumando a esto, para poder apreciar bien las evoluciones y no sobreponer varias evoluciones se ha optado por representar los resultados solo para cinco elementos. Variación de las concentraciones iniciales de OBayOCa Se realizan simulaciones con las concentraciones iniciales preestablecidas y con valores nulos de las mismas. Los valores de las concentraciones iniciales se encuentran reflejadas en el Apéndice A incluido a final de este trabajo. A continuación, se incluyen varias figuras que representan la evolución de la concentración de OBa y OCa para los cinco elementos mencionados anteriormente. a) Condiciones iniciales nulas b) Condiciones iniciales no nulas Figura 4.5 Evolución de la concentración de OCa en 5 elementos dispuestos en una dirección radial. El elemento 4025 es el interior y el 4153 el exterior, los demás son intermedios y se sitúan por orden de numeración. 38 Capítulo 4. Modelo local Figura 4.6 Evolución de la concentración de OCa en 5 elementos dispuestos en una dirección radial. El elemento 4025 es el interior y el 4153 el exterior, los demás son intermedios y se sitúan por orden de numeración. a) Condiciones iniciales nulas b) Condiciones iniciales no nulas Figura 4.7 Evolución de la concentración de OCa en 5 elementos dispuestos en una dirección radial. El elemento 4025 es el interior y el 4153 el exterior, los demás son intermedios y se sitúan por orden de numeración. Se ha representado inicialmente la evolución de OCa para concentraciones iniciales nulas en la Figura 4.6 para poder apreciar la evolución completa. A continuación, en la Figura 4.7.a se representa la misma evolución ampliada para facilitar la comparación con la Figura 4.7.b, en la que las que las condiciones iniciales no son nulas. En la Figura 4.5 y Figura 4.7 se puede apreciar como para distintas concentraciones iniciales los valores tienden a oscilar entre los mismos valores. Tanto para valores iniciales nulos de OCa y OBa como para los propuestos anteriormente, se observa un transitorio inicial y una vez estabilizado una evolución similar. Esto es interesante, ya que indica que los comportamientos tienden en el tiempo a ser similares para cualquier condición inicial. 4.2 Análisis de las simulaciones del Modelo local 39 Variación de los umbrales En este apartado se realizan varias simulaciones modificando el valor de los umbrales. En la Tabla 4.1 se indican los valores definidos para las simulaciones que se van a comparar. Tabla 4.1 Valores definidos de los umbrales para las simulaciones a comparar. Nombre Condiciones iniciales Umbral OCa Umbral OBa S1 Nulas 1,00 ·10−55.50 ·10−5 S2 No nulas 1.00 ·10−55.50 ·10−5 S3 Nulas 1.00 ·10−45.50 ·10−4 S4 No nulas 1.00 ·10−45.50 ·10−4 A continuación, se comparan la evolución de tejido óseo reabsorbido (RBV) y formado (FBV) para dos simulaciones que presentan umbrales iguales pero condiciones iniciales distintas. a) S1: RBV b) S2: RBV a) S1: FBV b) S2: FBV Figura 4.8 Evolución del RBV y FBV en 5 elementos dispuestos en una dirección radial. El elemento 4025 es el interior y el 4153 el exterior, los demás son intermedios y se sitúan por orden de numeración. 46 Capítulo 4. Modelo local FBV: Formed bone volume a) ˜ DTGF-β 10 b) ˜ DTGF-β 5 c) ˜ DTGF-β 2.5 Figura 4.14 Evolución del FBV en el elemento interior 4025 para distintos valores del parámetro ˜ DTGF-β. Finalmente, se analiza el volumen óseo generado. Lo primero a destacar es la presencia de las fases de quiescencia y formación buscadas. Se analiza los efectos que se producen en las mismas con la variación del parámetro estudiado. El mayor valor del parámetro que modela la velocidad de degradación del TGF- β provoca una menor formación de tejido óseo. Esto se debe a que se asume que el FBV es proporcional con la concentración de OBa y, como se ha comentado anteriormente, un mayor valor de la velocidad de degradación del TGF- β provoca una menor concentración de OBa . Por ello en la Figura 4.14.c no se generan ciclos de reabsorción. Al igual que con el RBV, para un valor mayor de la velocidad de degradación del TGF- β aumentan los ciclos (véase Figura 4.14) y las fases de quiescencia se acortan. También para los dos primeros valores tomados se observa uniformidad en las fases de formación y reposo (ver Figura 4.14.a y Figura 4.14.b). Conclusiones de la variación de la velocidad de degradación del TGF-β Las ideas obtenidas sobre la variación del parámetro ˜ DTGF-β a partir de los resultados se pueden resumir en: - Cambios apreciables en la producción de osteoblastos activos ( OBa ) aumentando su producción por ciclo al disminuir el valor de ˜ DTGF-β . Por otro lado, aumenta el número de ciclos de producción al aumentar su valor. 4.2 Análisis de las simulaciones del Modelo local 47 - Afecta levemente a la producción de osteoclastos activos ( OCa ) en cada ciclo. Pero reduce los ciclos presentes conforme se reduce la velocidad de degradación, disminuyendo la producción neta de OCa. - La concentración de TGF- β varía con la velocidad de degradación del mismo, aumentando el número de ciclos al reducir la velocidad de degradación. - La cantidad de RBV en cada ciclo es similar para los valores probados, disminuyéndose el número de ciclos con el menor valor del parámetro estudiado. - La cantidad de FBV se reduce en cada ciclo al aumentar la velocidad de degradación. Por el contrario, la cantidad de ciclos aumenta al aumentar ˜ DTGF-β. - Tanto para RBV y FBV los periodos de reposo o quiescencia acortan su duración en cada ciclo conforme el valor de ˜ DTGF-βes mayor. - En todas las evoluciones estudiadas se observa que para el mayor valor estudiado de la velocidad de degradación del TGF- β aparece cierta inestabilidad, lo que hace que se descarte ese valor del parámetro. 48 Capítulo 4. Modelo local 4.2.3 Variación velocidad Apoptosis OCa Para ver las consecuencias de esta modificación se ha seguido el mismo proceso que para la degradación del TGF- β : probar distintos valores de la velocidad de apoptosis. Se han empleado: AOCa/10 , AOCa/5 y AOCa/2.5 . De nuevo, se realizan las simulaciones y con los resultados se establece el valor de los umbrales y se repiten las simulaciones. Los umbrales definidos son: OBaTH =7.5·10−4 y OCaTH =1·10−3. Osteoblastos activos a) AOCa 10 b) AOCa 5 c) AOCa 2.5 Figura 4.15 Evolución de la concentración de OBa en dos elementos, 4025 interior y 4153 exterior, para distintos valores del parámetro ˜ DTGF-β. Para el estudio de esta concentración celular se ha representando la evolución para dos elementos, uno interior y otro exterior. La evolución de los otros tres elementos mencionados al comienzo del análisis abarca rangos intermedios entre los dos comportamientos representados. En primer lugar, se puede observar como los rangos abarcados para los dos elementos difieren al contrario de como ocurría para la variación de la velocidad del TGF-β(ver Figura 4.15). Para AOCa /10 se observa un aumento inicial muy notable (véase Figura 4.15.a) Para valores de la velocidad de apoptosis no tan pequeños este aumento se suaviza (ver Figura 4.15.b y Figura 4.15.c). También, se puede observar que los osteoblastos activos disminuyen al aumentar el valor de AOCa. Esto se debe a que a mayor valor de la apoptosis de los OCa menor cantidad TGF- β se liberará del tejido óseo. Y por lo comentado anteriormente, una menor concentración de TGF- β inhibe de manera indirecta la diferenciación de los OBaproduciéndose una menor cantidad de los mismos. 4.2 Análisis de las simulaciones del Modelo local 49 En los ciclos presentes en la evolución no se aprecia la uniformidad comentada para la variación de la velocidad de degradación del TGF-β. Osteoclastos activos a) AOCa 10 b) AOCa 5 c) AOCa 2.5 Figura 4.16 Evolución de la concentración de OCa en el elemento interior 4025 para distintos valores del parámetro ˜ DTGF-β. De nuevo, como ocurría en las anteriores simulaciones. La evolución de la concentración de OCa en los cinco elementos de estudio es similar, presentando cierto desfase entre ciclos. Por ello se representa únicamente un elemento. Se puede observar como al aumentar el valor de la velocidad de apoptosis de OCa aumentan las oscilaciones en la población de OCa (véase Figura 4.16). Esto se debe a que la apoptosis, o muerte celular programada, indica la cantidad de OCa "eliminada". Es coherente que al disminuir su valor aumente la concentración de OCa . Al reducir el valor de este parámetro se pretende aumentar la concentración promedio de OCa , así se superará con mayor frecuencia el umbral de reabsorción. Lo que supondrá fases activas de reabsorción más largas. La evolución para las tres simulaciones representada dista de la buscada debido a que se aprecian demasiados ciclos en periodos cortos de tiempo. Además, para los tres valores se aprecia inestabilidad en la evolución en distintos instantes temporales, lo cual no es deseado para la obtención de los ciclos de reabsorción porque presentarán una duración muy corta que no se asocia a la reflejada en la literatura. 50 Capítulo 4. Modelo local Concentración TGF-β a) AOCa 10 b) AOCa 5 c) AOCa 2.5 Figura 4.17 Evolución de la concentración de TGF- β en el elemento interior 4025 para distintos valores del parámetro ˜ DTGF-β. En la Figura 4.17 se puede observar mayores concentraciones para un menor valor de la apoptosis. Esto se debe a que un menor valor de la AOCa supone una mayor producción de OCa . Estos reabsorberán una mayor cantidad de volumen óseo liberando así una mayor cantidad de TGF- β . Matemáticamente lo comentado se refleja en la ecuación (4.2) y la ecuación (4.4). En la primera de ellas se observa que la concentración de TGF- β es proporcional a la concentración de OCa . En la segunda, en el primer término del segundo miembro se refleja de nuevo esa proporcionalidad con la concentración de OCa. Exceptuando las concentraciones producidas, la evolución de los ciclos es similar en todas las simulaciones, presentando demasiados ciclos en periodos cortos de tiempo. Tampoco se parecen a los ciclos tan interesantes que resultaban en las simulaciones anteriores en las que se variaba el parámetro ˜ DTGF−β. 4.2 Análisis de las simulaciones del Modelo local 51 RBV: Resorbed bone volume a) AOCa 10 b) AOCa 5 c) AOCa 2.5 Figura 4.18 Evolución del RBV en el elemento interior 4025para distintos valores del parámetro AOCa. En la Figura 4.18 se puede observar que conforme se aumenta el valor de la apoptosis osteoclástica se reduce la cantidad de volumen óseo reabsorbido. Esto se debe a que un aumento de la apoptosis disminuye la concentración de OCa . Siendo estas las células encargadas de realizar la reabsorción. La evolución presenta demasiados ciclos en periodos cortos de tiempo y no se observa uniformidad entre ellos. A esto se le suma los periodos de reposo casi inexistentes. En definitiva, la evolución presente dista de la interesante que resultaba al variar la velocidad de degradación del TGF-β. 52 Capítulo 4. Modelo local FBV: Formed bone volume a) AOCa 10 b) AOCa 5 c) AOCa 2.5 Figura 4.19 Evolución del FBV en el elemento interior 4025 para distintos valores del parámetro AOCa. Los resultados obtenidos para el volumen óseo no son tan interesantes como en el apartado anterior. Se puede observar que a mayor valor valor de la AOCa mayor cantidad de volumen óseo es formado (ver Figura 4.19). Esto se debe a que al reducir la apoptosis de los OCa aumenta la concentración de los mismos. Al aumentar estos, se libera una mayor concentración de TGF- β que inhibe de manera indirecta la diferenciación de los OBa . Siendo estas las células encargadas de la formación de tejido óseo y siendo FBV proporcional a dicha concentración en el modelo. Para el menor valor de la apoptosis de los OCa la evolución del FBV no presenta periodos de reposo (ver Figura 4.19.a). Mientras que para el valor intermedio si se producen algunos ciclos que apenas presentan entre ellos periodo de quiescencia (ver Figura 4.19.b). Por último, para el mayor valor de la apoptosis de los OCa se obtienen únicamente dos ciclos sin periodo de reposo entre ellos (ver Figura 4.19.c). Conclusiones de la variación de la velocidad de apoptosis de los OCa Los resultados obtenidos para la variación de la apoptosis de los OCa no han sido muy relevantes para el trabajo. Los cambios principales que se producen son en las concentraciones de las poblaciones celulares estudiadas, pero no alarga los ciclos en la evolución de las mismas, con lo que no afectará 4.2 Análisis de las simulaciones del Modelo local 53 periodos de quiescencia. Simplemente se reabsorberá o formará mayor o menor cantidad de tejido óseo en función del valor establecido para los umbrales. Con estas variaciones no se obtienen los ciclos interesantes observados con la variación de la degradación de la velocidad del TGF- β . En resumen, la variación de este parámetro no permite sacar ninguna conclusión que interese para el trabajo, pero sí para un mejor entendimiento del modelo. 4.2.4 Variación del factor regulador Πdam act Este parámetro se encuentra multiplicando a la concentración de Ot en la ecuación (3.12) para la obtención del RANKLe f f . Este se usa para el cálculo del RANKL en la ecuación (3.13). Al aumentar el valor del parámetro, aumenta la concentración de RANKLe f f y a su vez la producción de RANKL debido al daño. Siendo esto último la finalidad buscada en este apartado. Siguiendo la estructura de los apartados anteriores se estudian los cambios producidos para tres valores, en concreto: Πdam act ·102,Πdam act ·104yΠdam act ·107. Osteoclastos activos a) Πdam act ·102b) Πdam act ·104 c) Πdam act ·107 Figura 4.20 Evolución de los OCa en el elemento interior 4025 para distintos valores del parámetro Πdam act . En los resultados se obtienen variaciones significativas únicamente para la evolución de OCa . En la Figura 4.20 se puede observar como al aumentar el valor del parámetro los ciclos presentes se vuelven más regulares y uniformes (ver Figura 4.20.c), lo cual supone que los ciclos de reabsorción 54 Capítulo 4. Modelo local serás más estables. Realmente, entre la evolución presente en la Figura 4.20.a y Figura 4.20.b no se producen cambios relevantes, pero sí para el último valor probado. También es interesante destacar que para este último valor se satura la producción de RANKL debido al daño. Esto se refleja en la ecuación (3.13), si se multiplica por la concentración de RANKLe f f al numerador y denominador la ecuación queda: RANKL =  βRANKL +PRANKL βRANKL RANKLe f f +˜ DRANKL  ·A(4.6) siendo: A="1+OPG Kd,[RANKL−OPG] +RANK Kd,[RANKL−RANK] +ζ·Cp,den Kd,[RANKL−den]#−1 Atendiendo al denominador del segundo miembro de la ecuación (4.6) se observan dos términos. El primero de ellos es el cociente de una constante (velocidad de producción de RANKL, βRANKL ) entre la concentración de RANKLe f f . El segundo término es otra constante que modela la velocidad de degradación del RANKL ( ˜ DRANKL ), siendo esta de orden 10. La saturación ocurre debido a que al aumentar el valor de Πdam el cociente de la velocidad de producción de RANKL entre la concentración de RANKLe f f da como resultado un valor tres órdenes de magnitud inferior a la velocidad de degradación del RANKL. Por tanto, si se siguiese aumentando el valor de Πdam act no se produciría ningún cambio significativo. Finalmente, de estas simulaciones se obtiene que el aumento de la importancia del daño en la producción de RANKL provoca una evolución de la concentración de OCa con ciclos regulares, la cual resulta interesante para la obtención de las fases buscadas. 4.2.5 Variación conjunta de los parámetros estudiados Por último, se realizan varias simulaciones combinando variaciones de los parámetros estudiados. Se obtiene una simulación especialmente interesante, en ella los parámetros toman los valores: ˜ DTGF−β /10; AOCa /2.5; y Πdam act ·107 . A continuación, se representa la evolución de las variables de estado estudiadas para el mismo elemento interior que en los anteriores análisis. Figura 4.21 Evolución de la concentración de OBa en el elemento interior 4025 para ˜ DTGF−β /10, AOCa/2.5 y Πdam act ·107. 4.2 Análisis de las simulaciones del Modelo local 55 Figura 4.22 Evolución de la concentración de OCa en el elemento interior 4025 para ˜ DTGF−β /10, AOCa/2.5 y Πdam act ·107. Figura 4.23 Evolución de la concentración de TGF- β en el elemento interior 4025 para ˜ DTGF−β /10, AOCa/2.5 y Πdam act ·107. Para las dos concentraciones celulares y la biológica objeto de estudio, se pueden observar los ciclos buscados. Siendo estos uniformes al igual que su evolución a lo largo del tiempo (ver Figura 4.21, Figura 4.22 y Figura 4.23). Hay que destacar que inicialmente se produce un ciclo mayor para las tres concentraciones estudiadas, esto se debe a las condiciones iniciales de las concentraciones de OBa y OCa nulas impuestas. No es relevante por lo comentado anteriormente sobre el transitorio inicial que tiene lugar al comienzo de las simulaciones para cualquier valor de las concentraciones celulares. Al transcurrir días de simulación los valores de las variables de estudio tienden a ser similares. Por ello, no es relevante la información aportada por las gráficas para los primeros días de simulación. 5 Modelo no local Con las modificaciones anteriores se logra una mejora del modelo a nivel local, siendo el resultado más significativo la obtención de las fases de vida de una BMU para ciertas combinaciones de los parámetros. En este capítulo se aborda el estudio no local teniendo en cuenta el entorno cercano de cada punto de integración. Para ello se hace uso las vecindades definidas en el capítulo anterior y se estudiarán las concentraciones de OCa y las banderas de reabsorción. A estas se le impondrán más condiciones para su activación. En concreto, se pretende implementar las condiciones que deben darse en el entorno cercano a una BMU para que pueda esta activarse. Estas son: - La concentración de OCa en la BMU tiene que superar el umbral. Como se comentó anteriormente, no para todas las concentraciones de OCa se activa la reabsorción, es necesaria una concentración determinada para que ocurra el proceso. - La concentración de OCa tiene que ser la máxima de su entorno. Esto se debe a que una vez se activa una BMU los OCa de las demás BMUs presentes en el entorno cercano se desplazan para contribuir en la reabsorción llevada a cabo por la BMU activa.Por ello, la que presente una mayor concentración de OCa será la que se activará. El "suavizado" realizado en el capítulo anterior intenta simular esa redistribución producida, a parte de solventar los posibles problemas de inestabilidad y discontinuidades comentados. - Una BMU no puede activarse si ya hay otra activa en su entorno. Se debe a que en un entorno cercano del hueso no puede haber varias zonas de reabsorción activas. El primer paso para la obtención de este modelo no local es trasladarlas condiciones referidas anteriormente al modelo de elementos finitos. Las condiciones quedarían: - La concentración de OCaen el PI tiene que superar el umbral. - La concentración de OCaasociada a un PI tiene que ser la máxima de su vecindad. - Un PI no puede activarse si ya hay otro activo en su vecindad. Una vez obtenidos los resultados, se calcula la frecuencia de activación de las BMUs en el dominio de estudio. El cálculo de la misma se realiza obteniendo el número de veces que se activan los PI y dividiéndolo entre el volumen de estudio y entre el incremento de tiempo simulado. El valor obtenido se compara con la real, obtenido de la literatura. Esta frecuencia debe tomar un valor entre 0.05 y 0.005 #BMU mm3·día[9]. 63 64 Capítulo 5. Modelo no local 5.1 Implementación de las restricciones Para su implementación se parte de la modificación de los umbrales generada en el modelo local desarrollado en el Capítulo 4. En esta se definía una variable ( FOCa ) que toma valor unidad si se supera el umbral, o nulo si no se supera. De esta forma actúa como bandera de activación de la reabsorción. Dicha modificación implementa ya la primera restricción. El proceso seguido consiste en introducir el resto de restricciones mencionadas para la activación de la bandera de reabsorción. Para ello se hará uso de algunas variables definidas en el capítulo anterior. En concreto se hará uso de la matriz NEIGHB. Esta contiene en cada fila los puntos de Gauss vecinos y la distancia entre dichos puntos y el PI origen. En primer lugar se inicializa la bandera a cero y se introduce la condición del umbral. Si la concentración supera el umbral, la bandera toma valor unidad y comienza un bucle que recorre los PI vecinos. En dicho bucle se implementan las dos condiciones restantes en una, para ello se comprueba si la concentración de OCa es inferior a alguna de las concentraciones asociadas a los PI vecinos o alguna de sus banderas está activa (tiene valor unidad), en cuyo caso la bandera del PI origen se iguala a cero y se sale del bucle. En el caso de que ninguna de las dos últimas condiciones fuese cierta, el bucle se completaría hasta que no hubiese más PI en la vecindad tomando la bandera a la salida un valor unidad. La finalidad de activar la bandera del PIo antes de la entrada al bucle reside en optimizar el proceso, ya que si se incumple alguna de las condiciones se sale del bucle ahorrando tiempo de computación. Evidentemente, para la definición de la matriz NEIGHB y el número máximo de puntos de Gauss presentes en la vecindad son necesarias las funciones NEIGHBOURDHOOH1 yNEIGH- BOURDHOOD2 definidas en el Capítulo 4. Finalmente, para el cálculo de la frecuencia de activación se realiza el sumatorio de las veces que la bandera se activa en todos los PI del dominio de estudio y se divide entre el volumen estudiado (RVE) y entre el incremento de tiempo de simulación ( δt ). Este cálculo se realiza a partir de los datos obtenidos con Abaqus y el programa Matlab. La implementación es correcta, pero debido a la evolución de los OCa se ha encontrando una excepción para vecindades grandes en el planteamiento de las restricciones que aleja el comportamiento simulado del proceso de ROI real. La distribución de la concentración de OCa presenta un gradiente en dirección radial con el máximo localizado en el radio interior. Al producirse esta distribución, los elementos internos se activan pero los siguientes en ser analizados, en la dirección radial no, al no presentar la concentración máxima en la vecindad. Supóngase que en un determinado modelo hay 5 elementos en dirección radial, como se muestra en la Figura 5.1 en la que además se muestran las vecindades asociadas al primer y tercer elemento. Supóngase también por simplicidad que cada elemento consta de un único PI para así poder identificar elemento con PI. Con la evolución de OCa presente y el orden seguido por Abaqus para el cálculo, el elemento número 1 se activaría al superar el umbral, presentar la concentración máxima en su vecindad y no haber ningún elemento activo en la misma. A continuación, el elemento número 2 no se activa aunque supere el umbral porque el elemento 1 ya está activo. Tampoco presenta la concentración máxima en su vecindad (el elemento número 1 presenta una superior). Hasta aquí todo ocurre correctamente. Al llegar al elemento número 3, el umbral es superado y no hay ningún otro elemento activo en su zona pero la concentración de OCa presente en el elemento número 2 es superior, por tanto no se activaría. Los mismo ocurre para los elementos 4 y 5. La modificación funciona pero las condiciones impuestas no generan el comportamiento que presentaría una BMU en un dominio con esas condiciones. Lo que realmente tendría lugar sería la activación del elemento 3 y los OCa presentes en el elemento 2 contribuirían a la reabsorción en los elementos 1 y 3, como se ha comentado al comienzo de este capítulo. 5.2 Resultados y análisis de las simulaciones 65 Figura 5.1 Distribución de elementos en un radio del modelo con las vecindades de uno y tres. 5.2 Resultados y análisis de las simulaciones Con la problemática presente expuesta, se presentarán los resultados para varios valores del radio que define la vecindad. Es interesante recordar de nuevo las dimensiones características de los elementos comentadas en el capítulo anterior. En la dirección radial y axial los elementos presentan una longitud de 2 mm, en la dirección circunferencial la longitud aumenta conforme se alejan del eje del cilindro, presentando el elemento interior una longitud de 0.94 mm y el exterior de 2.94 mm (ver Figura 4.28). También es interesante reflejar que las distancias indicadas son entre nodos, no entre los PI en los que se realiza el estudio no local. Los elementos utilizados son lineales de tipo C3D8, estos presentan 8 nodos y 8 PI. Estos se localizan y numeran como refleja la Figura 5.2. Figura 5.2 Elemento C3D8 con sus nodos y puntos de integración numerados [44]. En un elemento lineal si las coordenadas de los nodos se sitúan en ±1 los PI en ±1/√3 . Es decir, los PI estarán separados una distancia menor a la que hay entre nodos. 66 Capítulo 5. Modelo no local Radio vecindad: R=5.0 mm Este radio de vecindad definido abarca dos elementos y medio, abarcando una cantidad de PI considerable en cada vecindad. En la Figura 5.3 se puede observar el gradiente radial presente en la concentración de OCaque fue mencionado anteriormente. Figura 5.3 Concentración de OCapara un radio de la vecindad igual a 5.0. A continuación, se muestran las zonas de activación para dos instantes determinados de la simulación. Como se puede comprobar en la Figura 5.4 únicamente se produce reabsorción en los elementos internos. Figura 5.4 Zonas de reabsorción para un radio de la vecindad igual a 5.0 mm. 5.2 Resultados y análisis de las simulaciones 67 Para esta simulación se obtienen los valores tomados por las banderas de reabsorción de los PI presentes en el dominio de estudio. Con estos se realiza el sumatorio de las veces que toman valor unidad ( veces que se activan) y se divide entre el volumen de estudio y el incremento de tiempo de simulación. Para la obtención del volumen se realiza el sumatorio de los volúmenes adscritos a los PI y el incremento de tiempo tendrá un valor de 100 días. Este cálculo es llevado a cabo con el programa Matlab. El resultado obtenido de la frecuencia de activación es: fact =0.0025#BMU mm3·día Este valor no se encuentra en el intervalo deseado. A pesar de ser inferior hay que tener en cuenta que únicamente se activan los elementos internos, si se lograse la activación coherente de los demás elementos este valor aumentaría, pudiéndose encontrar en el rango deseado. Si este valor no se encontrase en el intervalo deseado daría lugar a un nuevo análisis para ajustar los parámetros del modelo. En el siguiente capítulo se indicarán algunas líneas futuras que podrían solventar este problema. Radio vecindad: R=0.9 mm Con las dimensiones conocidas y la disposición de los PI en un elemento (ver Figura 5.2) se puede observar que el radio definido es tan pequeño que solo abarca algunos PI presentes en el mismo elemento. El uso de un radio pequeño simplemente es para verificar el funcionamiento del código planteado. Para ello se estudian los PI pertenecientes a un elemento interior, a partir de Abaqus se obtiene numeración de los nodos, su conectividad, y por tanto la numeración de los PI presentes en el elemento. A continuación se van a estudiar un elemento interior, el número 4025. Realizando el proceso comentado anteriormente se identifica la localización de los PI del elemento. En concreto, los PI 1 y 2 se encuentran en la dirección circunferencial. Al estar separados los nodos una longitud de 0.94 mm (Véase Figura 4.28), los PI estarán separados una distancia inferior a la del radio de la vecindad definido. Figura 5.5 Evolución de la bandera de reabsorción para dos PI pertenecientes a la misma vecindad. 68 Capítulo 5. Modelo no local En la Figura 5.5 se puede observar que nunca se activan en el mismo instante de tiempo, lo que nos indica un correcto funcionamiento del código generado y un comportamiento muy parecido a la realidad. A continuación se muestra la distribución de la concentración de los OCa en la Figura 5.6. Figura 5.6 Concentración de OCapara un radio de la vecindad igual a 0.9 mm. En este caso no ocurre esa distribución en forma de gradiente debido a que el radio de vecindad y suavizado es pequeño. Se pueden apreciar algunas discontinuidades, esto se debe al radio de la vecindad definido que es menor que el anterior y no consigue un suavizado adecuado. Las zonas de reabsorción se pueden observar en la Figura 5.7. Como se puede comprobar algunos elementos consecutivos están activos, lo cual es coherente con el radio definido pero no con el comportamiento real. Figura 5.7 Zonas en las que se produce reabsorción para un radio de vecindad igual a 0.9 mm. Para esta simulación se obtiene el siguiente valor de la frecuencia, cuyo resultado entra dentro del intervalo deseado, aproximándose la simulación al valor real. fact =0.0075#BMU mm3·día 5.2 Resultados y análisis de las simulaciones 69 Finalmente, se incluye una tabla con los distintos radios de vecindad probados y el valor de la frecuencia de activación que presentan. Tabla 5.1 Radios de vecindad probados y valor de la frecuencia que presentan. Rvecindad [mm] Tiempo [días] δt[días] fact #BMU mm3·día 0,5 400 100 0.0466 0,7 400 100 0.0412 0,8 400 100 0.0385 0,9 400 100 0.0075 1,0 400 100 0.0030 2,0 400 100 0.0026 5,0 400 100 0.0025 En la Tabla 5.1 se puede observar el valor límite del radio de vecindad que presenta una frecuencia de activación dentro del intervalo deseado es 0.9 mm, por ello se ha comentado esa simulación con mayor profundidad en el presente capítulo. También se ha comentado la simulación que presentan un radio de vecindad igual a 5 mm, ya que este se asocia a la distancia entre BMUs que permite activación de las mismas reflejada en la literatura. En este capítulo se ha desarrollado un modelo que realiza un estudio no local. En este estudio se han implementado las distintas condiciones que deben darse en el entorno de una BMU para su activación. Los resultados obtenidos han sido parcialmente satisfactorios, pero han permitido encontrar una excepción en el funcionamiento que permitirá desarrollar futuras investigaciones. Finalmente, se ha calculado la frecuencia de activación de las BMUs presentes en el dominio de estudio para comparar el valor obtenido con el reflejado en la literatura. 6 Conclusiones y líneas futuras En este capítulo se exponen la conclusiones obtenidas durante el desarrollo del trabajo y las líneas futuras que se pueden derivar del mismo. 6.1 Conclusiones A partir de los resultados obtenidos del modelo local se puede afirmar que las modificaciones introducidas en el código de la subrutina UMAT permiten la obtención de las fases presentes dentro del periodo de vida una BMU. Tanto para la reabsorción como para la formación, se obtienen etapas activas y estados de reposo o quiescencia, siendo este uno de los objetivos buscado y logrado. Con el estudio de sensibilidad realizado en el mismo capítulo, se obtiene conocimiento de los efectos que produce la variación de distintos parámetros, así como los valores que permiten ajustar mejor el comportamiento del modelo al proceso real con las modificaciones implementadas. Para el modelo desarrollado los valores estudiados que mejor adaptan las simulaciones al proceso real son: - Reducción de la velocidad de degradación del TGF-βinicial: ˜ DTGF−β=˜ Dinicial TGF−β/10 =0.2[día−1] - Reducción de la velocidad de apoptosis de los OCa: AOCa=Ainicial OCa/2.5=0.4[día−1] - Aumento del parámetro que modela la producción de RANKL en función del daño (Πdam act ): Πdam act =Πdam,inicial act ·107=1 Con el modelo no local se ha logrado una activación de los PI atendiendo a su dominio cercano, simulando las restricciones que se dan en el entorno de una BMU para su activación. El funcionamiento del modelo es correcto, pero gracias a la implementación se ha descubierto un caso particular en el que no funciona siendo el caso de un gradiente continuo en la concentración de OCa . Esto ha servido para plantear otras estrategias, que se enuncian en Líneas futuras (Apartado 6.2), que pueden conducir a resultados más cercanos a la realidad. Los resultados obtenidos para este modelo también sirven para resaltar la importancia de un estudio previo a la hora de abordar el problema y el desarrollo de la resolución. En este caso, un planteamiento incompleto de las condiciones de activación de las BMUs ha desembocado en unos resultados que solo son parcialmente satisfactorios. 71 78 Bibliografía [37] J. M. Garcia-Aznar, T. Rüberg, and M. Doblaré, “A bone remodelling model coupling microdamage growth and repair by 3D BMU-activity,” Biomechanics and modeling in mechanobiology, vol. 4, pp. 147–67, 12 2005. [38] C. Pattin, W. Caler, and D. Carter, “Cyclic mechanical property degradation during fatigue loading of cortical bone,” Journal of Biomechanics, vol. 29, no. 1, pp. 69–79, 1996. [Online]. Available: https://www.sciencedirect.com/science/article/pii/0021929094001561 [39] J. Martínez-Reina, J. García-Aznar, J. Domínguez, and M.Doblaré, “A bone remodelling model including the directional activity of BMUs,” Biomech Model Mechanobiol, vol. 8, pp. 111–127, 2009. [Online]. Available: https://doi.org/10.1007/s10237-008-0122-5 [40] J. D. Currey, “Tensile yield in compact bone is determined by strain, post-yield behaviour by mineral content,” Journal of Biomechanics, vol. 37, no. 4, pp. 549–556, 2004. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0021929003003294 [41] D. D. Marathe, A. Marathe, and D. E. Mager, “Integrated model for denosumab and ibandronate pharmacodynamics in postmenopausal women,” Biopharmaceutics Drug Disposition, vol. 32, no. 8, pp. 471–481, 2011. [Online]. Available: https://doi.org/10.1002/bdd.770 [42] J. L. Calvo-Gallego, P. Pivonka, J. M. García-Aznar, and J. Martínez-Reina, “A novel algorithm to resolve lack of convergence and checkerboard instability in bone adaptation simulations using non-local averaging,” International Journal for Numerical Methods in Biomedical Engineering, vol. 37, no. 2, p. e3419, 2021. [Online]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1002/cnm.3419 [43] Wikipedia. (2021) Entorno (matemática). [Online]. Available: https://es.wikipedia.org/wiki/ Entorno_(matem%C3%A1tica)#/media/Archivo:Neighborhood_illust1.svg [44] MIT. (2021) Eight-node brick element (c3d8 and f3d8). [Online]. Available: http: //web.mit.edu/calculix_v2.7/CalculiX/ccx_2.7/doc/ccx/node26.html