Numerical simulation of the behaviour of musculoskeletal tissue

==Resumen==

En este trabajo se muestra un modelo computacional tridimensional para la simulación del comportamiento mecánico de la unidad músculo tendón (UMT). El modelo ha sido formulado a través de una función densidad energía de deformación para materiales hiperelásticos que incorpora tanto el comportamiento pasivo de los tejidos conectivos como el activo, considerando para ambos casos direcciones de comportamiento preferencial no necesariamente coincidentes en el espacio. Las tensiones iniciales, presentes habitualmente en tejidos blandos, han sido implementadas mediante la definición de un gradiente de deformación para llevar a cabo la simulación de la UMT. El modelo se ha particularizado para un modelo de experimentación animal. La geometría utilizada se corresponde con la del músculo tibial anterior de rata reconstruido a partir de imagen médica procedente de resonancia magnética obtenida «''in-vivo'' ». También se estudia el efecto que sobre el comportamiento pasivo ejerce la fascia ajustando un modelo de comportamiento hiperelástico para la misma a partir de ensayos experimentales.

==Abstract==

This paper presents a three dimensional computational model to simulate the mechanical behavior of the muscle-tendon unit (MTU). The model has been formulated through a strain energy density function for hyperelastic materials incorporating both the passive and active behavior of the connective tissues, and taking into account preferential directions for both behaviours, not necessarily coincident. The initial stresses, usually present in soft tissues, have been considered when performing the MTU simulation. The geometry used corresponds to the rat tibialis anterior muscle reconstructed from «''in-vivo'' » magnetic resonance images. This paper also studies the effect of the fascia on the MTU passive behaviour and the setting of a hyperelastic model for this tissue by means of experimental tests.  

==Palabras clave==

Elementos finitos ; Tejido músculo-esquelético ; Comportamiento activo y pasivo

==Keywords==

Finite elements ; Musculoskeletal tissue ; Active and passive behavior

==1. Introducción==

El tejido músculo-esquelético es un tejido biológico blando responsable, fundamentalmente, del movimiento del cuerpo y de mantener la postura corporal. El cuerpo humano contiene aproximadamente 650 músculos, con formas diferentes en función del tipo de movimiento, representando aproximadamente el 40% del peso corporal [[#bib0005|[1]]] . Cada uno de estos músculos puede ser considerado como un conjunto de fibras musculares (elemento contractil) rodeadas, en diferentes niveles, por tejido conectivo (elemento elástico) que, de forma pasiva, es el responsable de absorber los alargamientos del músculo. Es bien sabido además que, las fibras musculares también contribuyen al comportamiento pasivo (elemento elástico en serie) [[#bib0010|[2]]]  and [[#bib0015|[3]]] . El tejido conectivo está constituido por fibras de colágeno y elastina embebidas en la sustancia fundamental y organizadas en tres niveles, el epimisio, perimisio y endomisio, que juegan un papel importante en la transmisión lateral de las fuerzas generadas por las fibras y, a través de los cuales, se distribuyen los vasos sanguíneos y ramificaciones nerviosas.      

Los recientes avances en la adquisición de imágenes, en la caracterización experimental de los tejidos biológicos y en el desarrollo de nuevos algoritmos de comportamiento han permitido el desarrollo de modelos 3D que permiten obtener el campo de tensiones y deformaciones en el interior de los tejidos vivos. En la mayoría de las ocasiones, la construcción de dichos modelos lleva asociada la adopción de hipótesis simplificativas, por ejemplo, las propiedades mecánicas se determinan mediante ensayos «''in-vitro'' »  [[#bib0020|[4]]]  and [[#bib0025|[5]]] , o se combinan parámetros de diferentes tejidos [[#bib0030|[6]]] , o bien se adoptan geometrías idealizadas  [[#bib0035|[7]]] , [[#bib0040|[8]]] , [[#bib0045|[9]]] , [[#bib0050|[10]]]  and [[#bib0055|[11]]] , lo que hace difícil validar los resultados obtenidos y por lo tanto poder extrapolar los resultados numéricos a estudios clínicos.      

Con el objetivo de minimizar las posibles causas de error atribuidas a las simplificaciones anteriores, en este trabajo se presenta el modelado del músculo tibial anterior (TA) de rata, partiendo de imagen procedente de resonancia magnética (RM). Se ha utilizado un modelo de comportamiento hiperelástico anisótropo capaz de reproducir la respuesta pasiva y activa del músculo partiendo de los resultados experimentales. También se han incorporado las deformaciones iniciales habituales en los tejidos blandos [[#bib0060|[12]]]  and [[#bib0065|[13]]] . El modelo parte de la caracterización experimental, tanto activa como pasiva de la UMT llevada a cabo por los autores en trabajos previos [[#bib0070|[14]]] , [[#bib0075|[15]]]  and [[#bib0080|[16]]] .      

Los resultados obtenidos en la simulación son validados con los resultados experimentales procedentes de ensayos uniaxiales a tracción «''in-vivo'' » de la UTM, en el caso del comportamiento pasivo  [[#bib0070|[14]]] , y de la estimulación eléctrica en condiciones isométricas, en el caso activo [[#bib0075|[15]]] .      

==2. Modelo matemático==

En este apartado se resumen las ecuaciones que definen el modelo de comportamiento hiperelástico anisótropo habitualmente empleado en la simulación de tejidos biológicos blandos, dentro siempre del dominio elástico. Dicho modelo estructural se plantea en el marco de la mecánica del continuo no lineal y teniendo en cuenta la posibilidad de ser incorporado en un código de elementos finitos.

===2.1. Definición de la deformación===

Se entiende por sólido tridimensional, ''Ω''<sub>''o''</sub> , a un subconjunto de <math display="inline">{\mathbb{R}}^3</math>  cuyos puntos se identifican mediante sus coordenadas en un sistema de referencia. Matemáticamente, podemos interpretar lo anterior a través de una función biunívoca de tres componentes '''''φ<sub>o</sub>'''''  aplicada sobre el sólido ''Ω''<sub>''o''</sub>  tal que, a cada punto '''''P'''''  ∈ ''Ω''<sub>''o''</sub> , adscribe tres valores que corresponden a las coordenadas de dicho punto '''''P'''''  en el sistema de referencia elegido, es decir:

<span id='eq0005'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\boldsymbol{\varphi }}_o:{\Omega }_o\rightarrow {\mathbb{R}}^3\quad {\boldsymbol{\varphi }}_o\equiv \boldsymbol{\mbox{X}};\quad {\varphi }_o^I=</math><math>X^I</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 1)
|}

A '''''φ<sub>o</sub>'''''  la denominaremos ''configuración inicial del sólido Ω''<sub>0</sub> , ''configuración de referencia''  o ''configuración indeformada'' . A lo largo del movimiento del sólido, la posición de cada uno de los puntos del mismo va variando. La configuración correspondiente define las coordenadas de los puntos del sólido en ese instante ''Ω''<sub>''t''</sub>  respecto a un sistema de referencia (que en general también puede variar con el tiempo) a través de una expresión similar a [[#eq0005|(1)]] , es decir:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\boldsymbol{\varphi }}_t:{\Omega }_t\rightarrow {\mathbb{R}}^3\quad {\boldsymbol{\varphi }}_t\equiv \boldsymbol{\mbox{x}};</math>
|-
|<math>{\boldsymbol{\varphi }}_t^I=x^i\quad con\quad t\in I\in \mathbb{R}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 2)
|}

A '''''φ'''''<sub>''t''</sub>  la denominaremos ''configuración actual del sólido Ω''<sub>''t''</sub>  o ''configuración deformada'' .      

Dado un movimiento <math display="inline">\boldsymbol{{\varphi }_t}:{\Omega }_t\rightarrow {\mathbb{R}}^3\quad C^1-</math><math></math>''regular'' , se define el gradiente de deformación '''F''' , como el campo tensorial sobre la configuración indeformada '''''φ'''''<sub>''o''</sub> (''Ω  '' ), <math display="inline">\boldsymbol{\mbox{F}}=</math><math>\frac{\partial \boldsymbol{\mbox{x}}}{\partial \boldsymbol{\mbox{X}}}</math> , siendo ''J''  ≡ ''det'''''F'''  > 0 el jacobiano de la transformación.      

Una de las características más relevantes de los tejidos biológicos blandos, en general, y del músculo en particular, es la ''incompresibilidad'' . Esta indica que el cuerpo al ser sometido a fuerzas o deformaciones no presenta un cambio de volumen, lo que se representa matemáticamente como que el jacobiano del gradiente de deformación es igual a la unidad, ''J''  ≅ 1. Siguiendo a [[#bib0085|[17]]] , definimos la siguiente descomposición multiplicativa del gradiente de deformación:

<span id='eq0015'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\mbox{F}}=J^{\frac{1}{3}}\overline{\boldsymbol{\mbox{F}}},\quad \overline{\boldsymbol{\mbox{F}}}=</math><math>J^{-\frac{1}{3}}\boldsymbol{\mbox{F}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 3)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\mbox{C}}=\boldsymbol{\mbox{F}}^T\boldsymbol{\mbox{F}},\quad \overline{\boldsymbol{\mbox{C}}}=</math><math>J^{-\frac{2}{3}}\boldsymbol{\mbox{C}}={\overline{\boldsymbol{\mbox{F}}}}^T\overline{\boldsymbol{\mbox{F}}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 4)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\mbox{b}}=\boldsymbol{\mbox{F}}\boldsymbol{\mbox{F}}^T,\quad \overline{\boldsymbol{\mbox{b}}}=</math><math>J^{-\frac{2}{3}}\boldsymbol{\mbox{b}}=\overline{\boldsymbol{\mbox{F}}}\quad {\overline{\boldsymbol{\mbox{F}}}}^T</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 5)
|}

el término <math display="inline">J^{\frac{1}{3}}\boldsymbol{\mbox{I}}</math>  representa la deformación volumétrica y <math display="inline">\overline{\boldsymbol{\mbox{F}}}</math>  la deformación desviadora. Denotamos por <math display="inline">\overline{\boldsymbol{\mbox{F}}}</math> , <math display="inline">\overline{\boldsymbol{\mbox{C}}}</math>  y <math display="inline">\overline{\boldsymbol{\mbox{b}}}</math>  al gradiente de deformación y los tensores de deformación de Cauchy-Green por la derecha e izquierda modificados, respectivamente.      

La anisotropía del material en un punto '''X'''  ∈ ''Ω''<sub>0</sub> , va asociada a la dirección de la fibras musculares, en el caso del comportamiento activo, y a la dirección de las fibras de colágeno en el caso del comportamiento pasivo. En los músculos fusiformes como es el caso del TA, la dirección de anisotropía coincide para ambos comportamientos, es decir, el comportamiento preferencial pasivo coincide con la dirección de las fibras musculares, pero no es así en músculos planos, como pueden ser los de la pared abdominal [[#bib0025|[5]]] . Por tanto, en una situación general, se puede partir de la hipótesis de que la dirección de anisotropía, asociada al comportamiento activo, se define por un vector unitario '''m'''<sub>0</sub> ('''X''' ), y la asociada al comportamiento pasivo mediante el vector '''n'''<sub>0</sub> ('''X''' ), con |'''m'''<sub>0</sub> | = 1, |'''n'''<sub>0</sub> | = 1. Definimos la base '''E'''<sub>''i''</sub>  = ('''m'''<sub>0</sub> , '''n'''<sub>0</sub> , '''l'''<sub>0</sub> ), donde el vector '''l'''<sub>0</sub>  = '''m'''<sub>0</sub>  ⊗ '''n'''<sub>0</sub>  está fuera del plano formado por los otros dos vectores.      

La deformación asociada a la actividad muscular (ver [[#fig0005|figura 1]] ) se puede modelar mediante un proceso ficticio de dos pasos [[#bib0090|[18]]] : el primero de ellos asociado al desplazamiento relativo de las proteínas actina y miosina y, el segundo, a una deformación elástica asociada a los puentes cruzados y a la titina. Esto puede ser expresado mediante una descomposición multiplicativa de <math display="inline">\overline{\boldsymbol{\mbox{F}}}</math>[[#bib0090|[18]]] :

<span id='fig0005'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_442579157-1-s2.0-S0213131512000193-gr1.jpg|center|357px|Esquema de la deformación pasiva y activa [18]. Sobre el punto X de la ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 1.

Esquema de la deformación pasiva y activa [[#bib0090|[18]]] . Sobre el punto '''X'''  de la configuración indeformada ''Ω''<sub>0</sub> , consideramos dos direcciones preferenciales definidas por '''m'''<sub>0</sub>  y '''n'''<sub>0</sub> . Esta configuración se transforma en la misma dirección por la contracción activa <math display="inline">{\overline{\boldsymbol{\mbox{F}}}}_a</math> , donde los nuevos vectores serán ''λ''<sub>1</sub>'''m'''<sub>0</sub> , ''λ''<sub>2</sub>'''n'''<sub>0</sub>  y ''λ''<sub>3</sub>'''l'''<sub>0</sub> . La deformación elástica transforma esta configuración intermedia activa a la configuración deformada ''Ω''<sub>''t''</sub>  por medio del gradiente de deformación <math display="inline">{\overline{\boldsymbol{\mbox{F}}}}_e</math> .                  

</span>
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\overline{\boldsymbol{\mbox{F}}}={\overline{\boldsymbol{\mbox{F}}}}_e{\overline{\boldsymbol{\mbox{F}}}}_a</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 6)
|}

donde <math display="inline">{\overline{\boldsymbol{\mbox{F}}}}_a</math>  representa la deformación asociada a la contracción muscular por el movimiento relativo de los filamentos de actina y miosina y <math display="inline">{\overline{\boldsymbol{\mbox{F}}}}_e</math>  representa la deformación elástica asociada a los puentes cruzados y a la titina. El tensor <math display="inline">{\overline{\boldsymbol{\mbox{F}}}}_a</math>  no corresponde a una transformación integrable, partes infinitesimales de ''Ω''<sub>0</sub>  se deforman de manera independiente y su unión puede dar lugar a configuraciones no compatibles. El tensor <math display="inline">{\overline{\boldsymbol{\mbox{F}}}}_e</math>  representa la deformación elástica y garantiza la compatibilidad de la nueva configuración.      

Si las fibras musculares se contraen en la dirección definida por el vector de comportamiento preferencial '''m'''<sub>0</sub>  y se supone que dicha dirección no cambia, el tensor de contracción activo vendrá dado por:

<span id='eq0035'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\overline{\boldsymbol{\mbox{F}}}}_a={\overline{\lambda }}_1\boldsymbol{\mbox{e}}_1\otimes \boldsymbol{\mbox{E}}_1+</math><math>{\overline{\lambda }}_2\boldsymbol{\mbox{e}}_2\otimes \boldsymbol{\mbox{E}}_2+</math><math>{\overline{\lambda }}_3\boldsymbol{\mbox{e}}_3\otimes \boldsymbol{\mbox{E}}_3</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 7)
|}

donde los vectores '''E'''<sub>''i''</sub>  forman una base dual a '''e'''<sub>''i''</sub>  definida como '''e'''<sub>''i''</sub>  · '''E'''<sub>''j''</sub>  = ''δ''<sub>''ij''</sub>  con ''δ''<sub>''ij''</sub>  la delta de Kronecker. <math display="inline">{\overline{\lambda }}_1=</math><math>{\overline{\lambda }}_a</math>  y el valor de <math display="inline">{\overline{\lambda }}_2=</math><math>{\overline{\lambda }}_3</math>  se determina al considerar el material incompresible.

<span id='eq0040'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>det{\overline{\boldsymbol{\mbox{F}}}}_a=1={\overline{\lambda }}_1{\overline{\lambda }}_2{\overline{\lambda }}_3</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 8)
|}

Sustituyendo los valores de <math display="inline">{\overline{\lambda }}_2</math>  y <math display="inline">{\overline{\lambda }}_3</math>  obtenidos al aplicar [[#eq0040|(8)]] , la ecuación [[#eq0035|(7)]]  se reduce a:

<span id='eq0045'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\overline{\boldsymbol{\mbox{F}}}}_a={\overline{\lambda }}_a\boldsymbol{\mbox{e}}_1\otimes \boldsymbol{\mbox{E}}_1+</math><math>{\overline{\lambda }}_a^{-1/2}\boldsymbol{\mbox{e}}_2\otimes \boldsymbol{\mbox{E}}_2+</math><math>{\overline{\lambda }}_a^{-1/2}\boldsymbol{\mbox{e}}_3\otimes \boldsymbol{\mbox{E}}_3</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 9)
|}

el tensor de contracción sólo depende del alargamiento de contracción <math display="inline">{\overline{\lambda }}_a</math> . Aplicando la regla de la cadena, la variación de <math display="inline">{\overline{\boldsymbol{\mbox{F}}}}_a</math>  con respecto al tiempo [[#eq0045|(9)]]  viene dada como:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\dot{\overline{\boldsymbol{\mbox{F}}}}}_a=\frac{\partial {\overline{\boldsymbol{\mbox{F}}}}_a}{\partial {\overline{\lambda }}_a}{\dot{\overline{\lambda }}}_a</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 10)
|}

donde <math display="inline">{\dot{\overline{\lambda }}}_a</math>  es la velocidad de contracción en la dirección de comportamiento preferencial que debe ser conocida y <math display="inline">\frac{\partial {\overline{\boldsymbol{\mbox{F}}}}_a}{\partial {\overline{\lambda }}_a}</math>  puede obtenerse como:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{{\boldsymbol{\partial \overline{\mbox{F}}}}_a}{\partial {\overline{\lambda }}_a}=</math><math>\boldsymbol{\mbox{e}}_1\otimes \boldsymbol{\mbox{E}}_1-</math><math>\frac{1}{2}{\overline{\lambda }}_a^{-3/2}\boldsymbol{\mbox{e}}_2\otimes \boldsymbol{\mbox{E}}_2-</math><math>\frac{1}{2}{\overline{\lambda }}_a^{-3/2}\boldsymbol{\mbox{e}}_3\otimes \boldsymbol{\mbox{E}}_3</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 11)
|}

