FISIOLOGIA E BIOMECÂNICA DA ATIVIDADE
MOTORA – AVALIAÇÃO E REABILITAÇÃO
Modelos Matemáticos de Avaliação
Dr. Luciano Luporini Menegaldo
E-mail [email protected]
Agrupamento de Sistemas de Controle
Divisão de Mecânica e Eletricidade
Instituto de Pesquisas Tecnológicas do
Estado de São Paulo - IPT
1
Tópicos a serem abordados:
1) O que é um modelo
2) Tipos de modelos matemáticos em
sistemas biomecânicos
3) Pesquisa em modelagem e simulação
de sistemas músculo-esqueléticos: um
exemplo
2
1) Modelos matemáticos
O que é um modelo?
‣




Uma representação intelectual de um
sistema real
Exemplos de modelos:
Representação dos átomos na química
Leis de Newton
Lei dos gases perfeitos
Fisiopatologia
3
Um bom modelo:
Funciona? Até que funciona, na
maior parte dos casos.
É um conhecimento verdadeiro? É.
Não é mentira.
É um conhecimento sempre
aperfeiçoável? Sem dúvida.
4
Modelo Matemático
É um tipo de modelo em que:
1) Propõe-se um sistema físico
equivalente ao sistema real
2) Do sistema físico equivalente se
estabelecem equações capazes de
descrever o comportamento desse
sistema
5
Exemplos:
1) Modelos da mecânica respiratória
2) Modelos da contração muscular
3) Modelos de carregamento de
ossos
6
Multiplicidade de modelos:
Sistemas reais admitem múltiplos
modelos
Modelos são bons ou ruins para fins
específicos
Exemplo: modelos de corpos rígidos para
controle de movimento e para análise
de tensões em ossos.
7
O que se pode fazer com um modelo?
Por que fazer modelos?
– quantificar relações e comportamentos
– projetar intervenções cirúrgicas
– projetar tratamentos fisioterápicos
(transferência de calor, exercícios etc.)
8
– projetar tratamentos farmacológicos
(quimioterapia, tratamentos de distúrbios
neurológicos etc.)
– projetar dispositivos tecnológicos
(respiradores artificiais, próteses
anatômicas e neurais, órgãos artificiais,
trajes etc.)
9
Tipos de modelos em sistemas
biomecânicos
• Modelos estáticos
– Modelos de distribuição de forças
– Modelos de elementos finitos
• Modelos cinemáticos
– Modelos para laboratório de marcha
(cálculo de ângulos articulares em função
das coordenadas dos marcadores
– Modelos geométricos da anatomia
10
• Modelos dinâmicos
– Dinâmica direta e dinâmica inversa
– Cálculo de momentos articulares
• Modelos lineares e não-lineares
• Modelos da mecânica muscular
11
Modelos estáticos
• Modelos de distribuição de forças
12
Modelos de elementos finitos
Geração da malha
13
Carregamento
seção transversal da tíbia
14
Análise
de tensões
15
Modelos cinemáticos
Modelos para laboratório de marcha (cálculo de ângulos
articulares em função das coordenadas dos marcadores
16
Modelos geométricos e
antropométricos
17
Modelos dinâmicos
•
Dinâmica direta e dinâmica inversa
•
Cálculo de momentos articulares
• Modelos lineares e não-lineares
• Modelos da mecânica muscular
18
2a. Parte: Pesquisa em modelagem de
sistemas biomecânicos
Biomecânica e controle da postura humana
Objetivo:
• Calcular sinais de excitação neuromuscular
• capazes de levantar o corpo humano
desde uma posição semi-agachada até
a postura ereta
• minimizando uma função de custo
19
Com que utilidade?
• Pesquisas básicas em teoria de controle
motor
• Estimar o efeito biomecânico e motor de
procedimentos cirúrgicos
• Determinar as estratégias ótimas de controle
motor que deveriam ser empregadas pelo
SNC depois de uma cirurgia, e sugerir
procedimentos de fisioterapia
• Determinar padrões ótimos de ativação para
Estimulação Elétrica Funcional (FES)
20
Modelo biomecânico
a) Sistema de múltiplos corpos rígidos
21
b) Modelo geométrico do membro inferior
Modelo de domínio público desenvolvido por Scott Delp (Univ.
Stanford), utilizado no SIMM (Musculographics Inc.), com 40
músculos e 5 articulações
22
Articulação
Quadril
Joelho
Tornozelo,
subtalar e
metatarsofalangeal
Características
Junta esférica com três graus de liberdade (Delp
et al., 1990)
Junta plana de um grau de liberdade, baseada em
Yamaguchi e Zajac (1989). Leva em conta o
ponto de contato variável da articulação tíbiofemoral, a cinemática patelo-femoral e o efeito do
aumento do braço de momento do tendão do
quadríceps (Delp et al., 1990)
Modeladas como juntas tipo pino (com um grau
de liberdade), com eixos localizados e orientados
segundo Inman (1976) com pequenas alterações
na orientação da articulação metatarso-falangeal
(Delp et al., 1990)
23
Referenciais posicionados e dimensionados segundo
acidentes anatômicos
Transformações cinemáticas: p / referencial inercial no quadril, as
coordenadas das origens e inserções musculares são expressas
através de uma seqüência de transformações (3 rots., 3 desloc. p/
cada articulação)
24
Funções cinemáticas
Exemplo: deslocamento dx da tíbia em função
do ângulo de flexão do joelho
-3
x 10
6
4
dx, ref. tibial
2
0
-2
-4
-6
-1 2 0
-1 0 0
-8 0
-6 0
-4 0
-2 0
0
â ng ulo d o jo e lho
25
–Parâmetros antropométricos e coordenadas de
origens e inserções determinados a partir de
vários cadáveres
• Calculo dos torques musculares
  
  r F
