Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques

1 2 3 ESCUELA TÉCNICA SUPERIOR DE INGENIERÍA INFORMÁTICA GRADO EN INGENIERÍA DE LA SALUD MODELADO DE CIRCUITOS NEURONALES EMPLEANDO TÉCNICAS DE PROGRAMACIÓN ORIENTADA A BLOQUES MODELLING OF NEURONAL CIRCUITS BY USING BLOCK-ORIENTED PROGRAMMING TECHNIQUES Realizado por Julia Perea Paizal Tutorizado por Francisco Javier Fernández de Cañete Rodríguez Departamento Ingeniería de Sistemas y Automática UNIVERSIDAD DE MÁLAGA MÁLAGA, Julio 2015 Fecha defensa: El Secretario del Tribunal 4 5 Resumen: El presente trabajo consiste en la implementación de un modelo del circuito completo del arco reflejo basado en la neurona biológica y en el estudio de su comportamiento dinámico. El circuito reflejo modular básico está compuesto por un receptor y su fibra aferente, una o más neuronas medulares intercalares, una motoneurona y un efector y su músculo. Las propiedades de cualquier reflejo son una función de las características de los receptores que inician el estímulo y de la naturaleza de los circuitos medulares y de las motoneuronas implicados. Para la implementación del circuito neuronal se ha optado por técnicas de modelado causal, en las que tienen más importancia las ecuaciones dinámicas que la estructura física del sistema. Para ello se ha utilizado la herramienta Simulink, que permite implementar la configuración utilizada basada en modelos compartimentales mediante la utilización de una estructura jerarquizada de diagramas de bloques. El modelado del circuito neuronal se realiza siguiendo un modelo matemático en términos de potenciales de membrana y potenciales de acción, implementando los cuatro elementos principales que componen una neurona: cuerpo celular o soma, dendritas, axón y terminales presinápticos, además del huso muscular y la placa motora, que actúan como transductores. Las ecuaciones matemáticas que definen el sistema se obtienen a partir de la analogía circuital que reproduce el tránsito de iones a través de la membrana. Construido el modelo se estudian los mecanismos de excitación, conducción y transmisión de estímulos en las distintas etapas que lo componen, el comportamiento en condiciones fisiológicas y patológicas, y los efectos que se producen tras la aplicación de tratamientos farmacológicos siguiendo una metodología de análisis cualitativa. Palabras claves: modelado de circuitos neuronales, modelado causal, técnicas de programación orientado a bloques, arco reflejo, modelo de Hodgking-Huxley. Abstract: This project consists in the implementation of a complete reflex arc model based in biological neurons and in the subsequent study of its dynamic behavior. The basic circuit of the reflex arc is composed by a receptor and its afferent fibers, one or several interneurons located in the spinal cord, a motoneuron and an effector connected to its muscle. The attributes of a reflex arc depend on the features of the sensory receptor that initiates the reflex and also on the physical nature of the medullar circuits and the motoneurons involved in it. The implementation of the neural circuit has been done using causal modeling techniques, which reflect rather the dynamics equations than the very structure of the system. The Simulink environment has been used for this purpose since it is based in a hierarchically block diagramming tool, which makes very efficient the implementation of compartmental models. The model of the neural circuit has been built following a mathematical approach based in membrane and action potentials. It includes the four major 6 components that embody a neuron: cell body, dendrite, axon and presynaptic terminals, as well as the muscle spindle and the motor junction, which behave as transductors. The mathematical equations that model the system are calculated from the circuit analogy that recreates the ion traffic across cellular membranes. Once the model is completely built, the excitation and conduction mechanisms, as well as the stimulus transmission, are studied. The system is analyzed under physiological and pathological conditions, studying the effects that some drugs produce in the model and using a qualitative analysis methodology to evaluate the obtained results. Keywords: neural circuit modelling, causal modelling, block-oriented programming techniques, reflex arc, Hodgking-Huxley model. 7 Índice INTRODUCCIÓN ............................................................................................................................. 9 DESARROLLO, PRUEBAS Y RESULTADOS ..................................................................................... 12 1. REVISIÓN DE LA FISIOLOGÍA NEURONAL ............................................................................ 12 1.1. Introducción ................................................................................................................ 12 1.2. Neuronas ..................................................................................................................... 12 1.3. Potenciales de membrana ........................................................................................... 13 1.4. Receptores sensoriales o transductores ..................................................................... 15 1.5. Aferencias segmentales que originan actividad refleja .............................................. 15 1.6. Efectores y su control .................................................................................................. 17 1.7. Reflejos modulares ...................................................................................................... 18 1.7.1. El reflejo miotático. ............................................................................................. 18 2. MODELADO DE NEURONA, RECEPTOR Y PLACA MOTORA. ................................................ 20 2.1. Modelo dendrítico ....................................................................................................... 21 2.2. Modelo axónico ........................................................................................................... 23 2.2.1. Modelo de hodgkin-huxley para los potenciales de acción ................................ 24 2.3. Modelo del soma ......................................................................................................... 26 2.4. Modelo sináptico ......................................................................................................... 26 2.5. Modelo del receptor ................................................................................................... 29 2.6. Modelo de la placa motora ......................................................................................... 32 3. HERRAMIENTAS DE MODELADO ORIENTADO A BLOQUES ................................................. 34 3.1. Introducción ................................................................................................................ 34 3.2. Simulink ....................................................................................................................... 34 4. IMPLEMENTACIÓN DE LA NEURONA UTILIZANDO SIMULINK ............................................ 36 4.1. Modelo dendrítico. ........................................................................................................... 36 4.2. Modelo axónico. ............................................................................................................... 38 4.3. Modelo del soma .............................................................................................................. 41 4.4. Modelo sináptico. ............................................................................................................. 42 4.5. Modelo del receptor sensorial. ........................................................................................ 43 4.6. Modelo de la placa motora .............................................................................................. 45 4.7. Implementación del modelo completo ............................................................................ 47 8 5. ESTUDIO BAJO CONDICIONES FISIOLÓGICAS, PATOLÓGICAS Y TRATAMIENTO FARMACOLÓGICO ....................................................................................................................... 50 5.1. Estudio bajo condiciones fisiológicas. ......................................................................... 50 5.2. Estudio bajo condiciones patológicas y tratamiento farmacológico .......................... 55 5.2.1. Enfermedades de la motoneurona ..................................................................... 55 5.2.2. Esclerosis múltiple. .............................................................................................. 67 CONCLUSIONES ........................................................................................................................... 72 Referencias bibliográficas. .......................................................................................................... 74 Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 9 INTRODUCCIÓN Motivación: El presente trabajo surge debido a la dificultad que supone el estudio de sistemas biológicos a partir de métodos experimentales. Los sistemas biológicos son inherentemente complejos y en los experimentos in situ que se realizan para su estudio a menudo surgen efectos no deseados debidos, entre otros factores, a condiciones estocásticas. Esto hace que sea muy difícil comprender su funcionamiento si no es invirtiendo un elevado número de recursos materiales y temporales. La utilización de modelos matemáticos que simulan el comportamiento de un sistema biológico hace posible que se puedan captar las características inherentes que lo definen. De esta forma se puede estudiar la respuesta del sistema biológico ante distintas señales externas e internas, siendo posible modificar sus condiciones iniciales para poder simular los efectos que provocan en él determinadas patologías o la aplicación de tratamientos farmacológicos. Las ventajas mencionadas anteriormente son de especial importancia en el campo de la neurociencia, caracterizado por la dificultad que supone la toma de datos experimentales que permiten evaluarlo ante distintas situaciones. Por esta razón, se ha considerado de gran interés el modelado de un circuito neuronal con un nivel de profundidad suficiente para permitir simular los efectos de ciertas patologías y tratamientos farmacológicos, además de entender mejor su naturaleza en condiciones de normalidad. Debido a la complejidad que suponen los circuitos neuronales, en los que una neurona media puede tener hasta 10000 sinapsis, se ha decidido estudiar el circuito neuronal más sencillo existente en el ser humano: el arco reflejo. Pese a su aparente simpleza, esta vía nerviosa es una herramienta indispensable para la supervivencia al permitir obtener respuestas en un periodo de tiempo mucho más reducido. Objetivos: Los objetivos de este trabajo consisten en la realización de un modelo del circuito del arco reflejo utilizando técnicas de modelado orientado a bloques. Las ecuaciones matemáticas que definen la transmisión de información a lo largo del circuito se obtienen empleando una analogía circuital. Para la elaboración del trabajo se ha optado por técnicas de modelado causal, en la que tienen más importancia las ecuaciones dinámicas que la estructura física del sistema que se está modelando. El modelado del circuito neuronal se realizará siguiendo un modelo matemático en términos de potenciales de membrana y potenciales de acción, implementando los cuatro elementos principales que componen una neurona: cuerpo celular o soma, dendritas, axón y terminales presinápticos. Construido el modelo, se pretenden estudiar los mecanismos de excitación, conducción y transmisión de estímulos en las distintas etapas que lo componen, su estudio tanto en condiciones fisiológicas como Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 16 Los transductores en nuestro caso van a ser los receptores sensoriales, concretamente las fibras de los husos musculares. Los husos musculares se encuentran en casi todos los músculos esqueléticos, y están particularmente concentrados en los músculos que ejercen un control motor preciso. El huso neuromuscular es una estructura en forma de huso de varios milímetros de largo y pocas micras de ancho. Tiene una vaina de tejido conectivo que envuelve entre 2 y 12 fibras musculares intrafusales. Cada extremo del uso está adosado a fibras musculares extrafusales. Existen dos tipos de fibras intrafusales: las fibras de la bolsa nuclear y las fibras de la cadena nuclear. Un huso típico se halla inervado por una fibra aferente Ia única y una sola fibra aferente del grupo II. Dichas fibras se ramifican dentro del huso para originar las terminaciones primarias y secundarias. Las fibras intrafusales, además de poseer innervación sensitiva a través de las motoneuronas α, poseen inervación motora del sistema nervioso central a través de las motoneuronas 𝛾, por lo que están bajo control centrífugo. Su estructura se muestra en la figura 1.3. Figura 1.3. A: Colocación de un huso muscular en paralelo con las fibras musculares extrafusales y un órgano tendinoso en serie con las fibras musculares. B: Estructura e inervación sensorial y motora de un huso muscular. C: Estructura e inervación de un órgano tendinoso. Imagen obtenida de Berne y Levy (2009). Berne y Levy Fisiología. Barcelona: Elsevier Los husos musculares responden a cambios en la longitud del músculo debido a que se sitúan en paralelo a las fibras extrafusales y, por tanto, también pueden estirarse o acortarse junto con las fibras extrafusales. Debido a que las fibras intrafusales, como todas las fibras musculares, muestran propiedades de tipo resorte, un cambio en su longitud puede cambiar la tensión a la que están sometidas, y este cambio se detecta mediante los mecanorreceptores de las aferencias del huso Ia y II. Las fibras del grupo Ia son sensibles tanto al grado de estiramiento del músculo como a su frecuencia, mientras que las fibras del grupo II responden principalmente al grado de Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 17 estiramiento. Por lo tanto, las fibras Ia tienen un comportamiento estático y dinámico mientras que las del grupo II sólo tienen un comportamiento estático. La inervación eferente de los husos musculares es extremadamente importante, ya que determina la sensibilidad del huso muscular al estiramiento. Cuando la porción extrafusal del músculo se contrae, el huso muscular se descarga mediante el acortamiento del músculo y la aferencia del huso muscular puede detener su descarga volviéndose insensible a cambios posteriores en la longitud muscular. Esta descarga del huso puede ser contrarrestada si las motoneuronas 𝛾 son estimuladas simultáneamente, provocando el acortamiento de las fibras musculares intrafusales del huso, junto con las fibras musculares extrafusales. Por tanto, el sistema de la motoneurona 𝛾 permite operar al huso muscular en un amplio rango de longitudes del músculo mientras retiene una alta sensibilidad a pequeños cambios de longitud. En resumen, el huso neuromuscular es un órgano receptor que posee aferencias dobles y al menos una sola eferencia. Los receptores responden al estiramiento muscular y proveen información al sistema nervioso con respecto a la longitud de las fibras nerviosas (terminaciones primarias y secundarias) y la velocidad de cambio en la longitud (terminaciones primarias). La sensibilidad de estos receptores se puede ajustar por medio del sistema nervioso central por vía de las motoneuronas 𝛾, y las inervaciones fusimotoras estática y dinámica permiten que las descargas fusimotoras estática y dinámica tengan un control propio para cada una. En la figura 1.4 se sintetiza el funcionamiento de las distintas fibras ante diversos estímulos. Figura 1.4. Respuestas de una terminación primaria (Ia) y una terminación secundaria (II) a cambios de longitud del músculo. Imagen obtenida de Berne y Levy (2009). Berne y Levy Fisiología. Barcelona: Elsevier 1.6. EFECTORES Y SU CONTROL Las salidas del extremo motor del sistema nervioso se expresan en última instancia mediante la contracción del músculo estriado, la contracción del músculo liso visceral y ocular y la secreción glandular. Los músculos reciben su inervación de una sola fuente: las motoneuronas de la médula espinal. Cada fibra muscular de un músculo recibe Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 18 inervación de una sola motoneurona y una sola motoneurona puede inervar varias fibras contenidas dentro de un mismo músculo. La motoneurona y las fibras que ésta inerva se conocen con el nombre de unidad motora, siendo la unidad motora la unidad básica para el movimiento. Las motoneuronas producen sólo una influencia excitatoria sobre el músculo, por lo que la activación de una moteneurona resulta sólo en la contracción de las fibras que inerva. La relajación o el estiramiento de un músculo sólo puede lograrse mediante una reducción en la descarga de esa motoneurona. Por lo tanto, la tensión generada en una unidad motora es función de la velocidad de descarga de la motoneurona y la tensión generada en un músculo es el resultado de la suma de las tensiones de todas las unidades motoras activas de ese músculo [3]. La estructura de la unidad motora se muestra en la figura 1.5. Figura 1.5. Representación de la estructura de la unidad motora, compuesta por una motoneurona y todas las fibras musculares que inerva. Imagen obtenida de Berne y Levy (2009). Berne y Levy Fisiología. Barcelona: Elsevier 1.7. REFLEJOS MODULARES Un circuito reflejo modular básico está compuesto por un receptor y su fibra aferente, una o más neuronas medulares intercalares (exceptuando el reflejo miotático), una motoneurona y su axón y un efector con su músculo. Las propiedades de cualquier reflejo son una función de las características de los receptores que inician los reflejos, la naturaleza de los circuitos medulares implicados y la distribución anatómica de las motoneuronas que participan. 1.7.1. El reflejo miotático. Este reflejo sirve para mantener la postura y es mediado por la médula espinal y originado por la activación de las terminaciones primarias Ia del huso neuromuscular. El segmento aferente del reflejo consiste en las terminaciones primarias del huso Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 19 neuromuscular y sus fibras aferentes correspondientes. Los circuitos medulares constan de una conexión excitatoria monisináptica con las motoneuronas homónimas (motoneuronas que inervan el músculo del cual se originó la fibra Ia). Por lo tanto, el reflejo es local por restringirse sus efectos al músculo del cual la fibra aferente activadora se origina. La acción refleja tiene una latencia corta. El circuito miotático incluye una inhibición mediada disinápticamente (con dos neuronas) conocida como inervación recíproca y que se dirige hacia los músculos antagonistas, provocando su relajación al mismo tiempo que la contracción de los agonistas. La estructura del reflejo miotático, que es el que se va a implementar en el modelo dinámico realizado, se muestra en la figura 1.6. Figura 1.6. Reflejo del arco patelar y la inhibición recíproca del músculo antagonista. Imagen obtenida de Saladin K.S. Anatomy & Physiology. McGraw-Hill Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 20 2. MODELADO DE NEURONA, RECEPTOR Y PLACA MOTORA. Un modelo es una representación de la realidad que implica cierto grado de aproximación. Los modelos matemáticos de un sistema fisiológico representan procesos fisiológicos en forma de ecuaciones matemáticas. La forma en la que se diseña el modelo y el grado de detalle que éste tiene están determinados principalmente por la aplicación que se le va a dar. Desde una perspectiva científica, los cuatro objetivos básicos para los que los modelos son desarrollados corresponden a cuatro categorías básicas: modelos descriptivos, interpretativos, predictivos y explicativos. Los modelos aplicados al campo de la neurociencia aportan potentes medios para investigar diferentes características del sistema neuronal, como pueden ser encontrar los efectos que provocan un mal funcionamiento de los canales dependientes de voltaje en determinadas patologías o la inhibición o excitación de ciertos canales iónicos provocados por agentes farmacológicos. Además, la simulación del modelo antes de realizar cualquier experimento biológico permite tener unos conocimientos previos a los resultados de los experimentos, siendo más fácil su evaluación. Como ya se ha comentado en el apartado anterior, la neurona puede dividirse en cuatro partes principales: cuerpo celular o soma, dendritas, axón y terminales presinápticos. El cuerpo celular o soma contiene el núcleo y es similar a la estructura de otras células. Las dendritas actúan como la entrada sináptica de la neurona previa excitatoria o inhibitoria. Tras la estimulación de las neuronas y sus dendritas, la permeabilidad de la membrana plasmática de la célula al sodio se intensifica y un potencial de acción se genera desde la dendrita al axón. Para modelar la neurona completa es necesario modelar cada una de las partes de forma independiente basándose en su mecanismo eléctrico para después poder conectarlas en un único sistema. En la figura 2.1 se muestra el mecanismo eléctrico de una neurona y su circuito equivalente que permite la implementación del modelado neuronal. Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 21 Figura 2.1. Esquema del mecanismo eléctrico de una neurona (a) y su correspondiente circuito equivalente (b), donde 𝑪𝒎representa la capacitancia eléctrica de la membrana y 𝑽𝒎 el potencial de membrana a través de esta capacitancia. Imagen obtenida de Xiao-Ling Z., A Mathematical Model of a Neuron with Synapses based on Physiology. 2.1. MODELO DENDRÍTICO Las dendritas pueden considerarse como delgados tubos envueltos de una membrana que se comporta como un buen aislante eléctrico en comparación con la resistencia que aporta el núcleo intracelular o el fluido extracelular. Para una longitud corta de la dendrita la corriente eléctrica tiende a fluir de forma paralela al eje del cilindro. Esta es la razón por la que sólo se considera una dimensión espacial (eje x), asumiendo que el voltaje a lo largo de la membrana es una función que depende del tiempo t y de la distancia x a lo largo del conductor. Matemáticamente, el modelo dendrítico compartimental es una aproximación discreta en diferencias finitas a la ecuación no lineal del modelo de cable (figura 2.2). En el modelo dendrítico compartimental (figura 2.3) se presupone que los segmentos dendríticos son isopotenciales y se representan en un solo compartimento de membrana, cada uno con una longitud y un diámetro predeterminados. Se asume que las propiedades eléctricas y citoplasmáticas de la membrana son pasivas (dependientes de voltaje). Cada compartimento tiene tres elementos eléctricos Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 22 pasivos: una fuente de fuerza electromotriz (emf), una resistencia y un condensador [4]. Figura 2.2. Modelo dendrítico de cable. Imagen obtenida de J.M. Bower & Beeman D. The Book of GENESIS Figura 2.3. Modelo dendrítico compartimental. Imagen obtenida de J.M. Bower & Beeman D. The Book of GENESIS Los compartimentos se conectan unos con otros mediante una resistencia axial y las diferencias en las propiedades físicas (diámetro, propiedades de membrana…) y el potencial ocurren entre compartimentos y no dentro de ellos. La entrada presináptica a la dendrita se modela como una fuente de corriente, de forma que la ecuación para el primer y último compartimento dendrítico quedan definidos por 𝐶𝑚𝑑𝑣1 𝑑𝑡 +𝑣𝑚1−𝑉𝑇𝐻 𝑅𝑒𝑞 +𝑣𝑚1−𝑣𝑚2 𝑅𝑎=𝑖𝑠 (2.1) 𝐶𝑚𝑑𝑣𝑚𝑛 𝑑𝑡 +𝑣𝑚𝑛−𝑉𝑇𝐻 𝑅𝐸𝑄 +𝑣𝑚𝑛−𝑣𝑚(𝑛−1) 𝑅𝑎=0 (2.2) siendo 𝑣𝑚1 el potencial de membrana del primer compartimento y 𝑣𝑚2 el del segundo. Los compartimentos dendríticos intermedios quedan definidos por la siguiente ecuación: 𝐶𝑚𝑑𝑣𝑚𝑛 𝑑𝑡 +𝑣𝑚𝑛−𝑉𝑇𝐻 𝑅𝐸𝑄 +𝑣𝑚𝑛−𝑣𝑚(𝑛−1) 𝑅𝑎+𝑣𝑚𝑛−𝑣𝑚(𝑛+3) 𝑅𝑎=0 (2.3) Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 23 La resistencia de membrana 𝑅𝑒𝑞, inversamente proporcional al área del cilindro, la capacitancia 𝐶𝑚, que es proporcional al área y la emf 𝑉𝑇𝐻, caracterizan cada compartimento, con 𝐶𝑚 expresada por 𝐶𝑚=𝜋·𝑙·𝑑·𝐶𝑀 (2.4) siendo 𝐶𝑀 la capacitancia específica de membrana en 𝐹/𝑚2 𝑦 𝑅𝑎 la resistencia axial en Ω/𝑚, que es proporcional a su longitud e inversamente proporcional a su sección transversal (𝑅𝑎=4𝑙𝑅𝐴 𝜋 𝑑2). El tiempo de subida (rise time) determina el retardo que se produce en la propagación del potencial a través de la dendrita. La mayor eficiencia se consigue con el mínimo número de compartimentos [10]. 2.2. MODELO AXÓNICO El modelado del axón que se va a utilizar está basado en el modelo de Hodgkin-Huxley (HH en adelante), consistente en cuatro ecuaciones no lineales que describen el potencial de membrana producido por cambios en la conductancia en el cono axónico. Además del modelo HH, existen un gran número de modelos que puede utilizarse para describir la transmisión en el axón [13] [26]:  Modelo Fitzhugh-Nagumo, en el que las variables voltaje y activación de sodio se representan como una sola variable dependiente del tiempo, y las variables correspondientes a la activación del potasio y la inactivación del sodio, quedan representadas por una sola variable. De esta forma se consigue simplificar mucho el modelo y conseguir un sistema autónomo al no depender las derivadas explícitamente del tiempo.  Modelo de Morris-Lecar. Muy usado en la neurociencia computacional, consiste en una ecuación para modelar el potencial de membrana con activación instantánea de la corriente de calcio y una ecuación adicional que describe la activación más lenta de la corriente de potasio.  Modelo Traub-Miles. Compuesto por un solo compartimento, incluye sólo el sodio, un rectificador para el retardo y una corriente de fuga. Contiene las mismas corrientes que el modelo de Hodgkin-Huxley, que modela a las neuronas de clase II, pero con los parámetros ligeramente cambiados, pues modela a las neuronas de clase I. La diferencia entre un tipo y otro de neuronas radica en el tiempo de latencia. En las neuronas de clase I su longitud es arbitraria y en las de clase II son finitas y cortas en comparación con el intervalo que hay entre una espiga de potencial y otra. Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 24 2.2.1. Modelo de Hodgkin-Huxley para los potenciales de acción Hodking y Huxley fueron los primeros en medir las propiedades de las corrientes iónicas en el axón del calamar gigante y en sintetizar los resultados en un modelo basado en las conductancias de los canales [12]. Este modelo utiliza tres corrientes iónicas diferentes: las provocadas por el flujo de iones de sodio, de potasio y una corriente de fuga compuesta principalmente por iones cloro. Los canales específicos dependientes de voltaje del sodio y del potasio controlan el flujo de iones a través de la membrana celular. La corriente de fuga modela el funcionamiento de otro tipo de canales que no han sido descritos de forma específica. La membrana de la célula neuronal es modelada como un condensador que separa la carga eléctrica. El flujo de corriente que atraviesa la membrana modifica el equilibrio eléctrico e inicia un potencial de acción. Los potenciales reversos 𝐸𝑁𝑎, 𝐸𝐾 y 𝐸𝐿existen debido al gradiente de concentración de iones a través de la membrana producido por la bomba sodio-potasio. La ecuación del balance de corrientes obtenida aplicando la ley de Kirchhoff es 𝐶𝑑𝑉𝑚 𝑑𝑡 = −𝐼𝑁𝑎 −𝐼𝐾−𝐼𝐿+𝐼 (2.5) donde 𝐼 corresponde con la corriente insertada en la célula. El modelo compartimental basado en una analogía circuital se muestra en la figura 2.4. Figura 2.4. Modelo compartimental de la membrana. Imagen obtenida de Jan Benda. Single Neuron Dynamics-Models Linking Theory and Experiment Las corrientes iónicas a través de los canales iónicos están determinadas por su apertura y su cierre. La variable del canal representa su grado de apertura y sigue una cinética de primer orden definida por las constantes dependientes de voltaje 𝛼(𝑉) y 𝛽(𝑉). Los canales de sodio son controlados por las variables de canal 𝑚 y ℎ, ambas dependientes del potencial de membrana. La variable de activación 𝑚 se incrementa con el potencial de membrana mientras que la variable de inactivación disminuye con él. La dinámica de los canales de sodio queda descrita por las siguientes ecuaciones: Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 25 𝐼𝑁𝑎 =𝑔𝑁𝑎𝑚3ℎ (𝑉−𝐸𝑁𝑎) (2.6) 𝑚󰇗 =𝛼𝑚(𝑉)(1−𝑚)−𝛽𝑚(𝑉)𝑚 (2.7) ℎ󰇗=𝛼ℎ(𝑉)(1−ℎ)−𝛽ℎ(𝑉)ℎ (2.8) 𝛼𝑚=2.5−0.1𝑢 exp(1−0.1𝑉)−1 (2.9) 𝛽𝑚=2.5−0.1𝑣 exp(2.5−0.1𝑉)−1 (2.10) 𝛼ℎ=0.07exp (−𝑉 20) (2.11) 𝛽ℎ=1 (exp(3−0.1𝑢)+1 (2.12) El flujo de corriente iónica a través de los canales de potasio está controlado por la variable de canal 𝑛, que crece a medida que se incrementa el potencial de membrana. La dinámica general de los canales de potasio queda descrita por las siguientes ecuaciones: 𝐼𝐾=𝑔𝐾𝑛4(𝐸𝐾−𝑉) (2.13) 𝑛󰇗 =𝛼𝑛(𝑉)(1−𝑛)−𝛽𝑛(𝑉)𝑛 (2.14) 𝛼𝑛=0.1−0.01𝑉 exp(1−0.1𝑉)−1 (2.15) 𝛽𝑛=0.125exp(−𝑉 80) (2.16) La corriente de fuga 𝐼𝐿 tiene un comportamiento lineal con el potencial de membrana que sigue la ecuación 𝐼𝐿=𝑔𝐿(𝐸𝐿−𝑉) (2.17) y las constantes que definen el comportamiento de los tres canales son 𝐸𝑁𝑎 =115𝑚𝑉 𝐸𝐾= −12𝑚𝑉 𝐸𝐿=10.6 𝑚𝑉 𝑔𝑁𝑎 =120𝑚𝑆/𝑐𝑚2 𝑔𝐾=36𝑚𝑆/𝑐𝑚2 𝑔𝐿=0.3 𝑚𝑆/𝑐𝑚2 Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 32 2.6. MODELO DE LA PLACA MOTORA La unidad motora está compuesta por la motoneurona alfa eferente procedente de la médula espinal. La depolarización de esta neurona es la que se traduce en una fuerza en el músculo receptor. En este proceso se distinguen tres fases fundamentales: la excitación de la motoneurona producida por el estímulo transmitido desde la neurona sensitiva, una activación dinámica que lleva a la neurona a su estado activo una vez que se ha producido la depolarización y la generación de la fuerza muscular [11]. Estas tres fases fundamentales se ilustran en la figura 2.10. La corriente de depolarización que circula por la motoneurona constituye la entrada del modelo y la salida es la fuerza muscular. Figura 2.10. Diagrama de bloques del modelo de activación muscular, compuesto por el bloque de excitación de la motoneurona, el bloque de adaptación y activación dinámicas y el bloque de la mecánica muscular. Imagen obtenida de Hannaford B. A Nonlinear Model of the Phasic Dynamics of Muscle Activation La fuerza isométrica producida por una corriente de entrada puede obtenerse a partir de dos ecuaciones diferenciales de primer orden obtenidas del modelo generalizado de Ding et al [8]. La primera relaciona el impulso de entrada con la concentración de calcio en filamentos musculares, mientras que la segunda relaciona la concentración de calcio con la fuerza muscular: Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 33 𝑑𝐶𝑛 𝑑𝑡 +𝐶𝑁 𝜏𝑐=𝛿(𝑡−𝑡𝑖) (2.29) 𝑑𝐹 𝑑𝑡 +𝐹 𝜏1+𝜏2 𝐶𝑁 𝑤 𝐶𝑁 𝑤+𝑘𝑚 𝑤 =𝐴 𝐶𝑁 𝑊 𝐶𝑁 𝑊+𝑘𝑚 𝑤 (2.30) La entrada del modelo está compuesta por un sumatorio de funciones delta, correspondientes a los picos que constituyen el impulso de corriente que llega al músculo. La salida del modelo es una fuerza isométrica, derivada a partir de un muelle lineal, un amortiguador y un motor colocados en serie. Estos elementos representan las componentes elásticas, conectivas y contráctiles durante la contracción muscular, respectivamente. Los coeficientes 𝜏𝑐,𝜏1,𝜏2,𝑘𝑚, 𝐴 𝑦 𝑤 son parámetros que varían en función de las características del sujeto. 𝐶𝑁 representa la concentración de calcio en los filamentos musculares y 𝐹 la fuerza isométrica. 𝐾𝑚 denota la sensibilidad para que se produzca una fuerza en función de 𝐶𝑁 y 𝐴 es un factor para escalarla. El segundo término de la segunda ecuación describe una disminución de la fuerza debido a las dos constantes de tiempo 𝜏1y 𝜏2. 𝜏1modela una disminución de la fuerza provocada por las componentes viscoelásticas de la estimulación del músculo cuando 𝐶𝑁 es pequeña, mientras que el parámetro 𝜏2 modela la disminución de la fuerza debido a las componentes viscoelásticas del músculo durante la estimulación. Los valores utilizados para los parámetros son los que se han empleado en el artículo [8], que han sido estimados a partir de datos experimentales. El modelo se comporta de forma que ante un tren de espigas de la motoneurona se produce una contracción con amplitud proporcional a la frecuencia del tren y una duración proporcional a la duración del tren [34] (figura 2.11). Figura 2.11. Representación de la fuerza en función del tiempo (B) para varios impulsos de entrada de distinta frecuencia (A). Imagen obtenida de Wilson et Al., A Predictive Model of the Isometric Force Response of the Locust Extensor Muscle. Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 34 3. HERRAMIENTAS DE MODELADO ORIENTADO A BLOQUES 3.1. INTRODUCCIÓN Una simulación es el proceso de diseñar un modelo de un sistema real o imaginario para poder realizar experimentos con ese modelo. El objetivo consiste en entender el comportamiento del sistema y vincular las observaciones realizadas con patrones entendibles. Para ello se hacen suposiciones basadas en relaciones matemáticas referidos el sistema que se está modelando. De esta forma se consigue un coste mucho más bajo, menos peligroso, más rápido y más práctico que probando de forma experimental con el sistema real. Existen dos técnicas de modelado básicas: causal y acausal. En las técnicas de modelado causal se le atribuye una mayor importancia a las ecuaciones dinámicas que a la estructura física del sistema que se está modelando. Estas técnicas funcionan muy bien con sistemas de control pero presentan algunas desventajas cuando representan ciertos sistemas físicos. Los sistemas físicos se expresan por medio de ecuaciones diferenciales que tienen que ser resueltas de forma simultánea y aproximada utilizando métodos numéricos. Las herramientas basadas en técnicas acausales se basan en la conexión de componentes físicos que constituyen el modelo. Además, las ecuaciones que definen el sistema se crean automáticamente por conectividad, permitiendo el modelado multidominio [17]. 3.2. SIMULINK La herramienta que se ha empleado para modelar el sistema siguiendo una metodología causal es Simulink, un entorno de programación visual que funciona sobre el entorno de programación Matlab, utilizando un nivel de abstracción mayor. Simulink se basa en diagramas de bloque para la simulación multidominio y el diseño de modelos (figura 3.1) Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 35 Figura 3.1. Estructura característica del entorno de trabajo de Simulink. Las herramientas de modelado orientado a bloques utilizan bloques conectados jerárquicamente. Las señales que se transmiten de un bloque a otro transfieren los valores de las variables individuales de la salida de un bloque a la entrada de otro bloque, obteniéndose un sistema de diferentes operadores matemáticos con entradas y salidas específicas. Por lo tanto, la conexión entre bloques refleja el procedimiento de cálculo más que la estructura de la realidad modelada. No obstante, Simulink permite el modelado de sistemas complejos a través del uso de bloques siguiendo una estructura jerarquizada desde un nivel general hasta diferentes niveles de detalle (figura 3.2). Figura 3.2. Representación de la estructura jerarquizada de Simulink. Bloques generales de la sinapsis y de la dendrita + soma de la motoneurona 𝜶 (A) y detalle de los subbloques que componen el bloque dendrita + soma en un nivel inferior (B). Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 36 4. IMPLEMENTACIÓN DE LA NEURONA UTILIZANDO SIMULINK La implementación del modelo neuronal utilizando una metodología causal va a realizarse haciendo uso de la herramienta Simulink. Para ello se van a modelar por separado los distintos componentes que conforman el circuito del arco reflejo, ya comentados en el apartado 2. En una primera etapa se modelarán la neurona aferente, compuesta por un axón, y la neurona eferente, que incluye el árbol dendrítico, el soma y el axón. Estas dos neuronas estarán conectadas por una sinapsis excitatoria. Una segunda vía conectará la neurona aferente con una interneurona inhibitoria que, a través de una sinapsis, transmitirá el impulso a la motoneurona 𝛼 del músculo flexor. El impulso inicial se obtiene a partir del modelado del receptor sensorial mientras que la fuerza final se obtiene a partir del modelado de la placa motora. El esquema de la estructura comentada se muestra en la figura 4.1. 4.1. MODELO DENDRÍTICO. Para la implementación del modelo dendrítico se han utilizado las ecuaciones mencionadas en el apartado 2.1. Cada terminal dendrítico se ha implementado con tres bloques compartimentales, siendo la entrada de cada uno la corriente del compartimento anterior. El primer compartimento del terminal dendrítico tiene como entrada la corriente postsináptica de la sinapsis anterior. El potencial de membrana generado en un compartimento constituye la entrada del compartimento siguiente y la Figura 4.1. Esquema del arco reflejo con las distintas estructuras que se van a modelar: receptor sensorial, neurona aferente, interneurona, neuronas eferentes y placa motora. Imagen obtenida de Saladin K.S. Anatomy & Physiology. McGraw-Hill Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 37 diferencia de potencial del último compartimento constituye la entrada del compartimento del soma. Los valores utilizados para las constantes son los empleados en el artículo [10], obtenidos a partir de datos experimentales. Se ha utilizado una fuerza electromotriz 𝑉𝑇𝐻 de -60 mV, una capacitancia 𝐶𝑚 de 0.45 𝜇𝐹, una resistencia 𝑅𝑒𝑞 de 3.1 𝑘Ω y una resistencia axial 𝑅𝑎 de 100 Ω, que es la que representa el retardo que supone el compartimento dendrítico en la transmisión de la señal. La implementación del modelo dendrítico completo se representa en la figura 4.2. Como se puede apreciar, el potencial de cada compartimento constituye la entrada del compartimento siguiente. Además, cada compartimento tiene como entrada el potencial generado en el compartimento siguiente, siendo en el último compartimento esta entrada el potencial producido por el soma. Cada uno de los compartimentos que constituyen el terminal dendrítico se implementan siguiendo las ecuaciones 2.1, 2.2 y 2.3, según corresponda. En la figura 4.3 se muestra un compartimento dendrítico intermedio, cuyas entradas son el potencial generado por el compartimento anterior (𝑉𝑚1) y el potencial generado por el compartimento siguiente (𝑉𝑚3). Figura 4.2. Implementación en Simulink del terminal dendrítico completo conectado al soma Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 38 Figura 4.3. Implementación en Simulink de un compartimento dendrítico intermedio. 4.2. MODELO AXÓNICO. El modelo axónico implementado puede dividirse en dos partes, una parte mielínica, que se representa con compartimentos pasivos de características muy similares a las empleadas en el modelo dendrítico, y una parte activa, que corresponde con los nódulos de Ranvier situados entre cada uno de los fragmentos mielínicos. La dinámica de los nódulos de Ranvier se modela utilizando las ecuaciones 2.6-2.17 del modelo de Hodgkin-Huxley descritas en el apartado 2.2.1. Cada axón tiene como entrada la corriente iónica que ha sido generada en el soma (𝐼𝑖𝑛) y la corriente de salida del compartimento axónico es la entrada hacia la placa motora o hacia la siguiente neurona, según corresponda. En el modelo realizado se ha simplificado la estructura real de un axón para reducir el gasto computacional total. Sólo se han incluido tres nódulos de Ranvier y los segmentos mielínicos situados entre ellos se han modelado por medio de tres segmentos pasivos cada uno. En un axón real la distancia entre los nódulos de Ranvier oscila entre los 0.4 y 0.8 mm y en el caso de la neurona sensitiva la longitud del axón es aproximadamente de medio metro, por lo que para la implementación de un modelo fiel a la realidad habría que haber incluido aproximadamente unos 1200 nódulos de Ranvier. No obstante, la simplificación realizada no influye demasiado en la transmisión del estímulo puesto que la atenuación producida por los nódulos es prácticamente imperceptible. Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 39 El compartimento axónico modelado se representa en la figura 4.4. Los bloques activos 1, 2 y 3 corresponden con los nódulos de Ranvier y los bloques pasivos con los segmentos mielínicos. El potencial de membrana 𝑉_𝑚𝑒𝑚𝑏𝑟𝑎𝑛𝑎 generado en el primer compartimento activo constituye la entrada del compartimento anterior, en este caso el soma. Del último compartimento pasivo se obtiene la entrada del siguiente bloque (𝐼_𝑝𝑙𝑎𝑐𝑎𝑚𝑜𝑡𝑜𝑟𝑎 en este caso). Figura 4.4. Implementación en Simulink del modelo axónico compuesto por compartimentos activos y pasivos. La estructura interna de un nódulo de Ranvier se muestra en la figura 4.5. Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 40 Figura 4.5. Implementación en Simulink de un compartimento activo, correspondiente al nódulo de Ranvier. Los tres bloques principales del modelo representan la corriente de fuga, la corriente iónica de los canales de sodio y la corriente iónica producida por los canales de potasio, que, junto con la intensidad 𝐼 de entrada procedente del soma, permiten obtener la corriente que circula por el condensador. Tras integrar y dividir entre 𝐶𝑚 esta corriente, se obtiene el potencial de membrana axónico. Puesto que las constantes utilizadas para la implementación del modelo axónico suponen un potencial de reposo de 0 𝑚𝑉, a este potencial se le suma el potencial de reposo de −65 𝑚𝑉 para que la salida axónica concuerde con las constantes empleadas en los compartimentos pasivos, cuyas constantes sí que suponen un potencial de reposo de −65𝑚𝑉. El potencial obtenido se divide entre la resistencia axial 𝑅𝑎 responsable de la atenuación en la transmisión del voltaje debido a la distancia. Los compartimentos de la corriente producida por los canales de sodio, la corriente producida por los canales de potasio y la corriente de fuga se representan en las figura 4.6, 4.7 y 4.8, respectivamente. En ellos se han implementado las variables de canal 𝑚,𝑛 y ℎ dependientes del potencial de membrana. Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 41 Figura 4.6. Compartimento para la obtención de la corriente de los canales de sodio Figura 4.7. Compartimento para la obtención de la corriente que circula a través de los canales de potasio Figura 4.8. Compartimento para la obtención de la corriente de fuga 4.3. MODELO DEL SOMA El soma se ha modelado utilizando un compartimento pasivo similar al empleado en el terminal dendrítico. El valor elegido para la resistencia axial es de 0.1 𝑘Ω en vez de 1 𝑘Ω, puesto que la transmisión del potencial sufre una atenuación menor en el soma debido, entre otros factores, a su menor tamaño en relación con el del árbol dendrítico. Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 48 su vez está conectado con el axón. El axón transmite el impulso a la placa motora, en la que se realiza la transducción en una fuerza de contracción. La tercera bifurcación del axón de la neurona sensitiva realiza la sinapsis con una neurona intermedia, que a su vez transmite el impulso a una motoneurona 𝛼 que se dirige al músculo flexor. La primera sinapsis se ha modelado como una sinapsis excitatoria. La interneurona situada en la médula espinal cuenta con un terminal dendrítico y un terminal axónico. Debido a que su tamaño es mucho más pequeño, el potencial de membrana sufre una atenuación mucho menor, por lo que se han incluido sólo tres compartimentos pasivos y uno activo en el modelado del axón en vez de los nueve pasivos y los tres activos utilizados en el modelado de las neuronas sensitivas y las motoneuronas. La interneurona permite modular y controlar de forma más eficiente la transmisión de la señal utilizando vías excitatorias o inhibitorias. Por último, esta interneurona realiza una sinapsis inhibitoria con la motoneurona alfa del músculo extensor, en la que el potencial de acción que se genera es nulo debido a la sinapsis inhibitoria. La estructura general del modelo implementado en Simulink se muestra en la figura 4.17. Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 49 Figura 4.17. Modelo completo del arco reflejo patelar implementado en Simulink Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 50 5. ESTUDIO BAJO CONDICIONES FISIOLÓGICAS, PATOLÓGICAS Y TRATAMIENTO FARMACOLÓGICO Una vez implementado el modelo completo se procede a analizar los resultados obtenidos en distintas condiciones, primero suponiendo un comportamiento fisiológico normal y posteriormente bajo situaciones patológicas para poder estudiar de forma cualitativa los efectos de diversos fármacos en la transmisión de la señal. El solver utilizado en la simulación es el método discreto ode3 (Bogacki-Sahmpine), con un tamaño de paso de 0.001. La duración de las simulaciones ejecutadas es de 500 𝑚𝑠. 5.1. ESTUDIO BAJO CONDICIONES FISIOLÓGICAS. Para estudiar el comportamiento en condiciones fisiológicas se introducen varios impulsos de entrada en los que la longitud del músculo varía desde un estiramiento mínimo, considerado como 0.90 𝐿𝑜, hasta un estiramiento máximo. El estiramiento mínimo que permite obtener una frecuencia de disparo adecuada es de 1.015 𝐿0, siendo 𝐿𝑜 la longitud inicial del fascículo muscular. A medida que se incrementa el estiramiento del músculo se obtiene una mayor frecuencia de disparo a la salida del bloque que modela el huso muscular. El estiramiento máximo que permite generar potenciales de acción corresponde con 1.17 𝐿0. Se han considerado unas velocidades de estiramiento muscular muy rápidas, del orden de 0.01 𝑚𝑚/𝑚𝑠. La intensidad que se introduce en el terminal axónico de la neurona sensitiva se obtiene multiplicando la frecuencia de salida del huso muscular por una ganancia de 0.35, puesto que la intensidad que circula a través de la neurona es el resultado del sumatorio ponderado de frecuencias axónicas, siendo 𝑘𝑛 la ganancia calculada para que el sistema exhiba la respuesta deseada. En la figura 5.1 se muestra el impuso de corriente introducido en el axón de la neurona sensitiva para los estiramientos mínimos y máximos. Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 51 Figura 5.1. Impulso de corriente obtenido para un estiramiento de 1.015 𝑳𝟎 (A) y 1.17 𝑳𝟎 (B) Si la amplitud y la duración de la intensidad son lo suficientemente elevados, ésta provoca que el potencial de membrana aumente, haciendo que las conductancias de los canales de sodio se incrementen. Como consecuencia, los iones positivos de sodio penetran en la célula e incrementan aún más el potencial de membrana. Si este feedback positivo es lo suficientemente grande, se producirá un potencial de acción. En nuestro caso, el umbral de iniciación del potencial de acción requiere aproximadamente una corriente de entrada de 10 𝜇𝐴 . Esta corriente genera un potencial de acción de unos 80 𝑚𝑉. Para corrientes de entrada inferiores no se produce potencial de acción, o este tiene una amplitud muy pequeña y no logra transmitirse hasta la última etapa del arco reflejo. Se observa que una corriente de Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 52 entrada constante (en nuestro caso de una duración de unos 150 𝑚𝑠) genera un tren de pulsos de potenciales de acción (figura 5.2) y que la frecuencia de estos potenciales está determinada por la intensidad de la corriente que los produce. Figura 5.2. Potencial de membrana generado en el axón de la neurona sensitiva para un estiramiento mínimo de 1.015 𝐋𝟎 (A) que produce una intensidad con una frecuencia reducida, y un estiramiento máximo de 1.17 𝐋𝟎 (B) que genera una intensidad de entrada de mayor frecuencia. La transmisión del potencial de acción generado a lo largo de las distintas etapas del circuito del arco reflejo produce una atenuación de la señal debido a las resistencias de membrana de cada bloque y a las resistencias axiales que conectan unos Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 53 compartimentos con otros. Esta es la razón por la que se observan menos potenciales a la salida. Además, se introduce un retardo en la transmisión de la señal provocado principalmente por la sinapsis, en la que los neurotransmisores requieren un tiempo para ser liberados y atravesar el espacio sináptico. El retardo total en la transmisión de la señal es de unos 20 𝑚𝑠, siendo despreciable el retardo introducido en la transmisión a lo largo del axón y del árbol dendrítico. Estos efectos pueden apreciarse en la figura 5.3, en la que se muestran los potenciales de acción generados al final del axón de la motoneurona 𝛼 del músculo extensor para el caso del estiramiento mínimo y máximo. Figura 5.3. Potencial de membrana en la última etapa del circuito del arco reflejo para un estiramiento mínimo de 1.015 𝐋𝟎 (A) y un estiramiento máximo de 1.17 𝐋𝟎 (B). A pesar de la atenuación introducida en el potencial de acción, se sigue observando que para un estiramiento máximo la frecuencia del potencial es mayor. Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 54 Finalmente, la corriente iónica que genera el potencial de membrana mostrado en la figura 5.3 se transmite a la placa motora induciendo una contracción del músculo extensor, de forma que cuanto mayor sea la frecuencia de entrada mayor será la fuerza de contracción generada en el músculo. En la figura 5.4 se muestra la evolución de esta fuerza de contracción normalizada frente al tiempo para los estiramientos mínimo y máximo de la entrada. Estos resultados concuerdan con los obtenidos en el modelo original, comentados en la figura 2.11. Figura 5.4. Representación de la fuerza normalizada de contracción frente al tiempo para un estiramiento mínimo de 1.015 𝐋𝟎 (A) y un estiramiento máximo de 1.17 𝐋𝟎 (B) El estiramiento máximo produce una mayor frecuencia en el impulso de corriente que llega al músculo, haciendo que su contracción sea mayor. Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 55 La contracción inducida en la motoneurona 𝛼 del músculo flexor es prácticamente nula puesto que la sinapsis inhibitoria que se produce entre la interneurona de la médula espinal y dicha motoneurona hace que el potencial de membrana generado tras la sinapsis se atenúe por completo. 5.2. ESTUDIO BAJO CONDICIONES PATOLÓGICAS Y TRATAMIENTO FARMACOLÓGICO 5.2.1. Enfermedades de la motoneurona Las enfermedades de la motoneurona son un grupo heterogéneo de enfermedades, algunas de ellas hereditarias, con afectación exclusiva o predominante de las neuronas motoras de la corteza cerebral, de los núcleos motores del tronco del encéfalo y/o de las astas anteriores de la médula espinal. En nuestro caso se van a estudiar algunas patologías pertenecientes a este grupo que afectan directamente a las motoneuronas inferiores, que son las que han sido implementadas en el modelo realizado. 5.2.1.1. Esclerosis lateral amiotrófica La esclerosis lateral amiotrófica (ALS en adelante) es una enfermedad degenerativa de tipo neuromuscular. Se produce cuando las motoneuronas van disminuyendo gradualmente su funcionamiento hasta que acaban muriendo, provocando una parálisis muscular progresiva de pronóstico mortal. En sus últimas etapas el paciente sufre una parálisis total acompañada de una exaltación de los reflejos tendinosos como consecuencia de una pérdida de los controles musculares inhibitorios. La ALS afecta tanto a las motoneuronas inferiores como a las motoneuronas superiores, que se encuentran situadas en la parte superior del cerebro ejerciendo cierto control sobre las motoneuronas inferiores y permitiendo que los movimientos musculares sean uniformes, dirigidos y con intensidad diversa. Cuando se pierden las neuronas motoras superiores y quedan las neuronas motoras inferiores, todavía pueden hacerse movimientos, pero se vuelven espasmódicos y menos precisos. Cuando la ALS afecta a las motoneuronas inferiores, los músculos se debilitan progresivamente cada vez más hasta llegar a su parálisis. Por lo general, en la ALS se observa una combinación de estos efectos debido a que tanto las neuronas superiores como las inferiores están muriendo. Las personas con ALS pueden tener músculos débiles y degenerados, acompañados por tensión (espasticidad). Las sacudidas y calambres musculares son comunes debido a los axones degenerados. Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 56 La evidencia que proporcionan enfermos de ALS apunta hacia una sobreabundancia de glutamato en el sistema nervioso. El glutamato es una sustancia química común en el sistema nervioso que las neuronas utilizan para enviar señales a otras neuronas. Si su concentración es excesiva induce la muerte de las células nerviosas que reciben la señal. Esto puede deberse al transporte inadecuado de glutamato lejos de las células nerviosas después de que haya terminado su función de envío de señales. Los experimentos sugieren que un defecto podría deberse también al exceso de producción o liberación de glutamato por las células transmisoras o que podría deberse a defectos en los receptores de glutamato en las células receptoras. Diversos estudios [38] sugieren que las propiedades fisiológicas de los canales dependientes de sodio se ven modificadas en las motoneuronas de pacientes con ALS, insinuando que una alteración de su funcionamiento en condiciones normales puede inducir esta patología [27]. Tratamiento con Riluzol Uno de los medicamentos más empleados para el tratamiento de la ALS es el Riluzol [19], cuyo principal mecanismo de acción se da por una inhibición en la liberación del glutamato debido a una inactivación de canales de sodio dependientes de voltaje en las terminaciones de glutamato de las neuronas. Riluzol bloquea pues, de manera reversible, los canales de sodio. Sabiendo que Riluzol puede ejercer su efecto sobre estos canales [32][38], se pueden realizar simulaciones para conocer el efecto que tiene sobre el sistema completo una modificación del número de canales de sodio activos. En este caso, se va a simular qué efecto tiene una disminución del potencial en los canales de sodio para reproducir en términos cualitativos el efecto del Riluzol en la contracción muscular. Se ha modificado la conductancia del sodio que tiene un habitualmente un valor de 120 𝑚𝑆/𝑐𝑚2 hasta un valor de 90 𝑚𝑆/𝑐𝑚2 para un estiramiento intermedio de 1.082 𝐿𝑜. En la figura 5.5 se observa que un bloqueo de los canales de sodio disminuye bruscamente la fuerza de contracción en el músculo efector. Estos resultados son coherentes con la forma de actuación del Riluzol [22] puesto que un bloqueo de los canales de sodio provoca una disminución en la concentración excesiva de glutamato y una disminución de la tensión muscular. Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 57 Figura 5.5. Fuerza de contracción muscular para un estiramiento promedio de 1.082 𝐋𝐨 con una conductancia estándar inicial en los canales de sodio de 120 𝐦𝐒/𝐜𝐦𝟐 (A) y una conductancia final disminuida de 90 𝐦𝐒/𝐜𝐦𝟐 (B), simulando la administración del fármaco Riluzol. El resultado obtenido en la simulación concuerda con las pruebas experimentales que se han realizado sobre este medicamento. En las figuras 5.6 y 5.7, correspondientes a un experimento realizado con sujetos reales [33], se muestra cómo después de la aplicación del Riluzol la fuerza de respuesta en el tobillo ante un estímulo aplicado disminuye considerablemente, mientras que utilizando un placebo ésta no sufre ninguna modificación. Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 64 Figura 5.12. Fuerza de contracción muscular para un estiramiento promedio de 1.082 𝐋𝐨 bajo condiciones fisiológicas normales (A) y bajo una inhibición presináptica (B), simulando la administración del medicamento tizanidina. Tratamiento con fenitoína La fenitoína, también llamada difenilhidantoína, es un antiepiléptico de uso común usado sobre todo en para evitar las convulsiones o los calambres. Actúa bloqueando la actividad neuronal no deseada mediante la reducción de la conductividad eléctrica entre las neuronas, bloqueando los canales de sodio sensibles al voltaje [19]. La fenitoína altera la conductancia de los canales de sodio, potasio y calcio, los potenciales de membrana y las concentraciones de ciertos aminoácidos, además de los neurotransmisores acetilcolina y ácido gamma-aminobutírico [14]. Se aplica en el tratamiento de enfermedades de la motoneurona. En primer lugar actúa disminuyendo el número de canales de sodio abiertos, afectando más en la duración en la que el canal está abierto que en la corriente que circula por ellos [14]. El efecto principal de la fenitoína en las corrientes de sodio sugiere que un importante Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 65 mecanismo farmacológico para conseguir una actividad anticonvulsionante consiste en reducir la duración en la que los canales de sodio están abiertos en lugar de centrarse en cerrarlos por completo [30]. Para simular los efectos de la fenitoína se han disminuido las conductancias de los canales de sodio y potasio hasta unos valores de 100 y 30 𝑚𝑆/𝑐𝑚2 en vez de los 120 y 36 𝑚𝑆/𝑐𝑚2 iniciales. La reducción de la actividad se traduce en una menor fuerza de contracción en el músculo efector (figura 5.13.A), obteniendo una fuerza máxima de contracción de 0.4 en vez de la fuerza normalizada de 0.9 que se registra en condiciones normales. En la figura 5.13.B se muestra el menor número de potenciales de membrana que se producen en el axón de la motoneurona 𝛼 Figura 5.13. Fuerza de contracción muscular para un estiramiento promedio de 1.082 𝐋𝐨 tras modificar las conductancias de sodio y de potasio (A) y potencial del membrana en el axón de la motoneurona 𝛂 bajo estas mismas condiciones (B) Estos resultados concuerdan con los obtenidos en estudios experimentales. En la figura 5.14 se muestran los resultados de un experimento realizado con ratas, en el que se aprecia cómo una administración del fármaco fenitoína (PHT) produce una Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 66 disminución en la corriente dependiente de voltaje. Este mismo efecto se observa en la figura 5.15, correspondiente a un estudio realizado con humanos, en el que los picos de corriente están mucho menos marcados tras la administración de fenitoína. Figura 5.14. Efecto de la fenitoína (PHT) en las corrientes dependientes de voltaje comparadas con una situación neutral de control (CTRL) y tras aplicarle un lavado a la muestra (wash). Imagen obtenida de Jeub et al., Effect of phenytoin on sodium and calcium currents in hippocampal CA1 neurons of phenytoin-resistant kindled rats. Figura 5.15. Comparación de los picos de corriente entre una situación de control (izda.) y la administración de fenitoína (dcha.) La fenitoína reduce los picos de corriente como resultado de un cierre de los canales de sodio. Imagen obtenida de Segal M. M. & Douglas A.F. Late Sodium Channel Openings Underlying Epileptiform Activity Are Preferentially Diminished by the Anticonvulsant Phenytoin. Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 67 5.2.2. Esclerosis múltiple. La esclerosis múltiple es una enfermedad autoinmunitaria causada por el daño en la vaina mielínica de la neurona. Cuando esta cubierta de los nervios se daña, los impulsos nerviosos disminuyen o se detienen. La disminución de la mielina causa además una reducción en la velocidad de propagación de la señal [27]. Diversos estudios sugieren que los síntomas permanentes pueden ser producidos por una lesión funcional del axón como consecuencia de un fallo de los mecanismos reparadores de la mielina, por una distrofia axonal debido al adelgazamiento en determinadas zonas del mismo, o bien por una pérdida crítica del número de axones [1]. La espasticidad es uno de los síntomas más comunes de la esclerosis múltiple y hay indicios de que puede estar producida por una desinhibición central del arco reflejo miotático. Las motoneuronas α se encuentran en un estado de hiperexcitabilidad como resultado de la reducción de la inhibición presináptica de las fibras del grupo Ia [1] y se activan por estímulos que normalmente no provocarían respuesta. Se va a simular en nuestro sistema cómo afectaría la esclerosis múltiple en la contracción del músculo extensor. Para ello se han asignado unos valores superiores a las constantes 𝑘𝑝 y 𝑘𝑖 que modelan nuestra sinapsis, implicando de esta forma un mayor tránsito de iones a través de los canales sinápticos. Este efecto simula una disminución de la inhibición sináptica al provocar un aumento de la excitación sináptica. En la figura 5.16 se aprecia cómo una mayor excitabilidad se traduce en una mayor contracción muscular, síntoma característico de la esclerosis múltiple. Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 68 Figura 5.16. Fuerza de contracción muscular para un estiramiento promedio de 1.082 𝐋𝐨 bajo condiciones fisiológicas normales (A) y bajo una inhibición presináptica (B) simulando los efectos de la esclerosis múltiple. El nivel de contracción muscular sufre un incremento considerable en la simulación de la inhibición presináptica. Los resultados obtenidos concuerdan con estudios realizados en sujetos reales. En la figura 5.17 se muestra cómo influye la espasticidad característica de la esclerosis múltiple en la fuerza muscular. Los pacientes que sufren espasticidad (CP) manifiestan para un ángulo determinado un mayor valor de fuerza muscular. Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 69 Figura 5.17. Curvas de la relación fuerza (Nm)-ángulo(°) entre pacientes que sufren espasticidad (línea continua) y pacientes sanos (línea punteada). Imagen obtenida de Peng Q. et al, Quantitative evaluations of ankle spasticity and stiffness in neurological disorders using manual spasticity evaluator. La disminución de la mielina propia de la esclerosis múltiple produce una importante atenuación de la señal puesto que disminuye el valor de la resistencia a través de la membrana. En condiciones fisiológicas normales el valor de la resistencia de la membrana se ve incrementado unas 200 veces gracias a la mielina, aumentando la constante de longitud que representa la distancia a la que se ha producido una atenuación del 63% de la diferencia de potencial. Se han estudiado los efectos que causa en nuestro sistema una disminución de la resistencia de membrana, modificando su valor inicial de 3.1 𝑘Ω a un valor de 0.4 𝑘Ω para simular la brusca reducción que ésta sufre en los pacientes con esclerosis múltiple. La atenuación producida en el potencial de membrana en su transmisión a través del axón de la motoneurona 𝛼 es evidente. En la figura 5.18.A se muestra el potencial de membrana al principio del axón de la neurona sensitiva y al final de este (figura 5.18.B) en condiciones fisiológicas normales. Se puede apreciar que la señal sufre cierta atenuación debido principalmente a la resistencia que ofrece la membrana celular al paso de la corriente. En la figura 5.18.C se muestra el potencial de membrana al final del axón para el caso de una resistencia interna modificada de 0.4 𝑘Ω. Se puede observar que la atenuación que sufre la señal es mucho mayor. Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 70 Figura 5.18. Atenuación sufrida en la transmisión del potencial de membrana a lo largo del axón de la motoneurona en condiciones fisiológicas normales (B) y con una resistencia de membrana atenuada (C). La figura A corresponde con el potencial de membrana todavía sin atenuar a la entrada del axón Estos resultados concuerdan con los efectos reales que produce la desmielinización en la atenuación de los potenciales de membrana debido a que se incrementa el umbral Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 71 de la corriente iónica necesaria para producir el potencial de acción. En la figura 5.19 se muestra cómo el umbral de la corriente iónica de una célula mielinizada (5.19.C) situado en 19.8·10−3 𝜇𝐴 se va incrementando para un grado intermedio de desmielinización del axón (5.19.B), alcanzando un umbral máximo de 75·10−3𝜇𝐴 cuando se produce un grado de desmielinización mayor. Figura 5.19. Estudio experimental de una desmielinización de las fibras nerviosas. La línea horizontal discontinua representa el valor del umbral mínimo de la corriente necesaria para producir un potencial de acción. C corresponde con una fibra nerviosa mielinizada en condiciones normales, B muestra un mayor grado de desmielinización y A muestra un estado más avanzado de la pérdida de mielina. Imagen obtenida de Stephanova D.I., Conduction Along Myelinated and Demyelinated Nerve Fibres During the Recovery Cycle: Model Investigations. Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 72 CONCLUSIONES Los objetivos principales del presente trabajo, consistentes en la implementación de un modelo circuital del arco reflejo completo considerando cada una de las distintas estructuras que componen la neurona, así como el receptor sensorial y el efector final, se han podido satisfacer. Una vez construido por completo el modelo, se ha podido evaluar su comportamiento de forma cualitativa bajo condiciones fisiológicas, patológicas y ante la administración de determinados fármacos. Para ello se han realizado simulaciones bajo el entorno de programación orientado a bloques SIMULINK modificando los parámetros correspondientes que se ven alterados en las diversas patologías consideradas. Los resultados obtenidos exhiben un comportamiento satisfactorio desde un punto de vista cualitativo, dada la complejidad de implementación del sistema del arco reflejo humano en cuanto al número de elementos que lo componen. El hecho de tener un modelo que simule el comportamiento de la transmisión de la señal a lo largo del circuito del arco reflejo y que además permita modificar parámetros, como pueden ser la conductancia de los canales de sodio y potasio, el grado de excitación e inhibición de las sinapsis entre dos neuronas o la resistencia al paso de corriente que ejerce la membrana celular en función de la cantidad de mielina que recubre el axón, supone grandes beneficios. Las implicaciones que conlleva tener un modelo que responda adecuadamente ante estos cambios son, entre otras, que se pueda conocer mejor cómo se propaga el impulso eléctrico a través de las distintas estructuras del circuito, qué factores pueden ser causantes de que esta señal se atenúe o se retarde o poder facilitar la labor del diseño de fármacos que actúen de forma más específica. Hay que reseñar que la implementación realizada no deja de ser un modelo muy simplificado de la realidad, por lo que presenta una serie de limitaciones. Las estructuras modeladas implican una gran complejidad y los mecanismos responsables del funcionamiento del arco reflejo no son conocidos en su totalidad. En la realización del modelo sólo se han tenido en cuenta los parámetros más generales que lo definen, haciendo posible la simulación del arco reflejo en determinadas situaciones bajo ciertas condiciones muy específicas. Por ello, en trabajos futuros se pretende incrementar la complejidad del modelo incluyendo más variables que permitan que su comportamiento se asemeje más a la realidad. Asimismo, haciendo uso de la misma estructura utilizada, se podrían implementar circuitos neuronales más complejos con múltiples sinapsis, como el circuito del reflejo extensor cruzado en el que está implicada la inervación recíproca. Modelado de circuitos neuronales empleando técnicas de programación orientada a bloques. 73 Los principios empleados en este trabajo para el modelado neuronal pueden aplicarse asimismo para modelar la actividad de un gran espectro de células, como pueden ser las fibras oculomotoras responsables de los movimientos sacádicos del ojo o las células cardíacas, por lo que en un futuro se pretende ampliar el campo de aplicación del modelo realizado. Otra línea futura de trabajo consistiría en implementar el mismo circuito pero utilizando bloques físicos, pudiendo incorporar las ventajas que conlleva el uso de una metodología acausal. De esta forma se podría conseguir que el modelo matemático se pueda realizar por interconexión de bloques físicos a través del uso de analogías eléctricas. Para la implementación del modelo siguiendo una metodología acausal se recurriría a herramientas de programación orientada a objetos como Simscape o Modelica.