Se supone que las fibras de colágeno se mueven de forma conjunta con la matriz extracelular, por lo que, el alargamiento <math display="inline">\overline{\lambda }</math>  es definido como la relación entre las longitudes de las fibras en las configuraciones deformada y de referencia, es decir:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\overline{\lambda }\boldsymbol{\mbox{n}}(\boldsymbol{\mbox{x}},t)=\overline{\boldsymbol{\mbox{F}}}(X,t)\boldsymbol{\mbox{n}}_0(\boldsymbol{\mbox{X}})</math>
|-
|<math>{\overline{\lambda }}^2=\boldsymbol{\mbox{n}}_0.{\overline{\boldsymbol{\mbox{F}}}}^T\overline{\boldsymbol{\mbox{F}}}.\boldsymbol{\mbox{n}}_0=\boldsymbol{\mbox{n}}_0.\overline{\boldsymbol{\mbox{C}}}.\boldsymbol{\mbox{n}}_0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 12)
|}

donde <math display="inline">\overline{\boldsymbol{\mbox{C}}}=</math><math>{\overline{\boldsymbol{\mbox{F}}}}^T\overline{\boldsymbol{\mbox{F}}}</math>  es el tensor de deformación de Cauchy-Green por la derecha modificado. Podemos definir <math display="inline">{\overline{\boldsymbol{\mbox{C}}}}_e=</math><math>{\overline{\boldsymbol{\mbox{F}}}}_e^T{\overline{\boldsymbol{\mbox{F}}}}_e=</math><math>{\overline{\boldsymbol{\mbox{F}}}}_a^{-T}\overline{\boldsymbol{\mbox{C}}}\quad {\overline{\boldsymbol{\mbox{F}}}}_a^{-1}</math>  como la medida de deformación de los puentes cruzados y la titina (no es una variable de estado ya que depende de <math display="inline">\overline{\boldsymbol{\mbox{C}}}</math>  y <math display="inline">{\overline{\lambda }}_a</math> ), '''m'''  y '''n'''  son los vectores unitarios que definen la dirección de comportamiento preferencial (asociada a las fibras musculares y de colágeno, respectivamente) en la configuración deformada ([[#fig0005|fig. 1]] ).      

===2.2. Respuesta tensional hiperelástica===

Para caracterizar un proceso isotermo en materiales reversibles sin disipación de energía, se suele postular la existencia de una única representación de la función densidad energía de deformación ''Ψ''[[#bib0095|[19]]] , en este caso, expresada de forma desacoplada como suma de la energía volumétrica y desviadora, para evitar los problemas de incompresibilidad, empleando en la hipótesis cinemática [[#eq0015|(3)]]  y siguiendo a [[#bib0100|[20]]] . A su vez, para el tejido muscular, la energía de deformación desviadora se puede expresar de forma aditiva como la energía pasiva almacenada en el tejido y asociada principalmente al colágeno y elástina, <math display="inline">{\overline{\Psi }}_p</math> , más la energía almacenada en la fibra muscular, que denominaremos energía activa, <math display="inline">{\overline{\Psi }}_a</math> . La energía almacenada dependerá de '''F'''  y de las variables de estado <math display="inline">{\overline{\lambda }}_a,\beta </math>  y de los tensores estructurales '''M''' , '''N'''  definidos como '''M'''  = '''m'''<sub>0</sub>  ⊗ '''m'''<sub>0</sub>  y '''N'''  = '''n'''<sub>0</sub>  ⊗ '''n'''<sub>0</sub> :

<span id='eq0065'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\Psi ={\Psi }_{vol}(J)+{\overline{\Psi }}_p(\overline{\boldsymbol{\mbox{C}}},\boldsymbol{\mbox{N}})+{\overline{\Psi }}_a({\overline{\boldsymbol{\mbox{C}}}}_e,{\overline{\lambda }}_a,\boldsymbol{\mbox{M}},\beta )</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 13)
|}

con ''Ψ''<sub>''vol''</sub> , <math display="inline">{\overline{\Psi }}_p</math>  y <math display="inline">{\overline{\Psi }}_a</math>  las partes volumétrica y desviadora asociadas a las componentes pasiva y activa, respectivamente, de la función densidad de energía de deformación. El tensor <math display="inline">{\overline{\boldsymbol{\mbox{C}}}}_e</math>  no es una variable de estado ya que depende de <math display="inline">\overline{\boldsymbol{\mbox{C}}}</math>  y <math display="inline">{\overline{\lambda }}_a</math> , pero es necesario para formular la energía libre. ''β''  representa el nivel de activación con el tiempo o, lo que es lo mismo, la variación de la fuerza activa generada con el tiempo.      

La función densidad energía de deformación ''Ψ''  puede expresarse en función de los invariantes definidos por:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\overline{I}}_1=tr\boldsymbol{\boldsymbol{\overline{C}}},\quad {\overline{I}}_2=\frac{1}{2}((tr\boldsymbol{\boldsymbol{\overline{C}}})^2-tr{\boldsymbol{\boldsymbol{\overline{C}}}}^2)</math>
|-
|<math>{\overline{I}}_4=\boldsymbol{\mbox{n}}_0.\boldsymbol{\boldsymbol{\overline{C}}}\boldsymbol{\mbox{n}}_0,\quad {\overline{I}}_5=\boldsymbol{\mbox{n}}_0.{\boldsymbol{\boldsymbol{\overline{C}}}}^2.\boldsymbol{\mbox{n}}_0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 14)
|}