26
– Relação entre as forças musculares e os
torques articulares
 r11 r12

   r21 r22
r
 31 r32
 F1 
 r1n  
 F2 
 r2 n  


 r3n  
 Fn 
– Comprimento do atuador músculo-tendíneo
calculado como distância entre origem e
inserção considerando os pontos de contorno
27
Ajuste de curvas de regressão múltipla para os braços
de momento e comprimento dos atuadores
- músculos agrupados segundo sua dependência das
mesmas coordenadas generalizadas
Grupos
1
2
3
4
5
6
7
8
9
10
11
12
13
14
n. de coord.
generalizadas
3
1
3
1
2
1
2
1
2
2
2
1
2
3
Músculos
gmed1, gmed3, gmin1, gmin3, addlong, tfl
gmed2, gmin2,
semimem, semiten, bifemlh, sar, gra
bifemsh, vasint, vaslat, ligpat
adbrev, amag1, pect, gmax1,
gmax2, gmax3
iliacus, psoas,
quadfem, gem,
peri
rf
medgas, latgas
sol
tibpost, tibiant, perbrev, perlong, pertert
flexdig, flexhal, extdig, exthal
28
– Geração automática em Matlab dos arquivos
de entrada para o SIMM: 20 pontos para cada
amplitude de movimento por coordenada
generalizada (para 3 coordenadas, 8000 pontos)
– Ajuste de curvas por mínimos quadrados
29
Equações de regressão propostas:
30
Seleção das equações com mínimo erro
Musculo
gmed1
gmed2
gmed3
gmin1
gmin2
gmin3
semimem
semiten
bifemlh
bifemsh
sar
addlong
addbrev
amag1
amag2
amag3
tfl
pect
LMT
2
1
2
2
1
2
2
2
1
1
1
2
4
4
4
4
2
4
r1
2
1
2
2
1
2
1
1
1
1
1
2
4
4
4
4
2
4
r2
2
x
1
2
x
1
2
2
2
x
2
1
1
1
1
4
2
4
r3
2
x
2
1
x
1
1
1
1
x
1
1
x
x
x
x
2
x
31
Comprimento do rectus femoris. Real (branco); ajustada
(cinza)
32
Braço de momento do rectus femoris (em relação ao
quadri)
33
c) Modelo da mecânica muscular
Kpe: elemento elástico paralelo
B: amortecimento
C: elemento contrátil
kT: rigidez do tendão
LT: comprimento do tendão
LST: comprimento do tendão
relaxado
: ângulo de empenamento

~ T ~ T ~ MT ~ M
(1a. eq.) F  k v  v cos 

34
Cálculo da rigidez no tendão
~T
k
35
Rigidez normalizada do tendão:
T

E   1   1,2GPa   1  37,5
~T
k   M   ~T   
  ~T   ~ T
  0   L s   32MPa   L s 
Ls
36
Cálculo de
~v
M
• relações f - l e v - l
37


~ PE ~ DE ~ CE ~ T
cos  F  F  F  F
Hipérbole de Hill:
lembrando que
~ M*
2
a  aF
M
~
ve  ~ M *
4F  a
*
~
M
 F M* 1   ~
ve
1 5
 
 

4  a
4  16
 a
