Activación discreta de BMUs en modelos de remodelación ósea interna

Equation Chapter 1 Section 1 Trabajo Fin de Máster Máster en Ingeniería Industrial Activación discreta de BMUs en modelos de remodelación ósea interna Autor: Elena Ramírez Solano Tutor: Javier Martínez Reina Departamento de Mecánica Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2017 Introducción 2 Trabajo de Fin de Máster Máster en Ingeniería Industrial Activación discreta de BMUs en modelos de remodelación ósea interna Autor: Elena Ramírez Solano Tutor: Javier Martínez Reina Profesor Contratado Doctor Departamento de Mecánica Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2017 Introducción 4 INDICE INDICE 5 1 Introducción 7 2 Biología ósea 9 2.1. INTRODUCCIÓN 9 2.2. COMPOSICIÓN DEL HUESO 9 2.2.1. Matriz ósea 9 2.2.2. Células óseas 12 2.3. TIPOS DE TEJIDO 13 2.3.1. Cortical y trabecular 13 2.3.2. Lamelar y fibroso 15 2.4. ESTRUCTURA DEL HUESO 16 2.5. REMODELACIÓN ÓSEA 18 2.5.1. Remodelación ósea externa 18 2.5.2. Remodelación ósea interna 19 2.6. PROPIEDADES MECÁNICAS 23 3 Modelos de remodelación ósea interna 27 3.1 INTRODUCCIÓN 27 3.2 MODELOS FENOMENOLÓGICOS ISÓTROPO DE STANFORD 27 3.3 MODELO MECANOBIOLÓGICO ISÓTROPO DE REMODELACIÓN ÓSEA BASADO EN LA ACTIVIDAD DE LAS BMUS 32 3.3.1 Introducción 32 3.3.2 Descripción matemática de la actividad de las BMUs 34 3.3.3 Crecimiento el daño 40 3.3.4 Activación de BMUs 42 3.3.5 Reabsorción y formación ósea 45 3.3.6 Proceso de mineralización 46 3.3.7 Degradación de las propiedades a fatiga con el contenido de mineral 49 3.4 MODELO MECANOBIOLÓGICO ANISÓTROPO DE REMODELACIÓN ÓSEA CON TRATAMIENTO DIFUSIVO DE LAS BMUS 51 3.4.1 Introducción 51 Introducción 6 3.4.2 Variables básicas del modelo 51 3.4.3 Actividad de las osteonas. Variacion del tensor MIL en el hueso cortical 56 3.4.4 Balance local direccional de hueso 60 3.4.5 Evolución del daño en el caso anisótropo 62 3.4.6 Avance de las BMUs 63 Reparto de la activiad de las BMUs 67 3.4.7 Propiedades elásticas del hueso 68 4 Modificación del modelo mecanobiológico isótropo 71 4.1 INTRODUCCIÓN 71 4.2 MODIFICACIONES PRINCIPALES 71 Evolución del daño 71 Constantes experimentales 72 Cálculo malla MEF 73 4.3 RESULTADOS 73 4.3.1 Daño localizado a tracción 73 4.3.2 Daño aleatorio a tracción 76 4.3.3 Daño homogéneo a flexión 78 5 Modificación del modelo mecanobiológico anisótropo 81 5.1 INTRODUCCIÓN 81 5.2 ACTIVIDAD DISCRETA DE LAS BMUS 82 5.3 DEFINICIÓN DE VARIABLES 86 5.3.1 CANDIDATOS 86 5.3.2 BMUPG 88 5.3.3 SENAL1 89 5.4 CONDICIONES DE ACTIVACIÓN DISCRETA 89 5.5 RESTRICCIONES Y PROBLEMÁTICA DEL MODELO 92 5.6 RESULTADOS 94 6 Conclusiones y desarrollos futuros 103 6.1 CONCLUSIONES 103 6.2 DESARROLLOS FUTUROS 104 Referencias 107 1 INTRODUCCIÓN En las últimas décadas la ingeniería ha ido tomando un papel de gran importancia en muchos ámbitos médicos. Con recursos ingenieriles es cada vez más frecuente dar respuesta a problemas relacionados con la salud y más usual emplear herramientas y planteamientos extraídos de la ingeniería para enfocar u optimizar distintos aspectos y cuestiones de la medicina. Un ejemplo de lo expuesto anteriormente y tema principal de este proyecto surge de la necesidad de conocer y prever el comportamiento del hueso y de sus propiedades mecánicas en el tiempo. Conocer estas variaciones incrementa de forma considerable las probabilidades de éxito en el diseño y la operación de una prótesis a un ser vivo. Esta prótesis tendrá una u otras características que dependerá de la evolución que vaya a seguir el hueso, evitando con estos conocimientos que este sea dañado. Para conocer el comportamiento del hueso es importante comprender todos los factores que afectan e influyen en la remodelación ósea interna, proceso mediante el cual el hueso va regenerando las grietas y el daño que surge en él tras estar sometido a diversas cargas o tras un tiempo de desuso. Son varios los modelos matemáticos creados en los últimos tiempos para plasmar este fenómeno. Estos modelos se implementan informáticamente en forma de lo que conocemos como subrutina, que no es más que un código informático que sirve de medio para resolver numéricamente el problema de la evolución del hueso. Una vez resuelto, con los modelos actuales se siguen observando diferencias entre la evolución de las propiedades del hueso que nos dicta la resolución numérica del modelo y aquellas que experimentalmente sabemos que tienen lugar. El presente proyecto surge como una propuesta de Trabajo de Fin de Máster con el objetivo de aproximar un poco más a la realidad los resultados aportados por los modelos de remodelación ósea intrena, eliminando algunas hipótesis establecidas que inicialmente simplificaban los cálculos pero que no representan de manera fiel el comportamiento del hueso, por lo que se distorsionan los resultados obtenidos de los modelos. El objetivo principal de este proyecto es discretizar la acción de las BMUs, el conjunto de células que efectúa el proceso de remodelación ósea. Es decir, si en un punto del hueso se requiere la puesta en marcha de la actividad de estas células esto se realizará por una única BMU y exclusivamente ahí, inhibiendo la posibilidad de que se activen y acudan fracciones de BMU de otras zonas cercanas, como se establecía en modeos previos. Este proyecto tiene un antecendente claro y es el Trabajo Fin de Grado llamado “Modelo de remodelación ósea interna discreto de BMUs” [64]. El objetivo de este trabajo era el mismo que el citado previamente, discretizar la actividad de las BMUs, pero sobre un modelo de remodelación ósea interna isótropo, es decir, que considera que el comportamiento del hueso es igual en todas las direcciones. En este proyecto, se hicieron grandes avances sobre la materia sin embargo los resultados Introducción 8 obtenidos no fueron del todo representativos respecto a lo que se pretendía obtener, por ello, en este trabajo se ha mejorado la rutina que plasmaba este modelo y depurado los errores de esta pudiendo obtenerse unos resultados más fieles para distintas casuísticas del estado del hueso. Por otra parte, considerar el hueso isótropo es en sí misma otra hipótesis que aleja algo más los resultados de la realidad que se sucede en el hueso, por ello el principal objtevo de este proyecto se ha ejecutado sobre un modelo de remodelación ósea interna que considera, por una parte que el hueso es anisótropo y su comportamiento varía según la dirección, y que, a diferencia del modelo usado en el TFG [64] las BMUs no son entidades estáticas, si no que a lo largo de su vida van desplazándose hacia zonas que requieran o vayan a requerir reparar su daño. Este modelo, es un modelo mucho más complejo, en el que intervienen numerosas variables y para el cual su programación se caracteriza por tener un nivel de dificultad elevado. Pero tras la ejecución de las modificicaciones se obtendrá un modelo mecanobiológico, que considera el hueso anisótropo, el movimiento de las BMUs, y que estas actúan de una manera discretizada, lo cuál se traduce en una versión más representativa que las rutinas que plasmaban esto previamente y por tanto un gran paso respecto a la previsión realista del comportamiento del hueso a largo plazo. A lo largo de este trabajo se hará en primer lugar (capítulo 2) una introducción a la biología ósea en la que se hablará sobre la composición del hueso y su estructura, los tipos de tejidos, el proceso de remodelación ósea y las propiedades mecánicas del hueso. A continuación (capítulo 3) se hará un repaso de algunos de los distintos modelos de remodelación ósea interna que se han expuesto a lo largo del tiempo y que han servido de antecedentes y puntos de partida del trabajo actual, pasando por los fenomenológicos, que no tienen en cuenta los procesos biológicos, y por los modelos mecanobiológicos que si tratan de incluir las dependencias de dichos procesos. Se hará especial hincapié en dos de los modelos mecanobiológicos: aquel que considera el hueso como un material isótropo y en el que la actividad delas BMUs está discretizada y aquel que considera el hueso anisótropo y la difusión de las BMUs. En el capítulo 4 se expondrá y se explicarán con detalle cuales han sido las modificaciones implementadas sobre el Modelo de Remodelación Ósea Interna discreto de BMUs para obtener resultados más certeros, así como se mostran dichos resultados para algunos casos de carga y daño del hueso. A lo largo del capítulo 5 se desarrollará cómo se ha modificado el modelo de remodelación ósea anisótropo que considera la difusión de las BMUs para discretizar la actividad de estas. Cuáles han sido las condiciones establecidas y las bases de dichas condiciones. Así como las variables que se han creado y las subrutinas de trabajo principales del modelo. Se mostrarán los resultados obtenidos tras las modificaciones y los problemas y limitaciones encontradas en el transcurso de la investigación. Finalmente, en el capítulo 6 se hará un resumen sobre las conclusiones finales a las que se han llegado. Además, se desglosarán las posibles líneas de trabajo a seguir en un futuro. 2 BIOLOGÍA ÓSEA 2.1. INTRODUCCIÓN El cuerpo humano está constituido por 206 huesos que forman parte del sistema musculoesquelético. Junto con los músculos y los ligamentos permiten la locomoción del cuerpo. Las funciones principales de los huesos pueden resumirse en las siguientes: a) Estructural: proporcionan estructura y forma al cuerpo, dando soporte y protegiendo los órganos vitales. b) Locomotora: con la interacción de huesos, músculos y articulaciones haciendo posible el movimiento. c) Reserva de fósforo y calcio y almacén de médula ósea: esta última se encuentra rellenando los poros y es la principal fuente de glóbulos rojos y blancos de nuestro organismo. Pese a las ideas preconcebidas que podamos tener por su apariencia estática los huesos son componentes muy dinámicos de nuestro cuerpo que están en continuo proceso de desgaste y regeneración. En este capítulo se hablará con detalle de este fenómeno, pasando antes por una breve descripción de las características básicas del tejido oseo: su composición, las células óseas y su función, los distintos tipos de tejido y su estructura. Tras hablar de de esto y del proceso de remodelación ósea se hará un recorrido por las distintas propiedades mecánicas de los huesos y sobre la influencia de algunos parámetros como la cantidad de mineral, la porosidad o la microestructura en ellas. 2.2. COMPOSICIÓN DEL HUESO El tejido óseo es el tejido conjuntivo que forma los huesos. Es un tejido duro constituido por una matriz extracelular mineralizada conocida como matriz ósea y por las células óseas. 2.2.1. Matriz ósea La matriz ósea se divide en tres componentes principales: una componente orgánica, otra inorgánica y agua. La componente inorgánica o mineral representa en torno a un 65% del total en peso, la componente orgánica un 20% y el agua el 10% restante [1]. Según el tipo de hueso que tengamos (en los apartados que siguen se explicará con más detalle la variedad que existe) las composiciones citadas antes pueden variar como muestra la figura 2.1 a continuación. Biología ósea 16 Figura 2.6 Osteona en la que las sucesivas lamelas alternan la dirección de las fibras de colágeno 90◦. Tomada de Martin et al. [2]. Figura 2.7 Lamelas sucesivas que cambian la dirección de las fibras de colágeno progresivamente con un aspecto helicoidal. Tomada de Martin et al. [2]. b) Fibroso, formado rápidamente y de una forma poco organizada. Menos rígido y resistente porque las fibras de colágeno y los cristales de mineral están dispuestos aleatoriamente. Para compensar esto contienen más mineral que el lamelar 2.4. ESTRUCTURA DEL HUESO Atendiendo a su forma y sus funciones existen tres tipos principales de huesos: huesos largos, huesos cortos y huesos planos. En los huesos largos una dimensión predomina respecto a las otras dos, en los huesos cortos todas las dimensiones son del mismo orden y en los huesos planos dos dimensiones predominan frente a la tercera. Los huesos cortos suelen encontrarse agrupados y permiten fundamentalmente el movimiento de las extremidades. Algunos ejemplos de huesos cortos serían el tarso, las vértebras, etc. Los huesos planos tienen como función proteger zonas más delicadas como en el caso del cráneo el cerebro, por lo tanto, apenas soportan carga. Otros ejemplos de huesos planos serían el omóplato y las crestas del ilion. Figura 2.8 Tipos de hueso. Los huesos largos como el fémur, la tibia y el húmero soportan normalmente cargas muy elevadas y de muchos tipos según el hueso: flexión-compresión, flexión-tracción o flexión predominante. Fue Borelli [7], en el siglo XVII, el que explicó por qué estas cargas podían llegar a ser tan elevadas y es debido a que gracias a la contracción de los músculos que tiran de los huesos es posible el movimiento de las extremidades, pero la contracción que producen estos es muy limitada y es necesario amplificarla para obtener movimientos normales. Por ello los huesos hacen de palanca, pero no con ganancia mecánica si no todo lo contrario y por ejemplo, para levantar algo con la mano el bíceps, que se inserta en el radio debe ejercer una fuerza unas 9 veces mayor que el peso a levantar. Los huesos largos se diferencian en varias zonas como se puede ver en la siguiente imagen. Transversalmente se observa una capa exterior llamada periosteo y una capa que recubre el tejido óseo interiormente, el endosteo. Longitudinalmente se diferencian tres zonas: la epífisis que son los dos extremos y la diáfisis que es toda la zona central. La diáfisis esta formada por hueso cortical mientras que la epífisis es principalmente hueso trabecular. Figura 2.9 Partes de un hueso largo: (1) epífisis, (2) placa epifisaria, (3) metáfisis, (4) diáfisis, (5) cavidad medular, (6) periostio, (7) cartílago articular. Adaptado de Fung [8]. Biología ósea 18 2.5. REMODELACIÓN ÓSEA Los huesos se encuentran en continua renovación. Son estructuras dinámicas que desde que nacemos hasta que morimos van cambiando su forma, o regenerándose mediante el fenómeno que se conoce como remodelación ósea. Existen dos tipos de remodelación ósea: la remodelación ósea externa y la remodelación ósea interna. Nuestro interés se centra fundamentalmente en la interna, aún así antes de profundizar sobre ella se darán unas nociones básicas sobre la remodelación ósea externa. 2.5.1. Remodelación ósea externa La remodelación ósea externa es aquella que se da cuando los huesos durante la infancia aumentan su tamaño y van esculpiendo la forma que tendrán cuando sean adultos y también cuando debido a algún trauma o fractura el hueso necesita reabsorberse en algunas zonas y formarse en otras. A diferencia de la remodelación ósea interna en la externa la actividad de los osteoclastos y los osteoblastos es necesaria pero de forma independiente, la velocidad de la ROE decrece con la edad mientras que la ROI ocurre durante toda la vida y la externa produce cambios de forma y tamaño pero la interna no afecta a esos aspectos. Algunos ejemplos de remodelación ósea externa podemos verlos en la figura 2.10: Figura 2.10 Ejemplos de ROE: (a) Reducción del diámetro de la metáfisis. (b) Aumento del diámetro de la diáfisis. (c) Modificación de la curvatura de a diáfisis. Tomada de Martin et al. [2]. 2.5.2. Remodelación ósea interna La remodelación ósea interna es un proceso biológico que variará según el trabajo y las cargas a las que esté sometido el hueso y según las condiciones mecánicas que afecten a sus diferentes partes. El hueso es una estructura óptima, es decir, consigue la máxima rigidez con el mínimo peso [9]. La función de la remodelación ósea interna es preservar esta forma de optimización y además reparar el daño que el hueso sufre a causa de la actividad diaria a la que se encuentra sometido. Esta regeneración o renovación del hueso tiene lugar gracias a una asociación de células que llamada BMU, Basic Multicellular Unit [10]. Las BMUs son conjuntos de osteoclastos y osteoblastos que de una forma secuencial y organizada que se explicará a continuación realizan todo el proceso de remodelación ósea. Una BMU está formada por unos 10 osteoclastos y por varios cientos de osteoblastos [2]. Tanto en el hueso cortical como en el hueso trabecular el proceso de remodelación es similar con la diferencia de que la forma de las BMUs es diferente en el cortical y que mientras que en el hueso trabecular el proceso tiene lugar solo en la superficie de la matriz ósea en el cortical puede llegar a atravesarla en su avance posterior. En las siguientes figuras se observa el avance y proceso de una BMU tanto en el hueso cortical como en el trabecular. Figura 2.11 BMUs de hueso esponjoso. Adaptada de Parfitt [11]. Biología ósea 20 Figura 2.12 Esquema de una BMU de hueso cortical. Se observa en varias secciones la formación del canal de Havers. Adaptada de Fridez [13]. Figura 2.13 Esquema de la progresión de un BMU de hueso cortical. Tomada de García [14]. El proceso de remodelación ósea se considera una secuencia invariable de fases ARF (activación, reabsorción y formación), [11], aunque realmente el proceso tiene más fases a saber: activación, reabsorción, inversión, formación, mineralización y fin de la adaptación (figura 2.14). Activación La activación es la fase más importante de la remodelación ósea. Durante ella se mandan las señales necesarias para que los osteoclastos precursores se diferencien y se fusionen, dando lugar a las grandes células multinucleadas que son los osteoclastos maduros. Los osteocitos son protagonistas de esta fase puesto que como ya se ha mencionado ellos "sienten" la carga mecánica a la que está sometido el hueso y según la teoría inhibitoria de Martin, [15], transmiten una señal bioquímica a las células de borde que impide que se active la actuación de las BMUs. En el momento en el que se genera alguna anomalía y los osteocitos dejan de transmitir esa señal inhibitoria a las células de borde se da comienzo a la remodelación ósea activando BMUs. Cuando se dice que los osteocitos "sienten" la carga mecánica se hace referencia a que estos son muy sensibles a los esfuerzos tangenciales. Por los canículos y las lagunas que conforman el entramado óseo y donde se alojan los osteocitos circulan nutrientes que vienen del torrente sanguíneo. La carga a la que se somete el hueso produce unas deformaciones que genera unas diferencias de presionen en los canículos, dando lugar así a un flujo de nutrientes que los osteocitos detectan como tensiones tangenciales. Si se da el caso de que el flujo se vea afectado, ya sea por desuso o por haber acumulado mucho daño, los osteocitos dejan de sentir las tensiones tangenciales puesto que el flujo se interrumpe y dejan de enviar la señal a las células de borde que inhibe la activación de BMUs. Figura 2.14 Fases de la actividad de una BMU de hueso trabecular: (0) Activación, (1) Reabsorción, (2) Inversión, (3) Comienzo de la formación, (4) Formación de osteoide, (5) Mineralización, (6) Fin de la adaptación. Adaptada de Fridez [13]. Biología ósea 22 Reabsorción ósea Cuando a las células de borde deja de llegarles la señal transmitida por los osteocitos en ese punto, llamado foco de la BMU, da comienzo la remodelación ósea. Las células de borde se encargan de segregar unas sustancias llamadas RANK-L que interaccionando con el RANK de los osteoclastos precursores permiten que estos se diferencien y se fusionen dando lugar a los osteoclastos maduros que forman parte de la BMU que actúa en ese punto. Las células de borde además en el momento previo al inicio de la reabsorción cambian su forma aplanada de manera que se encogen y aumentan su espesor dejando libre parte de la superficie ósea para que los osteoclastos puedan comenzar el proceso. Los osteoclastos reabsorben dando lugar a formas distintas según el tipo de hueso en el que se este produciendo la remodelación. Si es un hueso cortical se origina un túnel de sección aproximadamente circular mientras que si es un hueso trabecular el hueco es lo conocido como laguna de Howship. En unos estudios hechos por Jaworski y Lok, [16], se midió la velocidad de avance de las BMUs en costillas de perro obteniendo 39±14 μm/día. El frente de reabsorción tiene forma prácticamente elipsoidal, con un semieje mayor de 300 μm (en costillas de perro) y un semieje menor (el radio de la osteona) de 100 μm . Eso quiere decir que en un determinado punto de remodelación se tardan 300/39, aproximadamente 8 días, en abrir por completo la cavidad. A dicho lapso de tiempo se le denomina período de reabsorción, 𝑇𝑅. La velocidad de erosión de los osteoclastos en dirección perpendicular a la de avance de la BMU, sería entonces 100/8≃12μm/día. Respecto a cuanto dura el periodo de vida de los osteoclastos no es un tema que esté muy claro pero se cree que los osteoclastos continúan reabsorbiendo hueso y renovando sus núcleos durante todo el tiempo que está activa la BMU. Por lo tanto la apoptosis de los osteoclastos no llega hasta el fin de la actividad de las BMUs. Inversión La inversión es la fase durante la que tiene lugar la transición de actividad entre osteoclastos y osteoblastos. Este período en el ser humano dura unos 30 días. Formación ósea La fase de formación ósea comienza cuando los osteoblastos se sitúan sobre la zona que han reabsorbido los osteoclastos y depositan las primeras lamelas de osteoide. Los osteoblastos tienen su origen en las células mesenquimales, [17], y según Roberts et al. [18] los preosteoblastos se dividen por mitosis. Tras esta división y en función del nivel de tensiones una parte de los preosteoblastos se diferenciarán en osteoblastos. Estos volverán a dividirse por mitosis dando lugar a los osteoblastos definitivos. El proceso de diferenciación dura en torno a 60 horas. Si la BMU avanza a una velocidad de unos 40 μm/día, durante la diferenciación, el cono de reabsorción de la BMU ha Avanzado 100μm. Teniendo en cuenta el avance de la BMU durante otros procesos que tienen lugar antes de la formación el cono de reabsorción habrá avanzado entre 100−200μm. A diferencia de los osteoclastos se ha demostrado que los osteoblastos no avanzan con la BMU, si no que permanecen estacionarios en la sección de BMU, [19], rellenándola y cuando terminan se diferencian en osteocitos, células de borde o bien sufren apoptosis. Se calcula que durante el proceso de formación son necesarios unos 120 osteoblastos diarios para cada BMU. En el hueso cortical durante la formación ósea se va cerrando el hueco cilíndrico que los osteoclastos crearon y por ello la superficie en la que apoyan los osteoblastos es cada vez menor. Esto no sucede en el hueso trabecular puesto que en él la superficie de apoyo de los osteoblastos no tiene que disminuir necesariamente. Mineralización Durante la formación ósea lo que depositan los osteoblastos es osteoide, y este no tiene más que fase orgánica y agua. Por ello es necesario una posterior mineralización. Entre 10 y 20 días después de depositar el osteoide es cuando comienzan a aparecer los primeros cristales de hidroxiapatita desplazando parte del agua del osteoide. La mineralización se distingue en una fase primaria en la que se llega al 70% del contenido de mineral máximo, muy rápida y que dura unos pocos días y en una fase secundaria en la que se llega al 95% que puede durar entre 6 meses, según Parfitt, [20], y varios años, según Frost [21]. Durante esta fase la velocidad de deposición del mineral decrece exponencialmente. Fin de la adaptación Finalmente los osteoclastos sufren apoptosis y la BMU se detiene. Aún durante un cierto tiempo los osteoblastos seguiran rellenando las cavidades que correspondan de osteoide y luego sufrirán apoptosis, se convertirán en osteocitos o en células de borde. El osteoide depositado irá poco a poco mineralizándose aún cuando la BMU haya dejado de estar activa. 2.6. PROPIEDADES MECÁNICAS Las propiedades mecánicas del hueso, como era de esperar van a depender del tipo de hueso que se trate y por lo tanto de la microestructura que posea, la cual dependerá de las cargas que el hueso normalmente soporte. En este apartado se va a tratar la influencia de distintos parámetros en las propiedades mecánicas, exponiendo distintas correlaciones experimentales que representan dicha dependencia de estos parámetros en las propiedades mecánicas. Influencia del contenido de mineral Fueron Vose y Kubala, [23], los que observaron que el aumento de la fracción de ceniza en un pequeño porcentaje (del 63% al 71%) aumenta en casi 4 veces la resistencia del hueso. Biología ósea 24 Según Currey, [24], cuando se da un aumento excesivo de la fracción de ceniza el hueso se vuelve más frágil, sus propiedades a fatiga empeoran, es más débil a altas velocidades de deformación lo cual aumenta el riesgo de rotura. Se estima que la fracción de ceniza óptima oscile entre el 66-67%. Hernandez [22] obtuvo las siguientes correlaciones en humanos, que incluyen también la dependencia con la cantidad de material a través de la fracción volumétrica de hueso, 𝑣𝑏: 𝐸(𝑀𝑃𝑎)=84370 𝑣𝑏2,58· 𝛼2,74 𝜎𝑢𝑙𝑡(𝑀𝑃𝑎)=794,33𝑣𝑏1,92· 𝛼2,79 Influencia la porosidad La porosidad es un factor que afecta de forma considerable a las propiedades mecánicas de los huesos. Esta claro que cuanto mayor sea el volumen de poros en la estructura más débil será esta. Se suelen encontrar diferentes relaciones entre la porosidad y las propiedades mecánicas para el hueso cortical y para el trabecular. De hecho para un mismo valor de la porosidad el hueso cortical tendría mejores características que el hueso trabecular. Esto se debe a que la disposición del tejido no es la misma puesto que la estructura de las osteonas es major desde el punto de vista resistente. También se diferencian en el grado de mineralización debido a que el hueso trabecular tiene mayor tendencia remodelarse y por ello no cuenta con el tiempo suficiente para alcanzar un grado de mineralización alto. De esta forma el hueso trabecular tiene una rigidez y una resistencia menor que el trabecular. Beaupré et al., [25], obtuvieron experimentalmente unas expresiones del módulo elástico y del coeficiente de Poisson poniendo el límite de distinción entre hueso cortical y trabecular en una densidad aparente de 𝜌= 1,2 g/cm3. 𝐸={2014𝜌2,5 𝑠𝑖 𝜌≤1,2 𝑔/𝑐𝑚3 1763𝜌3,2 𝑠𝑖 𝜌≤1,2 𝑔/𝑐𝑚3 𝑣={0,2 𝑠𝑖 𝜌≤1,2 𝑔/𝑐𝑚3 0,32 𝑠𝑖 𝜌≤1,2 𝑔/𝑐𝑚3 Influencia de la microestructura Según el tipo de hueso la microestructura y su influencia en las propiedades mecánicas serán de una forma u otra (2.9a) (2.9b) ) (2.10) (2.11) En el hueso cortical el comportamiento de la capa de tejido es transversalmente isótropo y en la dirección longitudinal de mayor rigidez y resistencia. Esto se debe a que la dependencia principal de la microestructura es por la orientación de las osteonas que en los huesos largos lo están normalmente de forma longitudinal. También la orientación de las fibras de colágeno dentro de las osteonas. En las zonas donde predomina la tracción las fibras se orientarán longitudinalmente, mientras que en las de compresión predominante serán más frecuentes las de orientación transversal. En el hueso trabecular la microestructura está definida por la forma y orientación de los poros. El volumen de los mismos ya está recogido en la porosidad. Influencia del daño El nivel de daño es otro factor que afecta a las propiedades mecánicas del hueso. Cuando se habla de daño se hace referencia a una medida de las microgrietas que puede tener unhueso debido a la actividad diaria. El nivel de daño se considera un factor fundamental poner en marcha o no la remodelación ósea, lo cual es lógico puesto que si el hueso tiene más grietas será necesario remodelarlo con más apremio. Si el daño es elevado por lo tanto el hueso estará en continua remodelación evitando que se mineralice lo suficiente y haciendo esto que sea menos resistente. Hernandez et al. (2001) obtuvieron experimentalmente: 𝐸(𝑀𝑃𝑎)=84370𝑣𝑏2,58𝛼2,74(1−𝑑) (2.12) Modelos de remodelación ósea interna 32 32 Puesto que 𝑟󰇗 es un valor que depende del tiempo y de la superficie disponible en la que tener lugar la remodelación Martin, [36], definió una variable llamada superficie específica, 𝑆𝑣, que cuantifica en función de la porosidad local el valor de este parámetro, habiéndose encontrado un polinomio de quinto grado que parece ajustarse a las medidas experimentales. 𝑆𝑣=32,3𝑛−93,9𝑛2+134𝑛3−101𝑛4+28,8𝑛5 Según todo lo expuesto y recordando que este modelo considera que la matriz ósea está totalmente mineralizada, la variación de densidad aparente se define según la ecuación 3.11, teniendo en cuenta que k es el porcentaje de superficie disponible que está active para la remodelación y que normalmente se considera igual a la unidad. 𝜌󰇗=𝑘 𝑟 󰇗𝑆𝑣 𝜌 Por lo tanto la ley de variación de la densidad aparente en el tiempo sería la siguiente: 𝜌(𝑡+∆𝑡)=𝜌(𝑡)+𝜌󰇗(𝑡)∆𝑡 Obteniendo esta densidad en cada instante se determina el módulo de elasticidad y el coeficiente de Poisson siguiendo las correlaciones experimentales que se muestran a continuación. 𝐸={2014𝜌2,5 𝑠𝑖 𝜌≤1,2 𝑔/𝑐𝑚3 1763𝜌3,2 𝑠𝑖 𝜌≤1,2 𝑔/𝑐𝑚3 𝑣={0,2 𝑠𝑖 𝜌≤1,2 𝑔/𝑐𝑚3 0,32 𝑠𝑖 𝜌≤1,2 𝑔/𝑐𝑚3 Finalmente para obtener el tensor de comportamiento usamos la fórmula habitual. C=𝐸 1+𝑣(𝑣 1−2𝑣𝟏⨂𝟏+I) 3.3 MODELO MECANOBIOLÓGICO ISÓTROPO DE REMODELACIÓN ÓSEA BASADO EN LA ACTIVIDAD DE LAS BMUS 3.3.1 Introducción Los modelos mecanobiológicos de remodelación ósea interna, a diferencia de los fenomenológicos, tratan de tener en cuenta los procesos biológicos que participan en este fenómeno cuantificando su influencia para mejorar los resultados y ser más fieles de cara a representar el comportamiento real del hueso. Respecto al modelo fenomenológico explicado se hará más hincapié en este por ser uno de los que se usará en el resto del proyecto y aquel se pretende mejorar e introducir en él modificaciones. (3.10) (3.11) (3.12) (3.13) (3.14) (3.15) Este modelo parte de un antecesor que tiene en cuenta las diferentes células óseas, el concepto de BMU y la secuencia de actividad ARF. Algunas de las innovaciones que introduce García-Anar et al. [30] respecto a otros anteriores es que considera que el microdaño cambia las propiedades elásticas del hueso, que no solo dependen de la porosidad. Además, analiza los casos de desuso prolongado teniendo en cuenta que no se produce una pérdida total de hueso como proponen otros modelos, si no que antes de eso se alcanza una situación de equilibrio. Como ya se ha mencionado este modelo incluye las Antes de entrar a detallar el modelo es importante hacer un repaso por los parámetros definidos anteriormente para cuantificar la composición de la matriz ósea y definir algunos nuevos necesarios para comprender el modelo. En una muestra de hueso el volumen total, 𝑉𝑡, se define como la suma de los distintos volúmenes que lo conforman: el de los poros, 𝑉𝑣, más la suma de los distintos componentes que forman la matriz ósea, 𝑉𝑏, que son un volumen de tejido dañado, 𝑉𝑑, el volumen de la componente orgánica 𝑉𝑜 y el volumen del componente inorgánico 𝑉𝑚. 𝑉𝑡=𝑉𝑏+𝑉𝑣=𝑉𝑚+𝑉𝑜+𝑉𝑑+𝑉𝑣 La fracción voumétrica de hueso, 𝑣𝑏, se define como: 𝑣𝑏=𝑉𝑏 𝑉𝑡 Para medir el grado de mineralización usamos la fracción de ceniza, 𝛼: 𝛼= 𝑚𝑚 𝑚𝑚+𝑚𝑜 Para medir el nivel de daño al que está sometido el hueso hay dos enfoques. El primero consiste en medir la degradación de las propiedades mecánicas del hueso, como la rigidez. En la ecuación 3.19 se ve como se define la variable degradación, d según esta interpretación: 𝑑=1−𝐸 𝐸0 Siendo 𝐸 el módulo de Young y 𝐸𝑜 el módulo de Young del material virgen [31, 32, 33,34]. El segundo enfoque para la medida del daño consiste en cuantificar la densidad de microgrietas dentro del material definiendo una variable, ℎ, densidad de daño, como la fracción de volumen dañado por volumen total. ℎ=𝑉𝑑 𝑉𝑡 La relación entre ℎ y 𝑑 depende de muchos factores pero estudios experimentales han observado cierta linealidad en la cantidad de microgrietas en un hueso y la pérdida de propiedades mecánicas por lo que en este modelo se toma esta relación como lineal: ℎ=𝑘 𝑑 (3.16) (3.17) (3.18) (3.19) (3.20) (3.21) Modelos de remodelación ósea interna 34 34 donde 𝑘 es una constante a determinar experimentalmente. Como variable independiente hay dos opciones: usar la fracción volumétrica de hueso 𝑣𝑏, o bien la densidad aparente, relacionándose ambas de la siguiente manera: 𝜌=𝜌𝑡(𝑣𝑏−ℎ) donde 𝜌𝑡 es la densidad del tejido óseo y se relaciona con el grado de mineralización. Hernandez et al. [22, 35] hallaron la siguiente relación experimental: 𝜌𝑡(𝑔 𝑐𝑚3)=1,41+1,29𝛼 Los resultados son mejores cuando se usa la 𝑣𝑏 como variable independiente. La correlación entre 𝐸 y 𝑣𝑏 usada inicialmente era: 𝐸(𝑀𝑃𝑎)=84370𝑣𝑏2,58𝛼2,74 Posteriormente se ha introducido la influencia del daño usando la correlación expuesta al final del capítulo 2. 𝐸(𝑀𝑃𝑎)=84370𝑣𝑏2,58𝛼2,74(1−𝑑) El coeficiente de Poisson se considera constante pues su variación es poco importante (𝑣 = 0,3). En el resto de este capítulo se hará una breve descripción matemática sobre términos básicos del modelo y a continuación se explicará como afronta cada fase del proceso este modelo, que hipótesis plantea, cuales son sus fundamentos y las variables y expresiones matemáticas que rigen el fenómeno de remodelación ósea. Se hará además, un balance de masa ósea para explicar como varía la densidad aparente del hueso bajo ciertos estímulos. Partiendo de la expresión (3.22) y derivando con respecto al tiempo de forma que se obtenga: 𝜕𝜌 𝜕𝑡=𝜕𝜌𝑡 𝜕𝑡(𝑣𝑏−ℎ)+𝜌𝑡𝜕𝑣𝑏 𝜕𝑡−𝜌𝑡𝜕ℎ 𝜕𝑡=Π𝑏 Esta expresión sirve para determinar como evoluciona la masa ósea. El primer término describe cómo se mineraliza el osteoide, el segundo representa la evolución de la masa ósea por reabsorción y formación por parte de las BMUs y el tercero representa el efecto que sobre la masa ósea tiene la evolución del daño. 3.3.2 Descripción matemática de la actividad de las BMUs Como ya se ha explicado anteriormente las BMUs, están formadas por osteoblastos y osteoclastos que se encargan de formar y reabsorber hueso respectivamente a unas velocidades 𝑣󰇗𝐹 y 𝑣󰇗𝑅 . De tal forma la velocidad a la que cambia la fracción volumétrica se rige por la siguiente ecuación: 𝑣󰇗𝑏=𝑣󰇗𝐹−𝑣󰇗𝑅 En una situación de equilibrio la velocidad de cambio de la fracción volumétrica es igual a cero porque (3.22) (3.23) (3.24) (3.25) (3.26) (3.27) 𝑣󰇗𝐹 es igual a 𝑣󰇗𝑅 pero eso no quiere decir que no esté teniendo lugar la remodelación ósea, el tejido está en continua renovación. El valor de estas velocidades va a dependía en modelos previos del número de BMUs que se encontraran activas. Sin embargo, tras la discretización no es necesaria esta consideración puesto que será una la BMU que actúe en cada momento y no hay que cuantificarla. Una vez activada la BMUs actúa a una velocidad 𝑣𝐵𝑀𝑈, recorren una longitud 𝐿𝐵𝑀𝑈 antes de morir de forma que 𝐿𝐵𝑀𝑈 = 𝜎𝐿 𝑣𝐵𝑀𝑈 siendo 𝜎𝐿 la vida de la BMU. Las velocidades de reabsorción y de formación pueden expresarse de coincidirán con 𝐴󰇗(𝑥). Siendo 𝑥 la dirección de avance de la BMU y 𝐴󰇗(𝑥) la velocidad con la que varía el volumen por unidad de longitud en 𝑥 debido a la acción de las BMUs. Esta variación de la sección del hueso toma valores distintos para la zona de reabsorción y de formación Debido a que el número de osteoblastos y osteoclastos que actúan no es constante y va variando en el tiempo el valor de 𝐴󰇗𝐹(𝑥)y el de 𝐴󰇗𝑅(𝑥) tampoco lo será y en función de la fase en la que se esté tendrán unos valores u otros. Inicialmente se asumían valores constantes de las variaciones de secciones de la siguiente forma: 𝐴󰇗(𝑥)= { 𝐴𝐵𝑀𝑈 𝑇𝑅𝑓𝑐 en la zona de reabsorción 0 en la zona de inversión 𝐴𝐵𝑀𝑈 𝑇𝐹𝑓𝑏 en la zona de formacion 𝑇𝑅 y 𝑇𝐹 son los tiempos que duran, en una sección las fases de reabsorción y formación respectivamente. 𝑓𝑏 y 𝑓𝑐 son variables cuyo valor oscila entre 0 y 1 y representan el nivel de actividad de los osteoblastos y los osteoclastos. 𝐴𝐵𝑀𝑈 es el área de la sección de una BMU, distan para el hueso trabecular y para el cortical. 𝐴𝐵𝑀𝑈= { (𝑑𝑜2 4−𝑑ℎ 2 4)𝜋 𝜋4𝑑𝑒𝑑𝐵𝑀𝑈 (3.28) (3.29) Modelos de remodelación ósea interna 36 36 Figura 3.6 Forma simplificada del volumen remodelado por una BMU, en hueso trabecular (izquierda) y en hueso cortical (derecha). Fueron Martínez-Reina et al. [37] quienes obtuvieron unas expresiones de 𝐴󰇗𝑅(𝑡′′) y 𝐴󰇗𝐹(𝑡′′) para el hueso cortical más fieles a la realidad y que inicialmente se aplicaron al modelo anisótropo pero que se usan en este proyecto por ser perfectamente válidas para el modelo isótropo. Estas nuevas expresiones tienen en cuenta todas las fases durante la reabsorción y formación, sin olvidar que no todos los osteoblastos ni los osteoclastos trabajan simultáneamente. Figura 3.7 Fases de la actividad de una BMU en el hueso cortical. 𝐴󰇗𝑅(𝑥)= { 𝜋𝑑𝑜2 4(𝑡′′ 𝑇𝑅)2 Fase 1 𝜋𝑑𝑜2 4 Fase 2-4 𝜋𝑑𝑜2 4[1−𝑡′′−𝜎𝐿 𝑇𝑅]2 Fase 5 𝐴󰇗𝐹(𝑥)= { 𝜋𝑑𝑜2 4(1−𝑓12(𝑡′′)) Fase 3 𝜋(𝑑𝑜2−𝑑𝐻 2) 4 Fase 4-6 𝜋𝑑𝑜2 4(𝑓12(𝑡′′)−𝑑𝐻 2 𝑑𝑜2) Fase 7 donde las funciones 𝑓1 y 𝑓2 vienen dadas por las siguientes expresiones: 𝑓1(𝑡′′)=1−(1−𝑑𝐻 𝑑𝑜)·𝑡′′−𝑇𝐼−𝑇𝑅 𝑇𝐹 (3.30) (3.31) (3.32) (3.33) Modelos de remodelación ósea interna 38 38 𝑓2(𝑡′′)=1−(1−𝑑𝐻 𝑑𝑜)·𝑡′′−𝑇𝐼−𝑇𝑅−𝜎𝐿 𝑇𝐹 𝑇𝐼 representa la duración de la fase de inversión. Combinando las distintas ecuaciones obtendremos el valor de los volúmenes de reabsorción y de formación en cada momento de la vida de la BMU según su edad que hace variar 𝐴𝐹 y 𝐴𝑅 : 𝑣𝐹=𝐴𝐹 𝑣𝐵𝑀𝑈 𝑓𝑏 𝑣𝑅=(1−𝑝) 𝐴𝑅 𝑣𝐵𝑀𝑈𝑓 𝑐 Donde el factor (1−𝑝) se introduce para tener en cuenta el hecho de que las BMUs reabsorben un tejido poroso con una porosidad de valor 𝑝. Por otra parte, otro factor que controla la cantidad de hueso remodelado es el balance local de hueso (focal bone balance) [22]. Esta variable registra el balance de actividad entre osteoclastos y osteoblastos, y se define: 𝑓𝑏𝑏 = (𝑓𝑏/𝑓𝑐). Es función del estímulo mecánico 𝜉 [38, 39, 40, 43], el cual ya se ha mencionado en apartados anteriores, y de la edad [2, 44]. Entre 𝑓𝑏𝑏 y 𝜉 se establece una relación que se puede observar en la figura 3.8, similar a la que ya se vió cuando se definió el estímulo mecánico en el modelo isótropo de Stanford. La dependencia con la edad resulta difícil de establecer. Según García-Aznar [30] el hueso de una persona adulta tiene un balance local de hueso neutro, es decir, 𝑓𝑏𝑏,𝑚𝑎𝑥 = 𝑓𝑏𝑏,𝑚𝑖𝑛 =1,0 , mientras que en un hueso en crecimiento o en un hueso viejo 𝑓𝑏𝑏,𝑚𝑎𝑥 >1 y 𝑓𝑏𝑏,𝑚𝑖𝑛 <1 pudiendo producirse cambios significativos en la masa ósea. Figura 3.8 Balance local de hueso en función del estímulo mecánico. (3.34a) (3.34a) El estímulo mecánico que influye en el balance local de hueso se define como un escalar que depende del nivel de deformaciones y del número de ciclos que corresponden a cada caso de carga 𝑖. Esta propuesta hecha por Mikic y Carter, [41], es análoga a la del modelo isótropo de Stanford, salvo que en esta se relaciona directamente con las deformaciones en lugar de con el nivel de tensiones. En la ecuación 3.35 se plasma la definición del estímulo mecánico de este modelo, teniendo en cuenta que 𝑛𝑖 es el número de ciclos que actúa el caso de carga 𝑖, 𝜀𝑖 la deformación efectiva correspondiente a ese caso de carga y 𝑚 el exponente que mide la importancia del número de ciclos y del nivel de deformaciones y cuyo valor típico se ha determinado experimentalmente y es 𝑚 = 4. 𝜉=(∑𝑛𝑖 𝜀𝑖𝑚 𝑖)1/𝑚 La deformación efectiva para cada caso de carga, 𝜀𝑖, se define en función de la densidad de energía de deformación que provoca el caso de carga correspondiente 𝑈𝑖 y del módulo de Young del tejido óseo, 𝐸. 𝜀𝑖=√2𝑈𝑖 𝐸 En la figura 3.8, se observa que en una zona entorno al estímulo de equilibro el balance local de hueso toma valor unitario. Esto se debe a que para ese estímulo la actividad de las BMUs no tiene efecto neto sobre el volumen de hueso a largo plazo. Un valor típico del estímulo de equilibrio es 𝜉∗= 2,5 x 10−3𝑐𝑖𝑐𝑙𝑜𝑠1 4𝑚𝑚/𝑚𝑚. Éste es el valor del estímulo en el hueso cortical que soporta una tensión de 50 [25] a razón de 10000 ciclos diarios. Sin embargo, este valor parece ser aplicable sólo a huesos cuya misión sea la de soportar carga, como el fémur. Huesos menos cargados, como el cráneo y otros huesos planos, cuya misión es principalmente de protección, no soportan valores tan altos. La diferencia de comportamiento entre huesos largos y huesos planos se obtiene en el modelo considerando que el valor de 𝜉∗ no es una constante característica de un hueso, sino que varía con el tiempo, adaptándose al valor del estímulo externo aplicado, según el principio de acomodación celular propuesto por Turner [42]. Éste sostiene que el tejido se acomoda a la deformación que soporta. De esta forma la deformación de referencia tiende hacia la deformación normal de uso: 𝑑𝜀𝑜 𝑑𝑡=𝜑(𝜀−𝜀𝑜) Adaptando la expresión al estímulo mecánico para este modelo, la ecuación queda totalmente análoga: 𝑑𝜉∗ 𝑑𝑡=𝜑(𝜉−𝜉∗) (3.35) (3.36) (3.37) (3.38) Modelos de remodelación ósea interna 40 40 𝜑 es el parámetro que controla la velocidad de adaptación y que toma en el modelo un valor de 0.002días-1. Con esta teoría un estímulo muy bajo, como el que sufren los huesos planos, no produce reabsorción indefinida, sólo hasta que el estímulo de equilibrio se adapte a esa carga tan baja que produce un estímulo exterior pequeño. Así se explica el hecho de que en huesos no muy cargados, como el cráneo, no se pierda masa ósea como consecuencia de un desuso prolongado. 3.3.3 Crecimiento el daño Las cargas diarias a las que están sometidos los huesos provocan que se acumule un daño que acaba dando lugar al crecimiento de grietas. En este proceso no hay aportación de hueso, Π𝑏, ni varía la densidad del tejido, 𝜕𝜌 𝜕𝑡=0, con lo que del balance de masa de la ecuación (3.26) se deduce que: 𝜕ℎ 𝜕𝑡=𝜕𝑣𝑏 𝜕𝑡 Esta ecuación establece que el aumento del daño por acumulación de microgrietas desemboca en un aumento de la cantidad de tejido. Para comprender esto es importante recordar que el hueso dañado forma parte del tejido óseo y no de los huecos. La ecuación (3.40) indica que al abrirse las grietas, el volumen asociado expande el tejido completo, disminuyendo la proporción de volumen ocupada por los poros. Para saber como se relacionan el daño acumulado en el hueso con las cargas a las que éste está sometido hay que estudiar el comportamiento a fatiga del hueso. Este ha sido estudiado de manera exhaustiva por diversos autores [32, 33, 31, 45, 46]. Como se explicó antes, una de las maneras de medir el daño acumulado por fatiga es a través de la degradación del módulo de Young con la variable 𝑑. Análogamente a la ecuación 3.21 la variación de esta variable y la del volumen dañado cumplen la siguiente relación lineal: 𝜕ℎ 𝜕𝑡=𝑘𝜕𝑑 𝜕𝑡 Lemaitre y Chaboche expresaron el crecimiento de daño por fatiga de forma genérica [47]. Usaron expresiones del tipo 𝜕𝑑 𝜕𝑡=𝐺(𝑑,𝜀) donde 𝑛 es el número de ciclos. Por otro lado, la vida a fatiga, 𝑛𝐹 , se suele aproximar por expresiones del tipo: (3.39) (3.40) (3.41) 𝑛𝐹=𝐶𝜀𝛿 La evolución del daño y la expresión de la vida a fatiga siguen leyes distintas en función si las cargas a las que el hueso está sometido son de tracción o de compresión. García-Aznar et al. [30] proponen las siguientes leyes de evolución del daño para cada caso: 𝜕𝑑𝑐 𝜕𝑛=𝐺𝑐(𝑑𝑐,𝜀)=𝐶1 𝛾1𝑒𝛾1𝑑𝑐𝜀𝛿1 en compresión 𝜕𝑑𝑡 𝜕𝑛=𝐺𝑡(𝑑𝑡,𝜀)=𝐶2 𝐶3𝛾2(1−𝑑𝑡)1−𝛾2𝑒−𝐶3(1−𝑑𝑡)𝛾2𝜀𝛿2 en tracción donde 𝑑𝑐 y 𝑑𝑡 son los parámetros de degradación del módulo de Young en compresión y tracción, respectivamente; 𝑁 es el número de ciclos y 𝐶1,𝐶2,𝐶3, 𝛾1, 𝛾2 , 𝛿1 y 𝛿2 son parámetros del modelo. La forma particular de las funciones anteriores fue elegida por García-Aznar et al. para ajustar los resultados experimentales obtenidos por Pattin et al. [33]. La integración de las ecuaciones (3.43) conduce a las siguientes expresiones: 𝑑𝑐=−1 𝛾1[ln(1−𝐶1𝜀𝛿1𝑛)] 𝑑𝑡=−1− √1 𝐶3ln(𝑒𝐶3−𝐶2𝜀𝛿2𝑛) 𝛾2 que tienen la forma representada en la figura 3.9. Figura 3.9 Evolución del daño en el tejido óseo, en tracción (izquierda) y compresión (derecha). (3.42) (3.43a) (3.43b) (3.44b) (3.44a) Modelos de remodelación ósea interna 48 48 Martinez-Reina, [55], propuso una nueva formulación para simular el proceso de mineralización. A continuación, se va exponer ésta empezando por la ecuación que representa la cantidad de mineral en un instante 𝑡+𝛥𝑡, la cual en términos de 𝑣𝑚 en lugar de 𝛼 quedaría de la siguiente forma: 𝑣𝑚(𝑡+𝛥𝑡)=𝑣𝑚(𝑡)𝑉𝑏(𝑡)−𝑣𝑚(𝑡)𝛥𝑉𝑟(𝑡)+𝛥𝑉𝑚𝑖𝑛(𝑡+𝛥𝑡) 𝛥𝑉𝑏(𝑡+𝛥𝑡) El primer término representa el volumen de mineral presente en el instante 𝑡. Con el segundo término se resta el mineral que es eliminado en la matriz ósea durante la fase de reabsorción, siendo 𝑉𝑟(𝑡) el tejido eliminado desde t hasta t+Δt. El último término representa el volumen de mineral que va apareciendo debido a la mineralización. 𝛥𝑉𝑚𝑖𝑛(𝑡+𝛥𝑡) viene dado por la siguiente expresión: 𝛥𝑉𝑚𝑖𝑛(𝑡+𝛥𝑡)=∫ ∇𝑓(𝑡+𝛥𝑡)𝑑𝑣𝑚(𝜏) 𝑑𝜏 = ∞ 0∫ ∇𝑓(𝑡,𝜏)𝑑𝑣𝑚(𝜏) 𝑑𝜏 ∞ 0 donde ∇𝑓(𝑡,𝜏) es el volumen de tejido que se formó hace 𝜏 días y que está presente en el instante 𝑡. En la integral se debe tener en cuenta la mineralización de todo el tejido formado. Sin embargo y puesto que el hueso está en continua regeneración, aquel tejido que formó hace mucho tiempo ya debe haber sido reabsorbido, por lo que no será mineralizado y es necesario sacarlo del integrando (∇𝑓(𝑡,𝜏)→ 0 cuando 𝜏→∞). Por este motivo se pueden reducir los límites de integración de entre 0 y ∞ a entre 0 y 𝑡𝑅, tiempo de residencia del tejido en la matriz ósea. Para calcular el valor de esta integral se aproximó por la adición de incrementos diarios de la siguiente forma: 𝛥𝑉𝑚𝑖𝑛(𝑡+𝛥𝑡)≅∑∇𝑓(𝑡,𝑖)[𝑣𝑚(𝑖)−𝑣𝑚(𝑖−1)] 𝑡𝑅 𝑖=0 Para calcular o al menos estimar ∇𝑓(𝑡,𝑖) se supone que el tejido formado en los días previos se encontraba uniformemente distribuido en el volumen de hueso. Así el volumen reabsorbido de un cierto ∇𝑓(𝑡,𝑖) es proporcional a ∇𝑓(𝑡,𝑖)/𝑉𝑏(𝑡) y claro está, al volumen reabsorbido 𝛥𝑉𝑟 (𝑡) . Se define entonces un algoritmo recursivo para calcular el volumen de tejido dañado, ∇𝑓(𝑡+1,𝑖+1) que observamos en la ecuación 3.61 y representado en la figura 3.12. Respecto al mencionado tiempo de residencia, tR, se considera variable y depende aquí de 𝛥𝑉𝑟 y 𝛥𝑉𝑓 , que en un estado de equilibrio de remodelación deberán ser iguales. Para disminuir el coste computacional de las simulaciones la cola de la figura 3.12 se trunca en el día 𝑖 = 𝑡𝐶(𝑡𝐶 < 𝑡𝑅). El volumen de tejido formado hace 𝑡𝐶 se supone que se mineraliza instantáneamente hasta el máximo contenido mineral. Esto es aproximado, pero si 𝑡𝐶 es demasiado largo, la cantidad de tejido envuelta (3.56) (3.57) (3.55) es tan pequeña que el efecto de la aproximación es despreciable. ∇𝑓(𝑡+1,𝑖+1)=∇𝑓(𝑡,𝑖)−𝛥𝑉𝑟(𝑡)∇𝑓(𝑡,𝑖) 𝑉𝑏(𝑡) Figura 3.12 Algoritmo empleado para calcular ∇𝑓(𝑡+1,𝑖+1). Se usa una cola (FIFO) con la ley de recurrencia expresada en la ecuación 3.61. 𝑉𝑟 y 𝑉𝑓 son los volúmenes de tejido reabsorbido y formado el día 𝑡, respectivamente. 3.3.7 Degradación de las propiedades a fatiga con el contenido de mineral El contenido de mineral aporta rigidez al hueso pero también ocasiona que este se vuelva más frágil. Conforme el hueso se mineraliza, la deformación de rotura disminuye como observamos en la figura 3.13. Las referencias de esta sección pueden consultarse en el artículo de Martínez-Reina, [55]. Los tests de fatiga predicen la evolución del daño según el nivel de deformaciones o tensiones y según el número de ciclos. Esta relación viene dada con una expresión del tipo (teniendo en cuenta que 𝑖 puede ser t o c en función de si hablamos de compresión o tracción): 𝐾𝑖 𝜀𝛿𝑖 (3.58) (3.59) Modelos de remodelación ósea interna 50 50 Figura 3.13 Relación entre la tensíon última de rotura y la concentración de calcio (masa de calcio por masa seca). Se establecen una serie de hipótesis para poder representar este fenómeno por el cual las propiedades a fatiga disminuyen con el contenido mineral y a causa de la carencia de estudios sobre este tema. • La forma de las curvas 𝑑−𝑁, expuestas anteriormente tanto para los casos de compresión como de tracción se mantienen sin incluir el contenido de mineral. • Sólo la vida a fatiga se ve afectada por el contenido de mineral realizando un cambio en la constante 𝐾𝑡 que hemos visto en las ecuaciones anteriores. El valor de 𝛿𝑡 no se cambia. Estos cambios modifican la curva 𝑑−𝑁 puesto que 𝐶𝑡1 depende de 𝐾𝑡. De esta forma incrementando 𝐾𝑡 se obtiene una vida a fatiga mayor. • Se ha establecido una vida de 107 ciclos para el límite de fatiga. Este límite de fatiga se da respecto a una fracción de la deformación del límite de rotura que depende del material mediante una constante 𝛽, 𝜀𝑢/𝛽. Así y teniendo en cuenta la ecuación 3.62 podemos despejar el valor de 𝐾𝑡: 𝐾𝑡([𝐶𝑎])=107 𝜀𝑢 ([𝐶𝑎])𝛿𝑡 𝛽 • En base a los resultados experimentales que se muestran en la figura 3.13 la función por la que se obtiene 𝜀𝑢=𝜀𝑢 ([𝐶𝑎]) es la siguiente: log𝜀𝑢=25,452−11,341log[𝐶𝑎] (3.61) (3.60) 3.4 MODELO MECANOBIOLÓGICO ANISÓTROPO DE REMODELACIÓN ÓSEA CON TRATAMIENTO DIFUSIVO DE LAS BMUS 3.4.1 Introducción El modelo que se explica a continuación surge tras la propuesta hecha por Martinez Reina [3] en su tesis, al desarrollar un modelo en el que el comportamiento del hueso fuera anisótropo y tras las aportaciones que hizo I. Reina [12] a este modelo considerando el movimiento de las BMUs a lo largo de su vida. El modelo inicial se trata de una continuación del modelo isótropo de BMUs propuesto por García- Aznar et al. [30], el cual es el antecesor del modelo de remodelación ósea discreto de BMUs que se ha explicado en el apartado anterior. Este modelo predecía los cambios del hueso teniendo en cuenta aspectos biológicos significativos. Sin embargo, no contemplaba la direccionalidad en el comportamiento de este. Las bases del modelo que sí las contempla se explicarán a continuación. Este modelo anisótropo citado, sin embargo, no incluía el movimieto de las BMUs, las cuales una vez activadas permanecían estáticas en el punto de activación. Tras las aportaciones de I. Reina [12] esto se subsana, considerando que la BMU avanza a lo largo del transcurso de su vida, en una dirección concreta y repartiendo un efecto remodelatorio en su entorno. Así se generan flujos de BMUs que pueden remodelar más allá de su entorno de origen. De esta manera surge un modelo mecanobiologico que considera tanto la direccionalidad del hueso como el movimiento de las BMUs a trevés del hueso. Cabe destacar que en todo lo referido a mineralización del hueso y a la degradación que sufren las propiedades de fatiga debido a la cantidad de calcio en el hueso son totalmente análogas a las explicadas en los apartados 3.3.6 y 3.3.7, ya que el planteamiento es el mismo en ambos modelos. Por ello, y para no caer en la redundacia no se ha repetido la misma información en esta sección. 3.4.2 Variables básicas del modelo En este apartado se repasan las variables generale básicas del modelo. Gran parte de ellas son simiares a las vistas en el modelo isótropo, pero realizando su extensión al caso anisótropo convirtiéndolas en tensores. En el resto de subsecciones se verán más variables del modelo asociadas a áreas más concretas de este, como el crecimiento del daño y el tratamiento difusivo de las BMUs. Porosidad y tensor MIL Para cuantificar la anisotropía del entramado óseo se definen nuevas variables que representen estas modificaciones. Modelos de remodelación ósea interna 52 52 Al igua que 𝑣𝑏 establece la cantidad de tejido que hay en un punto, el tensor 𝐋, describe como se reparte ese tejido en el espacio, proporcionando así una medida de la anisotropía e influyendo en las propiedades del hueso. El tensor MIL, L, se define como aquél cuya proyección en una determinada dirección, representa la longitud del poro promedio contenido en la unidad de volumen medida en esa dirección. Se supone que el poro promedio contenido en la unidad de volumen tiene forma ortoédrica. Figura 3.14 Esquema de una estructura anisótropa de hueso trabecular y una línea de muestreo para la evaluación de la MIL. Tomada de Martin et al. [2]. La porosidad y el tensor MIL se relacionan mediante la siguiente expresión: 𝑝=det (L) Esto se debe a que debido a la hipótesis de poro ortoédrico los autovectores del tensor L serás las normales a las caras de este y los autovalores las dimensiones. El tensor M es un tensor paralelo al tensor L con autovalores proporcionales, pero de traza unitaria. M=1 𝑡𝑟(𝐋)L Este tensor solo contiene información del grado de anisotropía del entramado óseo, que afecta directamente a las propiedades elásticas del hueso. Mientras que el tensor L se ve afectado directamente por la actividad de las BMUs, ya que la formación y reabsorción son fenómenos que modifican la longitud de los poros y por tanto las componentes del tensor L. De esta manera es necesario el uso simultáneo de los dos tensores, L y M. (3.63) (3.62) La velocidad con la que varía la porosidad se obtiene a partir de la velocidad a la que varía el tensor L mediante la siguiente expresión: 𝑝󰇗=𝑝 𝑡𝑟(𝐋−1𝐋) Número de BMUs A diferencia del modelo isótropo expicado previamente, en este modelo las BMUs no se activan de manera discreta, si no que fracciones de BMUs se activan de manera continua en el hueso. Por lo que es necesario una variable que cuantifique el número de BMUs activas que avanzan en una dirección determinada por unidad de volumen. Puesto que además en función de la dirección de avance de estas el entramado óseo variará las características de su anisotropía de una manera u otra. Esta variable será el tensor 𝐍, que se define como aquel cuya proyección en una determinada dirección, 𝒏𝑇𝐍 𝒏 proporciona en número de BMUs activas que avanzan endicha dirección n por unidad de volumen. De manera análoga el tensor 𝐍󰇗 se define como aquel cuya proyección en una determinada dirección 𝒏𝑇𝐍󰇗 𝒏 proporcióna el número de BMUs activas por unidad de tiempo y de volumen que avanzan en dicha dirección. Tensor estímulo y dirección de avance Cuando un hueso está sometido a varios estados de carga durante la actividad normal del individuo, las direcciones principales de deformación irán cambiando con cada actividad. Hay que tener en cuenta que el tiempo medio de variación de las deformaciones durante la actividad normal es mucho menor que el tiempo que tardan en activarse las BMUs. Por tanto, dichas BMUs no pueden seguir las direcciones principales de deformación de manera instantánea. Se supondrá en ese caso que el hueso recopila información de todo lo que ocurre diariamente y las BMUs se activan en función de esa historia diaria más reciente. Esta hipótesis requiere la definición de una deformación acumulada, análoga al estímulo mecánico defnido por Carter et al. [26]. Este estímulo mecánico ya se usó en el modelo isótropo de BMUs y tiene en cuenta las distintas actividades que soporta el hueso, así como el número de ciclos de cada actividad. El tensor de deformación acumulada se define como la superposición ponderada de los diferentes estados de deformación que producen las actividades diarias a las que se encuentra sometido el hueso: 𝜀=(∑𝑛𝑗 𝜀𝑗𝑚 𝑁 𝑗=1 )1/𝑚=(∑∑(𝑛𝑗 𝜀𝑖𝑗 𝑚 3 𝑖=1 ) e𝑖𝑗(𝑡)⊗e𝑖𝑗(𝑡) 𝑁 𝑗=1 )1/𝑚 (3.64) (3.64) Modelos de remodelación ósea interna 54 54 Siendo 𝑛𝑗 el número de ciclos correspondientes al caso de carga 𝑗, N el número de actividades distintas que soporta el hueso, 𝑚 el exponente introducido por Carter et al. [26], 𝜀𝑖𝑗 la deformación principal 𝑖 correspondiente alcaso de carga 𝑗 y e𝑖𝑗 la dirección principal de deformación 𝑖 correspondiente al caso de carga 𝑗. Por otra parte, se define el tensor estímulo mecánico como un tensor paralelo al tensor de deformación acumulada cuyos valores representan una energía de deformación asociada a la dirección principal correspondiente de este último tensor. 𝝃=∑𝜉𝑖e𝑖(𝑡)⊗e𝑖(𝑡) 3 𝑖=1 Estos autovalores, 𝜉𝑖, se expresan como: 𝜉𝑖=12𝜎𝑖𝑃𝜀𝑖=12(e𝑖𝑇𝝈 e𝑖)𝜀𝑖 Donde 𝜎𝑖𝑃 es la proyección el tensor de tensiones en la dirección principal de deformación acumulada e𝑖 y el autovalor 𝜉𝑖 es la energía de deformación asociada a dicha dirección, tomando como deformación la acumulada 𝜀𝑖. Pese a lo comentado, el estímulo mecánico y por tanto las deformaciones sufridas en el hueso no serán el único factor que determine la dirección a seguir por las BMUs, si no que será una combinación de estos y el daño del hueso. La BMU seguirá en cada instante una dirección 𝐞, que salde de combinar la dirección de estímulo mecánico promedio máximo, 𝐞𝑚𝑎𝑥, y la dirección de máximo gradiente del daño ∇𝐝. La influencia del daño dependerá de laintensidad de este y la distancia a la que se encuentre la BMU del foco dañado. La relación entre estos parámetros es la siguiente: 𝐞=𝐾 ∇𝐝 ||∇𝐝||+( 1−𝐾)𝐞𝑚𝑎𝑥 ||𝐞𝑚𝑎𝑥|| Siendo 𝐾 un factor de ponderación que varía entre 0 y 1 que permite cuantificar el nivel de importancia del daño en la trayectoria de la BMU. Si valiera 0 la dirección de avance vendría regida exclusivamente por el estímulo mecánico promedio máximo, y si valiera 1 evolucionará solo en función del gradiente del daño. Las BMUs son atraídas por los osteocitors apoptóticos presentes en la zona daña y esta apostosis está controlada por el equilibrio de dos moléculas: el Bax y el Bcl-2. La diferencia entre las concentraciones (3.65) (3.66) (3.67) de estas molécuas, que llamamos señal neta, determina si los osteoclastos se deben dirigir a la zona dañada. Por tanto, el factor de ponderación 𝐾 dependerá directamente de esta señal neta, 𝑆𝑛𝑒𝑡𝑎. Se puede comprobar, en función de estas dependencias, que la relación es la siguiente: Figura 3.15: Influencia del daño en función del valor de la señal neta [12]. Siendo 𝑑0 un umbral antes del cual la influencia del daño es nula y 𝑑𝑡ℎ un nivel que da paso a una zona de influencia máxima del daño. El tensor de daño considerado en este modelo, d, es análogo a la degradación, d, del modelo isótropo. La proyección del tensor de daño en una determinada dirección está relacionada con la pérdida de rigidez del material en dicha dirección. El tensor de daño está relacionado con el volumen de tejido dañado, h, como ya se discutió en el modelo isótropo, a través de una constante k. Ahora, sin embargo, se ha de tener en cuenta la degradación en todas las direcciones, por medio de la traza del tensor d. Así: ℎ=𝑘3𝑡𝑟( 𝐝) Frecuencia de activación Para cada dirección de avance de BMUs se define una frecuenci de activación cuyo cálculo se hace en dospartes. En primer lugar, se estima la frecuencia de activación total, o frecuencia de activación isótropa y en segundo lugar se repartendichas BMUs según las direcciones en las que progresan. La frecuencia de activación isótropa 𝑓𝑜𝑟 𝑖𝑠𝑜 es el número de BMUs que se activan por unidad de tiempo y por unidad de superficie disponible para la activación. 𝑓𝑜𝑟 𝑖𝑠𝑜=𝑓𝑏𝑖𝑜 ( 1− 𝜉𝑖𝑠𝑜 𝜉𝑖𝑠𝑜+𝑐(1−𝑑𝑖𝑠𝑜)𝑎 (3.67) (3.68) Modelos de remodelación ósea interna 56 56 donde 𝑐 y 𝑎 son constantes del modelo; 𝑓𝑏𝑖𝑜 es un parámetroo que tiene en cuenta diversos factores biológicos que influyen en la activación de BMUs como por ejemplo la edad; la variable 𝑑𝑖𝑠𝑜 es la componente octaédrica del tensor de daño (𝑑𝑖𝑠𝑜=𝑡𝑟(𝐝) 3) que es una medida del nivel de daño global y el estimulo 𝜉𝑖𝑠𝑜=𝑡𝑟(𝝃) 3 es la componente octaédrica del tensor de estímulo. Según la expresión anterior, la cantidad de BMUs que se activen dependerá de la energía de deformación total que esté soportando el hueso y del nivel de daño global, sin distinguir la dirección de este daño. Una vez que se conoce el número total de BMUs que se activan por unidad de tiempo y por unidad de superficie específica, se hace un reparto de las BMUs según su dirección de avance. Para ello se define primero una señal inhibitoria direccional y después el tensor frecuencia de activación. Para cada dirección principal de estímulo 𝑖, hay una señal inhibitoria, 𝑠𝑖, que depende del estímulo y del daño asociados a esa dirección. El tensor señal inhibitoria se define como un tensor paralelo al tensor estímulo y cuyos autovalores son: 𝑠𝑖(𝜉𝑖,𝑑e 𝑖)= 𝜉𝑖 𝜉𝑖+𝑐(1−𝑑e 𝑖)𝑎 donde 𝜉𝑖 es el autovalor del tensor de esímulo correspondiente a la dirección 𝑖 y 𝑑e 𝑖=e𝑖𝑇𝐝 e𝑖 es la proyección del tensor de daño en dicha dirección. Las constantes 𝑐 y 𝑎 son las mismas de la ecuación anterior. El tensor frecuencia de activación es un tensor también paralelo al tensor de estímulo: 𝐟𝑜𝑟=∑𝑓𝑜𝑟 𝑖(e𝑖𝑇𝝈e𝑖) 3 𝑖=1 cuyos autovalores dependen de la señal inhibioria en la dirección correspondiente, 𝑠𝑖, y de la frecuencia de activación isótropa, 𝑓𝑜𝑟 𝑖𝑠𝑜, 𝑓𝑜𝑟 𝑖= 𝑓𝑜𝑟 𝑖𝑠𝑜 1−𝑠𝑖 ∑ (1−𝑠𝑖) 3𝑖=1 3.4.3 Actividad de las osteonas. Variacion del tensor MIL en el hueso cortical En este modelo, se trabaja únicamente con el hueso cortical, puesto que las distintas fórmulas que se emplean para el área de paso de la BMU, (3.30) y (3.31), son únicamente válidas para huesos poco porosos. Además, la influencia de la superficie específica, 𝑆𝑣 es diferente para hueso trabecular y cortical, por lo que es necesario tratarlos de manera independiente. La variable 𝑓𝑜𝑟 𝑖, vista previamente proporciona el número de BMUs activas por unidad de tiempo, unidad de volumen y unidad de superficie disponible, por lo que: (3.69) (3.70) (3.71) 𝑁󰇗𝑖= 𝑓𝑜𝑟 𝑖 𝑆𝑣 En la siguiente figura se esquematiza el progreso de una BMU de hueso cortical. Las osteonas abren nuevos poros en la matriz ósea durante la fase de reabsorció y después esos poros se ven llenando durante la fase de formación. El avance de esta BMU porvoca la variación de a longitud media de los poros en las direcciones perpendiculares a la de avance, manteniendo así la isotropía transversal, característica del hueso cortical. Una BMU que avanza en una determinada dirección 𝑖 y está enfase de reaborción provocará una variación del tensor L tal que: 𝑑𝐋=𝑑𝑙( ê𝑗⊗ê𝑗+ê𝑘⊗ ê𝑘) Figura 3.16: Progreso de la BMU en hueso cortical [12] 𝑑𝐋 puede relacionarse con el cambio de porosidad que se experimente en el hueso de la siguiente manera: 𝑑𝑝=𝑝 𝑡𝑟 (𝐋−𝟏( ê𝑗⊗ê𝑗+ê𝑘⊗ ê𝑘))𝑑𝑙 Cada BMU, durante su vida, reabsorbe y forma un volumen igual a 𝑉𝐵𝑀𝑈=𝜋 4 𝐿𝐵𝑀𝑈(𝑑𝑂 2−𝑑𝐻 2). Como ya vimos, 𝐿𝐵𝑀𝑈= 𝜎𝐿 𝑣𝐵𝑀𝑈, siendo la longitud recorrida por la BMU. La variación de porosidad en la reabsorción y la formación debida a una cantidad concreta de BMUs activadas entre 𝑡′ y 𝑡′+𝑑𝑡′es: 𝑣󰇗𝑅=𝑝󰇗=𝑁󰇗𝑖(𝑡′)𝐴𝑅(𝑡′)𝑣𝐵𝑀𝑈 𝑓𝑐𝑖𝑑𝑡′ (3.72) (3.72) (3.73) (3.74a) Modelos de remodelación ósea interna 64 64 científicas resultado de distintas investigaciones muestran que las BMUs se mueven y realizan giros y trayectorias más complicadas en función de parámetros como el daño o el estímulo mecánico. Se ha demostrado además [56] que el microdaño que repara una BMU es mucho mayor que el que teóricamente podría reparar si se considera que avanza en línea recta. Por ello no considerar la difusión de las BMUs reduce la fiabilidad de los resultados de las simulaciones de estos modelos. Como ya se explicó en el apartado de las variables básicas del modelo,las BMUs siguen una trayectoria dada por una combinación etre el estímulo mecánico promedio máximo y el gradiente del daño a través de la siguiente expresión: 𝐞=𝐾 ∇𝐝 ||∇𝐝||+( 1−𝐾)𝐞𝑚𝑎𝑥 ||𝐞𝑚𝑎𝑥|| La idea cuando las BMUs, más bien el número de BMUs promediado diariamente, sigan esta trayectroria es que repartan su actividad remodelatoria dentro de un volumen de influencia definido por un radio de influencia, 𝑅𝑖𝑛𝑓. El efecto ded la remodelación dentro de la zona de influencia no tiene por qué ser el mismo en todos los puntos, si no que decrecerá conforme aumentela distanciaa centro de la zona, allá donde se encuentre la BMU. A la hora de considerar la trayectoria de las BMUs y de procesar las direcciones que va tomando en cada paso de iteración que tiene lugar en una simulación en la que se implemente el modelo, es muy importante tener en cuenta que hay tres puntos de dicha trayectoria que es importante ir almacenando y cuyo papel en este modelo es fundamental: frente guí, frente de reabsorción y frente de formació. El frente guía es el punto más adelantado dela BMU que va marcado el camino a seguir según el daño y el estímulo mecánico., como hemos visto antes. El frente de reabsorción va siguiendo al frente guía repartiendo su actividad de remover hueso conforme avanza. El frente de formación va más atrasado que el de formación y va cubriendo con hueso nuevo las zonas reabsorbidas. En definitiva, el direccionamiento de la BMU se produce en el frente guía y los otros le seguirán con un retraso varible. La implementación informática del frente de reabsorción considerandolo con forma cónica, que es la forma de la cavidad de reabsorción conforme se avanza, tiene un alto coste computacional por lo que se simplifica a considerar la posición del frente de reabsoricón como el centro geométrico de dicho cilindro, que irá variando conforme aumenta el tamaño de la cavidad. (3.94) Figura 3.19: Esquema con las diferentes etapas del proceso de reabsorción [12] De igual modo la evolución del frente de formación vendrá dada por su centro geométrico y por la velocidad de formación en cada instante. Figura 3.20: Esquema con las diferentes etapas del proceso de formación [12] Para poder llevar a cabo ese traspaso de posiciones de frente guía a frente de reabsorción y a frente de formación es necesario almecenar una parte de la trayectoria seguida por el frente guía en un vector de posiciones, para cada conjunto de BMUs activas, durante todos los días de activaión. Lo cual, para una sección de hueso simulada en un programa de elementos finitos que emplee una malla moderadamente fina y con un numero adecuado de puntos de integración se traduce un coste computacional muy elevado. Sin embargo, esto es necesario para poder ir transmitiendo la información al frente de formación y de Modelos de remodelación ósea interna 66 66 reabsorción de qué posiciones deben alcanzar para repartir su efecto remodelatorio. Figura 3.21: Distancia recorrida por los frentes [12] El vector de posiciones que recoge todas las últimas tomadas por el frente guía se llama SENDA. Cada grupo de BMUs tendrá asignada una SENDA, por lo que como decíamos antes es necesario almacenar gran cantidad de información. Por ello, el tamaño de SENDA se ajusta de manera que recoja lo estrictamente necesario. En la figura siguiente, se observa que la distancia máxima de posiciones a almacenar será (𝑇𝑅+𝑇𝐼+𝑇𝐹/3) 𝑣 𝐵𝑀𝑈 : Figura 3.22: Esquema con distancias entre los frentes de formación y reabsorción La rutina que actualiza la variable SENDA se ejecuta cada 𝑇𝑆𝐸𝑁𝐷𝐴 días. Como la velocidad de avance de la BMU se considera constante, la distancia existente entre posiciones consecutivas del vector es 𝑇𝑆𝐸𝑁𝐷𝐴𝑣𝐵𝑀𝑈, por ello, el número de posiciones que compone SENDA es: 𝑁º 𝑃𝑜𝑠𝑖𝑐𝑖𝑜𝑛𝑒𝑠=(𝑇𝑅+𝑇𝐼+𝑇𝐹 3)𝑣 𝐵𝑀𝑈 𝑇𝑆𝐸𝑁𝐷𝐴𝑣𝐵𝑀𝑈 =(𝑇𝑅+𝑇𝐼+𝑇𝐹 3) 𝑇𝑆𝐸𝑁𝐷𝐴 Para reducir el coste computacional 𝑇𝑆𝐸𝑁𝐷𝐴tiende a ser mayor que uno, y se actualiza SENDA cada más tiempo. Cuando se actualice SENDA se determinan las posiciones de los frentes: el frente guía es la primera posición y el de reabsorción y formación dependen de cuánto hace que se activó el grupo de BMUs y en que fase se encuentra. Reparto de la activiad de las BMUs El efecto de las BMUs según este modelo se reparte entre los puntos de Gauss de un entorno cercano. Este reparto no es igualitario, y para llevarlo a cabo se definen unos pesos en función de a distancia al foco. 𝑃𝑖∗=(𝑅𝑖𝑛𝑓−𝑑𝑖𝑠𝑡𝑖) 𝑅𝑖𝑛𝑓 Siendo 𝑅𝑖𝑛𝑓 el radio de influencia de la actividad de las BMUs y 𝑑𝑖𝑠𝑡𝑖 la distancia desde el frente hasta el punto de Gauss 𝑖. Si se normaliza esta función para que el total valga 1 la expresión queda así: 𝑃𝑖=𝑃𝑖∗ ∑𝑃𝑖∗ 𝑁𝑃𝐺 𝑖=1 =(𝑅𝑖𝑛𝑓−𝑑𝑖𝑠𝑡𝑖) 𝑁𝑃𝐺 𝑅𝑖𝑛𝑓−∑𝑑𝑖𝑠𝑡𝑖 𝑁𝑃𝐺 𝑖=1 Donde 𝑁𝑃𝐺 es el número de puntos de Gauss que hay dentro de la zona de influencia en la que se realiza el reparto. Según esto, en cada punto de Gauss habrá un efecto remodelatorio, pudiendo ser tanto de reabsorción como de formación en función del conjunto de BMUs que esté repartiendo su efecto ahí y cuando se haya originado este. En la siguiente imagen se muestra un esquema del algoritmo realizado por I. Reina [12] en su proyecto: (3.95) (3.96) (3.97) Modelos de remodelación ósea interna 68 68 Figura 3.23: Esquema del algoritmo [12] 3.4.7 Propiedades elásticas del hueso Como ya se ha comentado, las propiedades elásticas de un trozo de hueso dependen de la distribución de los poros, a través del tensor L, del nivel de daño con el tensor d y de las propiedades elásticas del tejido óseo. A continuación, se establece como se relacionan estas variables tensoriales y dichas propiedades del hueso y como se definen estas últimas. Estas propiedades dependen principalmente de la porosidad y de la composición del tejido óseo, concretamente de la fracción de ceniza, 𝛼, que varía durante la mineralización. En el caso isótropo estas variables se relacionaban con la correlación de Hernandez et al. [22,35] 𝐸(𝑀𝑃𝑎)=84370 𝑣𝑏2,58· 𝛼2,74 El proceso de mineralización ya fue descrito en el modelo isótropo y aquí se supondrá idéntica evolución para el contenido mineral del tejido. La única diferencia entre el caso isótropo y el anisótropo estriba en la distinta forma que ambos modelos tienen de evaluar: ∆𝑣𝑅, ∆𝑣𝐹 y ℎ(𝑡). A continuación, se incorpora el efecto de la porosidad. La influencia de la cantidad de material y de la orientación de los poros se tiene en cuenta de forma independiente. De esta forma lo hicieron Zysset et al. [62]. Estos autores relacionaron las propiedades elásticas del hueso con la densidad aparente y las componentes del tensor MIL, normalizado con traza unitaria. Obtuvieron experimentalmente las siguientes correlaciones para las constantes elásticas: 𝐸𝑖=17607 𝜌3,2𝑚𝑖3,2 𝐸𝑖 𝑣𝑖𝑗=45800 𝜌2,7𝑚𝑖2,4𝑚𝑗2,4 𝐺𝑖𝑗=7799 𝜌3,3𝑚𝑖1,6𝑚𝑗1,6 En las cuales se ve que está claramente desacoplada la dependenciade la porosidad y la anisotropía. Se ha visto que la dependencia de las propiedades elásticas con la porosidad se tiene en cuenta de manera indirectaa, por medio de la fracción volumétrica de hueso. La dependencia con la orientación del entramado óseo se establece con el tensor MIL normalizado de traza unitaria, como hacen Zysset et al. [61, 62]. Los exponentes a los que se elevan sus autovalores, 𝑚𝑖, se han obtenido mediante el ajuste por mínimos cuadrados de resultados experimentales. Sin embargo, hay que tener en cuenta que la interpretación del tensor MIL que se ha hecho en este trabajo es distinta de la que hacen Zysset et al. [62]. Mientras que en este proyecto el MIL representa la longitud promedio de poros, para Zysset et al. representa la longitud promedio de las trabéculas. Ambas están relacionadas y es posible encontrar una correlación experimental con esa otra interpretación de 𝑚𝑖, aunque evidentemente con distintos exponentes y constantes. Por esa razón, se ha optado por mantener los exponentes teóricos de estas expresiones obtenidos del modelo Mecánica del Daño. Así, las componentes del tensor de flexibilidad de un hueso con cierta porosidad, pero sin daño serían: 𝐷 𝑖𝑖𝑖𝑖=1 𝐸𝑖=1 84370 𝑣𝑏2,58· 𝛼2,74 1 (3𝑚𝑖)4 𝐷 𝑖𝑖𝑗𝑗=−𝑣𝑖𝑗 𝐸𝑖=− 𝑣 84370 𝑣𝑏2,58· 𝛼2,74 1 (3𝑚𝑖)2(3𝑚𝑗)2 (3.100) (3.98) (3.99) Modelos de remodelación ósea interna 70 70 𝐷 𝑖𝑗𝑖𝑗=1 2𝐺𝑖𝑗=1+𝑣 2·84370 𝑣𝑏2,58· 𝛼2,74 1 (3𝑚𝑖)2(3𝑚𝑗)2 donde 𝑣 es el coeficinte de Poisson del tejido óseo igual a 0,3. Las propiedades previas están dadas en ejes de ortopría y coinciden con los ejes principales del tensor MIL en el hueso sin daño microestructural. Estas escuaciones reproducen la correlación de Hernandez et al. [22,35] si se considera que el tensor MIL está normalizado con traza unitaria y se particulariza para un material isótropo en el que 𝑚𝑖= 1/3 ∀𝑖. A continuación, se incorpora efecto del daño microestructural al tensor de comportamiento. Para tenerlo en cuenta se sigue la teorí de dañoo anisótropo de Cordebois y Sideroff [63], de la misma forma que hacen Doblaré y García [43] en el modelo de remodelación ósea anisótropo basado en la Mecánica del Daño, con la diferencia de que allí la variable de daño medía la porosidad y no el daño microestructural. Se define en primer lugar la tensión efectiva 𝝈 en función del tensor de daño y la tensión real 𝝈, de tal forma que conserva la simetría, 𝝈=(1−𝒅)−12 𝝈 (1−𝒅)−12 La ecuación constitutiva deldaño establece que la tensión efectiva es aquella que aplicada al material sin daño daría lugar a una energía de deformación igual a la que existe realmente en el material dañado, sujeto a las tensiones reales. Esto es: 𝑈(𝝈,𝒅)=𝑈(𝝈,0) Si se supone un comportamiento elástico lineal, la ecuación anterior conduce al siguiente tensor de fexibilidad del material poroso y dañado, D: 𝐃=𝐑−1𝐑−1𝐃  𝐑−1𝐑−1 donde 𝐃  es el tensor de flexibilidad del material poroso, pero sin daño, cuyas componentes ya hemos definido antes y 𝐑=(1−𝐝)1 2 Con esta expresión se considera la influencia del daño en las propiedades elásticas del material. Sustituyendo en las anteriores, se obtienen las constantes elásticas del hueso real teniendo en cuenta dañoy porosidad. (3.101) (3.102) (3.103) 4 MODIFICACIÓN DEL MODELO MECANOBIOLÓGICO ISÓTROPO 4.1 INTRODUCCIÓN Con el modelo de remodelación ósea isótropo discreto de BMUs se hicieron avances significativos en la discretización de la actividad de las BMUs. Rescindir esta hipótesis del modelo que desvirtuaba los resultados al no representar un fenómeno biológico es un avance de cara a encontrar una rutina más correcta que simule el comportamiento real del hueso. Sin embargo, tras las modificaciones pertinentes para efectuar el cambio en la rutina no todos los resultados fueron los esperados. Se consiguió la activación discreta de BMUs y la evolución de estas en el tiempo, pero el efecto que esto producía en algunas variables básicas como la evolución del daño o la fracción volumétrica de hueso, 𝑣𝑏, no concordaba con lo esperado. Es por esto que se ha analizado la rutina, detectando erratas o mejorando aproximaciones para poder mejorar los resultados. 4.2 MODIFICACIONES PRINCIPALES Evolución del daño En primer lugar se modificó la forma de calcular el daño, puesto que las correlaciones anteriores planteadas por García-Aznar et al. [30] no daban a priori buenos resultados. Por lo que aunque no se cambiaron las correlaciones en sí, si cambió el uso de estas. Por un lado, se eliminó la acumulación de daño por compresion, ya que se considera que Cierra las posibles fisuras del hueso. Modificación del modelo mecanobiológico isótropo 72 72 Para el daño acumulado por tracción se modificó el método de cálculo de la siguiente manera: en cada step, de actividad, se saca el número de ciclos N que habrían sido necesarios para tener el daño que se tiene en ese momento. Este cálculo se hace mediante la expresión que relaciona el daño y N, para cada valor de la deformación efectiva. La forma de esta expresión viene dada en la siguiente imagen. 𝑑𝑡=−1− √1 𝐶3ln(𝑒𝐶3−𝐶2𝜀𝛿2𝑛) 𝛾2 Figura 4.1 Evolución del daño en el tejido óseo, en tracción Una vez calculado n, se suman los ciclos aplicados ese día teniendo un nuevo número de ciclos total. Con este nt se vuelve a usar la correlación obteniendo ahora el daño correspondiente a ese número de ciclos, que será el daño final de ese step. Constantes experimentales Por otra parte, se modificó los valores de las constantes experimentales que se usan para la corrección de las propiedades a fatiga según el valor de la fracción de ceniza. Estos valores fueron sacados del artículo de Martínez-Reina, [55]. (4.1) Parámetro Descripción Valor 𝑡𝑚𝑙𝑡 Tiempo de retraso de mineralzación 12 días 𝑡𝑝𝑟𝑖𝑚 Duración de la fase primaria 10 días 𝑣𝑝𝑟𝑖𝑚 Valor de 𝑣𝑚 tras la fase primaria 0,121 𝑣𝑚𝑎𝑥 Máximo volumen específico de mineral 0,442 𝜅 Tasa de la fase secundaria 0,0005 Cálculo malla MEF En último lugar se solventaron errores cometidos durante el cálculo de los elementos necesarios para la malla. Debido a que este modelo depende fuertemente de la malla este error tuvo que ser resuelto tras un proceso iterativo de pruebas. Como se sabe el tamaño de los elementos está influenciado directamente por dos factores: en primer lugar, por el tamaño de la cavidad de reabsorción, ya que al considerar las BMUs estáticas el elemento debe ser como mínimo del tamaño de esta cavidad y en segundo lugar del valor del Radio de influencia. Este radio delimita una zona en la cual se incluyen todos los elementos en los que si ya hay una BMU activa en el central no se activarán otras en ellos. Por tanto, los elementos deben ser suficientemente pequeños como para que en esa zona de influencia entren varios. El error cometido en los cálculos ofrecía una mayor libertad y no hacía incompatibles estas dos restricciones encontrando un tamaño de los elementos para los que se cumplían ambas. Sin embargo, al observar el fallo y corregirlo, se hacía incompatibles, por lo que se tuvieron que variar una serie de parámetros como la velocidad de avance de la BMU y el factor biológico 𝑓𝐵𝐼𝑂 para poder encontrar un valor de los elementos aceptable. El factor 𝑓𝐵𝐼𝑂 representa el máximo número de BMUs que se van a activar por unidad de volumen y tiempo. 4.3 RESULTADOS A continuación se muestran resultados del modelo que se obtuvieron para distintos tipos de daño inicial y distintos casos de carga. 4.3.1 Daño localizado a tracción Para este primer caso se hicieron simulaciones en las que se consideraban dañados inicialmente solo una serie de elementos del centro de la malla. Se aplicaron unas cargas moderadas a tracción y se compararon los resultados con el modelo anterior en que el que la activación de las BMUs no era discreta. Las condiciones de contorno correspondientes a este caso, figura 4.1, consisten en un empotramiento en una de las caras de los extremos. La cara opuesta se carga en la dirección 𝑥 con distintos valores, generando tensiones de tracción para analizar distintos casos de activación. Modificación del modelo mecanobiológico isótropo 80 80 También podemos ver cómo sería la distribución de volumen reabsorbido tanto para el modelo nuevo, donde para un caso de desuso se observa como el volumen reabsorbido tiene lugar en la línea neutra, como para el modelo anterior, donde se observa de forma más distribuida el volumen reabsorbido. Figura 4.14: Volumen reabsorbido para un caso de desuso en el modelo nuevo Figura 4.15: Volumen reabsorbido para un caso de desuso en el modelo antiguo 5 MODIFICACIÓN DEL MODELO MECANOBIOLÓGICO ANISÓTROPO 5.1 INTRODUCCIÓN En este capítulo se recogen las principales aportaciones de este trabajo de Fin de Master a modelos mecanobiológicos de remodelación ósea interna. La importancia de estas mejoras se basa en la necesidad de acercar a la realidad los resultados de la evolución en el tiempo de las propiedades mecánicas del hueso. Para poner en pie los modelos matemáticos que representan el comportamiento del hueso es necesario un soporte informático, una herramienta que permita implementar el complejo algoritmo que son estos modelos. Para los modelos explicados en este trabajo la herramienta usada es Abaqus, el cual permite, simular utilizando el método de los elementos finitos el comportamiento del hueso a largo plazo. Sin embargo, es necesario programar además, con algún lenguaje de programación, en este caso FORTRAN (Formula Translating System) las ecuaciones que según se ha estudiado rigen ese comportamiento. El conjunto de ecuaciones, variables, condiciones y en general todo lo que conforma la programación del modelo se llaman subrutina. El modo de proceder en esta combinación de programación más simulador de elementos finitos se observa en la figura 5.1.: Figura 5.1: Funionamiento de una subrrutina con Abaqus como programa principal Modificación del modelo mecanobiológico anisótropo 82 82 De forma más específica: partiendo de unos valores iniciales de daño, fracción de ceniza, fracción volumétrica y de todas las conocidas como variables de estado de la subrutina, de unas condiciones de contorno, de unas cargas aplicadas al hueso, de una definición de la malla a usar en la resolución del método de elementos finitos y de un valor inicial de las propiedades mecánicas del hueso, Abaqus calcula la tensión y la deformación en cada punto de integración de la malla. Estos valores los recibe la subrutina, conocida como UMAT, en la cual está implementado el modelo de remodelación ósea y realiza los cálculos necesarios para obtener la variación en las propiedades mecánicas. Con esto actualizamos los tensores de tensión, deformación y todo lo necesario y así, podemos repetir el proceso el número de veces que queramos que dure la simulación. A diferencia de la rutina asociada al modelo isótropo sobre el que se habló en el capítulo 3 y en el 4 la programación de este modelo adquiere un nivel de dificultad muy elevado: considerar la anisotropía del hueso, la difusión de las BMUs y todos los mecanismos físicos y biológicos tenidos en cuenta en este caso carga la subrutina de una complejidad mayor que dificulta su comprensión previa y por tanto su modificación. Esto es debido a que cualquier modificación se aclopa en muchasfases y partes del algoritmo dificultando su resolución airosa. Mientras que a nivel conceptual discretizar la actividad de las BMUs es intuitivo y relativamente sencillo de comprender, incluir en las líneas de programación de este modelo las modificaciones necesarias para ello sobre lo ya programado en el algortimo ha sido un trabajo dificultoso. En este capítulo se explicará inicialmente cuales son las bases de estas modificaciones, qué sentido biológico tienen y en qué consisten. A continuación, se desarrollarán las funciones de las variables más importantes y representativas definidas en el modelo. Despuésse explicará en que consisten las condiciones y circunstancias que se tienen que dar para que se activen BMUs en el hueso. Finalmente se hará un repaso por las distintas limitaciones y problemática encontrada en el modelo, así como por los diferentes resultados obtenidos tras el desarrollo programado. 5.2 ACTIVIDAD DISCRETA DE LAS BMUS Cuando en un punto del hueso el daño es demasiado elevado se activa en ese lugar una BMU para repararlo y este hecho inhibe la posibilidad de que se activen más BMUs alrededor de esta hasta que muera. Sin embargo, el modelo usado aquí no era eso lo que realizaba, en lugar de activar una BMU en el punto dañado activaba distintas fracciones de BMU de los puntos de alrededor que acudían a reparar el mismo daño. Modelo previo Supongamos que tenemos una sección de hueso sometida a una serie de cargas por otro lado con una malla definida de forma que un elemento tenga ocho puntos de integración de Gauss, como se muestra en la figura 5.2: Figura 5.2: Modelo de una sección de hueso sometida a cargas. Si al inicio de una simulación, día 0, en distintos puntos, aparece un daño lo suficientemente grande como para hacer desaparecer la señal inhibitoria y necesitar remodelarse, el modelo anterior estaba implementado de forma que se activaban distintas fracciones de BMUs en los puntos de alrededor a efectuar la remodelación, cómo se puede ver en la figura 5.3. Fracciones de BMU se activaban en distintos puntos de integración teniendo únicamente en cuenta para ello el valor de la señal inhibitoria, sin considerar en ningún caso si había más BMUs activas en la zona. Figura 5.3: Modelo de una sección de hueso sometida a cargas en la que se activan fracciones de BMU (día 0) Estas fracciones de BMUs iban siguiendo una dirección en base al daño y al estímulo mecánico del entorno para reparar el daño cercano, según la trayectoria que se muestra en la figura 5.4. a seguir por Modificación del modelo mecanobiológico anisótropo 84 84 las BMUs activadas el día 0 una vez comienza el día 1. Además de seguir una cierta trayectoria las BMUs ya activas, el día 1 volvían a activarse BMUs donde la señal inhibitoria fuera baja. Fracciones de BMUs se activaban donde fuera necesario y se movían en la dirección pertinente. De esta manera, día tras día se activaban fracciones de BMUs en un mismo punto y se desplazaban repartiendo su efecto remodelatorio allá donde se movieran. Los elementos naranjas de la figura 5.4 representan las fracciones de BMUs activadas ese mismo día. Figura 5.4: Modelo de una sección de hueso sometida a cargas en la que se activan fracciones de BMU y avanzan las ya activadas (día 1) Modelo modificado Sin embargo, y en contraposición del modelo planteado anteriormente, las evidencias obtenidas de la investigación concluyen que si en una zona el daño es elevado y debe remodelarse el hueso se activara una única BMU ahí que reabsorberá y formará hueso para reparar el daño y además se desplazará reparando el daño de zonas cercanas en las que sea necesario. Por ello, en una zona del entorno del punto de activación inical, delimitada por un radio de influencia 𝑅𝑖𝑛𝑓, no se activarán más BMUs, puesto que la ya existente reparará todo el daño de dicha zona durante el transcurso de su vida. Es decir, si debido a las necesidades remodelatorias del hueso, se activa una BMU el día 0 en un determinado punto, se establece una zona de influencia alrededor de este definido por 𝑅𝑖𝑛𝑓 como se muestra en le figura 5.5. Dentro de este volumen no se podrán activar más BMUs mientras la actividad existente y desarrollada por esa BMU activada no cese. Figura 5.5: Modelo nuevo de una sección de hueso sometida a cargas en la que se activa una BMU (día 0). Tras la activación de una BMU, el día 1 tomará una dirección al igual que dictaba el modelo previo en función del daño y el estímulo mecánico del entorno, de manera que acuda a la zona con más necesidad de remodelación. Mientras, en la misma zona delimitada previamente como se observa en la figura 5.6. seguirá inhibiendose la activación de otras BMUs. Figura 5.6: Modelo nuevo de una sección de hueso sometida a cargas en la que se activa una BMU que avanza (día 1) Modificación del modelo mecanobiológico anisótropo 86 86 Esto se ha implementado en el algoritmo del modelo estableciendo una serie de condiciones de activación que permitan o inhiban la posibilidad de activar BMUs. Con el fin de simplificar la programación el bloque principal para la discretización de la subrutina se divide en dos partes: activación de las BMUs y reparto del efecto de las BMUs. En activación se imponen una serie de condiciones que se explicaran en secciones posteriores. Si se cumplen las cuatro condiciones se activará una BMU en el punto de Gauss en el que se esté. Si alguna de las condiciones no se cumple ese día en ese punto de Gauss no se activará ninguna BMU. El bloque del reparto de las BMUs se ha mantenido en esencia como en el modelo previo, adaptándolo a la discretización. Se divide en dos bucles: el bucle de reabsorción y el bucle de formación. En cada punto de gauss se analiza qué BMUs activas están lo suficientemente cerca como para causar un efecto de reabsorción de hueso. Se suman los efectos de todas las BMUs que contribuyan a esto, y análogamente se hace el de formación. Obteniendo las contibuciones de las BMUs activas que está en la zona de influencia del punto de Gauss analizado podemos calcular, las variaciones, de la porosidad y demás variable sdel hueso como se explicó en capítulos anteriores. Figura 5.7: Posibilidades de actividad relacionada con las BMUs: activación o reparto 5.3 DEFINICIÓN DE VARIABLES Las variables más importantes empleadas en este modelo en cuanto a la discretización de la actividad de las BMUs son CANDIDATOS, BMUPG y SENAL1 5.3.1 CANDIDATOS CANDIDATOS es una matriz de datos que ya estaba definida en el modelo anterior, pero cuya función en los objetivos de este es fundamental por lo que es importante destacarla en esta sección. Es una variable global que toma valores al principio de la rutina y permanecen constants a lo largo de las simulaciones. Su misión es recoger la siguiente información para cada punto de Gauss de cada elemento: qué puntos de Gauss de qué elementos se encuentran a menos de una distancia 𝑅𝑖𝑛𝑓 del punto de integración en cuestión. En la siguiente imagen, figura 5.8, vemos un ejemplo de como se iría creando candidatos para el punto de Gauss j del elemento i. Se calcularía la distancia a cada uno de los puntos de la malla, en este ejemplo sería la distancia al punto de gauss l del element j. Si esa distancia es menor que 𝑅𝑖𝑛𝑓 entonces se incluye a ese puntos en la posición adecuada de los candidatos correspondientes al punto 𝑗 del elemento 𝑖. Si no es menor se pasa al siguiente punto. Figura 5.8: Procedimiento a seguir para crear Candidatos. Esta variable se usa para saber cuando inhibir la posibilidad de que se activen BMUs en una zona y funciones más específicas que se explicarán en el apartado de Condiciones de activación discreta. Modificación del modelo mecanobiológico anisótropo 88 88 5.3.2 BMUPG BMUPG es una de las variables básicas del modelo creadas para poder registrar la actividad de las BMUs activadas. Esta variable se compone de un vector de valores por cada punto de Gauss de cada elemento. Es decir, una serie de valores por cada punto de integración del modelo. El vector de datos se compone de la siguiente información: edad de la BMU, posición del frente guía, posición del frente de reabsorción y posición del frente de formación. Si en un punto de integración se dan las cuatro condiciones de activación discreta la componente 1 de su BMUPG se pondrá a 1 (edad de la BMU) y se actualizarán las componentes del frente guía y del frente de reabsorción., pasando lo mismo con las componentes del frente de formación cuando esta fase de la remodelación ósea comience a tener lugar. Si en el punto de integración a analizar, no se dan las circunstancias de activación, todas las componentes irán a 0. Si por el contrario ya hay una BMU activa asociada a ese punto de Gauss, su edad se aumentará una unidad y las direcciones que toman los distintos frentes se actualizaran como corresponda. Esta variable nos permite controlar las BMUs que se activan, su posición y cómo evolucionan. De esta manera la hora de decidir si activar o no BMUs podemos saber si en el entorno ahí observando la variable BMUPG de los CANDIDATOS del punto de integración a analizar. A continuación, se muestra una tabla ejemplificativa con las distintas componentes de las que constaría la variable BMUPG. ELEMENTO PUNTO DE GAUSS INFORMACIÓN 1 1 Edad BMU Componente X FRENTE GUÍA Componente Y FRENTE GUÍA Componente Z FRENTE GUÍA Componente X FRENTE REABSORCIÓN Componente Y FRENTE REABSORCIÓN Componente Z FRENTE REABSORCIÓN Componente X FRENTE FORMACIÓN Componente Y FRENTE FORMACIÓN Componente Z FRENTE FORMACIÓN 2 … 3 … 4 … 5 … 6 … 7 … 8 … 2 1 … 2 … 3 … Figura 5.9: Estructura de la variable BMUPG 5.3.3 SENAL1 Esta variable, al igual que la anterior proporciona información de los puntos de integración de la malla. La información almacenada es el valor de la señal inhibitoria que emiten los osteocitos, calculada en base al daño y al estímulo mecánico. Disponer de esta información nos sirve para hacer una comparativa entre qué puntos de una cierta zona requieren más actividad remodeatoria porque tengan una señal inhibitoria menor de manera que podamos tener en cuenta este factor a la hora de imponer las condiciones de activación. SENAL1 es un vector de tantas componentes como (𝑁º 𝑑𝑒 𝑒𝑙𝑒𝑚𝑒𝑛𝑡𝑜𝑠 𝑥 𝑁º 𝑑𝑒 𝑝𝑢𝑛𝑡𝑜𝑠 𝑑𝑒 𝐺𝑎𝑢𝑠𝑠) haya. 5.4 CONDICIONES DE ACTIVACIÓN DISCRETA Para que se efectuara la activación correcta de una BMU cuando fuera necesario y diera comienzo al proceso de remodelación ósea se definieron cuatro condiciones de activación. Las cuatro condiciones que se implantaron para que se pudiera activar una BMU un punto de integración son: • No hubiera BMUs activas en el punto. • La señal inhibitoria del punto fuera menor que un valor límite llamado señal minima. Modificación del modelo mecanobiológico anisótropo 96 96 Figura 5.14: Evolución de la edad de las BMUs con 𝑑0=0,25 Al no ser una carga muy elevada se observa como el número de BMUs originadas es reducido. En la siguiente figura se ve como evoluciona la variable que representa la edad de las BMUs para cada punto de Gauss. En la mayoría de los puntos la edad de la BMU se mantiene a cero, y en algunos casos al activarse la BMU esta variable evoluciona de manera proporcional como en la imagen: Figura 5.15: Evolución de la edad de las BMUs para cada punto de Gauss Día 1 Día 10 Día 25 Día 35 Si obtenemos imágenes de la fracción volumétrica de hueso como se observa en la siguiente figura, puede comprobarse como para las zonas en las que se activan BMUs la fraccion volumétrica de hueso comienza a variar debido a la actividad remodelatoria de los osteoclastos: Figura 5.16: Evolución de la fracción volumétrica de hueso con 𝑑0=0,25 Carga moderada y 𝒅𝟎=𝟎,𝟎𝟓 Para casos en los que el daño inicial del hueso es considerablemente menor y la carga se mantiene como en el caso anterior se activan muchas menos BMUs y su efecto en la fracción volumétrica de hueso es más imperceptible: La distribución de tensiones de Von Mises toma unos valores más reducidos al inicio y al final de la simulación y al ser una carga no muy elevada y no estar apenas dañado el hueso hay muy poca variación, como se ve en la Figura 5.14: Día 1 Día 10 Día 25 Día 35 Modificación del modelo mecanobiológico anisótropo 98 98 Figura 5.17: Distribución de tension de Von Mises a comienzo y al final de la simulación. La activación de BMUs es por tanto mucho menos necesaria, ya que el porcentaje de tejido dañado que es necesario reparar es mucho menor. De esta manera, en la figura 5.18 se observa como solo se activan algunas BMUs en la zona de aplicación decarga y en la linea neutral de la malla. Figura 5.18: Evolución de la edad de las BMUs con 𝑑0=0,05 De igual manera, si representamos mediante una gráfica, figura 5.19, cómo evoluciona la edad de las BMUs activas, se observa como se activan un mínimo número de estas que reparan el poco daño encontrado. Día 1 Día 10 Día 25 Día 35 Figura 5.19: Evolución de la edad de las BMUs para cada punto de Gauss La fracción volumétrica de hueso se verá afectada en zonas muy concretas donde se activen BMUs. Como muestra la figura 5.20 estas variaciones de tejido óseo ocasionadas por los osteoclastos, será poco perceptible en la diáfisis. Figura 5.20: Evolución de la fracción volumétrica de hueso con 𝑑0=0,05 Día 1 Día 35 Día 25 Día 10 Modificación del modelo mecanobiológico anisótropo 100 100 Carga elevada y 𝒅𝟎=𝟎,𝟐𝟓 Si por el contrario se aumenta el valor de la carga aplicada, adquiriendo valores de 8 𝑁 en cada nodo del extremo opuesto al empotramiento se puede apreciar como aumenta el número de BMUs activadas. Además, las tensiones adoptan valores mucho mayores, dándose la siguiente distribución de tensiones: Figura 5.21: Distribución de tension de Von Mises a comienzo y al final de la simulación. Para este caso, las BMUs se activan de manera más abundante como se puede ver en la Figura 5.22, debido a la necesidad de reparar el daño ocasionado por el alto valor de la carga aplicada. Figura 5.22: Evolución de la edad de las BMUs con 𝑑0=0,25 y carga elevada Día 1 Día 25 Día 35 Día 10 Representando la edad de las BMUs en una gráfica que recoja elcomportamiento de esta variable para todos los puntos a lo largo del tiempo, de la misma manera que en los casos anteriores, se aprecia como el número de BMUs activadas es mucho mayor y mucho más seguido. Figura 5.23: Evolución de la edad de las BMUs para cada punto de Gauss Según estas premisas la fracción volumétrica de hueso evolucionará consecuentemente según se puede ver en la figura 5.24, variando en todos esos puntos en los que se activan BMUs, de manera que se reducirá su valor por la acción de los osteoclastos, para aumentar posteriormente cuando los osteoblastos comiencen la formación de hueso. Día 1 Día 10 Modificación del modelo mecanobiológico anisótropo 102 102 Figura 5.24: Evolución de la fracción volumétrica de hueso con 𝑑0=0,25 y carga elevada Día 25 Día 35 6 CONCLUSIONES Y DESARROLLOS FUTUROS 6.1 CONCLUSIONES Puede establecerse que los objetivos de este proyecto eran fundamentalmente dos: • Pulir la rutina del modelo isótropo obteniendo resultados coherentes tras la discretización de la actividad de las BMUs. • Discretizar la actividad de las BMUs en un modelo de remodelación ósea que considera el hueso un material anisótropo y la difusión de las BMUs Respecto al primer objetivo, y a la vista de los resultados se puede concluir que se ha cumplido. La evolución de las BMUs una vez discretizada su actividad muestra unos resultados coherentes con la evolución real del hueso, sucediéndose las distintas fases de reabsorción, inversión y formación que pueden verse reflejadas en la evolución de la fracción volumétrica de hueso. Se han acotado los problemas que originaban la pérdida de veracidad de los resultados obtenidos en trabajos previos, permitiendo esto introducir distintos cambios en el modelo. Entre estos cambios están el cálculo del daño permitiendo solo la acumulación de este por tensiones de tracción, la modificación de las constantes experimentales, basándose en investigaciones nuevas y la modificación de la malla de elementos finitos usada para las simulaciones. Tras estos cambios, se ha conseguido implementar como algoritmo un modelo mecanobiológico isótropo que considera la activación discreta de las BMUs de forma que realizando simulaciones con elementos finitos se puede ver la evolución del hueso de manera más realista que antes. Sin embargo, las limitaciones de este modelo siguen siendo destacables. No considerar la anisotropía del hueso, ni la difusión de las BMUs son premisas que distan mucho de la realidad. Además, la implementación del modelo está hecha de manera que cualquier resultado o posible extrapolación del mismo a distintos casos o formas de hueso no es viable. Esto se debe a que el correcto funcionamiento del modelo depende de la malla definida. Dependecia que sería interesante erradicar en un posible desarrollo futuro. Conclusiones y desarrollos futuros 104 104 Sobre la discretización en el modelo anisótropo de difusión de BMUs se ha implementado un algortimo, que permite eliminar una hipótesis simplificativa del modelo previo, según la cual las BMUs se activaban de manera continua, aspecto que la experiencia nos dice que no sucede y desvirtua los resultados. Las modificaciones implementadas en el algoritmo del modelo se basan en una serie de condiciones de activación dicreta que respondena a las siguientes preguntas a la hora de realizar las simulaciones en todos los puntos de integración de la malla del hueso: ¿hay ya una BMU activa en el punto?, ¿la señal inhibitoria del punto es menor que el parámetro 𝑠𝑚𝑖𝑛?, ¿es la señal inhibitoria del punto la de menor valor de la zona de influencia? ¿hay alguna BMU activa en la zona de influencia con su frente de reabsorción suficientemente lejos? Según las modificaciones llevadas a cabo, si no hay BMU en el punto, si su señal inhibitoria es menor que 𝑠𝑚𝑖𝑛 y que toda la zona de influencia y si no hay BMUs activas en la zona de influencia con su frente de reabsorción cerca se activará una única BMU en el punto. El tratamiento sobre las BMUs ha conseguido llevarse a cabo con éxito, consiguiendo la activación discreta y evolución deseada. Las BMUs se activan y evolucionan conforme a las evidencias de la experiencia que se han obtenido por lo que supone un avance eiminar estas hipótesis. Sin embargo, a la hora de realizar simulaciones con Abaqus y observar como se han visto afectadas el resto de variables no comprueba que no se han conseguido los resultados esperados. Esto se debe a que adapatar la rutina previa a las modificaciones implementadas no es inmediato ni sencillo. Cada variable se ve afectada en muchas líneas de código por cualquier variación. La rutina requiere ser pulida y repasar todas las posibles incongruencias que pueda haber entre el modelo previo y las modificaciones implementadas. Tras esto debería poderse obtener simulaciones largas que representen la evolución del hueso de manera que pueda preverse con cierta fiabilidad el comportamiento del hueso y la variación de sus propiedades mecánicas a lo largo del tiempo. 6.2 DESARROLLOS FUTUROS Como ya se ha explicado ambos modelos requieren cambios, modificaciones y avances que permitan dar un paso más en la remodelación ósea y su implementación en softwares de simulaciones. Sería interesante seguir mejorando, puliendo los errores cometidos en los dos modelos y solventando las limitaciones encontradas que restringían resultados. Alguna de las posibles mejoras a seguir serían: 1. Extender el modelo isótropo a elementos con varios puntos de Gauss. 2. Modificar el modelo isótropo generlizando su aplicabilidad de manera que no sea dependiente de la malla de elementos finitos usada. 3. Considerar en el modelo isótropo el movimiento de las BMUs teniendo en cuenta su avance y giro aplicando las investigaciones hechas. 4. Mejorar la eficiencia computacional de los modelos de manera que las simulaciones no ralenticen tanto el proceso de investigación. 5. Pulir la rutina del modelo anisótropo consiguiendo obtener simulaciones prolongadas en el tiempo con resultados admisibles según la investigación. 6. Analizar con detalle los distintos casos de porosidad del hueso de manera que pueda incluirse en ambos modelos la aplicación al hueso trabecular. 7. Incluir el efecto de las poblaciones celulares para el desarrollo de las BMUs en sus osteoclastos y osteoblastos.