<math display="inline">{\overline{I}}_1</math>  y <math display="inline">{\overline{I}}_2</math>  son el primer y segundo invariante de <math display="inline">\overline{\boldsymbol{\mbox{C}}}</math>  y los pseudo-invariantes <math display="inline">{\overline{I}}_4,{\overline{I}}_5</math>  caracterizan la anisotropía asociada al comportamiento pasivo. <math display="inline">{\overline{I}}_4</math>  tiene un claro significado físico ya que se define como el alargamiento de la fibras de colágeno al cuadrado e <math display="inline">{\overline{I}}_5</math>  está asociado a la deformación transversal de las fibras de colágeno. De forma análoga se definirán los pseudo-invariantes asociados a las fibras musculares, donde el tensor de deformación <math display="inline">\boldsymbol{\boldsymbol{\overline{C}}}</math>  se ha sustituido por <math display="inline">{\boldsymbol{\boldsymbol{\overline{C}}}}_e</math> :

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\overline{J}}_4=\boldsymbol{\mbox{m}}_0.{\boldsymbol{\boldsymbol{\overline{C}}}}_e\boldsymbol{\mbox{m}}_0,\quad {\overline{J}}_5=</math><math>\boldsymbol{\mbox{m}}_0.{\boldsymbol{\boldsymbol{\overline{C}}}}_e^2.\boldsymbol{\mbox{m}}_0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 15)
|}

Es norma habitual en biomecánica omitir la dependencia de la función ''Ψ  ''  de los invariantes <math display="inline">{\overline{I}}_5,{\overline{J}}_5</math> , como consecuencia de la fuerte correlación de <math display="inline">{\overline{I}}_5</math>  con <math display="inline">{\overline{I}}_4</math>  y de <math display="inline">{\overline{J}}_5</math>  con <math display="inline">{\overline{J}}_4</math>[[#bib0105|[21]]] . Con lo cual se consigue disminuir el número de parámetros del material y, por lo tanto, facilitar el ajuste a los resultados experimentales. Adicionalmente, los invariantes <math display="inline">{\overline{I}}_4</math>  y <math display="inline">{\overline{J}}_4</math>  poseen una clara interpretación física. En lo que sigue, se planteará la formulación únicamente considerando los invariantes <math display="inline">{\overline{I}}_4</math>  y <math display="inline">{\overline{J}}_4</math> , con lo que la ecuación [[#eq0065|(13)]]  vendrá dada como:

<span id='eq0080'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\Psi ={\Psi }_{vol}(J)+{\overline{\Psi }}_p({\overline{I}}_1,{\overline{I}}_2,{\overline{I}}_4)+</math><math>{\overline{\Psi }}_a({\overline{J}}_4,{\overline{\lambda }}_a,\beta )</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 16)
|}

La respuesta en tensiones se puede obtener a partir de la ''Ψ''  teniendo en cuenta la desigualdad de Clausius-Planck:

<span id='eq0085'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>D_{int}=-\dot{\Psi }+\frac{1}{2}\boldsymbol{\mbox{S}}:\dot{\boldsymbol{\mbox{C}}}+</math><math>S_a{\dot{\overline{\lambda }}}_a+S_c\dot{\beta }\geq 0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 17)
|}

siendo ''S''<sub>''a''</sub>  la tensión generada por el golpe de fuerza en los puentes cruzados y ''S''<sub>''c''</sub>  por la variación del impulso en el tiempo. Teniendo en cuenta las ecuaciones [[#eq0080|(16)]]  y [[#eq0085|(17)]]  se obtiene:

<span id='eq0090'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left(\frac{1}{2}\boldsymbol{\mbox{S}}-\frac{\partial \Psi }{\partial \boldsymbol{\mbox{C}}}\right):\dot{\boldsymbol{\mbox{C}}}-\frac{\partial \Psi }{\partial \boldsymbol{\mbox{C}}_e}:{\boldsymbol{\dot{\mbox{C}}}}_e</math>
|-
|<math>+\left(S_a-\frac{\partial \Psi }{\partial {\overline{\lambda }}_a}\right){\dot{\overline{\lambda }}}_a+\left(S_c-\frac{\partial \Psi }{\partial \beta }\right)\dot{\beta }\geq 0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 18)
|}

Si suponemos que el primer término de la ecuación [[#eq0090|(18)]]  no depende de <math display="inline">{\dot{\overline{\lambda }}}_a</math>  y <math display="inline">\dot{\beta }</math> , el segundo tensor de tensiones de Piola-Kirchhoff vendrá dado por:

<span id='eq0095'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\mbox{S}}=2\frac{\partial \Psi }{\partial \boldsymbol{\mbox{C}}}</math>
|-
|<math>=\boldsymbol{\mbox{S}}_{vol}+{\boldsymbol{\overline{\mbox{S}}}}_p+{\boldsymbol{\overline{\mbox{S}}}}_a</math>
|-
|<math>=Jp\boldsymbol{\mbox{C}}^{-1}+J^{-\frac{2}{3}}(I-1/3\boldsymbol{\mbox{C}}^{-1}\otimes \boldsymbol{\mbox{C}}):({\tilde{\boldsymbol{\mbox{S}}}}_p)</math>
|-
|<math>+{\overline{\boldsymbol{\mbox{F}}}}_a^{-1}(J^{-\frac{2}{3}}(I-1/3\boldsymbol{\mbox{C}}_e^{-1}\otimes \boldsymbol{\mbox{C}}_e):({\tilde{\boldsymbol{\mbox{S}}}}_a)){\overline{\boldsymbol{\mbox{F}}}}_a^{-T}</math>
|-
|<math>=Jp\boldsymbol{\mbox{C}}^{-1}+J^{-\frac{2}{3}}DEV[{\tilde{\boldsymbol{\mbox{S}}}}_p]+{\overline{\boldsymbol{\mbox{F}}}}_a^{-1}(J^{-\frac{2}{3}}{DEV}_e[{\tilde{\boldsymbol{\mbox{S}}}}_a]){\overline{\boldsymbol{\mbox{F}}}}_a^{-T}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 19)
|}

con:

<span id='eq0100'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>p=\frac{d{\Psi }_{vol}(J)}{dJ}</math>
|-
|<math>{\tilde{\boldsymbol{\mbox{S}}}}_p=2\frac{\partial \overline{\Psi }({\overline{I}}_1,{\overline{I}}_2,{\overline{I}}_4)}{\partial \overline{\boldsymbol{\mbox{C}}}}\quad {\tilde{\boldsymbol{\mbox{S}}}}_a=2\frac{\partial \overline{\Psi }({\overline{J}}_4,{\overline{\lambda }}_a,\beta )}{\partial {\overline{\boldsymbol{\mbox{C}}}}_e}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 20)
|}

siendo ''p  ''  la presión hidrostática, <math display="inline">{\tilde{\boldsymbol{\mbox{S}}}}_p</math>  y <math display="inline">{\tilde{\boldsymbol{\mbox{S}}}}_a</math>  los tensores de tensiones de Piola-Kirchhoff modificados correspondientes a la parte pasiva y activa respectivamente y ''DEV''  el operador desviador en descripción material definido como:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>DEV[\cdot ]=(\cdot )-1/3(\boldsymbol{\mbox{C}}:(\cdot ))\boldsymbol{\mbox{C}}^{-1}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 21)
|}

El valor del segundo tensor de Piola-Kirchhoff, '''S''' , en función de los invariantes <math display="inline">{\overline{I}}_1</math> , <math display="inline">{\overline{I}}_2,{\overline{I}}_4</math>  y <math display="inline">{\overline{J}}_4</math>  vendrá dado por:

<span id='eq0110'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\mbox{S}}=Jp\boldsymbol{\mbox{C}}^{-1}</math>
|-
|<math>+2\left(\frac{\partial {\overline{\Psi }}_p}{\partial {\overline{I}}_1}+{\overline{I}}_1\frac{\partial {\overline{\Psi }}_p}{\partial {\overline{I}}_2}\right)\boldsymbol{\mbox{I}}-\frac{\partial {\overline{\Psi }}_p}{\partial {\overline{I}}_2}\boldsymbol{\mbox{C}}+{\overline{I}}_4\frac{\partial {\overline{\Psi }}_p}{\partial {\overline{I}}_4}\boldsymbol{\mbox{n}}_0\otimes \boldsymbol{\mbox{n}}_0</math>
|-
|<math>-\frac{1}{3}\left(\frac{\partial {\overline{\Psi }}_p}{\partial {\overline{I}}_1}{\overline{I}}_1+2\frac{\partial {\overline{\Psi }}_p}{\partial {\overline{I}}_2}{\overline{I}}_2+\frac{\partial {\overline{\Psi }}_p}{\partial {\overline{I}}_4}{\overline{I}}_4\right)\boldsymbol{\mbox{C}}^{-1}</math>
|-
|<math>+2{\overline{\boldsymbol{\mbox{F}}}}_a^{-1}\left[\left({\overline{J}}_4\frac{\partial {\overline{\Psi }}_a}{\partial {\overline{J}}_4}\right)\boldsymbol{\mbox{m}}_0\otimes \boldsymbol{\mbox{m}}_0-\frac{1}{3}\left(\frac{\partial {\overline{\Psi }}_a}{\partial {\overline{J}}_4}{\overline{J}}_4\right)\boldsymbol{\mbox{C}}_e^{-1}\right]{\overline{\boldsymbol{\mbox{F}}}}_a^{-T}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 22)
|}

El tensor de tensiones de Cauchy '''''σ'''''  es 1/''J''  veces el empuje (''push-forward'' ) de '''S'''  ('''''σ'''''  = ''J''<sup>−1</sup>'''''φ'''''<sub>*</sub> ('''''S''''' )), o en notación indicial, ''σ''<sub>''ij''</sub>  = ''J''<sup>−1</sup>''F''<sub>''iI''</sub>''F''<sub>''jJ''</sub>''S''<sub>''IJ''</sub> . Utilizando la relación entre el operador ''DEV'' [•], definido en la descripción material, y el operador ''dev'' [•], definido en la espacial  [[#bib0110|[22]]] :

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>J^{-2/3}\boldsymbol{\mbox{F}}(DEV[\bullet ])\boldsymbol{\mbox{F}}^T=</math><math>dev[\overline{\boldsymbol{\mbox{F}}}(\bullet ){\overline{\boldsymbol{\mbox{F}}}}^T]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 23)
|}