M
M
~
~
ve   v
~
M
a fl  aF
M
~
ve  ~ M
4 F  afl
2
38
M),
Isolando as forças nos elemento contrátil (~
F
de rigidez em paralelo e de amortecimento
~v
M
e substituindo em (*), calcula-se
através da solução de uma eq. algébrica de 2o.
grau

MT
~
v é imposto pelo movimento
39
2a. equação: dinâmica da ativação
da dt  (u  a )( k 1u  k 2 )
onde Tact=1/(k1+k2) e Tdeac=1/k2.
40
Equações dinâmicas
 x4



 x5

x

6

 x 1  



 
 x 7n 
2






r
r

r
x

 11 12
 x2  
1n 
4





x




7

n

1

1
2
 x   M ( x , x , x ) [D] r r  r 
 [C( x1 , x 2 , x 3 )] x 5   g( x1 , x 2 , x 3 )  
1
2
3
21
22
2n 


 3  


 x2 






r
r

r

 x4  
3 n 

 6
  31 32

x
 72n 
 x  


 5    a 1  f1 (u1 , a1 )

 x 6  


a

f
(
u
,
a
)
2
2
2
2

 


x
 7  

 x 8  


a

f
(
u
,
a
)

  n n n n


~
~
~
~

  T
MT
T
T
T
M

F

g
(
a
,
L
,
F
,
k
,
L
,
F
,...)
1
1
1
s
0
 x


~MT ~ T ~ T ~T M
 7  2 n   F T  g (a , L
, F , k , L s , F0 ,...)
 2

2
2




~MT ~ T ~ T ~T M
 F T  g (a , L

, F , k , L s , F0 ,...)
n
n
41
 n

2) Problema de controle ótimo
-
Malha aberta / malha fechada
-
Controle ótimo
42
Controle da postura em malha aberta
utilizando controle ótimo
Objetivo: levar o modelo proposto de uma
condição inicial do agachamento até a postura
ereta, minimizando uma função de custo
Os controles obtidos correspondem às
excitações de cada músculo ao longo do
tempo de simulação
Principal vantagem do controle ótimo: solução
do problema da redundância de atuadores
43
b
 



min
f ( u, )  g o ( , x ( b))   l o ( t, x , u ) dt 
OCP ( u ,  )Lm , 2 [ a , b ]x n 
a


x  h( t, x, u )
equações dinâmicas
x (a )  
condições iniciais
t  [a , b ]
j
j
u min
( t )  u j ( t )  u max
( t ),
j  1,..., m
vínc. de controle
g ei (, x (b))  0,   q ei
vínculos de desigualdade
g ee (, x (b))  0,   q ee
vínculos de igualdade
x(t) é n x 1, u(t) é m x 1, l e g são escalares
44
Solução do Problema de controle ótimo
Utilização de Algoritmos de controle ótimo
baseados na Teoria das aproximações
consistentes (RIOTS)
Diversos problemas numéricos precisaram
ser resolvidos
45
6. Alguns resultados
Modelo com 10 atuadores musculares nãolineares
Hipóteses simplificadoras:
1. Contração isométrica
2. Relação força x comprimento constante
3. Seleção de 10 grupos musculares,
agrupando os músculos de função e
morfologia semelhantes, eliminando
músculos com r muito pequeno
46
47
Músculo
grupos
LM0

