Repositorio Institucional de Documentos

Proyecto Fin de Carrera Ingeniería Industrial Modelo computacional de la migración celular en matrices 3D basado en estímulos químicos y mecánicos Autor Javier Santos Antón Director Carlos Borau Zamora Ponente José Manuel García Aznar Departamento Ingeniería Mecánica Escuela de Ingeniería y Arquitectura (EINA) Convocatoria Junio Curso Académico 2011/2012 1 Agradecimientos Quiero que estas primeras líneas sirvan como muestra de agradecimiento a todas las personas que me han ayudado y apoyado durante la realización de este Proyecto. En primer lugar y de una forma muy especial, a Carlos Borau y José Manuel García, director y ponente de este Proyecto, por la gran ayuda y atención recibida, sin la cual la realización del mismo no hubiese sido posible. De igual modo me gustaría agradecer al resto de componentes del Área de Mecánica de Medios Continuos y Teoría de Estructuras por haberme permitido trabajar junto a ellos durante este tiempo. Por último, me gustaría acordarme de todas aquellas personas que de una forma u otra están y han estado a mi lado hasta llegar a donde he llegado, porque sin su ayuda no hubiese sido posible. Especialmente de mis padres por su cariño, paciencia y comprensión en todo momento. 2 3 Resumen El presente proyecto tiene como objetivo principal, el desarrollo de un modelo computacional 3D de la migración de una célula individual embebida en una matriz de gel. Con ello se pretende contribuir a la comprensión de los procesos subyacentes en este proceso, así como los diferentes estímulos que afectan al proceso general de la migración. Por otra parte se plantea ofrecer una plataforma de simulación de nuevos experimentos que contribuya a la obtención de unos mejores resultados y permita determinar de una forma predictiva cuáles serán las condiciones de ensayo más optimas. Este proyecto se enmarca dentro de un proyecto de investigación sobre el modelado computacional multiescala de la migración celular. En concreto el modelo aquí desarrollado ofrece un enfoque a escala celular, para lo que se lleva a cabo una discretización de la célula mediante vóxels. Esto posibilita el estudio de características referentes a la morfología de la célula como el factor de forma, el área de contacto célula matriz o la distribución de tensiones en la superficie celular. Al mismo tiempo, también permite el estudio de características propias de un modelo de migración a escala macro como son las velocidades o la distancias recorridas por la célula en el proceso migratorio. El proceso de migración celular se determina a través de la aparición/desaparición de los vóxels que constituyen la célula, al actualizarse la forma de la célula en el transcurso del tiempo se computa la migración como el cambio de situación del centroide geométrico del conjunto de elementos que la constituyen. Las funciones que rigen la aparición/desaparición de vóxels recogen los estímulos que condicionan la migración celular en matrices 3D. En este modelo, de los múltiples estímulos que pueden condicionar la migración celular se han tenido en cuenta los factores mecánicos (mechanotaxis), químicos (chemotaxis) y de flujo. La mechanotaxis se implementa mediante un modelo de mechano-sensing basado en un trabajo anterior[1], según el cual la célula contrae su citoesqueleto y ejerce fuerzas sobre la matriz extracelular (ECM) en función de las propiedades de la misma. Los estímulos químicos y de flujo, se obtienen mediante la simulación de un dispositivo microfluídico utilizando para ello las condiciones de un experimento real [2]. Los resultados obtenidos son postprocesados ofreciéndose diferentes modos de visualización entre los que destacan representaciones animadas de la célula migrando a través de la ECM. Estas animaciones muestran a la vez el estado tensional, así como la morfología de la célula en cada instante. 4 5 Indice Resumen ........................................................................................................................................ 3 Lista de Figuras .............................................................................................................................. 7 Lista de Tablas ............................................................................................................................... 9 1. Introducción ........................................................................................................................ 11 1.1 Motivación .................................................................................................................. 11 1.2 Origen del proyecto ..................................................................................................... 11 1.3 Migración celular: estado del arte .............................................................................. 12 1.3.1 Experimentación ................................................................................................. 12 1.3.2 Modelado ............................................................................................................ 13 2. Modelo ................................................................................................................................ 15 2.1 Modelo previo. ............................................................................................................ 15 2.1.1 Mecanismo de mecano-sensing .......................................................................... 15 2.1.2 Adaptación del citoesqueleto ............................................................................. 17 2.1.3 Migración ............................................................................................................ 18 2.2 Modelo actual ............................................................................................................. 19 2.2.1 Mechano-sensing ................................................................................................ 20 2.2.2 Factores fluido-químicos ..................................................................................... 21 2.2.3 Funciones de probabilidad .................................................................................. 22 2.2.4 Degradación ECM ................................................................................................ 25 2.3 Implementación numérica .......................................................................................... 25 2.3.1 Análisis fluido-químico ........................................................................................ 26 2.3.2 Análisis mecánico ................................................................................................ 27 2.3.3 Cálculo de probabilidades .................................................................................. 28 2.3.4 Actualización forma célula .................................................................................. 30 2.3.5 Técnicas de visualización ..................................................................................... 31 2.3.6 Software utilizado ............................................................................................... 31 3. Resultados ........................................................................................................................... 33 3.1 Simulación del microdispositivo .................................................................................. 33 3.2 Simulación fluido-química a nivel celular .................................................................... 35 3.3 Simulación mecánica ................................................................................................... 38 3.4 Simulación completa ................................................................................................... 43 3.4.1 Caso (i) ................................................................................................................. 43 6 3.4.2 Caso (ii) ................................................................................................................ 44 3.4.3 Resultados ........................................................................................................... 44 4. Discusión ............................................................................................................................. 49 5. Conclusiones y Trabajo Futuro ............................................................................................ 53 6. Referencias .......................................................................................................................... 55 Anexo I: Estado del arte. Experimentación migración celular. ................................................... 59 Anexo II: Estado del arte. Tipos de modelos existentes para simular la migración celular. ....... 63 Anexo III: Análisis paramétrico de las funciones de probabilidad. .............................................. 65 Anexo IV: Validación mecánica del modelo computacional. ...................................................... 71 Anexo V: Validación fluido-química del modelo computacional. ................................................ 77 Anexo VI: Determinación factor de forma de la célula. .............................................................. 81 7 ListadeFiguras Figura. 1: Posibles estímulos condicionantes de la migración celular. ....................................... 14 Figura. 2: Esquema de las diferentes partes de la célula involucradas en el proceso de mechano-sensing y su modelado. ............................................................................................... 16 Figura. 3: Fuerza ejercida por el sistema actina-miosina (AM) en función del solapamiento entre filamentos ( c ε ). ................................................................................................................. 16 Figura. 4: Representación esquemática reorientación citoesqueleto celular ............................ 18 Figura. 5: Esquema iterativo del modelo de migración celular en cada incremento de tiempo.19 Figura. 6: Dispositivo microfluídico utilizado para los ensayos. .................................................. 27 Figura. 7: Posición y dimensiones iniciales de la célula. ............................................................. 28 Figura. 8: Esquema de la conectividad considerada-6 utilizado en la simulación. ..................... 29 Figura. 9: Esquema de la migración en función de la dirección del gradiente de tensión. ........ 29 Figura. 10: Técnicas de visualización. .......................................................................................... 31 Figura. 11: Software utilizado en la simulación, postprocesado y visualización de resultados. . 32 Figura. 12: Simulación fluídica del dispositivo. ........................................................................... 34 Figura. 13: Caída de presión en el dispositivo microfluídico. ...................................................... 35 Figura. 14: Evolución de la célula para la simulación fluido-química a nivel celular. ................. 35 Figura. 15: Volumen de control considerado par simulación fluido-química a nivel celular. ..... 36 Figura. 16: Líneas de corriente y factor autocrino en la simulación fluido-química a nivel celular. ..................................................................................................................................................... 37 Figura. 17: Distribución de presiones causadas por el flujo intersticial en la superficie celular. 37 Figura. 18: Simulación de difusión de un factor químico a nivel celular..................................... 38 Figura. 19: Gradiente de rigidez del gel (ECM) a lo largo de la dirección X. ............................... 39 Figura. 20: Velocidad y distancia recorrida por la célula tras 16 y tras 100 steps. ..................... 40 Figura. 21: Área proyectada y factor de forma para simulación mecánica. ............................... 41 Figura. 22: Mapa de tensiones en la célula y desplazamientos en la matriz extracelular. ......... 42 Figura. 23: Degradación de la matriz extracelular. ..................................................................... 43 Figura. 24: Esquema de los diferentes estímulos considerados en la simulación. Caso (i). ....... 44 Figura. 25: Esquema de los diferentes estímulos considerados en la simulación. Caso (ii). ...... 44 Figura. 26: Trayectoria de la migración celular bajo la influencia de diferentes estímulos........ 45 Figura. 27: Resultados simulación basada en estímulos mecánicos y químicos. ........................ 46 Figura. 28: Célula adherida al extremo del dominio. .................................................................. 47 Figura. 29: Representación esquemática de los distintos tipos de dispositivos microfluídicos. 60 Figura. 30: Imágenes capturadas mediante 3D TFM durante la migración celular. ................... 61 Figura. 31: Variación de las funciones de probabilidad en función de los estímulos considerados. .............................................................................................................................. 66 Figura. 32: Variación de las funciones de probabilidad con el ángulo θ. .................................... 67 Figura. 33: Variación de las funciones de probabilidad con la constante de sensibilidad σ λ ..... 68 Figura. 34: Variación de las funciones de probabilidad con K⁰ ................................................... 69 Figura. 35: Velocidad media y velocidad efectiva para simulación............................................. 72 Figura. 36: Distancia total y distancia efectiva recorrida por la célula. ...................................... 72 Figura. 37: Factor de forma de la célula para diferentes valores de rigidez. .............................. 73 1. Introducción 14 Figura. 1: Posibles estímulos condicionantes de la migración celular. El comportamiento de las células viene determinado por diversidad de factores que pueden dividirse estímulos químicos (haptotaxis y chemotaxis), mecánicos (mechanotaxis, tensotaxis y durotaxis) y topográficos (contact guidance y curvature) [17]. Existe un gran debate sobre el papel de cada uno de los factores según los diferentes tipos de células y procesos. Una misma célula puede comportarse de manera diferente en función de las condiciones que la rodean, e incluso un leve cambio puede producir resultados contradictorios. Es por ello que cada vez son más los modelos que se desarrollan para intentar arrojar algo de luz, y aportar información valiosa para mejorar el diseño de nuevos experimentos. Existen diferentes modelos de migración, tanto para células individuales, como para poblaciones. Al desarrollarse a la par que los experimentos, existen muchos más modelos bidimensionales, aunque cada día son más los de migración 3D [18-21]. Una explicación más detallada sobre los tipos de modelos existentes y sus principales características se recoge en el Anexo II. Como se puede observar en este breve resumen, aunque existen multitud de aproximaciones, todas se suelen centrar en un solo aspecto concreto de la migración celular. El paso natural, por lo tanto, es tratar de combinar cada una de esas aportaciones desarrollando un modelo multidisciplinar capaz de capturar a diferentes escalas el comportamiento celular y los procesos más relevantes de la migración. 15 2. Modelo Este proyecto, como ya se ha comentado anteriormente, está integrado en el desarrollo de un proyecto de investigación sobre el modelado computacional multiescala de la migración celular. En dicho proyecto se desarrolla un modelo de migración celular (a escala macro) de elementos finitos en el cual, la célula se simula como un solo elemento. Su movimiento queda definido por un equilibrio de fuerzas celulares cuya relación con la matriz extracelular (ECM) se establece a través del proceso de mecano-sensing, como se explicará brevemente a continuación. Por otro lado, se simula la red de actina que compone el citoesqueleto celular y se estudia el fenómeno del mecano-sensing desde la perspectiva micro. El modelo desarrollado en el presente proyecto, corresponde también a una aproximación macro-escala de una única célula que incorpora no sólo estímulos mecánicos (mechano-sensing), sino también, químicos y de flujo. Las diferencias con el modelo anterior son que el presente no se trata de un modelo basado en fuerzas y, además, se adopta un enfoque más centrado en la célula. Para ello la célula es modelada como un conjunto de elementos, cada uno de los cuales presenta comportamientos validados con el modelo micro. Con todo ello, el modelo es capaz de obtener resultados relativos al enfoque macro como la velocidad de migración, y otros referentes a nivel celular, como la morfología de la célula, la deformación sufrida por el citoesqueleto o la distribución de tensiones en el cuerpo celular. 2.1 Modelo previo. El modelo propuesto en este trabajo se desarrolla a partir de otro anterior [1] basado en elementos finitos y en el fenómeno del mecano-sensing, el cual será referido como modelo “macro” de aquí en adelante. En ese modelo, la célula, embebida en una matriz elástica 3D, se discretiza como un único elemento, no teniendo en cuenta por lo tanto ningún aspecto morfológico de la célula ni ningún estímulo químico, pero manteniendo la interacción mecánica entre matriz y célula. El proceso de migración se divide en 3 procesos diferentes: mechano-sensing, adaptación del citoesqueleto y movimiento celular. A pesar de que estos tres procesos ocurren en la realidad de forma simultanea, en el modelo macro, por simplicidad, ocurren de manera independiente aunque interrelacionada, como se verá a continuación. 2.1.1 Mecanismo de mecano-sensing El citoesqueleto celular es una estructura mecánica que juega un papel fundamental en la locomoción celular [22]. El comportamiento del citoesqueleto celular es altamente dinámico y complejo, y sus características reológicas [23, 24] y contráctiles [25] se han estudiado intensivamente. Un modelo simplificado del citoesqueleto, basado en una estructura de muelles y un actuador, se utiliza para simular el comportamiento mecánico de la célula 2. Modelo 16 (Figura. 2). Los elementos celulares considerados por este modelo son: los filamentos de actina (simulados mediante un muelle de constante act K ), el sistema actina-miosina (simulado mediante un actuador contráctil) y por último la rigidez pasiva de la célula, debida en su mayor parte a los microtúbulos (simulada mediante un muelle de constante elástica pas K ). El citoesqueleto celular se une a la matriz extracelular por medio de adhesiones focales e integrinas transmembrana, que a modo de simplificación se consideran completamente rígidas. Figura. 2: Esquema de las diferentes partes de la célula involucradas en el proceso de mechano-sensing y su modelado. El citoesqueleto se simplifica a una estructura de muelles que simulan la resistencia de los filamentos internos ( act K , referente a la rigidez de los filamentos de actina y pas K rigidez pasiva de la célula producida por microtúbulos), más un actuador ( AM ) que simula el carácter contráctil del sistema actina-miosina. La deformación de la célula depende por tanto de la fuerza ejercida por el actuador y del valor de rigidez con que se simulan las diferentes partes de la célula[1]. La tensión activa generada por el actuador contráctil de acto-miosina ( c p ) se relaciona con el solapamiento entre los filamentos de actina-miosina ( c ε ) (Figura. 3). Figura. 3: Fuerza ejercida por el sistema actina-miosina (AM) en función del solapamiento entre filamentos ( c ε ). Zona 1: Comportamiento pasivo, el AM es sometido a compresión por cargas externas. Zona 2: Contracción, el AM ejerce fuerzas de para detectar las propiedades del medio. Zona 3: Tensión, el AM es traccionado por cargas externas pero todavía tiene capacidad de comprimirse para detectar las propiedades del medio. Zona 4: Comportamiento pasivo, las cargas de tracción son tan elevadas que el AM no puede ejercer una fuerza suficiente para contrarrestarlas. 2.1 Modelo Previo 17 Así pues la tensión que la célula ejerce sobre la matriz ( cell p ) es dependiente de la deformación celular ( cell ε ), que a su vez viene determinada por las condiciones de contorno. Por tanto, cell p se puede expresar en función de los parámetros mecánicos como: ( ) ( ) ( ) ( ) min max min min max min max max max max max max max max ( ) pas cell cell pas cell pas cell cell act act cell cell pas cell pas cell act cell act pas cell cell K K p K p K K p pK p K p K K p K ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε <   − + < < −  = − + < < −  <   ec. 2-1 La célula se considera un cuerpo esférico de forma y volumen constante. Además tiene un comportamiento anisótropo, por tanto la tensión ejercida por la célula tiene dos componentes, una componente isótropa dependiente de la tensión calculada y de la deformación volumétrica (la célula se contrae por igual en todas las direcciones) y una segunda componente anisótropa que depende de la dirección de polarización de la célula, cuya obtención se explicará en la siguiente sección. Por lo tanto, la tensión ejercida por la célula se puede expresar como: ( ) ( ) cell cell cell pol pol cell cell p p ε θ = ⊗ + σd d I ec. 2-2 La resolución del problema implica la resolución de la ecuación de equilibrio global tanto para la célula como para el sustrato: ( ) 0 en cell cell ∇ = σ Ω ec. 2-3 ( ) en subs ext subs ∇ = σfΩ ec. 2-4 2.1.2 Adaptación del citoesqueleto Observaciones experimentales recientes sugieren que el citoesqueleto celular forma fibras de tensión (stress-fibers)[26] orientadas preferentemente paralelas con la dirección principal de tensión[17]. En este modelo esto corresponde con la dirección de menor deformación y por tanto se considera que la célula se alinea gradualmente con ella mediante la siguiente expresión: ( ) ( ) ( ) pol i i pol pol pol d t k t dt ε ∆ = ∆ = × × ∆ d d d d d ec. 2-5 2. Modelo 18 Figura. 4: Representación esquemática reorientación citoesqueleto celular El citoesqueleto celular se reorienta para cada incremento del tiempo en función del vector de polarización anterior( i pol d ) y el vector de variación ( pol ∆d ) [1]. donde pol d es la dirección actual de polarización de la célula, ε d la dirección principal de mínima deformación, α el ángulo que forman y κ la constante de reorientación ( -1 min ). Como se observa en la Figura. 4, la dirección de polarización se actualiza en cada intervalo de tiempo alineándose progresivamente con la dirección de deformación mínima, obteniéndose así la nueva dirección de polarización ( 1i pol + d ). 2.1.3 Migración La última fase del modelo consiste en la resolución del equilibrio de todas las fuerzas implicadas en el movimiento celular. 0 trac drag prot ++= F F F ec. 2-6 Las fuerzas que participan en este proceso son: fuerza de protrusión ( prot F ), calculada como un vector aleatorio y de módulo proporcional a ( ) cell cell p ε , fuerza de arrastre ( drag F ), proporcional a la velocidad y la viscosidad efectiva η de la ECM y la fuerza de tracción ( trac F ) ejercida por la célula, proporcional al área donde esta se aplique ( a a ). Las fuerzas de tracción en la parte frontal y en la parte posterior tienen diferentes valores debido a la diferente densidad de adhesiones focales, ( F β y B β ). Por tanto, estas fuerzas pueden expresarse como: drag cell η = F υ υυ υ ec. 2-7 ( ) trac F B cell a pol p a β β = −F d ec. 2-8 y la velocidad puede ser obtenida como: ( ) ( ) 1 cell prot F B cell a pol p a β β η − = + −F d υ υυ υ ec. 2-9 2.2 Modelo actual 19 2.2 Modelo actual En el nuevo modelo, a diferencia del anterior, el proceso de migración se analiza a nivel celular, en él, la célula pasa de ser discretizada como un único elemento a ser discretizada como un conjunto de vóxels. El vóxel es la unidad mínima procesable de una matriz tridimensional, o en otras palabras, cada una de las unidades cúbicas que compondrá el cuerpo celular y que posibilitará principalmente el estudio de su morfología. Estos vóxels están embebidos en una matriz elástica 3D y aparecen y desaparecen con una probabilidad que viene dada en función de diversos parámetros, lo cual hace avanzar a la célula así como cambiar su forma. Como se ha descrito anteriormente, el proceso de migración celular está influenciado por multitud de factores. Este modelo incorpora factores mecánicos, químicos y de flujo. Se tienen en cuenta por tanto: las propiedades mecano-químicas de la ECM (gradientes de rigidez, concentraciones químicas, direcciones de flujo) así como su geometría. Se trata de un modelo iterativo (Figura. 5) que resuelve el equilibrio mecánico entre célula y matriz en cada incremento de tiempo. Por otro lado, el problema de flujo y difusión química en un dispositivo microfluídico se considera estacionario, con lo que las condiciones fluido-químicas son constantes durante las simulaciones del proyecto actual. Las hipótesis tomadas respecto a cada uno de los factores, son explicadas en detalle en la sección correspondiente. Figura. 5: Esquema iterativo del modelo de migración celular en cada incremento de tiempo. La velocidad de la célula no es resuelta mediante ecuaciones de equilibrio, sino que viene dada intrínsecamente por medio de la aparición y desaparición de vóxels, dependientes a su vez de las funciones de probabilidad. Estas funciones (detalladas en los siguientes apartados) dependen de los diferentes estímulos comentados (mecánicos, químicos y de flujo) así como del incremento temporal. Los elementos que componen la célula son actualizados en cada incremento de tiempo, determinando así la forma, el área, el volumen, el contorno etc. 2. Modelo 20 A la hora de realizar el modelo se ha tratado, en la medida de lo posible, ceñirse lo máximo al proceso real de migración celular en el que se suceden temporalmente los cinco procesos ya mencionados en el apartado 1.3.2 [4]: protrusión de pseudópodos en la parte frontal, formación de contactos focales, proteólisis focalizada, contracción del sistema actinamiosina y desprendimiento del extremo trasero y consecuente avance. Cada uno de estos procesos se incorpora al modelo bien como un proceso en sí o bien como combinación de alguno de ellos. La protrusión de pseudópodos se implementa en dos fases: la primera incluye la formación de filopodios, estructuras de carácter explorativo (simuladas como elementos en las partes más externas del contorno en varias direcciones). La segunda es la formación de lamelipodios (estructuras más amplias que extienden la membrana celular) y se corresponde en el modelo con la aparición de elementos en la parte frontal de la célula. La formación de contactos focales no se modela específicamente, ya que se considera que los elementos perimetrales de la célula están perfectamente adheridos a la ECM. La proteólisis focalizada se asocia con la degradación que experimenta la matriz como consecuencia del paso de la célula, y se implementa como se explica en el apartado 2.2.4. Por último, el desprendimiento de la parte trasera se corresponde con la desaparición de vóxels, que unido a la aparición en la parte delantera producen la migración de la célula. 2.2.1 Mechano-sensing Las células tienen la capacidad de detectar las propiedades mecánicas del medio en que se encuentran y actuar en consecuencia. Tal como se ha explicado anteriormente, este proceso se denomina mechano-sensing. Al igual que ocurría en el modelo previo, el citoesqueleto se modela mediante una estructura consistente en varios muelles y un actuador (Figura. 2), con la salvedad de que en este caso la célula deja de ser un único elemento y se convierte en un conjunto de vóxels, cada uno de los cuales se comporta según el modelo anteriormente descrito. A diferencia del enfoque macro, ya no se tiene en cuenta la dirección de polarización de la célula completa, ya que cada elemento se evalúa de forma individual. Por tanto, en este caso, la ecuación que describe el comportamiento tensional de cada elemento que conforma la célula solo consta del término isótropo. ( ) cell cell cell p θ = I σ σσ σ ec. 2-10 Igualmente, se debe satisfacer la ecuación de equilibrio global mecánico planteada anteriormente ( ec. 2-3 y ec. 2-4). La resolución de este problema mecánico se realizara mediante un software de elementos finitos (Abaqus) como se detalla más adelante en el apartado 2.3.2. La tensión de equilibrio alcanzada por cada elemento de la célula, servirá, entre otros factores, para decidir la aparición de un nuevo elemento o la desaparición del mismo. 2.2 Modelo actual 21 2.2.2 Factores fluido-químicos El proceso de migración celular, además de estar condicionado por los estímulos mecánicos resultantes de la interacción de la célula con el medio, está influenciado por estímulos químicos. Como se ha comentado anteriormente, existen dos modos diferentes de migración celular en función de si los factores químicos presentes en el medio se encuentran diluidos (chemotaxis) o son insolubles y se encuentran dispersos por la ECM (hatpotaxis). Este modelo, por simplicidad, se centra exclusivamente en la chemotaxis. Dentro de la misma, se distingue entre: quimiotaxis autologa (autologue chemotaxis), cuando la propia célula segrega el factor químico que estimula el movimiento celular, o quimiotaxis paracrina (parachrine chemotaxis), término utilizado para denominar la respuesta de la célula frente a un factor químico esterno a ella [27]. Mediante experimentos in-vitro, se ha constatado también la frecuente interacción entre la migración producida por factores químicos y la producida por la presencia de fluido intersticial [28], debido al arrastre de la concentraciones químicas que altera los gradientes percibidos por la célula. En el presente modelo se incorpora la dependencia de la migración respecto a factores fluido-químicos, incorporando la citada dependencia a las funciones de probabilidad para la aparición/desaparición de vóxels. Se precisa por tanto conocer el valor tanto de la dirección del flujo, como de la concentración de factores químicos en cada uno de los elementos que rodean la célula. Para calcular el valor de las líneas y direcciones de flujo, así como del valor de concentración de los factores químicos se utilizará la ley de Brinkman para materiales porosos [29]: 2 0 P µ µκ ∇ − −∇ = ν ν ν νν ν ν ν ec. 2-11 donde µ designa a la viscosidad del fluido, ν νν ν la velocidad del fluido, κ hace referencia a la permeabilidad y P la presión. La determinación de la difusión de un factor químico se obtiene a través de la ley de Fick: D c =− ∇ j ec. 2-12 donde j designa al flujo molar del factor químico difundido, D hace referencia al coeficiente de difusión y c a la concentración de factor químico. Así mismo, se incorporan las ecuaciones de convección-difusión para el transporte de especies diluidas: ( ) D c c R ∇⋅ − ∇ + ⋅∇ = v ec. 2-13 D c c =− ⋅∇ + ⋅ j v ec. 2-14 2. Modelo 22 donde D vuelve a ser la difusividad del medio, c la concentración de la sustancia que difunde, v la velocidad del fluido, R designa la tasa de consumo/secreción de la sustancia difundida y j , al igual que en el caso anterior, el flujo molar del factor difundido. 2.2.3 Funciones de probabilidad La migración de la célula, está basada en la aparición y desaparición de los vóxels que la constituyen. En cada incremento de tiempo, ciertos elementos que eran parte de la célula dejan de serlo y algunos que anteriormente correspondían a la matriz extracelular pasan a ser ocupados por la célula, de modo que el cuerpo de la célula se va desplazando progresivamente durante el transcurso de la simulación. La aparición y desaparición de vóxels viene gobernada por dos funciones de probabilidad en las que se consideran los estímulos mecánicos y fluido-químicos mencionados anteriormente. A cada uno de ellos le corresponde uno o varios parámetros como se detallará a continuación. Las probabilidades de aparición/desaparición de vóxels son respectivamente: ( ) 0 0 max 1 C C F F k d d d d dt p p p e σ σ σ σ λ λ λ λ ∆ ∆ + + + + + + + + − + + + + + +   = + −     ec. 2-15 ( ) 0 0 max 1 C C F F k d d d d dt p p p e σ σ σ σ λ λ λ λ ∆ ∆ − − − − − − − − − + + + − − −   = + −     ec. 2-16 donde , d σ + − , , d σ ∆ + − , , C d + − , , F d + − son parámetros variables en función de las condiciones que afectan a cada vóxel, y k y λ son parámetros constantes ajustados mediante ensayos numéricos (Anexo III). Sus valores y significados pueden ser consultados en las Tabla 1. Designación Significado Valor Unidades 0 p + , 0 p − Probabilidades mínimas de aparición/desaparición 0.1, 0.1 [-] max p + , max p − Probabilidades máximas de aparición/desaparición 0.8, 0.4 [-] 0 k + , 0 k − Ratio de aparición/desaparición 0.4, 0.4 [ 1 min − ] σ λ + , σ λ − Constante de sensibilidad de aparición desaparición relativo a la tensión máxima 0.0035, 0.0035 [-] σ λ ∆ + , σ λ ∆ − Constante de sensibilidad de aparición/desaparición para gradientes de tensión 0.004, 0.004 [-] C λ + , C λ − Constante de sensibilidad de aparición/desaparición frente a un gradiente concentración factor químico 0.3, 0.3 [-] F λ + , F λ − Constante de sensibilidad de aparición/desaparición frente al flujo 0.004, 0.004 [-] dt Incremento de tiempo 5 [ min ] d Parámetro de degradación 0.01 [ - ] Tabla 1: Significado de las constantes presentes en las funciones de probabilidad 2.2 Modelo actual 23 Estas funciones se aplican en cada uno de los elementos del contorno de la célula en función de las condiciones mecano-fluido-químicas. Con los parámetros de dichas funciones, se reflejan los estímulos considerados en este modelo para la migración celular: la magnitud de la tensión celular, el gradiente de tensión, la concentración de factores químicos y la dirección de flujo. La tensión celular se introduce en las funciones a través de d σ + y d σ ∆ + , cuyo valor para la aparición de elementos viene dado por: max cos - 2 2 3 0 2 2 cell d σ σ π π θ θ σ π π θ +  < <   = ≤ ≤   cos - 2 2 3 0 2 2 d σ π π θ θ π π θ ∆ +  < <   = ≤ ≤   ec. 2-17 y para la desaparición por: max 0 - 2 2 3 sin 2 2 cell d σ π π θ σ π π θ θ σ −  < <   = ≤ ≤   0 - 2 2 3 sin 2 2 d σ π π θ π π θ θ ∆ −  < <   = ≤ ≤   ec. 2-18 donde max σ es la tensión máxima generada por la célula (determinada según los parámetros del modelo de mechano-sensing) y θ el ángulo que forma cada elemento con la dirección del gradiente de tensión. Estas definiciones, permiten por un lado enfatizar en la direccionalidad según el gradiente mecánico, así como obtener formas celulares alargadas en la dirección de movimiento, simulando la polarización mostrada por las células en experimentos in-vitro. Los factores químicos que pueden condicionar la migración incluyen a través de C d . Este parámetro compara el valor de la concentración de factor químico en cada elemento con los elementos adyacentes. Por tanto: max cell C C dC ∆ = ec. 2-19 donde cell C ∆ es la variación de la concentración de factor químico en cada elemento considerado y max C el valor máximo que puede tomar la concentración de un determinado factor químico. La dependencia respecto al fluido intersticial se introduce a través de F d . Con ello se tiene en cuenta la dirección del flujo como posible condicionante de la migración celular. En condiciones normales, se ha detectado experimentalmente [2] que poblaciones pequeñas de células migran siguiendo las líneas de flujo, Sin embargo, se ha visto que la razón para ello es el arrastre aguas abajo de factores químicos segregados por la propia célula. Al bloquear ciertos receptores de la membrana, la célula cambia completamente su comportamiento, migrando 2. Modelo 30 probabilidad de aparición de elementos y la mitad trasera de desaparición, con lo cual se espera que la célula avance preferentemente alineada con el factor correspondiente. Por lo tanto, utilizando las funciones descritas hasta ahora, se asigna a cada uno de los elementos adyacentes a la célula un valor de probabilidad. Sin embargo, pueden surgir conflictos como consecuencia de la geometría (un mismo elemento puede aparecer por estar en contacto con diferentes elementos del contorno). Para solventarlos, se computan todas las probabilidades al mismo tiempo, y se elige aquella que sea mayor para cada posible elemento. Cabe recordar, que aparte de las probabilidades condicionadas por factores, también se incluye una probabilidad de carácter estocástico. El cálculo es análogo al del ejemplo anterior, con la única diferencia de que la dirección “preferente” es obtenida aleatoriamente de una distribución equiprobable. Este factor, al igual que los anteriores, puede ser activado o desactivado en función del carácter de la simulación a realizar. Para evitar posibles problemas derivados del excesivo decrecimiento del volumen de la célula o que su tamaño llegue a adquirir valores demasiado grandes, se acota el tamaño celular superior e inferiormente. El tamaño mínimo que la célula puede alcanzar es el de la célula de partida (515 elementos), y el máximo es 1,15 veces su volumen (592 elementos). Este valor se ha definido de manera arbitraria debido a la falta de datos experimentales en 3D en literatura. 2.3.4 Actualización forma célula Una vez se ha asignado un valor probabilístico de aparición/desaparición a los elementos susceptibles de ser añadidos/eliminados, y de nuevo mediante un Script de MATLAB, se generan valores aleatorios para cada elemento y se computa cuales de ellos finalmente aparecerán/desaparecerán. A continuación se actualiza la forma de la célula almacenando diferentes parámetros (volumen, área de adhesión, centroide geométrico, factor de forma etc.) así como los listados de los elementos que son célula, cuáles pertenecen al contorno, etc. Esta nueva geometría, es introducida de nuevo como entrada para el próximo incremento de tiempo, por lo que el problema mecánico es calculado de nuevo con la nueva distribución. Todo ello, además, es almacenado en ficheros para ser post-procesado permitiendo una visualización 3D de los cambios acontecidos en cada incremento de tiempo. Por último, un script se encarga de simular la formación de filopodios. Para ello se calcula el elemento más alejado del centroide de la célula en cada dirección (y cada sentido de la misma) y se añade un elemento extra con una sola cara en contacto con el contorno. Este elemento “exploratorio”, gracias al proceso de mechano-sensing, permitirá a la célula un cálculo más preciso de la dirección de tensión máxima en el próximo incremento. 2.3 Implementación numérica 31 2.3.5 Técnicas de visualización La visualización de los resultados obtenidos se realiza a través de diversas técnicas en función de la naturaleza de los mismos. Por un lado todas las variables de interés (tensiones, desplazamientos, trayectoria, factor de forma, área de adhesión, elementos que conforman la célula etc.) son almacenadas en ficheros de texto en formatos entendibles o bien por ParaView o bien por MATLAB, lo que permite su análisis de forma directa. Por otro lado y por medio de scripts, se automatiza el postprocesado de datos, haciendo que el análisis sea mucho más fácil y rápido. Se hace especial hincapié en enfatizar que la migración simulada es un proceso tridimensional, y para ello los videos son la mejor herramienta para visualizar los resultados. Con tal fin, se ha desarrollado una rutina que permite obtener las coordenadas de los elementos del contorno y transformarlas a un formato legible por el software VMD. Con ello se consigue generar una superficie celular más suavizada y se obtienen las diferentes capturas que compondrán la correspondiente película. Figura. 10: Técnicas de visualización. Las coordenadas del contorno y del centroide celular, son procesadas por un script de Matlab y transformadas a un formato entendible por VMD para una visualización 3D de mayor calidad. 2.3.6 Software utilizado A lo largo de la implementación numérica se ha citado el software utilizado en cada parte de la simulación. Dada la variedad de software utilizado, se introduce esta sección a modo de síntesis (Figura. 11). El código principal con el que se lleva a cabo la computación de la migración celular se ha implementado mediante diferentes scripts de MATLAB. La computación del proceso de migración, se basa en diferentes estímulos, tanto mecánicos como químicos. Los estímulos mecánicos provienen de resolver el problema mecánico (ecuaciones de equilibrio global tanto en la célula como en el sustrato, ec. 2-3 y ec. 2-4) para lo cual se utiliza el software ABAQUS. Los estímulos fluido-químicos provienen de la simulación del microdispositivo. Para ello se resuelven las ecuaciones de convección difusión para el transporte de sustancias químicas, así como la ecuación de Brinkman para flujo a través de medios porosos (ec. 2-11) a través del software COMSOL MULTIPHYSICS. 2. Modelo 32 La visualización de los resultados obtenidos tanto de la computación de la migración como de la resolución del problema mecánico se realiza mediante ParaView y VDM, ambos software libre. Para la simulación fluido-química a nivel celular, es necesario llevar a cabo un postprocesado de los resultados obtenidos. Mediante puntos del contorno se reconstruye por completo una superficie suavizada de la célula. Esto se hace mediante el software Rhinoceros, indicado especialmente para el tratamiento de imágenes. Figura. 11: Software utilizado en la simulación, postprocesado y visualización de resultados. El modelo computacional utilizado en este trabajo, se implementa principalmente a través de diferentes scripts de MATLAB, que componen el código principal. La computación se realiza en base a ciertos estímulos tanto mecánicos como fluido químicos, obtenidos de la realización de simulaciones independientes mediante el software ABAQUS y COMSOL respectivamente. Para ciertas simulaciones específicas, los resultados son procesados mediante un software de tratamiento de imágenes, Rhinoceros. Por último la visualización de resultados y el postproceso se produce a través del software libre ParaView y VMD. 33 3. Resultados 3.1 Simulación del microdispositivo Como se comenta en el apartado 1.1, el efecto de una única célula sobre el estudio del flujo intersticial en un dispositivo completo se considera mínimo y por tanto se asume como despreciable. Además, se busca resolver el dispositivo en situación estacionaria, ya que es una vez desarrollado el flujo cuando se realizan las medidas experimentales. Por lo tanto y dado que las condiciones fluídicas adquieren un valor constante en el transcurso temporal, basta con realizar una simulación del dispositivo de forma previa a la simulación de la migración. Bajo las condiciones expuestas en el apartado 2.3.1 se obtienen el valor de la velocidad de flujo, las líneas de flujo, caída de presión y concentración de sustancias químicas en todos los puntos del dispositivo. Los resultados muestran como la velocidad aumenta en las zonas de menor sección de paso, viéndose afectada en gran parte la parte central de la zona de gel donde la velocidad aumenta en la región existente entre los postes. El las zonas situadas tras los postes la velocidad decrece debido al “efecto sombra” ejercido por los postes de sujeción. Es en la zona situada en la parte central del gel en la que más debe centrarse el análisis, ya que de los valores aquí obtenidos dependerán los resultados a obtener de la simulación del proceso de migración. La velocidad máxima obtenida es 2.93μm/seg , un valor muy similar al obtenido en [2] (3μm/seg) tanto experimentalmente como en una simulación 2D realizada con las mismas condiciones y geometría. Las líneas de corriente en la región central, presentan una orientación prácticamente alineada con la dirección X, como se observa en la (Figura. 12A). El hecho de llevar a cabo una simulación tridimensional del dispositivo posibilita el estudio de la variación de la velocidad y líneas de corriente con la profundidad. Como se observa en los resultados el efecto de los límites superiores e inferiores, es claramente reconocible por una menor velocidad de flujo en la cercanía del contorno (hasta hacerse nula en las paredes) y mayor velocidad en la parte central de las secciones de paso (Figura. 12B). De nuevo los resultados de mayor interés son los pertenecientes a la región central de la zona de gel, en donde se estudiará la migración celular utilizando el modelo desarrollado. 3. Resultados 34 A B Figura. 12: Simulación fluídica del dispositivo. Se simula el dispositivo bajo la influencia de un flujo intersticial estacionario obteniéndose (A) el perfil de velocidades en todo el dispositivo y líneas de flujo y perfil de velocidades en la zona de gel, donde se simula la migración celular, (B) diferentes planos transversales al dispositivo en los que se muestran la variación de la velocidad en esta dirección y evidencian la tridimensionalidad de la simulación. Con respecto a la presión se observa que la mayor caída se produce en el medio poroso (gel), mientras que la pérdida de carga en los canales en estado estacionario es prácticamente nula. La caída de presión en la región central de la zona de gel adquiere un valor prácticamente lineal con X, más concretamente cae desde 25Pa hasta adquirir un valor de 15Pa (Figura. 13). De igual modo la concentración química en estado estacionario presenta un perfil eminentemente lineal, como se muestra en el apartado siguiente. 3.1 Simulación del microdispositivo 35 Figura. 13: Caída de presión en el dispositivo microfluídico. La caída de presión producida en los canales por pérdidas de carga es totalmente despreciable respecto a la producida en la zona de gel (material poroso). En la región central de la zona de gel, la caída de presión adquiere un valor lineal con la coordenada x. 3.2 Simulación fluido-química a nivel celular Una vez resuelto el problema global, para valorar el efecto que la presencia de una célula tiene sobre el flujo intersticial y viceversa, y sobre la posible presencia de un factor químico autocrino secretado por la célula, se lleva a cabo una simulación fluido-química a nivel celular. Partiendo de una célula generada en una simulación anterior y elegida por su morfología polarizada, se malla su superficie en función de los vóxels que la constituyen utilizando un script de Matlab combinado con el software VMD. Posteriormente se lleva a cabo un renderizado de la malla mediante el Software Rhinoceros (Figura. 14). De esta forma se obtiene un volumen que puede ser importado por COMSOL para someterse a la simulación. Para ello se considera la célula embebida en un volumen control cúbico en cuyas fronteras se imponen condiciones de contorno obtenidas en la simulación del dispositivo microfluídico (Figura. 15). A B C Figura. 14: Evolución de la célula para la simulación fluido-química a nivel celular. La célula se discretiza en el modelo como (A) un conjunto de vóxels que aparecen y desaparecen durante el transcurso de la simulación de acuerdo con las funciones de probabilidad. (B) Mediante VMD (software de dinámica molecular) se realiza una representación aproximada de la membrana y su núcleo y (C) mediante Rhinoceros se adapta la malla para para simular en COMSOL la interacción entre la célula, el fluido y un factor autocrino segregado por la propia célula. 3. Resultados 36 Figura. 15: Volumen de control considerado par simulación fluido-química a nivel celular. Concretamente se establecen los perfiles de velocidad en la cara frontal y trasera del dominio, considerando en el resto de caras condiciones periódicas. Así mismo, se establece que la velocidad en la superficie celular es 0, y que dicha superficie mantiene de forma constante una concentración de factor autocrino. El resto de condiciones y parámetros utilizados para llevar a cabo la simulación se recogen el la siguiente tabla: Condición Descripción Valor Unidades i v Velocidad entrada 2.96e -6 [m/s] o v Velocidad salida 2.96e -6 [m/s] ,cell surf v Velocidad superficie celular 0 [m/s] D Coeficiente de difusión 10 -9 [m²/s] i c Concentración entrada 1 [mol/m 3 ] o c Concentración salida 0 [mol/m 3 ] cell c Concentración superficie celular 1 [mol/m 3 ] Tabla 3: Condiciones y parámetros en la simulación fluido-química a nivel celular. A la vista de los resultados obtenidos, se observa como la presencia de la célula apenas altera la distribución de las líneas de corriente y que al alejarse ligeramente de la superficie celular su efecto es despreciable incluso en un dominio de estudio tan pequeño, siendo las líneas prácticamente paralelas a la dirección del eje X (Figura. 16A). El factor químico autocrino previamente comentado, es segregado a través de la superficie celular de forma constante. Como puede observarse en la Figura. 16B, el factor químico se ve arrastrado por la corriente a la que se ve sometida la célula, provocando una concentración media mayor aguas abajo de la célula. 3.2 Simulación fluido-química a nivel celular 37 A B Figura. 16: Líneas de corriente y factor autocrino en la simulación fluido-química a nivel celular. (A) La interacción sólido-célula provoca que las líneas de flujo se vean alteradas por su presencia, sin embargo, su efecto es tan pequeño que puede considerarse despreciable a nivel del dispositivo completo. (B) En el caso de que la célula segregue un factor químico, la presencia de flujo intersticial provoca que este se vea arrastrado en el sentido de la corriente, creando una acumulación de concentración mayor aguas abajo. Al igual que la célula ofrece una ligera influencia sobre las líneas de corriente, el fluido interacciona con la superficie celular generando tensiones en función de la exposición de la célula al flujo. La zona que recibe mayor presión es aquella en que la célula es atacada por el flujo, viéndose reducido su valor al avanzar en la dirección del flujo (Figura. 17). Figura. 17: Distribución de presiones causadas por el flujo intersticial en la superficie celular. Las presiones producidas por el flujo sobre la superficie celular son mayores en la zona de la célula qué es atacada directamente por el flujo. Por otro lado, se lleva a cabo una simulación de un factor químico presente en el medio, cuyos resultados se muestran en la Figura. 18. Como se puede observar la presencia de la célula no ejerce apenas influencia sobre la difusión del factor químico y la concentración de este varía únicamente linealmente con la dirección X. 3. Resultados 38 Figura. 18: Simulación de difusión de un factor químico a nivel celular. Al llevarse a cabo la simulación de difusión de un factor paracrino con carácter estacionario, se obtiene una distribución lineal con la dirección X del gradiente de concentración, no ejerciendo la presencia de la célula apenas influencia. Estos resultados, en concordancia con otras simulaciones y experimentos recientes [31], confirman las hipótesis utilizadas en la definición de las funciones de probabilidad, aunque como se ha comentado anteriormente, se ha simplificado por el exceso de esfuerzo computacional que supondría una simulación de estas características en cada incremento temporal. 3.3 Simulación mecánica Tras llevar a cabo el análisis y la validación del modelo mecánico (Anexo IV) se está en condiciones de realizar una simulación completa de la migración considerando el comportamiento mecánico de la célula dentro del dispositivo de ensayo. Para esta simulación, se considera que la célula migra exclusivamente debido a la interacción mecánica célula-matriz (mechano-sensing) presentada en el apartado 2.2.1 . Por lo tanto, tanto los factores fluidoquímicos como el aleatorio no son considerados en este análisis. En este estudio se han considerado geles de igual geometría pero con diferentes rigideces para poder comparar los resultados obtenidos. En primer lugar se utiliza una matriz cuyo módulo elástico varía linealmente con la coordenada X, (valor 0,001kPa para valor 0 de abcisa y 100kPa en el extremo contrario) (Figura. 19). En segundo lugar se utiliza el mismo tipo de matriz, lineal en x con valor máximo de rigidez 100kPa, pero se considera la matriz de un material degradable. Por último se emplea una matriz, de material no degradable, cuyo módulo elástico toma un valor bajo (0.1kPa) constante al inicio, aumenta exponencialmente en la zona media y se mantiene en un valor alto constante en la zona final (100kPa) (Figura. 19). Los cálculos realizados simulan un tiempo de 500 minutos (8 horas y 20 minutos) implementadas en 100 incrementos temporales de 5 minutos. 3.3 Simulación mecánica 39 Figura. 19: Gradiente de rigidez del gel (ECM) a lo largo de la dirección X. Para cada una de las matrices en las que se representa la variación de su módulo elástico en función de la coordenada x, así como el valor de la rigidez en el punto donde se posiciona la célula previo a la simulación A continuación se exponen, para los tres casos considerados, los resultados obtenidos en relación a: la trayectoria (posición del centroide en cada incremento temporal), el camino recorrido, la velocidad de migración, el factor de forma, el área de adhesión y el mapa de tensiones y desplazamientos. Se analizan por una parte los resultados al comienzo de la simulación y por otra los resultados finales. Con este análisis se pretende remarcar el efecto inmediato de un cambio de rigidez tan brusco (exponencial) y cómo sus efectos se suavizan una vez se alcanza la saturación. Tras los primeros 16 incrementos (punto donde se alcanza dicha saturación) se calcula: la distancia total recorrida (teniendo en cuenta todos los movimientos), la distancia efectiva (distancia entre posición final e inicial), la velocidad media (utilizando la distancia total) y la velocidad efectiva (distancia efectiva recorrida en el tiempo considerado). En todos los casos la distancia total recorrida es muy similar [25μm (Figura. 20A), siendo ligeramente mayor para el caso que considera degradación de la matriz, 27.2μm. Sin embargo, la distancia efectiva es claramente mayor en el caso con gradiente exponencial (20.9μm frente a 6.9μm y 3.4μm). Así mismo, se observa como, a pesar de que la velocidad media adquiere un valor muy similar para los tres casos considerados ( [0.32μm/min), la velocidad efectiva para las matrices de gradiente lineal es bastante menor que la obtenida para la matriz con gradiente de rigidez exponencial (Figura. 20C), concretamente 0.04 y 0.08μm/min frente a 0.25μm/min. Para las matrices de gradiente lineal, a pesar de que la dirección de migración sigue preferentemente la dirección del gradiente (eje X), tiene un marcado carácter aleatorio y por tanto el avance efectivo es mucho menor. Este comportamiento se explica debido a que las dimensiones de la célula son tan pequeñas que apenas existe diferencia entre los valores de rigidez entre el frente y la parte trasera de la célula. Además, para estos valores de rigidez la fuerza ejercida por la célula ha alcanzado casi su valor de saturación, por lo que para variaciones pequeñas de la rigidez, la fuerza ejercida por la célula casi no varía. Por lo tanto, aunque las probabilidades de aparición son altas debido a las altas tensiones, el direccionamiento no está claramente definido y por lo tanto el avance es muy lento. Por otro lado, para la matriz con gradiente de rigidez exponencial la diferencia de valores entre la rigidez de la ECM en la parte frontal y la posterior son bastante pronunciados durante los 3. Resultados 46 A B C D Figura. 27: Resultados simulación basada en estímulos mecánicos y químicos. En todas las figuras se representan al mismo tiempo los resultados correspondientes a los casos (i) y (ii) de simulación, de forma que se puedan comparar los valores obtenidos para uno y otro caso. A pesar de migrar la célula en diferentes direcciones por los estímulos considerado en cada uno de los ejes coordenados, la distancia recorrida (A), la velocidad de migración (B), el área de proyección (C) y el factor de forma (D) adquieren valores similares. Los valores correspondientes al factor de forma en los últimos incrementos de la simulación son debidos a que la célula ha alcanzado los límites de la región de gel elegida para llevar a cabo la simulación mecánica y al no ser posible el avance en la dirección X, el crecimiento se produce en el plano XY. El área de la célula (Figura. 27C), como ya se comenta anteriormente en el apartado 1.1, aumenta en los primeros incrementos temporales adquiriendo tras ellos un valor medio prácticamente constante y acotado de aproximadamente 13000μm² para ambos casos. El factor de forma de la célula (Figura. 27D) adquiere valores alternos debido al carácter cambiante de la morfología celular, comprendidos en su gran mayoría entre 1 y 2. Sin embargo, para el caso (ii) se observa como al final de la simulación el factor de forma aumenta hasta valores superiores. Este resultado se explica por el hecho de que debido a las grandes velocidades obtenidas, para este momento la célula ha alcanzado los límites del fragmento de gel considerado para la simulación mecánica. Al no poder avanzar hacia adelante, la célula se extiende en el plano XY, aumentando así su factor de forma. Aunque esto se debe dejar fuera del análisis por tratarse de un caso extremo, es interesante observar lo que ocurriría al llegar a una superficie plana de material infinitamente rígido para caracterizar la respuesta de la célula ante tal situación (por ejemplo un pilar del microdispositivo). Al no poder deformar el sustrato, la célula ejerce una tensión alta, sobre todo en los puntos en contacto con la superficie (Figura. 3.4 Simulación completa 47 28). Además, al no poder añadir elementos por la parte delantera, comenzaría a extenderse lateralmente, hasta alcanzar con el suficiente tiempo una forma de disco adherida a la pared, lo cual ha sido observado en experimentos. Figura. 28: Célula adherida al extremo del dominio. Debido a las altas velocidades obtenidas en el caso(ii), la célula alcanza el contorno del gel (cuyos desplazamientos están restringidos). Al no poder deformar la pared, la célula queda adherida, ejerciendo altas tensiones y expandiéndose lateralmente. 48 49 4. Discusión De entre todos los tipos de modelos de migración celular existentes, no existe hasta la fecha ningún modelo de simulación que permita al mismo tiempo ofrecer una perspectiva general del proceso de migración (trayectoria, velocidad) y por otro lado obtener información cuantitativa acerca de aspectos relevantes a nivel celular como el área de adhesión, el factor de forma o las tensiones generadas. Por lo tanto, con este modelo se pretende cubrir en cierto modo esta carencia, planteándolo de tal manera que pueda ser fácilmente manipulado y perfeccionado en el futuro. Este modelo simula la migración celular mediante la aparición y desaparición de vóxels, que en el transcurso del tiempo van actualizando la forma celular. La aparición y desaparición de vóxels ocurre de acuerdo a unas funciones de probabilidad que asignan a cada elemento un valor en base a varios estímulos. De los múltiples posibles estímulos que condicionan la migración celular, únicamente algunos de ellos han sido incorporados, concretamente de carácter mecánico, químico y de flujo. En primer lugar, siguiendo el orden seguido para llevar a cabo la simulación, se discuten los resultados obtenidos de la simulación del microdispositivo, así como de la simulación fluido-química a nivel celular. Más adelante se discutirán los resultados obtenidos de la simulación mecánica y por último los resultantes de la simulación completa (teniendo en cuenta todos los factores). Los resultados obtenidos de la simulación del microdispositivo en estado estacionario, muestran el perfil de velocidades y la presión en todos los puntos. Como se observa el la Figura. 12, los valores de velocidad obtenidos para la zona de gel se sitúan en el rango comprendido entre 1 y 1.5μm/seg, valores consistentes con los utilizados en ensayos sobre la influencia del flujo intersticial en la migración celular [31]. Por otra parte si se observa el valor de la velocidad se ve considerablemente influida por la geometría del dispositivo, observándose los valores máximos de velocidad allí donde las secciones de paso son menores o en las partes, en las que debido a la geometría se producen cambios en las direcciones del flujo. Referente a la caída de presión, según se muestra en los resultados y acorde con las condiciones de contorno adoptadas (40Pa a la entrada y 0 a la salida), se produce casi en su totalidad en la región de gel. Estos resultados, se tomarán como datos de entrada para la simulación de migración completa. Por otra parte, con los resultados obtenidos de la simulación fluido-química a nivel celular se valida una de las simplificaciones más importantes asumidas en el modelo. Como se observa en la Figura. 16A, la dirección de las líneas de corriente apenas se ven alteradas por la presencia de la célula aun considerando las dimensiones minúsculas del dominio, por lo tanto, la aproximación de considerar las líneas de corriente constantes durante toda la simulación (o asumir estado estacionario para el flujo) es válida. De igual modo, se valida la simplificación aplicada a las funciones de probabilidad según la activación/inhibición de receptores en concordancia con [31]. En la Figura. 16B se muestra como el factor autocrino es arrastrado por la corriente de flujo intersticial, lo cual haría migrar a la célula aguas abajo. También se observa 4. Discusión 50 como las tensiones en la parte delantera de la célula son mayores (Figura. 17) al estar su superficie más expuesta al flujo. Inhibiendo los receptores del factor autocrino, la diferencia de presiones explicaría la tendencia de las células a migrar contracorriente[2]. El proceso de mechano-sensing ha sido heredado de un modelo previo [1], el cual determina que la fuerza ejercida por la célula crece con la rigidez hasta alcanzar la saturación, lo cual ha sido demostrado experimentalmente [32, 33]. En el apartado 1.1 se ha comprobado como esto también se cumple en el modelo actual. La Figura. 22 muestra como la tensión de la célula es mayor para mayor rigidez, y como los desplazamientos sufridos por la matriz disminuyen cuanto mayor es el módulo elástico. Además, de forma consistente a lo obtenido mediante experimentos [34], la velocidad de migración celular también aumenta conforme aumenta la rigidez (Figura. 20) hasta alcanzar un valor constante para valores elevados de la rigidez (para los cuales la fuerzas ejercida por la célula están próximas a la saturación,). Los resultados obtenidos para la velocidad de migración celular media y efectiva para matrices con gradientes de rigidez lineales y exponencial presentan diferencias significativas. Si se centra el análisis en los primeros 16 incrementos temporales se observa como la velocidad efectiva para gradiente de rigidez exponencial (0.25μm/min) es muy superior a la obtenida para gradiente de rigidez lineal (0.8μm/min en el caso máximo). Esto se explica debido a que la variación en los valores de rigidez entre la región posterior y anterior de la célula es mucho mayor que en el caso lineal, y por tanto la migración tiene un marcado carácter direccional en la dirección de este gradiente. Estas diferencias se reducen conforme transcurre el tiempo de simulación, ya que en las matrices de gradiente de rigidez lineal, la célula va migrando progresivamente hacia posiciones más rígidas y por tanto presentará mayor velocidad de migración. En cualquier caso, tanto los valores de velocidad como de longitud recorrida por las células al migrar son consistentes tanto con los valores obtenidos en la simulaciones de validación del modelo mecánico como el químico. Además, estos valores se encuentran dentro de los rangos experimentales (ej.: la velocidad de los fibroblastos oscila alrededor de 0.2-0.7μ/min) [11, 35]. Es preciso destacar sin embargo, que existe una gran disparidad de resultados en función del ensayo y del tipo de célula considerado. En referencia al área proyectada, los resultados obtenidos por el modelo (Figura. 21 y Figura. 27), son consistentes con los obtenidos en las simulaciones de validación (Figura. 37) pues para todos los casos se obtienen valores dentro del rango. Aunque bien es cierto que este parámetro está condicionado por la restricción de tamaño impuesta a la célula, es interesante ver como el valor medio se mantiene para todas las simulaciones. Sin embargo, los valores obtenidos son hasta diez veces superiores a los logrados mediante ensayos experimentales [36]. Esta diferencia es fundamentalmente debida a la diferencia en los métodos de medida utilizados en los ensayos y en la computación. En los ensayos se suele determinar el área de adhesión determinando el contorno a partir de imágenes 2D y realizando estimaciones. Otras técnicas estiman el área de contacto celular a partir de tintes fluorescentes. Estos tintes se vuelven visibles únicamente en zonas de la célula activas (ej. Adhesiones focales), por lo que la medida predice el área sobre la que la célula ejerce fuerzas. En esta simulación, sin embargo, se considera toda el área de contacto célula-matriz, cuyo valor se calcula como la suma de 4. Discusión 51 todas las caras libres de los elementos del contorno. Por otro lado, al discretizar el volumen con vóxels, se sobrestima el área que tendría un volumen suavizado de dimensiones similares. Comparando los resultados de velocidad media y efectiva obtenidos en la simulación exclusivamente mecánica con los obtenidos en la simulación completa, se observa como la diferencia entre las velocidades media y efectiva es mucho mayor en el caso mecánico. Para los casos con gradiente de rigidez lineal, los valores de velocidad media y efectiva son [0.4μm/min y 0.05μm/min respectivamente, sin embargo, para los casos que incluyen además estímulos químicos y flujo, la velocidad media es de [0.37μm/min y la efectiva [0.26μm/min. Este resultado evidencia que la migración debida a estímulos fluidoquímicos presenta un marcado carácter direccional, mientras que la debida exclusivamente a estímulos mecánicos es más estocástica. Esta aleatoriedad, viene en parte dada por las pequeñas diferencias de rigidez en la ECM (salvo en casos extremos como el exponencial), y por otro lado por el modo en el que el gradiente de tensión es computado. Cuando dos o más elementos presentan tensiones altas muy similares (algo habitual en casos normales), el elemento que determinará la dirección a seguir se escoge al azar. Considerando la complejidad del proceso, a la hora de desarrollar el modelo ha sido necesario asumir ciertas simplificaciones. Por ejemplo, para resolver el problema mecánico que da como resultado la tensión en cada elemento de la célula se considera que la matriz de “Colágeno I” tiene un comportamiento elástico lineal, cuando en realidad se trata de un material visco-elástico. Aunque se considera que esta aproximación es suficientemente válida, sería interesante incorporar este comportamiento si se buscasen unos resultados detallados y más realistas de la interacción célula-matriz. Por otra parte, siguiendo con los estímulos mecánicos, solo se consideran las fuerzas generadas por el comportamiento activo de la célula, dejándose a un lado por ejemplo la fuerza de arrastre, de vital importancia cuando se considera el comportamiento viscoso del material que constituye la ECM. En este modelo, la velocidad de migración siempre aumenta con la rigidez (mayor rigidez implica mayor tensión y por lo tanto mayores probabilidades de aparición/desaparición). Sin embargo, para valores muy altos de rigidez la velocidad de migración experimenta en realidad una reducción al aumentar arrástrela resistencia viscosa al avance sobre la superficie celular [1]. Aunque por simplicidad esto no se ha tenido en cuenta, sería fácilmente incorporable al modelo. Bastaría con añadir un parámetro con peso negativo y proporcional a la rigidez a las funciones de probabilidad, de manera que su efecto fuese bajo para sustratos blandos y alto para sustratos rígidos, consiguiendo así un óptimo de velocidad a rigideces intermedias[1, 37]. En cuanto a los estímulos fluido-químicos, en este modelo se incluyen en las funciones de probabilidad de una forma sencilla debido al alto coste computacional que supondría el realizar los cálculos en cada incremento temporal. Como se explica anteriormente, las células en presencia de flujo intersticial pueden segregar un factor químico autocrino que es arrastrado por el flujo y en última instancia provoca que la célula migre en la dirección de la corriente. Sin embargo, si se bloquean ciertos receptores, las células dejan de verse afectadas por este factor químico, migrando en dirección contraria. Según se recoge en [2, 31], la presión ejercida por el flujo sobre la célula es la causante de este comportamiento: en ausencia de 4. Discusión 52 estímulos químicos, la célula migra en dirección al gradiente de presión que siente sobre su cuerpo, es decir, contra la dirección del flujo que es el causante de esta presión. Con el fin de observar los efectos del flujo intersticial, tanto en la presión ejercida sobre la célula como sobre el arrastre del factor autocrino, se ha realizado un estudio en detalle de la interacción fluido-célula. Los resultados de tal simulación confirman como la parte de la célula que sufre las mayores tensiones por efecto del fluido es la más expuesta al flujo. Los resultados obtenidos en la validación del modelo fluido-químico (Anexo II) muestran que, los estímulos fluido-químicos, pese a ser implementados de una forma sencilla, presentan tendencias similares a las obtenidas en [2]. Por otra parte, la implementación de la migración celular bajo la influencia de un flujo intersticial tiene únicamente en cuenta la dirección de dicho flujo, pero no la magnitud de la velocidad. Esto es importante para determinar el arrastre del factor autocrino, y su efecto ha sido estudiado experimentalmente, Además, las líneas de corriente se consideran constantes durante todo el transcurso de la simulación (estado estacionario) sin considerarse la influencia que pudiera tener la forma celular. Aunque en cierto modo realizar una simulación fluidoquímica en cada incremento temporal arrojaría resultados mucho más precisos, como se ha expresado anteriormente su coste computacional sería demasiado elevado y por tanto se prefiere simplificar el modelo. La degradación del gel por parte de la célula se ha implementado mediante la reducción del módulo elástico del material que constituye la matriz. Los resultados obtenidos no muestran variación aparente con respecto a los casos en los que la degradación no está activada. Los estudios muestran sin embargo, que para ciertos tipos de matrices, el proceso de degradación es de vital importancia en el proceso de migración [6], dejando las células a su paso un canal en la matriz, por el que otras células migran con mayor facilidad. Para nuestro caso, al estudiarse la migración de una célula, este fenómeno no tiene demasiado valor. A pesar de ello, se trata de un fenómeno altamente complejo que todavía se está dilucidando en la actualidad, por lo que merece la pena tenerlo en cuenta para versiones futuras. Para ello sería necesario incorporar un modelo completo dedicado al proceso de degradación. El modelo desarrollado presenta la complejidad de que la velocidad, distancia y sentido de la migración no se determinan por medio de la resolución numérica de una serie de ecuaciones, sino que está basado en funciones de probabilidad que determinan la aparición/desaparición de vóxels. A pesar de que esto permite la activación/desactivación o incluso la incorporación de nuevos factores de manera sencilla, también conlleva un complejo proceso de ajuste. Por lo tanto, al existir esta dependencia del azar, son necesarias numerosas pruebas hasta poder concluir que los resultados son consistentes. 53 5. Conclusiones y Trabajo Futuro Como se ha venido haciendo hincapié a lo largo de toda la memoria, la migración celular es un proceso complejo en el que influyen numerosos factores de diferente naturaleza (estímulos mecánicos, químicos, de flujo, topográficos etc.) y su comprensión es de vital importancia para el avance en ingeniería tisular y en el tratamiento de numerosas enfermedades. En este trabajo, se ha desarrollado un modelo computacional 3D basado en probabilidades capaz de simular la migración de una célula en función de factores mecánicos, químicos y de flujo, creando así una herramienta que permite analizar de una vez multitud de características como la interacción célula-matriz (tensión y deformación generada, degradación), la velocidad de migración, el camino recorrido, la forma adquirida por la célula etc. cuyo estudio sería de otra forma muy costoso, o requeriría de un gran número de experimentos. Para su desarrollo se ha utilizado software diverso y se ha requerido el uso de un clúster de computación. Al estar el modelo basado en probabilidades, y por tanto depender hasta cierto punto del azar, presenta una fuerte dependencia de los parámetros utilizados. Este, aparentemente, hándicap es sin embargo una de las mayores virtudes del mismo. Y es que, sin establecer directamente una ley de comportamiento, el modelo es capaz de capturar el comportamiento emergente a nivel macro de una célula migrando en un ambiente 3D. En este punto reside el verdadero potencial de esta aproximación, y en su mejora debe centrarse el desarrollo futuro. Con este trabajo se ha completado la primera versión de un modelo cuyo potencial se extenderá en líneas de trabajo futuras como: - Introducción de nuevos factores influyentes en la migración (ej. estímulos eléctricos o presión): nuevos avances y experimentos son realizados día a día, por lo que si con el modelo se pretende ofrecer una plataforma de ayuda, debe estar en permanente actualización. - Modelado en detalle de factores simplificados (ej. degradación ECM o química de la célula): las simplificaciones consideradas en la primera versión del modelo, deben revisarse en función de la orientación que se le dé en un futuro. - Mejora del modelo de mechano-sensing introduciendo una dependencia temporal en función de la activación de los motores moleculares: esta línea de trabajo, actualmente ya en funcionamiento, pretende mejorar el modelo teórico de mechanosensing, introduciendo un nuevo componente al esquema celular que simule la activación paulatina de los motores moleculares (encargados de realizar la contracción). Con ello, el modelo será capaz de estudiar por ejemplo, la evolución temporal de la tensión antes de que la célula migre de su posición. - Simulado específico de adhesiones focales en ciertos puntos del contorno celular: la formación de adhesiones focales es un tema crucial para comprender la fuerza 5. Conclusiones y Trabajo Futuro 54 generada por las células. Debido a su fuerte componente bioquímico, se trata de un problema de gran complejidad, por lo que el seguimiento de esta línea conllevaría una orientación del modelo muy focalizada a este fenómeno. - Incorporación de más células a la simulación: aunque la simulación de células individuales es de gran interés y sirve para el estudio de factores de manera aislada, la incorporación de poblaciones es un paso necesario y obligatorio si se pretende introducir el modelo en el campo de regeneración de tejidos. Además de los cambios pertinentes en la programación, sería necesario el desarrollo y la incorporación a los presentes códigos de un modelo de interacción célula-célula. 55 6. Referencias 1. Borau, C., R.D. Kamm, and J.M. García-Aznar, Mechano-sensing and cell migration: a 3D model approach. Physical Biology, 2011. 8(6): p. 066008. 2. Polacheck, W.J., J.L. Charest, and R.D. Kamm, Interstitial flow influences direction of tumor cell migration through competing mechanisms. Proceedings of the National Academy of Sciences of the United States of America, 2011. 108(27): p. 11115-11120. 3. Lauffenburger, D.A. and A.F. Horwitz, Cell migration: A physically integrated molecular process. Cell, 1996. 84(3): p. 359-369. 4. Friedl, P. and K. Wolf, Tumour-cell invasion and migration: diversity and escape mechanisms. Nat Rev Cancer, 2003. 3(5): p. 362-74. 5. Yamaguchi, H., J. Wyckoff, and J. Condeelis, Cell migration in tumors. Current Opinion In Cell Biology, 2005. 17(5): p. 559-564. 6. Ehrbar, M., et al., Elucidating the Role of Matrix Stiffness in 3D Cell Migration and Remodeling. Biophysical Journal, 2011. 100(2): p. 284-293. 7. Tabata, Y., Biomaterial technology for tissue engineering applications. J R Soc Interface, 2009. 6 Suppl 3: p. S311-24. 8. Moreo, P., J.M. Garcia-Aznar, and M. Doblare, Modeling mechanosensing and its effect on the migration and proliferation of adherent cells. Acta Biomaterialia, 2008. 4(3): p. 613-621. 9. Simpson, M.J., et al., Migration of breast cancer cells: Understanding the roles of volume exclusion and cell-to-cell adhesion. Physical Review E, 2010. 82(4): p. 041901. 10. Fraley, S.I., et al., A distinctive role for focal adhesion proteins in three-dimensional cell motility. Nat Cell Biol, 2010. 12(6): p. 598-604. 11. Hakkinen, K.M., et al., Direct Comparisons of the Morphology, Migration, Cell Adhesions, and Actin Cytoskeleton of Fibroblasts in Four Different Three-Dimensional Extracellular Matrices. Tissue Engineering Part A, 2011. 17(5-6): p. 713-724. 12. Chung, S., et al., Microfluidic Platforms for Studies of Angiogenesis, Cell Migration, and Cell-Cell Interactions. Annals of Biomedical Engineering, 2010. 38(3): p. 1164-1177. 13. Aznavoorian, S., et al., Signal transduction for chemotaxis and haptotaxis by matrix molecules in tumor cells. J Cell Biol, 1990. 110(4): p. 1427-38. 14. Sanz-Herrera, J.A., et al., On the effect of substrate curvature on cell mechanics. Biomaterials, 2009. 30(34): p. 6674-86. 15. Bischofs, I.B. and U.S. Schwarz, Cell organization in soft media due to active mechanosensing. Proceedings of the National Academy of Sciences of the United States of America, 2003. 100(16): p. 9274-9279.