se obtiene la descripción espacial de [[#eq0110|(22)]] , es decir, el tensor de Cauchy. Para obtener la componente asociada al comportamiento activo, el empuje se ha de realizar con el tensor <math display="inline">{\overline{\boldsymbol{\mbox{F}}}}_e</math>  en lugar del tensor <math display="inline">\overline{\boldsymbol{\mbox{F}}}</math>  utilizado para la componente volumétrica y pasiva. Operando, se obtiene:

<span id='eq0120'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\sigma }=p\boldsymbol{1}</math>
|-
|<math>+\frac{2}{J}\left(\frac{\partial {\overline{\Psi }}_p}{\partial {\overline{I}}_1}+{\overline{I}}_1\frac{\partial {\overline{\Psi }}_p}{\partial {\overline{I}}_2}\right)\overline{\boldsymbol{\mbox{b}}}-\frac{\partial {\overline{\Psi }}_p}{\partial {\overline{I}}_2}{\overline{\boldsymbol{\mbox{b}}}}^2+{\overline{I}}_4\frac{\partial {\overline{\Psi }}_p}{\partial {\overline{I}}_4}\boldsymbol{\mbox{n}}\otimes \boldsymbol{\mbox{n}}</math>
|-
|<math>-\frac{1}{3}\left(\frac{\partial {\overline{\Psi }}_p}{\partial {\overline{I}}_1}{\overline{I}}_1+2\frac{\partial {\overline{\Psi }}_p}{\partial {\overline{I}}_2}{\overline{I}}_2+\frac{\partial {\overline{\Psi }}_p}{\partial {\overline{I}}_4}{\overline{I}}_4\right)\boldsymbol{1}</math>
|-
|<math>+{\overline{J}}_4\frac{\partial {\overline{\Psi }}_a}{\partial {\overline{J}}_4}\boldsymbol{\mbox{m}}\otimes \boldsymbol{\mbox{m}}-\frac{1}{3}\frac{\partial {\overline{\Psi }}_a}{\partial {\overline{J}}_4}{\overline{J}}_4\boldsymbol{1}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 24)
|}

con '''1'''  el tensor identidad de segundo orden, '''m'''  ⊗ '''m'''  y '''n'''  ⊗ '''n'''  los tensores estructurales asociados a las fibras musculares y de colágeno, respectivamente, en la configuración espacial.      

===2.3. Tensor de comportamiento elástico===

Para obtener la solución numérica del problema no lineal empleando una técnica iterativa tipo Newton, es necesario obtener la linealización de las ecuaciones constitutivas [[#bib0115|[23]]] . Conocido el segundo tensor de tensiones de Piola-Kirchhoff, '''S''' , en un punto, su variación con respecto al tensor de Cauchy-Green por la derecha '''C'''  puede escribirse de la forma:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>d\boldsymbol{\mbox{S}}=\mathbb{C}:\frac{1}{2}d\boldsymbol{\mbox{C}},\quad \mathbb{C}=</math><math>2\frac{\partial \boldsymbol{\mbox{S}}(\boldsymbol{\mbox{C}})}{\partial \boldsymbol{\mbox{C}}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 25)
|}

con <math display="inline">\mathbb{C}</math>''el tensor elástico''  en la configuración material  [[#bib0120|[24]]] . Partiendo de la ecuación [[#eq0095|(19)]]  se obtiene la contribución volumétrica y desviadora del tensor, esta última expresada como la suma de la componente pasiva y activa:

<span id='eq0130'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mathbb{C}={\mathbb{C}}_{vol}+{\overline{\mathbb{C}}}_p+</math><math>{\overline{\mathbb{C}}}_a=2\frac{\partial \boldsymbol{\mbox{S}}_{vol}}{\partial \boldsymbol{\mbox{C}}}+</math><math>2\frac{\partial {\boldsymbol{\overline{\mbox{S}}}}_p}{\partial \boldsymbol{\mbox{C}}}+</math><math>2\frac{\partial {\boldsymbol{\overline{\mbox{S}}}}_a}{\partial \boldsymbol{\mbox{C}}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 26)
|}

donde <math display="inline">{\mathbb{C}}_{vol}</math> , <math display="inline">{\overline{\mathbb{C}}}_p</math>  y <math display="inline">{\overline{\mathbb{C}}}_a</math>  pueden escribirse de la forma siguiente:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\mathbb{C}}_{vol}=2\boldsymbol{\mbox{C}}^{-1}\otimes \left(p\frac{\partial J}{\partial \boldsymbol{\mbox{C}}}+J\frac{\partial p}{\partial \boldsymbol{\mbox{C}}}+2Jp\frac{\partial \boldsymbol{\mbox{C}}^{-1}}{\partial \boldsymbol{\mbox{C}}}\right)</math>
|-
|<math>=J\tilde{p}\boldsymbol{\mbox{C}}^{-1}\otimes \boldsymbol{\mbox{C}}^{-1}-2Jp\boldsymbol{\mbox{C}}^{-1}\odot \boldsymbol{\mbox{C}}^{-1}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 27)
|}

Por conveniencia, se introduce la función escalar <math display="inline">\tilde{p}</math> , definida por:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\tilde{p}=p+J\frac{dp}{dJ}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 28)
|}

con la ecuación constitutiva para ''p''  dada por  [[#eq0100|(20)]] .

<span id='eq0145'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\overline{\mathbb{C}}}_p=-\frac{4}{3}J^{-\frac{4}{3}}\left(\frac{\partial {\overline{\Psi }}_p}{\partial \overline{\boldsymbol{\mbox{C}}}}\otimes {\overline{\boldsymbol{\mbox{C}}}}^{-1}+{\overline{\boldsymbol{\mbox{C}}}}^{-1}\otimes \frac{\partial {\overline{\Psi }}_p}{\partial \overline{\boldsymbol{\mbox{C}}}}\right)</math>
|-
|<math>+\frac{4}{3}J^{-\frac{4}{3}}\left(\frac{\partial {\overline{\Psi }}_p}{\partial \overline{\boldsymbol{\mbox{C}}}}:\overline{\boldsymbol{\mbox{C}}}\right)\left(I_C^{-1}-\frac{1}{3}{\overline{\boldsymbol{\mbox{C}}}}^{-1}\otimes {\overline{\boldsymbol{\mbox{C}}}}^{-1}\right)</math>
|-
|<math>+J^{-\frac{4}{3}}{\overline{\mathbb{C}}}_{\overline{w}}^p</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 29)
|}

con:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>(I_{C^{-1}})_{IJKL}=\frac{\partial C_{IJ}^{-1}}{\partial C_{KL}}=-(\boldsymbol{\mbox{C}}^{-1}\odot \boldsymbol{\mbox{C}}^{-1})_{IJKL}=</math>
|-
|<math>-\frac{1}{2}(C_{IK}^{-1}C_{JL}^{-1}+C_{IL}^{-1}C_{JK}^{-1})</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 30)
|}

y

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\overline{\mathbb{C}}}_{\overline{w}}^p=4\frac{{\partial }^2{\overline{\Psi }}_p}{\partial \overline{\boldsymbol{\mbox{C}}}\partial \overline{\boldsymbol{\mbox{C}}}}-\frac{4}{3}\left(\frac{{\partial }^2{\overline{\Psi }}_p}{\partial \overline{\boldsymbol{\mbox{C}}}\partial \overline{\boldsymbol{\mbox{C}}}}:\overline{\boldsymbol{\mbox{C}}}\right)\otimes {\overline{\boldsymbol{\mbox{C}}}}^{-1}</math>
|-
|<math>+{\overline{\boldsymbol{\mbox{C}}}}^{-1}\otimes \left(\frac{{\partial }^2{\overline{\Psi }}_p}{\partial \overline{\boldsymbol{\mbox{C}}}\partial \overline{\boldsymbol{\mbox{C}}}}:\overline{\boldsymbol{\mbox{C}}}\right)</math>
|-
|<math>+\frac{4}{9}\left(\overline{\boldsymbol{\mbox{C}}}:\frac{{\partial }^2{\overline{\Psi }}_p}{\partial \overline{\boldsymbol{\mbox{C}}}\partial \overline{\boldsymbol{\mbox{C}}}}:\overline{\boldsymbol{\mbox{C}}}\right){\overline{\boldsymbol{\mbox{C}}}}^{-1}\otimes {\overline{\boldsymbol{\mbox{C}}}}^{-1}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | 
|}

Para calcular la componente desviadora del tensor elástico asociada al comportamiento activo, se parte de la hipótesis de que las actualizaciones de <math display="inline">{\overline{\lambda }}_a</math>  se calcularán de forma explícita en cada incremento, por lo tanto <math display="inline">{\boldsymbol{\overline{\mbox{S}}}}_a</math>  no dependerá de <math display="inline">{\overline{\boldsymbol{\mbox{F}}}}_a</math> . Operando se obtiene:

<span id='eq0160'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\overline{\mathbb{C}}}_a=</math>
|-
|<math>{\overline{\boldsymbol{\mbox{F}}}}_a^{-1}\odot {\overline{\boldsymbol{\mbox{F}}}}_a^{-1}-\frac{4}{3}J^{-\frac{4}{3}}\left(\frac{\partial {\overline{\Psi }}_a}{\partial {\overline{\boldsymbol{\mbox{C}}}}_e}\otimes {\overline{\boldsymbol{\mbox{C}}}}_e^{-1}+{\overline{\boldsymbol{\mbox{C}}}}_e^{-1}\otimes {\frac{\partial {\overline{\Psi }}_a}{\partial \overline{\boldsymbol{\mbox{C}}}}}_e\right)</math>
|-
|<math>+\frac{4}{3}J^{-\frac{4}{3}}\left(\frac{\partial {\overline{\Psi }}_a}{\partial {\overline{\boldsymbol{\mbox{C}}}}_e}:{\overline{\boldsymbol{\mbox{C}}}}_e\right)\left(I_{\boldsymbol{C}_e^{-1}}-\frac{1}{3}{{\overline{\boldsymbol{\mbox{C}}}}_e}^{-1}\otimes {\overline{\boldsymbol{\mbox{C}}}}_e^{-1}\right)</math>
|-
|<math>+J^{-\frac{4}{3}}{\overline{\boldsymbol{\mbox{C}}}}_{e\overline{w}}^a{\overline{\boldsymbol{\mbox{F}}}}_a^{-T}\odot {\overline{\boldsymbol{\mbox{F}}}}_a^{-T}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 31)
|}

con:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\overline{\mathbb{C}}}_{e\overline{w}}^a=4\frac{{\partial }^2{\overline{\Psi }}_a}{\partial \boldsymbol{{\overline{\mbox{C}}}_\mbox{e}}\partial {\overline{\boldsymbol{\mbox{C}}}}_e}-\frac{4}{3}\left(\frac{{\partial }^2{\overline{\Psi }}_a}{\partial {\overline{\boldsymbol{\mbox{C}}}}_e\partial {\overline{\boldsymbol{\mbox{C}}}}_e}:{\overline{\boldsymbol{\mbox{C}}}}_e\right)\otimes {\overline{\boldsymbol{\mbox{C}}}}_e^{-1}</math>
|-
|<math>+{\overline{\boldsymbol{\mbox{C}}}}_e^{-1}\otimes \left(\frac{{\partial }^2{\overline{\Psi }}_a}{\partial {\overline{\boldsymbol{\mbox{C}}}}_e\partial {\overline{\boldsymbol{\mbox{C}}}}_e}:{\overline{\boldsymbol{\mbox{C}}}}_e\right)</math>
|-
|<math>+\frac{4}{9}\left({\overline{\boldsymbol{\mbox{C}}}}_e:\frac{{\partial }^2{\overline{\Psi }}_a}{\partial {\overline{\boldsymbol{\mbox{C}}}}_e\partial {\overline{\boldsymbol{\mbox{C}}}}_e}:{\overline{\boldsymbol{\mbox{C}}}}_e\right){\overline{\boldsymbol{\mbox{C}}}}_e^{-1}\otimes {\overline{\boldsymbol{\mbox{C}}}}_e^{-1}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 32)
|}