LST
gmed1
gmed2
gmed3
gm 1
semimem
semiten
bifemlh
gm2
bifemsh
gm3
gmax1
gmax2
gmax3
gm4
iliacus
psoas
gm5
rf
gm6
vasmed
vasint
vaslat
gm7
medgas
latgas
gm8
sol
tibpost
gm9
tibant
perlong
extdig
exthal
gm10
1
1
1
0.0535
0.0845
0.0643
0.0681
0.0800
0.2010
0.1090
0.1108
0.1730
0.1730
0.1420
0.1470
0.1440
0.1448
0.1000
0.1040
0.1018
0.0840
0.0840
0.0890
0.0870
0.0840
0.0857
0.0450
0.0640
0.0507
0.0300
0.0310
0.0301
0.0980
0.0490
0.1020
0.1110
0.1006
8.0
0.0
19.0
9.9344
15.0
5.0
0.0
7.5027
23.0
23.0
5.0
0.0
5.0
3.0026
7.0
8.0
7.4496
5.0
5.0
5.0
3.0
5.0
4.3788
17.0
8.0
14.3097
25.0
12.0
23.6436
5.0
10.0
8.0
6.0
6.0362
0.0780
0.0530
0.0530
0.0595
0.3590
0.2620
0.3410
0.3363
0.1000
0.1000
0.1250
0.1270
0.1450
0.1333
0.0900
0.1300
0.1080
0.3460
0.3460
0.1260
0.1360
0.1570
0.1407
0.4080
0.3850
0.4011
0.2860
0.3100
0.2885
0.2230
0.3450
0.3450
0.3050
0.2649
2
2
2
3
4
4
4
5
5
6
7
7
7
8
8
9
9
10
10
10
10
F0M
546
382
435
1363
1030
328
717
2075
402
402
382
546
368
1296
429
371
800
779
779
1294
1365
1871
4530
1113
488
1601
2839
1270
4109
603
754
341
108
1052
Fr0M1
-6.825
-8.595
-10.875
-26.295
-56.65
-21.32
-46.605
-124.58
-15.28
-27.30
-25.76
-68.34
15.44
12.61
28.05
36.61
36.61
F0M
r2
F0Mr3
-41.20
-13.78
-39.435
-94.415
-20.10
-20.10
22.98
22.98
40.89
40.95
57.06
139.90
-16.69
-6.10
-22.79
-43.96
-19.76
-63.72
-109.02
-12.70
-121.72
25.69
6.03
13.64
4.54
43.87
LMT
0.12
0.13
0.11
0.1191
0.42
0.47
0.45
0.4398
0.21
0.21
0.20
0.21
0.24
0.2191
0.20
0.26
0.2270
0.45
0.45
0.23
0.25
0.26
0.2464
0.41
0.41
0.41
0.30
0.35
0.3052
0.30
0.40
0.44
0.40
0.3539
48
trajetória resp13-12-3
3
2
x, xp (rad)
1
0
-1
perna
coxa
tronco
vel. perna
vel. coxa
vel. tronco
-2
-3
0
0.05
0.1
0.15
0.2
t(s)
0.25
0.3
0.35
0.4
49
controle ótimo resp13-12-3
1
0.8
0.6
0.4
u(t)
0.2
0
-0.2
-0.4
-0.6
-0.8
0
0.05
0.1
0.15
0.2
t(s)
0.25
gmed
sm,st,bflh
bfsh
gmax
ilipsoas
rf
vas
gas
sol,tip
ta,pl,exd,exh
0.3
0.35
0.4
50
ativações resp13-12-3
0.7
gmed
sm,st,bflh
bfsh
gmax
ilipsoas
rf
vas
gas
sol,tip
ta,pl,exd,exh
0.6
0.5
a(t)
0.4
0.3
0.2
0.1
0
0
0.05
0.1
0.15
0.2
t(s)
0.25
0.3
0.35
0.4
51
forças resp13-12-3
2500
gmed
sm,st,bflh
bfsh
gmax
ilipsoas
rf
vas
gas
sol,tip
ta,pl,exd,exh
2000
F(N)
1500
1000
500
0
0
0.05
0.1
0.15
0.2
t(s)
0.25
0.3
0.35
0.4
52
53
Resultados com tf = 1.0s e
gravidade variável
- Gravidade: 2 m/s2  4 m/s2  6 m/s2  8
m/s2  9 m/s2  9.81 m/s2
- Tempo total de 1 s de simulação: +- 30 dias
de CPU e 490 MB RAM Pentium 600 MHz
54
trajetória resp8-1-1
2
perna
coxa
tronco
vel. perna
vel. coxa
vel. tronco
1.5
x, xp (rad)
1
0.5
0
-0.5
-1
-1.5
0
0.2
0.4
0.6
0.8
1
t(s)
55
controle ótimo resp8-1-1
0.5
gmed
sm,st,bflh
bfsh
gmax
ilipsoas
rf
vas
gas
sol,tip
ta,pl,exd,exh
0.45
0.4
0.35
u(t)
0.3
0.25
0.2
0.15
0.1
0.05
0
0
0.2
0.4
0.6
0.8
1
t(s)
56
ativações resp8-1-1
0.5
gmed
sm,st,bflh
bfsh
gmax
ilipsoas
rf
vas
gas
sol,tip
ta,pl,exd,exh
0.45
0.4
0.35
a(t)
0.3
0.25
0.2
0.15
0.1
0.05
0
0
0.2
0.4
0.6
0.8
1
t(s)
57
forças resp8-1-1
2500
gmed
sm,st,bflh
bfsh
gmax
ilipsoas
rf
vas
gas
sol,tip
ta,pl,exd,exh
2000
F(N)
1500
1000
500
0
0
0.2
0.4
0.6
0.8
1
t(s)
58
59
-Padrões de excitação semelhantes ao
caso anterior
-Níveis de u(t), a(t) e F(t) mais baixos
-Velocidades máximas inferiores
-Queda no início do movimento mais
pronunciada
-Oscilação maior do tronco
60
Trabalhos futuros
1. Novas funções de custo como,
por exemplo, a maximização da
altura do centro de massa.
2. Realizar estudos do movimento
de levantar de uma cadeira e da
marcha
61
3. Introduzir no modelo biomecânico
expressões de momento passivo
gerado por ligamentos do joelho
4. Formulação de protocolos para
levantamento de parâmetros
antropométricos individuais
62
3) Projetos em andamento
1. Comprovação experimental dos
resultados através de laboratório de
marcha e análise de padrões EMG
(Temático FAPESP / IOT)
2. Nova versão do RIOTS: CAOS
(Consistent Approximations Optimal
control Solver)
63