como puede observarse el valor de [[#eq0160|(31)]]  es idéntico a [[#eq0145|(29)]]  reemplazando <math display="inline">\overline{\boldsymbol{\mbox{C}}}</math>  por su componente elástica <math display="inline">{\overline{\boldsymbol{\mbox{C}}}}_e</math>  y posteriormente realizar la operación tirón (''pull-back'' ).      

El tensor elástico en la configuración espacial, denotado por <math display="inline">\mathbb{C}</math> , se define como el empuje de <math display="inline">\mathbb{C}</math>  escalado por el factor ''J''<sup>−1</sup> , esto es:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mathbb{C}=J^{-1}{\boldsymbol{\varphi }}_\ast (\mathbb{C}),\quad {\mathbb{C}}_{abcd}=</math><math>J^{-1}F_{aA}F_{bB}F_{cC}F_{dD}{\mathbb{C}}_{ABCD}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 33)
|}

Operando, se obtiene la expresión volumétrica y desviadora en la configuración deformada:

<span id='eq0175'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mathbb{C}={\mathbb{C}}_{vol}+{\mathbb{C}}_p+{\mathbb{C}}_a</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 34)
|}

siendo:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\mathbb{C}}_{vol}=(\tilde{p}\boldsymbol{1}\otimes \boldsymbol{1}-2pI),</math>
|-
|<math>{\overline{\mathbb{C}}}_p=\frac{2}{3}tr({\boldsymbol{\tilde{\sigma }}}_p)\mathbb{P}-\frac{2}{3}(\boldsymbol{1}\otimes {\boldsymbol{\overline{\sigma }}}_p+{\boldsymbol{\overline{\sigma }}}_p\otimes \boldsymbol{1})+{\mathbb{C}}_w^p</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 35)
|}

con:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\boldsymbol{\tilde{\sigma }}}_p=\frac{2}{J}\left[\left(\frac{\partial {\overline{\Psi }}_p}{\partial {\overline{I}}_1}+{\overline{I}}_1\frac{\partial {\overline{\Psi }}_p}{\partial {\overline{I}}_2}\right)\overline{\boldsymbol{\mbox{b}}}-\frac{\partial {\overline{\Psi }}_p}{\partial {\overline{I}}_2}{\overline{\boldsymbol{\mbox{b}}}}^2+{\overline{I}}_4\frac{\partial {\overline{\Psi }}_p}{\partial {\overline{I}}_4}\boldsymbol{\mbox{n}}\otimes \boldsymbol{\mbox{n}}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 36)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\boldsymbol{\overline{\sigma }}}_p=\frac{2}{J}\left(\frac{\partial {\overline{\Psi }}_p}{\partial {\overline{I}}_1}+{\overline{I}}_1\frac{\partial {\overline{\Psi }}_p}{\partial {\overline{I}}_2}\right)\overline{\boldsymbol{\mbox{b}}}-\frac{\partial {\overline{\Psi }}_p}{\partial {\overline{I}}_2}{\overline{\boldsymbol{\mbox{b}}}}^2+{\overline{I}}_4\frac{\partial {\overline{\Psi }}_p}{\partial {\overline{I}}_4}\boldsymbol{\mbox{n}}\otimes \boldsymbol{\mbox{n}}</math>
|-
|<math>-\frac{1}{3}\left(\frac{\partial {\overline{\Psi }}_p}{\partial {\overline{I}}_1}{\overline{I}}_1+2\frac{\partial {\overline{\Psi }}_p}{\partial {\overline{I}}_2}{\overline{I}}_2+\frac{\partial {\overline{\Psi }}_p}{\partial {\overline{I}}_4}{\overline{I}}_4\right)\boldsymbol{1}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 37)
|}

y <math display="inline">\mathbb{P}</math>  el tensor proyección:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mathbb{P}=I-\frac{1}{3}\boldsymbol{1}\otimes \boldsymbol{1}.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 38)
|}

Para calcular la parte desviadora de la componente activa en la configuración espacial es necesario realizar el empuje con el tensor <math display="inline">{\overline{\boldsymbol{\mbox{F}}}}_e</math> :

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\overline{\mathbb{C}}}_a=\frac{2}{3}tr({\boldsymbol{\tilde{\sigma }}}_a)\mathbb{P}-</math><math>\frac{2}{3}(\boldsymbol{1}\otimes {\boldsymbol{\overline{\sigma }}}_a+</math><math>{\boldsymbol{\overline{\sigma }}}_a\otimes \boldsymbol{1})+</math><math>{\mathbb{C}}_w^a</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 39)
|}

con:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\boldsymbol{\tilde{\sigma }}}_a=\frac{2}{J}\left[{\overline{J}}_4\frac{\partial {\overline{\Psi }}_a}{\partial {\overline{J}}_4}\boldsymbol{\mbox{m}}\otimes \boldsymbol{\mbox{m}}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 40)
|}

y

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\boldsymbol{\overline{\sigma }}}_a=\frac{2}{J}\left[{\overline{J}}_4\frac{\partial {\overline{\Psi }}_a}{\partial {\overline{J}}_4}\boldsymbol{\mbox{m}}\otimes \boldsymbol{\mbox{m}}-\right. </math><math>\left. \frac{1}{3}\frac{\partial {\overline{\Psi }}_a}{\partial {\overline{J}}_4}{\overline{J}}_4\boldsymbol{1}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 41)
|}

<math display="inline">{\overline{\mathbb{C}}}_{\overline{w}}^i</math> , ''i''  = ''p'' , ''a  '' , es el empuje de <math display="inline">{\overline{\mathbb{C}}}_{\overline{w}}^i</math>  y su valor también se puede expresar mediante:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\overline{\mathbb{C}}}_w^i=\mathbb{P}:{\overline{\mathbb{C}}}^i:\mathbb{P}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 42)
|}

Operando se puede llegar a la expresión del tensor constitutivo desviador <math display="inline">\overline{\mathbb{C}}</math>  expresado en función de los invariantes <math display="inline">{\overline{I}}_1,{\overline{I}}_2,{\overline{I}}_4</math>  y <math display="inline">{\overline{J}}_4</math> .      

===2.4. Imposición de deformaciones iniciales===

El músculo, al igual que la mayoría de los tejidos blandos, se encuentra sometido a deformaciones iniciales en la configuración de referencia. Para su incorporación en la formulación hiperelástica se sigue la metodología propuesta por [[#bib0125|[25]]] . Se supone que el tensor gradiente de deformación correspondiente al estado deformado '''F'''<sub>''r''</sub>  se puede expresar mediante una descomposición multiplicativa '''F'''<sub>''r''</sub>  = '''FF'''<sub>0</sub> , donde, '''F'''<sub>0</sub>  representa el gradiente de deformación correspondiente a las deformaciones iniciales y '''F'''  es el gradiente de deformación correspondiente a la aplicación de las cargas desde la configuración inicial ''Ω''<sub>0</sub>  (ver [[#fig0010|figura 2]] ).

<span id='fig0010'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_442579157-1-s2.0-S0213131512000193-gr2.jpg|center|282px|Descomposición multiplicativa del gradiente de deformación. Ωsf, estado libre de ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 2.

Descomposición multiplicativa del gradiente de deformación. ''Ω''<sub>''sf''</sub> , estado libre de tensiones (una vez diseccionado el tendón); ''Ω''<sub>0</sub> , estado de referencia o configuración libre de cargas (configuración fisiológica) y ''Ω''<sub>''t''</sub> , estado deformado actual.                  

</span>
|}

Como '''F'''<sub>0</sub>  no es fácil de evaluar, en este caso se supone que '''F'''<sub>0</sub>  corresponde al alargamiento longifudinal ''λ''<sub>0</sub>  en la dirección de las fibras musculares '''m'''<sub>0</sub>  en el estado de referencia ''Ω''<sub>0</sub>[[#bib0125|[25]]] . En un sistema local de coordenadas (<sup>*</sup> ) donde la dirección de las fibras musculares '''m'''<sub>0</sub>  está alineada con el eje ''X''<sub>1</sub> , <math display="inline">\boldsymbol{\mbox{F}}_0^\ast </math>  puede ser expresado por:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left[\boldsymbol{\mbox{F}}_0^\ast \right]=\left[\begin{array}{ccc}
{\lambda }_0 & 0 & 0\\
0 & \frac{1}{\sqrt{{\lambda }_0}} & 0\\
0 & 0 & \frac{1}{\sqrt{{\lambda }_0}}
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 43)
|}

y trasformado al sistema global mediante el tensor de rotación <math display="inline">\boldsymbol{\mbox{F}}_0=</math><math>\boldsymbol{\mbox{R}}\boldsymbol{\mbox{F}}_0^\ast \boldsymbol{\mbox{R}}^T</math>[[#bib0065|[13]]] .      

===2.5. Extensión al elemento membrana===

A la hora de abordar el modelado mediante elementos finitos de distintas unidades musculares, aparece la dificultad de que determinados tejidos presentan un espesor despreciable frente a los tejidos circundantes. Una posible alternativa sería despreciar su contribución mecánica, lo cual es factible siempre que su rigidez sea inferior a la de los otros tejidos. Sin embargo, no es posible despreciar la contribución de las estructuras que aportan rigidez, como es el caso de la fascia. En este caso se ha optado por modelar el tejido fascial mediante elementos membrana, por lo que las tensiones normales en la dirección del espesor se consideran despreciables. Por lo tanto, las seis componentes del tensor de tensiones considerado en la formulación 3D (ecuaciones [[#eq0110|(22)]]  y [[#eq0120|(24)]] ) se han de reducir en la formulación de membrana, por lo que la condición de tensiones nulas en el espesor del elemento ha de ser incorporada si se trabaja con la ley constitutiva planteada (ecuaciones [[#eq0130|(26)]]  y [[#eq0175|(34)]] ) para materiales hiperelásticos anisótropos. Este ha sido el esquema seguido en este trabajo, para el cual se ha empleado el algoritmo propuesto por [[#bib0130|[26]]] , que ha sido implementado en el código Abaqus [[#bib0135|[27]]]  mediante una subrutina UMAT.      

Para el modelado del comportamiento de la membrana, se utiliza un sistema de coordenadas local definido por tres vectores. ''e''<sub>1</sub>  y ''e''<sub>2</sub>  están incluidos en el plano de la membrana y ''e''<sub>3</sub>  es perpendicular a dicho plano. Estos vectores giran como lo hace el sólido rígido y las deformaciones se expresan en el sistema local. El gradiente de deformación '''F'''  en notación matricial, para el elemento membrana, viene dado por:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\mbox{F}}=\left(\begin{array}{ccc}
F_{11} & F_{12} & 0\\
F_{21} & F_{22} & 0\\
0 & 0 & F_{33}\\
 &  & 
\end{array}\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 44)
|}

por lo que las componentes ''C''<sub>33</sub>  y ''b''<sub>33</sub>  del tensor de deformaciones de Cauchy Green por la derecha y la izquierda no son nulas, así como las componentes ''S''<sub>33</sub>  y ''σ''<sub>33</sub> .      

Para comenzar con el algoritmo, la ley constitutiva 3D se expresa agrupando los tensores de tensión y deformación en los términos no nulos '''S'''<sub>''m''</sub>  = (''S''<sub>11</sub> , ''S''<sub>22</sub> , ''S''<sub>12</sub> , ''S''<sub>13</sub> , ''S''<sub>23</sub> )<sup>''T''</sup>  y los que tendrían que ser nulos '''S'''<sub>''z''</sub>  = (''S''<sub>33</sub> ) = 0.

<span id='eq0230'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left(\begin{array}{c}
d\boldsymbol{\mbox{S}}_m\\
d\boldsymbol{\mbox{S}}_z\\

\end{array}\right)=\left(\begin{array}{cc}
{\mathbb{C}}_{mm} & {\mathbb{C}}_{mz}\\
{\mathbb{C}}_{zm} & {\mathbb{C}}_{zz}\\
 & 
\end{array}\right)\left(\begin{array}{c}
d\boldsymbol{\mbox{C}}_m\\
d\boldsymbol{\mbox{C}}_z\\

\end{array}\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 45)
|}

Para desarrollar el algoritmo de imposición de tensión plana de forma local, se parte de la ecuación [[#eq0230|(45)]]  y se considera que '''S'''<sub>''z''</sub>  = 0 y '''C'''<sub>''z''</sub>  es la componente desconocida del tensor de deformación. Expresando la tensión en serie de Taylor:

<span id='eq0235'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\mbox{S}}_z^{\left(i+1\right)}=\boldsymbol{\mbox{S}}_z^{\left(i\right)}+</math><math>\frac{\partial \boldsymbol{\mbox{S}}_z^{\left(i\right)}}{\partial \boldsymbol{\mbox{C}}_z^{\left(i\right)}}\Delta \boldsymbol{\mbox{C}}_z+</math><math>\ldots \doteq 0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 46)
|}

donde el superíndice ''i''  indica el número de la iteración local. En las iteraciones sucesivas se modifica la componente en z del tensor de deformación de Cauchy Green por la derecha, '''C'''<sub>''z''</sub>  hasta alcanzar la condición de '''S'''<sub>''z''</sub>  = 0. Despreciando los términos de alto orden en la ecuación [[#eq0235|(46)]]  se obtiene <math display="inline">{\mathbb{C}}_{zz}^{\left(i\right)}=</math><math>\frac{\partial \boldsymbol{\mbox{S}}_z^{\left(i\right)}}{\partial \boldsymbol{\mbox{C}}_z^{\left(i\right)}}</math>  y con ello la deformación incremental:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\Delta \boldsymbol{\mbox{C}}_z=-[{\mathbb{C}}_{zz}^{\left(i\right)}]^{-1}\boldsymbol{\mbox{S}}_z^{\left(i\right)}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 47)
|}

con lo que la deformación en la siguiente iteración viene dada por:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\mbox{C}}_z^{\left(i+1\right)}=\boldsymbol{\mbox{C}}_z^{\left(i\right)}+</math><math>\Delta \boldsymbol{\mbox{C}}_z</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 48)
|}

La matriz de rigidez tangente, utilizada en cada iteración en el algoritmo de Newton-Raphson, podría estar asociada con la variación de '''S'''<sub>''m''</sub>  con respecto a '''C'''<sub>''m''</sub>  pero dependerá del estado de deformación completo '''C''' . Para obtener la matriz de rigidez tangente consistente, con la condición de tensión impuesta por la formulación del elemento, la matriz de rigidez tangente se ha de condensar. En este caso concreto, la condición es que la tensión sea nula en el espesor de la membrana. Si ''d'''''C'''<sub>''z''</sub>  = 0, de la segunda ecuación de [[#eq0230|(45)]]  se obtiene:

<span id='eq0250'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>d\boldsymbol{\mbox{C}}_z=-{\mathbb{C}}_{zz}^{-1}{\mathbb{C}}_{zm}d\boldsymbol{\mbox{C}}_m</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 49)
|}

Insertando [[#eq0250|(49)]]  en la primera ecuación de [[#eq0230|(45)]]  se obtiene

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>d\boldsymbol{\mbox{S}}_m={\mathbb{C}}_{psc}d\boldsymbol{\mbox{C}}_m</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 50)
|}

con:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\mathbb{C}}_{psc}=[{\mathbb{C}}_{mm}-{\mathbb{C}}_{mz}{\mathbb{C}}_{zz}^{-1}{\mathbb{C}}_{zm}]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 51)
|}

donde <math display="inline">{\mathbb{C}}_{psc}</math>  es la matriz de rigidez tangente correspondiente al elemento membrana.      

==3. Ejemplos==

Es esta sección se presentan diferentes resultados obtenidos de la simulación del músculo TA de rata que muestran tanto su respuesta pasiva como activa. Con el objetivo de validar el modelo numérico propuesto en la sección anterior, estas simulaciones se comparan con ensayos experimentales llevados a cabo sobre el mismo músculo.

===3.1. Modelo geométrico del músculo TA===

La geometría del músculo TA se reconstruyó a partir de la imagen axial de resonancia magnética (RM) de una rata adulta Wistar de 215 ± 15 g de peso tal y como se describe en [[#bib0080|[16]]] . Dicha geometría se discretizó con elementos hexaédricos utilizando el software ABAQUS [[#bib0135|[27]]] .      

Es necesario tener en cuenta que la geometría del músculo se ha obtenido «''in-vivo'' », es decir, en la configuración fisiológica, por lo que está sujeto a tensiones iniciales que es necesario introducir en el modelo numérico. En este caso se introduce un gradiente de deformación inicial, '''F'''<sub>0</sub> , definido por el acortamiento ''λ''<sub>0</sub>  = 1.2035 obtenido experimentalmente al diseccionar el tendón en la zona distal (4.54 ± 1.325 mm) [[#bib0080|[16]]] .      

===3.2. Función densidad energía de deformación===

Para caracterizar el comportamiento pasivo de la UMT, se utilizó la siguiente función densidad energía de deformación:

<span id='eq0265'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\overline{\Psi }}_p=c_1({\overline{I}}_1-3)+{\overline{\Psi }}_{pf}</math>
|-
|<math>{\overline{\Psi }}_{pf}=0\quad {\overline{I}}_4<{\overline{I}}_{4_0}</math>
|-
|<math>{\overline{\Psi }}_{pf}=\frac{c_3}{c_4}(e^{c_4({\overline{I}}_4-{\overline{I}}_{4_0})}-c_4({\overline{I}}_4-{\overline{I}}_{4_0})-1)\quad {\overline{I}}_4>{\overline{I}}_{4_0}\quad y\quad {\overline{I}}_4<{\overline{I}}_{4_{ref}}</math>
|-
|<math>{\overline{\Psi }}_{pf}=c_5\sqrt{{\overline{I}}_4}+\frac{1}{2}c_6ln({\overline{I}}_4)+c_7\quad {\overline{I}}_4>{\overline{I}}_{4_{ref}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 52)
|}

donde <math display="inline">{\overline{I}}_1</math>  representa el primer invariante del tensor de deformación de Cauchy-Green modificado, <math display="inline">{\overline{I}}_4>{\overline{I}}_{4_0}</math>  caracteriza la respuesta mecánica en la dirección de las fibras de colágeno, <math display="inline">{\overline{I}}_{4_{ref}}</math>  alargamiento de transición entre la zona lineal y exponencial. ''c''<sub>1</sub>  > 0, ''c''<sub>3</sub>  > 0, ''c''<sub>5</sub>  > 0 y ''c''<sub>6</sub>  < 0 son parámetros con unidades de tensión, siendo ''c''<sub>3</sub>  > 0 una escala de la parte exponencial, ''c''<sub>4</sub>  > 0 es un parámetro adimensional y ''c''<sub>7</sub>  es un parámetro con unidades de energía de deformación. Nótese que ''c''<sub>5</sub>  > 0, ''c''<sub>6</sub>  > 0 y ''c''<sub>7</sub>  > 0 no son parámetros independientes, sino que vienen determinados por las condiciones de continuidad de la energía, de la tensión y de la primera derivada de la tensión (en términos globales del tensor elástico tangente).      

Para caracterizar el comportamiento pasivo de la fascia se utilizó la función isótropa de Yeoh:

<span id='eq0270'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\overline{\Psi }}_p=c_{10}({\overline{I}}_1-3)+c_{20}({\overline{I}}_1-</math><math>3)^2+c_{30}({\overline{I}}_1-3)^3</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 53)
|}

donde ''c''<sub>10</sub>  > 0, ''c''<sub>20</sub>  > 0 y ''c''<sub>30</sub>  > 0 son parámetros con unidades de tensión. En las ecuaciones [[#eq0265|(52)]]  y [[#eq0270|(53)]]  solo se expresa la parte isocórica.      

Para caracterizar el comportamiento activo, proponemos la función <math display="inline">{\overline{\Psi }}_a</math>  expresada como una serie de funciones adimensionales que escalan la magnitud de la tensión isométrica máxima ''σ''<sub>0</sub>[[#bib0075|[15]]] . Estas funciones se obtienen de la relación tensión-alargamiento <math display="inline">f_1({\overline{\lambda }}_a)</math>  asociada al solapamiento efectivo entre los filamentos de miosina y actina en el elemento contráctil, de la relación tensión-velocidad <math display="inline">f_2({\dot{\overline{\lambda }}}_a)</math> , de la deformación asociada a los puentes cruzados <math display="inline">f_3({\overline{J}}_4)</math>  y del nivel de activación ''f'' (''f''<sub>''r''</sub> , ''t'' ):

<span id='eq0275'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\overline{\Psi }}_a={\sigma }_0f_1({\overline{\lambda }}_a)f_2({\dot{\overline{\lambda }}}_a)f_3({\overline{J}}_4)\beta (f_r,t)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 54)
|}

La función <math display="inline">f_1({\overline{\lambda }}_a)</math>  ha sido formulada por diferentes autores [[#bib0140|[28]]]  and [[#bib0050|[10]]] , en este trabajo utilizamos la función propuesta por [[#bib0075|[15]]] :

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>f({\overline{\lambda }}_a)={\int }_{0,5}^{1,5}e^{\frac{-({\overline{\lambda }}_a-{\lambda }_{opt})^2}{2(1-\alpha )^2}}d{\overline{\lambda }}_a</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 55)
|}

siendo ''λ''<sub>''opt''</sub>  el alargamiento óptimo del músculo y ''α  ''  corresponde a la curvatura de la función. La relación tensión-velocidad <math display="inline">f_2({\dot{\overline{\lambda }}}_a)</math>  tomará el valor unidad en el caso de contracciones isométricas, para el caso en el que varíe la longitud del músculo durante la contracción, se adopta la propuesta por [[#bib0050|[10]]] . La deformación asociada a los puentes cruzados:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>f_3({\overline{J}}_4)=\frac{2}{3}\left(2{\overline{J}}_4^{3/2}-\right. </math><math>\left. \frac{3}{2}{\overline{J}}_4+1\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 56)
|}

Finalmente la función de activación ''β'' (''f''<sub>''r''</sub> , ''t'' ) puede expresarse como  [[#bib0075|[15]]] :

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\beta (f_r,t)=\sum_{i=1}^n\left[\left(1-{re}^{(f_rT^{{'}})/c}\right)\left(P^{{'}}\frac{t-t_i}{T^{{'}}}e^{1-(\frac{t-t_i}{T^{{'}}})}\right)\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 57)
|}

siendo ''T'' ′ el tiempo de contracción aparente para todo el músculo generado en un pulso y ''P'' ′ es la amplitud aparente de la tensión generada en el mismo pulso, ''t''<sub>''i''</sub>  es el intervalo de tiempo entre estímulos que corresponde a 1/''f''<sub>''r''</sub> , ''r''  determina la tasa entre la tensión en un pulso y un tétano, ''c''  es la tasa de incremento de fuerza conforme aumenta la frecuencia y ''n''  es el número de estímulos.      

===3.3. Simulación ensayo uniaxial===

Para validar el modelo de comportamiento pasivo se ha reproducido numéricamente un ensayo de tracción uniaxial sobre la UMT de dos formas, considerando la inclusión de la fascia a través de elementos membrana y sin incluirla. Después, se comparan los resultados numéricos obtenidos con los experimentales [[#bib0070|[14]]] . En ambos casos, las fibras se consideran en la dirección fusiforme y el espesor de la fascia, medido previamente en ensayos experimentales [[#bib0145|[29]]] , se toma igual a 0.3 mm. Las constantes de la función densidad energía de deformación para el músculo, tendón y fascia se recogen en las tablas [[#tbl0005|1]]  y [[#tbl0010|2]] , respectivamente. Es necesario destacar aquí que, si bien los parámetros que caracterizan el comportamiento pasivo de músculo y tendón son bien conocidos de un trabajo previo [[#bib0070|[14]]] , aquellos que caracterizan la fascia han sido determinados ajustando la respuesta numérica a los ensayos experimentales. Para reproducir el ensayo se fijan los desplazamientos de la zona proximal del músculo TA y se aplica un desplazamiento impuesto de 10 mm en la zona distal (nodos correspondientes al tendón). Como el ensayo se realiza una vez que se ha producido la disección del tendón, se ha de partir en la simulación de la configuración libre de tensiones iniciales (''Ω''<sub>''sf''</sub> ), por lo que inicialmente se ha de obtener dicha configuración liberando las tensiones iniciales.

<span id='tbl0005'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabla 1.

Parámetros que caracterizan el comportamiento pasivo del músculo y tendón. Constantes en MPa [[#bib0070|[14]]] .                  

|-

! 
! ''c''<sub>1</sub>
! ''c''<sub>3</sub>
! ''c''<sub>4</sub>
! ''c''<sub>5</sub>
! ''c''<sub>6</sub>
! ''c''<sub>7</sub>
! <math display="inline">{\overline{I}}_{4_0}</math>
! <math display="inline">{\overline{I}}_{4_{ref}}</math>
|-

| Músculo
| 0,001
| 0,053915
| 0,782802
| 5,742780
| −9,035709
| −4,875961
| 1,254400
| 3,189796
|-

| Tendón
| 0,01
| 0,054292
| 6,860021
| 57,738461
| −66,243575
| −57,078409
| 1,0
| 1,44
|}

<span id='tbl0010'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabla 2.

Parámetros que caracterizan el comportamiento pasivo de la fascia. Constantes en MPa.

|-

! 
! ''D''
! ''C''<sub>10</sub>
! ''C''<sub>20</sub>
! ''C''<sub>30</sub>
|-

| Fascia
| 0,01
| 0,008
| 0,05
| 1,0
|}

La [[#fig0015|figura 3]] a presenta el modelo de elemenos finitos para el músculo TA, donde se ha indicado una región central del vientre muscular sobre la que se visualizan los campos de tensiones máximas principales (TMP). Las TMP que se muestran en la [[#fig0015|figura 3]] b se corresponden con la simulación del modelo sin incluir la fascia mientras que la [[#fig0015|figura 3]] c incluye este tejido. Para establecer una comparación directa entre ambas imágenes, los elementos correspondientes a la fascia han sido suprimidos en la [[#fig0015|figura 3]] .c. Se observa que la distribución de TMP evidencia valores mayores en la [[#fig0015|figura 3]] b, donde no está considerada la fascia.

<span id='fig0015'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_442579157-1-s2.0-S0213131512000193-gr3.jpg|center|360px|(a) Modelo de elementos finitos del TA. En azul se indica la parte central del ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 3.

(a) Modelo de elementos finitos del TA. En azul se indica la parte central del músculo, donde se visualizan los campos de TMP. (b) Campo de TMP en el músculo cuando no se incluye la fascia. (c) Campo de TMP en el músculo cuando se considera la fascia. La fascia ha sido retirada para visualizar solamente el músculo y poder comparar con (b). (d) Fuerza vs. alargamiento obtenido de las simulaciones del ensayo uniaxial, con y sin fascia, comparadas con los resultados experimentales [[#bib0070|[14]]] .                  

</span>
|}

Las relaciones fuerza-desplazamiento obtenidas en las simulaciones numéricas se muestran en la [[#fig0015|figura 3]] .d donde, a su vez, se comparan con las curvas experimentales. Aquí se aprecia la diferencia de rigidez entre ambas situaciones, evidenciando una rigidez notablemente mayor para el caso en el que se incluye la fascia.      

===3.4. Simulación de la contracción===

Se han llevado a cabo 2 grupos diferentes de simulaciones para estudiar la contracción isométrica (la longitud del músculo es constante) y la contracción concéntrica (la longitud del músculo disminuye) a velocidad constante (isotónica). La [[#tbl0015|tabla 3]]  recoge los parámetros utilizados para las simulaciones.

<span id='tbl0015'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabla 3.

Parámetros que caracterizan el comportamiento activo del músculo [[#bib0075|[15]]] .                  

|-

! ''σ''<sub>0</sub>  (MPa)                                                    
! ''λ''<sub>''opt''</sub>
! ''α''
! ''T'' ′ (s)                                                    
! ''P'' ′ (N)                                                    
! ''t''<sub>''i''</sub>  (s)                                                    
! ''r''
! ''c''
|-

| 0,8
| 1
| 0,83616
| 0,04
| 0,11
| 0,01667
| 1,0535
| 1,1245
|}

En primer lugar, para estudiar la contracción isométrica, se fijan ambos extremos de la malla manteniendo fijos los grados de libertad de los nodos en dichas regiones. Puesto que en una contracción isométrica, la unidad mínima funcional del músculo (sarcómero) no altera su longitud inicial, no existe desplazamiento relativo entre los filamentos de actina y miosina. Es por ello que, en la descomposición multiplicativa del gradiente de deformación <math display="inline">{\overline{\boldsymbol{\mbox{F}}}}_a=</math><math>\boldsymbol{1}</math> , por lo tanto <math display="inline">\overline{\boldsymbol{\mbox{F}}}=</math><math>{\overline{\boldsymbol{\mbox{F}}}}_e</math> . El valor de la función <math display="inline">f({\dot{\overline{\lambda }}}_a)</math>  en la ecuación [[#eq0275|(54)]]  al no existir dependencia con la velocidad se toma como <math display="inline">f({\dot{\overline{\lambda }}}_a)=</math><math>1</math> .      

La [[#fig0020|figura 4]]  muestra la fuerza de reacción total registrada en los nodos de la malla del extremo superior para cuatro frecuencias del estímulo. El efecto que la frecuencia del estímulo ejerce sobre la deformación y la TMP del tejido se observa en la [[#fig0025|figura 5]] .

<span id='fig0020'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_442579157-1-s2.0-S0213131512000193-gr4.jpg|center|338px|Fuerza de reacción isométrica para cuatro frecuencias de estímulo.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 4.

Fuerza de reacción isométrica para cuatro frecuencias de estímulo.

</span>
|}

<span id='fig0025'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_442579157-1-s2.0-S0213131512000193-gr5.jpg|center|320px|(a) Configuración inicial de la UMT. TMP alcanzada para el estímulo a (b) 10Hz. ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 5.

(a) Configuración inicial de la UMT. TMP alcanzada para el estímulo a (b) 10Hz. (c) 30Hz. (d) 60Hz. (e) 90Hz. Nota: tensiones por encima de 2, 0 MPa se alcanzan en la región del tendón.

</span>
|}

Por otro lado, para llevar a cabo la simulación de una contracción concéntrica, se liberan las condiciones de contorno establecidas en el caso anterior para los nodos del extremo superior de la malla. De esta manera, el músculo acortará su longitud inicial libremente.

Bajo estas condiciones, y realizando sucesivos incrementos de tiempo en un análisis estático, es necesario definir la velocidad de la contracción. En las simulaciones llevadas a cabo en este trabajo se ha considerado que la contracción se realiza a velocidad constante (contracción isotónica) tomando un valor <math display="inline">{\dot{\overline{\lambda }}}_a=</math><math>-5</math>  [1/s] [[#bib0170|[34]]]  por ello, <math display="inline">{\overline{\boldsymbol{\mbox{F}}}}_a</math>  se supone conocida y, en la dirección de la fibra muscular tendrá la siguiente forma, con <math display="inline">{\overline{\lambda }}_a^{n+1}=</math><math>{\overline{\lambda }}_a^n+{\dot{\overline{\lambda }}}_a\vartriangle t</math> :

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left[\boldsymbol{\mbox{F}}_a\right]=\left[\begin{array}{ccc}
{\overline{\lambda }}_a & 0 & 0\\
0 & \frac{1}{\sqrt{\overline{\lambda }a}} & 0\\
0 & 0 & \frac{1}{\sqrt{\overline{\lambda }a}}
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 58)
|}

En esta situación, el valor de <math display="inline">f({\dot{\overline{\lambda }}}_a)</math>  se ha tomado como <math display="inline">f({\dot{\overline{\lambda }}}_a)=</math><math>0.5</math>[[#bib0050|[10]]] .      

La [[#fig0030|figura 6]]  representa el desplazamiento del extremo superior de la UMT para las cuatro diferentes frecuencias consideradas. La [[#fig0035|figura 7]]  representa, junto con la configuración inicial del tejido, la deformada y el campo de desplazamientos en la dirección longitudinal para los 4 estímulos. Se observa cómo para la frecuencia mayor se obtiene el mayor acortamiento del músculo.

<span id='fig0030'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_442579157-1-s2.0-S0213131512000193-gr6.jpg|center|334px|Desplazamiento del extremo superior de la UMT para 4 frecuencias de estímulo.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 6.

Desplazamiento del extremo superior de la UMT para 4 frecuencias de estímulo.

</span>
|}

<span id='fig0035'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_442579157-1-s2.0-S0213131512000193-gr7.jpg|center|376px|Campo de desplazamientos en la dirección longitudinal del músculo para una ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 7.

Campo de desplazamientos en la dirección longitudinal del músculo para una contracción concéntrica a 4 diferentes frecuencias del estímulo. (a) Configuración indeformada. (b) 10Hz. (c) 30Hz. (d) 60Hz. (e) 90Hz.

</span>
|}

==4. Conclusiones==

En este trabajo se ha presentado y validado un modelo numérico de comportamiento 3D que permite simular la respuesta mecánica del tejido músculo-esquelético de mamíferos. Esta respuesta está fuertemente influenciada por la geometría y la arquitectura de las fibras musculares [[#bib0150|[30]]] . Por ello, en este trabajo se ha modelado el músculo TA partiendo de imagen de RM considerando una distribución de dichas fibras siguiendo la geometría fusiforme del tejido. En la mayoría de los trabajos de la literatura se parte de geometrías de músculos procedentes de cadáveres [[#bib0155|[31]]]  and [[#bib0050|[10]]]  o bien utilizan geometrías idealizadas [[#bib0045|[9]]] , [[#bib0140|[28]]] , [[#bib0050|[10]]]  and [[#bib0160|[32]]] .      

El modelo numérico propuesto se ha particularizado para el músculo TA de rata. Para ello, las constantes que definen el comportamiento pasivo y activo se han tomado de trabajos experimentales previos desarrollados por los autores [[#bib0070|[14]]]  and [[#bib0075|[15]]] . El comportamiento pasivo se ha validado comparando los resultados numéricos con los experimentales cuando a la UMT se le solicita a un estado uniaxial de tracción ([[#fig0015|fig. 3]] d). Este procedimiento ha permitido, a través de las curvas fuerza-desplazamiento obtenidas experimentalmente, ajustar un modelo de comportamiento isótropo para la fascia que envuelve el tejido muscular. De esta forma se ha conseguido estimar las constantes que definen su comportamiento indirectamente a través de la simulación computacional. Este hecho permite caracterizar de manera precisa el tejido evitando llevar a cabo ensayos experimentales de extrema dificultad técnica. La [[#fig0015|figura 3]]  muestra además el efecto que ejerce la fascia sobre la distribución de tensiones en el vientre muscular. Este tejido, debido a que presenta una rigidez mayor, disminuye la tensión que soporta el músculo para un mismo nivel de extensión pasiva [[#bib0165|[33]]] .      

En el caso de la simulación del comportamiento activo, la fuerza de reacción obtenida para la contracción isométrica ([[#fig0020|figura 4]] ) coincide con la observada experimentalmente en [[#bib0075|[15]]]  para las frecuencias de estímulo de 30, 60 y 90 Hz. Una limitación importante del modelo bajo contracción isométrica consiste en haber considerado <math display="inline">{\overline{\boldsymbol{\mbox{F}}}}_a=</math><math>\boldsymbol{1}</math> . Si bien es cierto que para este tipo de contracción no hay desplazamiento de los extremos de la UMT, las fibras musculares se acortan hasta alcanzar una posición de equilibrio en la dirección formada por los dos extremos fijos. En el caso de contracciones concéntricas los resultados presentan la misma tendencia obtenida por otros autores para geometrías diferentes [[#bib0045|[9]]] , [[#bib0050|[10]]]  and [[#bib0170|[34]]] . La simplificación llevada a cabo de no considerar la masa y efectos inerciales del tejido así como la ausencia de cargas exteriores condiciona que la simulación se acerque a la realidad del comportamiento de la UMT.

==Agradecimientos==

Esta investigación ha sido financiada por el Ministerio de Economía y Competitividad de España mediante el proyecto DPI2011-27939-C02-01 y el Instituto de Salud Carlos III (ISCIII) mediante la iniciativa CIBER. Los autores también agradecen al Ministerio de Ciencia e Innovación por la beca (BES-2009-021515) concedida a B. Hernández-Gascón.

==References==

<ol style='list-style-type: none;margin-left: 0px;'><li><span id='bib0005'></span>
[[#bib0005|[1]]] B.R. MacIntosh, P.F. Gardiner, A.J. McComas; Skeletal Muscle Form and Function;  (2nd edition)Human Kinetics (2006)</li>
<li><span id='bib0010'></span>
[[#bib0010|[2]]] A.V. Hill; The heat of shortening and the dynamic constants of muscle; P. ROY. SOC. B-BIOL. SCI., Vol. 126 (1938), pp. 136–195</li>
<li><span id='bib0015'></span>
[[#bib0015|[3]]] A. Magid, D.J. Law; Myofibrils bear most of the resting tension in frog skeletal muscle; Science, 230 (4731) (1985), pp. 1280–1282</li>
<li><span id='bib0020'></span>
[[#bib0020|[4]]] J.A.C. Martins, E.B. Pires, R. Salvado, P.B. Dinis; Numerical model of passive and active behavior of skeletal muscles; Comput. Methods Appl. Mech. Eng., 151 (1998), pp. 419–433</li>
<li><span id='bib0025'></span>
[[#bib0025|[5]]] B. Hernández-Gascón, E. Peña, G. Pascual, M. Rodríguez, B. Calvo, M. Doblaré, J.M. Bellón; Mechanical and histological characterization of the abdominal muscle. a previous step to modelling hernia surgery; J. Mech. Behav. Biomed. Mater., 4 (3) (2011), pp. 392–404</li>
<li><span id='bib0030'></span>
[[#bib0030|[6]]] J.W. Fernandez, M.L. Buist, D.P. Nickerson, P.J. Hunter; Modelling the passive and nerve activated response of the rectus femoris muscle to a flexion loading: a finite element framework; Med. Eng. Phys., 27 (10) (2005), pp. 862–870</li>
<li><span id='bib0035'></span>
[[#bib0035|[7]]] T. Johansson, P. Meier, R. Blickhan; A finite-element model for the mechanical analysis of skeletal muscles; J. Theor. Biol., 206 (1) (2000), pp. 131–149</li>
<li><span id='bib0040'></span>
[[#bib0040|[8]]] C.A. Yucesoy, B.H.F.J.M. Koopman, P.A. Huijing, H.J. Grootenboer; Three-dimensional finite element modeling of skeletal muscle using a two-domain approach: linked fiber-matrix mesh model; J. Biomech., 35 (9) (2002), pp. 1253–1262</li>
<li><span id='bib0045'></span>
[[#bib0045|[9]]] C.P. Tsui, C.Y. Tang, C.P. Leung, K.W. Cheng, Y.F. Ng, D.H. Chow, C.K. Li; Active finite element analysis of skeletal muscle-tendon complex during isometric, shortening and lengthening contraction; Bio-Med. Mater. Eng., 14 (2004), pp. 271–279</li>
<li><span id='bib0050'></span>
[[#bib0050|[10]]] M. Böl, S. Reese; Micromechanical modelling of skeletal muscles based on the finite element method; Comput. Methods Biomech. Biomed. Eng., 11 (2008), pp. 489–504</li>
<li><span id='bib0055'></span>
[[#bib0055|[11]]] S.-W. Chi, J. Hodgson, J.-S. Chen, V. Reggie Edgerton, D.D. Shin, R.A. Roiz, S. Sinha; Finite element modeling reveals complex strain mechanics in the aponeuroses of contracting skeletal muscle; J. Biomech., 43 (7) (2010), pp. 1243–1250</li>
<li><span id='bib0060'></span>
[[#bib0060|[12]]] J.C. Gardiner, J.A. Weiss, T.D. Rosenberg; Strain in the human medial collateral ligament during valgus lading of the knee; Clin. Orthop. Relat. Res., 391 (2001), pp. 266–274</li>
<li><span id='bib0065'></span>
[[#bib0065|[13]]] E. Peña, M.A. Martinez, B. Calvo, M. Doblaré; On the numerical treatment of initial strains in biological soft tissues; Int. J. Numer. Methods Engin., 68 (8) (2006), pp. 836–860</li>
<li><span id='bib0070'></span>
[[#bib0070|[14]]] B. Calvo, A. Ramŕez, A. Alonso, J. Grasa, F. Soteras, R. Osta, M.J. Muñoz; Passive nonlinear elastic behavior of skeletal muscle: Experimental results and model formulation; J. Biomech., 43 (2) (2010), pp. 318–325</li>
<li><span id='bib0075'></span>
[[#bib0075|[15]]] A. Ramírez, J. Grasa, A. Alonso, F. Soteras, R. Osta, M. Muñoz, B. Calvo; Active response of skeletal muscle: In vivo experimental results and model formulation; J. Theor. Biol., 267 (4) (2010), pp. 546–553</li>
<li><span id='bib0080'></span>
[[#bib0080|[16]]] J. Grasa, A. Ramírez, R. Osta, M.J. Muñoz, F. Soteras, B. Calvo; A 3D active-passive numerical skeletal muscle model incorporating initial tissue strains. validation with experimental results on rat tibialis anterior muscle; Biomech. Model. Mechanobio., 10 (2011), pp. 779–787</li>
<li><span id='bib0085'></span>
[[#bib0085|[17]]] P.J. Flory; Thermodynamic relations for high elastic materials; Trans. Faraday Soc., 57 (1961), pp. 829–838</li>
<li><span id='bib0090'></span>
[[#bib0090|[18]]] J. Staalhand, A. Klarbring, G.A. Holzapfel; A mechanochemical 3d continuum model for smooth muscle contraction under finite strains; J. Theor. Biol., 268 (1) (2011), pp. 120–130</li>
<li><span id='bib0095'></span>
[[#bib0095|[19]]] J.C. Simo, R.L. Taylor; Quasi-incompresible finite elasticity in principal stretches. continuum basis and numerical algorithms; Comput. Methods Appl. Mech. Eng., 85 (1991), pp. 273–310</li>
<li><span id='bib0100'></span>
[[#bib0100|[20]]] A. Spencer; Continuum Physics. Theory of Invariants; Academic Press, New York (1954)</li>
<li><span id='bib0105'></span>
[[#bib0105|[21]]] G.A. Holzapfel; Nonlinear Solid Mechanics; A continuum approach for engineering, Wiley (2000)</li>
<li><span id='bib0110'></span>
[[#bib0110|[22]]] J.C. Simo, T.J.R. Hughes; Computational Inelasticity; Springer-Verlag, New York (1998)</li>
<li><span id='bib0115'></span>
[[#bib0115|[23]]] T.J.R. Hughes, K.S. Pister; Consistent linearization in mechanics of solids and structures; Comput. Struct., 8 (1978), pp. 391–397</li>
<li><span id='bib0120'></span>
[[#bib0120|[24]]] R.W. Ogden; Non-linear Elastic Deformations; Dover, New York (1996)</li>
<li><span id='bib0125'></span>
[[#bib0125|[25]]] J.C. Gardiner, J.A. Weiss; Subject-specific finite element analysis of the human medial collateral ligament during valgus knee loading.; J. Orthop. Res., 21 (6) (2003), pp. 1098–1106</li>
<li><span id='bib0130'></span>
[[#bib0130|[26]]] S. Klinkel, S. Govindjee; Using finite strain 3d-material models in beam and shell elements; Eng. Comput., 19 (8) (2002), pp. 902–921</li>
<li><span id='bib0135'></span>
[[#bib0135|[27]]] Hibbit, Karlsson, Sorensen, Abaqus users guide, v. 6.5, HKS inc. Pawtucket, RI, USA., 2005.</li>
<li><span id='bib0140'></span>
[[#bib0140|[28]]] B.B. Blemker, P.M. Pinsky, S.L. Delp; A 3d model of muscle reveals the causes of nonuniform strains in the biceps brachii; Ann. Biomed. Eng., 38 (2005), pp. 657–665</li>
<li><span id='bib0145'></span>
[[#bib0145|[29]]] A. Ramírez, Modelado y simulación del tejido músculo-esquelético. validación experimental con el músculo tibial anterior de rata, Ph.D. thesis (2011).</li>
<li><span id='bib0150'></span>
[[#bib0150|[30]]] B.B. Blemker, S.L. Delp; Three-dimensional representation of complex muscle architectures and geometries; Ann. Biomed. Eng., 33 (2005), pp. 661–673</li>
<li><span id='bib0155'></span>
[[#bib0155|[31]]] D. d’Aulignac, J.A. Martins, E.B. Pires, T. Mascarenhas, R.M. Jorge; A shell finite element model of the pelvic floor muscles; Comput. Methods Biomech. Biomed. Eng., 8 (5) (2005), pp. 339–347</li>
<li><span id='bib0160'></span>
[[#bib0160|[32]]] C.Y. Tang, G. Zhang, C.P. Tsui; A 3d skeletal muscle model coupled with active contraction of muscle fibres and hyperelastic behaviour.; J. Biomech., 42 (7) (2009), pp. 865–872</li>
<li><span id='bib0165'></span>
[[#bib0165|[33]]] J.M. Rijkelijkhuizen, H.J.M. Meijer, G.C. Baan, P.A. Huijing; Myofascial force transmission also occurs between antagonistic muscles located within opposite compartments of the rat lower hind limb; J. Electromyogr. Kinesiol., 17 (6) (2007), pp. 690–697</li>
<li><span id='bib0170'></span>
[[#bib0170|[34]]] P. Meier, R. Blickhan, Skeletal muscle mechanics: from mechanisms to function, Ch. FEM-Simulation of Skeletal Muscle: The Influence of Inertia During Activation and Deactivation, 1st ed. New York, 2000, pp. 207-223.</li>
</ol>