Modelação da Metástase Óssea utilizando Derivadas ... · fracionárias em equações diferenciais, a análise do corpotamento dinâmico do remodelamento do osso na presença ou

$Page 1: Modelação da Metástase Óssea utilizando Derivadas ... · fracionárias em equações diferenciais, a análise do corpotamento dinâmico do remodelamento do osso na presença ou$
Modelação da Metástase Óssea utilizando Derivadas

Fracionárias

Luiz Filipe Christ

Dissertação para obtenção do Grau de Mestre em

Engenharia Mecânica

Orientadores: Prof. Duarte Pedro Mata de Oliveira Valério

Prof. Susana de Almeida Mendes Vinga Martins

Júri

Presidente: Prof. João Rogério Caldas Pinto

Orientador: Prof. Duarte Pedro Mata de Oliveira Valério

Vogal: Prof. António Joaquim dos Santos Romão Serralheiro

Julho de 2015

ii

Agradecimentos

Aos meus orientadores, por toda a paciência, dedicação, atenção e ajuda ao longo de todos esses

meses. Aos meus pais e amigos, pelo suporte nos momentos difíceis, e a todos que conspiraram a

favor dessa dissertação.

Por fim, um agradecimento especial ao hospital Santa Maria e ao Rui Coelho, por fornecerem

material e apoio fundamentais para a elaboração dessa dissertação, e ao projecto CancerSys, o qual

esse trabalho está inserido.

iii

Abstract

The adult skeleton is a highly specialized organ that undergoes constant remodeling over time. This

process is performed throught two highly specialized cells: osteoclasts, which are responsible for bone

resorption, and osteoblasts, for bone rebuilding. This dynamic is responsible for leaving the body in

homeostasis.

The bone metastasis disease is characterized by affecting the dynamic of these cells. As the tumor

develops, there is an increase in the population of osteoclasts and osteoblasts, a reduction in, so that

there is a decrease in pacient’s bone density.

For this reason, the study of this dynamic model is essential to the development of better therapies for

the disease. The work presented here performs, through implementation of fractional derivatives in

differential equations, the analysis of the dynamic bone remodeling behavior in the presence or

absence of tumor and treatment when considered a discretized single point of the bone. Later, these

cases are extended to one spatial dimension.

The resulting dynamic behavior has the characteristic expected from the literature, where higher

fractional orders leads to a more fast and oscillatory system whereas smaller orders have a more

damped behavior.

Keywords: Bone Remodeling, Bone Metástasis, Fractional Derivatives, Dynamic System.

Resumo

O esqueleto adulto é um órgão altamente especializado que sofre constante remodelamento ao longo

do tempo. Tal processo é realizado através de duas células altamente especializadas: os

osteoclastos, responsáveis por destruir o osso velho, e os osteoblastos, que reconstroem o osso.

Essa dinâmica é responsável por deixar o corpo em homeostase.

A metástase óssea caracteriza-se por afetar a dinâmica dessas células. A medida que o tumor se

desenvolve ocorre um aumento na população de osteoclastos e uma redução na de osteoblastos, de

modo a haver uma diminuição na densidade óssea do paciente.

Por essa razão, o estudo desse modelo dinâmico é fundamental para o desenvolvimento de melhores

terapias para a doença. O trabalho apresentado aqui realiza, através da implementação de derivadas

fracionárias em equações diferenciais, a análise do corpotamento dinâmico do remodelamento do

osso na presença ou ausência de tumor e tratamento para o caso onde considera-se um ponto

discreto do osso. Em seguida esses casos são estendidos para uma dimensão espacial.

O comportamento dinâmico resultante apresenta as características esperadas pela literatura, onde

ordens fracionárias mais elevadas levam o sistema a um comportamento mais oscilatório e rápido

enquanto que menores apresentam um comportamento mais amortecido.

Palavras-chave: Remodelamento do Osso, Metástase óssea, Derivadas Fracionárias, Sistema

Dinâmico.

iv

Índice

Agradecimentos ........................................................................................................................................ ii

Abstract.................................................................................................................................................... iii

Resumo ................................................................................................................................................... iii

Índice ....................................................................................................................................................... iv

Índice de figuras ..................................................................................................................................... vii

Índice de tabelas ..................................................................................................................................... ix

Abreviaturas ............................................................................................................................................. x

1. Introdução ........................................................................................................................................1

1.1. Considerações Iniciais ............................................................................................................1

1.2. Proposta e objetivos ................................................................................................................2

1.3. Estrutura da dissertação .........................................................................................................3

2. Cálculo fracionário ...........................................................................................................................4

2.1. Função Gama..........................................................................................................................4

2.2. Cálculo de ordens inteiras .......................................................................................................5

2.3. Cálculo de ordens fracionárias................................................................................................6

2.3.1. Definições ........................................................................................................................ 6

2.3.2. Operador não-local .......................................................................................................... 7

2.3.3. Operadores lineares ........................................................................................................ 7

2.3.4. Relações entre as diferenes definições .......................................................................... 7

2.3.5. Lei dos expoentes............................................................................................................ 7

2.3.6. Princípio da Memória Curta (“Short Memory”) ................................................................ 8

2.4. Transformada de Laplace no cálculo fracionário ....................................................................8

2.4.1. Riemann-Liouville e Grünwald-Leitnikoff ......................................................................... 8

2.4.2. Caputo ............................................................................................................................. 8

3. Funções de transferência ................................................................................................................9

3.1. Definições ................................................................................................................................9

3.1.1. Função de transferência SISO ........................................................................................ 9

3.1.2. Função de transferência MIMO ....................................................................................... 9

3.2. Estabilidade .............................................................................................................................9

3.3. Resposta no domínio do tempo e da frequência ................................................................. 10

3.3.1. Resposta em frequência de uma função de transferência genérica ............................. 10

3.3.2. Estudo de 𝟏/(𝒔 𝒂⁄ )𝟐𝜶 + 𝟐𝜻(𝒔 𝒂⁄ )𝜶 + 𝟏 ........................................................................... 10

4. Método aplicado ........................................................................................................................... 13

4.1. Matrizes triangulares ............................................................................................................ 13

4.2. Produto Kronecker de matrizes............................................................................................ 14

4.3. Eliminadoras......................................................................................................................... 15

4.4. Derivadas fracionárias com matrizes triangulares superiores (inferiores) ........................... 15

4.4.1. Derivada fracionária esquerda ...................................................................................... 15

4.4.2. Derivada fracionária direita ............................................................................................ 16

v

4.5. Solução numérica de equações diferenciais ....................................................................... 17

4.6. Solução numérica de equações diferenciais parciais .......................................................... 18

4.6.1. Introdução ...................................................................................................................... 18

4.6.2. Metodologia ................................................................................................................... 18

4.7. Método do passo largo (“Large Step”) ................................................................................. 20

5. Dinâmica do Osso ........................................................................................................................ 22

5.1. Regeneração do osso .......................................................................................................... 22

5.2. Formação e morte dos osteoclastos e osteoblastos............................................................ 22

5.3. Remodelamento do osso ..................................................................................................... 23

5.4. Tumor e seus efeitos ............................................................................................................ 24

5.5. Tratamento ........................................................................................................................... 24

5.6. Terapias direcionadas .......................................................................................................... 25

6. A farmacocinética e a farmacodinâmica ....................................................................................... 26

6.1. Farmacocinética ................................................................................................................... 26

6.1.1. Introdução ...................................................................................................................... 26

6.1.2. Modelos ......................................................................................................................... 26

6.1.3. Aplicações ..................................................................................................................... 28

6.2. Farmacodinâmica ................................................................................................................. 29

6.2.1. Introdução ...................................................................................................................... 29

6.2.2. Modelamento ................................................................................................................. 30

6.2.3. Efeito de múltiplos fármacos ......................................................................................... 31

7. Estudo dos modelos dinâmicos .................................................................................................... 32

7.1. Caso sem tumor e sem tratamento 0D ................................................................................ 32

7.2. Caso com tumor e sem tratamento 0D ................................................................................ 35

7.3. Caso com tumor e com tratamento de inibidores de proteassoma 0D................................ 37

7.4. Caso com tumor e com tratamento utilizando PK/PD 0D .................................................... 39





8. Aplicação de derivadas fracionárias ............................................................................................. 47









9. Conclusão ..................................................................................................................................... 74

vi

9.1. Trabalho futuro ..................................................................................................................... 75

Referências ........................................................................................................................................... 77

Anexo ..................................................................................................................................................... 81

vii

Índice de figuras

Figura 1. Interação das células na remodelação do osso (extraído de [3]). ........................................... 1

Figura 2. (a) Mapa de estabilidade para uma função de transferência de α=q menor que 1 e (b) para

α=q maior que 1. (extraído de [31]). ...................................................................................................... 10

Figura 3. Comportamento da frequência (extraído de [43]). ................................................................. 11

Figura 4. Resposta do sistema 𝐺(𝑠) = 𝑠2𝛼 + 3𝑠𝛼 + 1 a um degrau unitário para diversas ordens de

derivada. ................................................................................................................................................ 12

Figura 5. Discretização para derivadas de ordem inteira. (extraído de [36]). ....................................... 18

Figura 6. Discretização do tempo para derivadas de ordem fracionárias (extraído de [36]). ............... 19

Figura 7. Discretização no tempo e no espaço para derivadas de ordem fracionárias bilaterais

(extraído de [36]). .................................................................................................................................. 19

Figura 8. Modelo de um compartimento (extraído de [11]). .................................................................. 27

Figura 9. Modelo de dois compartimentos (extraído de [11]). ............................................................... 28

Figura 10. Alteração da concentração de um farmaco no plasma no modelo de dois

compapartimentos (extraído de [11]). ................................................................................................... 28

Figura 11. Modelo farmacocinético para a administração de uma dose oral (extraído de [9]). ............ 29

Figura 12. Modelo farmacocinético para múltiplas doses (extraído de [9]). ......................................... 29

Figura 13. Modelo farmacodinâmico de um farmaco (extraído de [9]). ................................................ 30

Figura 14. Dinâmica das populações provocadas por um distúrbio na população de osteoclastos. ... 34

Figura 15. Dinâmica das populações com remodelação constante. .................................................... 34

Figura 16. Dinâmica das populações na presença de um tumor. ......................................................... 37

Figura 17. Dinâmica das populações na presença de um tumor e tratamento. ................................... 38

Figura 18. Modelo de tratamento com o uso de PK/PD. ....................................................................... 41

Figura 19. Distribuição inicial da população de osteoclastos (extraído de [3]). .................................... 42

Figura 20. Comportamento dinâmico para o caso com uma dimensão espacial. ................................ 42

Figura 21. Comportamento dinâmico na presença de tumor para uma dimensão espacial. ............... 44

Figura 22. Comportamento dinâmico na presença de tumor e tratamento com uma dimensão

espacial.................................................................................................................................................. 45

Figura 23. Modelo de tratamento com o uso de PK/PD com uma dimensão espacial. ........................ 46

Figura 24. Dinâmica das populações para α=1.3. ................................................................................ 49






Figura 30. Comparação do comportamento de diversas ordens fracionárias para um mesmo

parâmetro de equações ......................................................................................................................... 51

Figura 31. Sistema não apresenta comportamento cíclico quando α=0.8. ........................................... 52

Figura 32. Comportamento das populações para α=1.3 na presença do tumor. ................................. 53


viii





Figura 38. Gráficos comparativos entre diversas ordens fracionárias. Os três primeiros correspondem

a α≤0.95 e os três últimos para 1≤α≤1.3. .............................................................................................. 56

Figura 39. Comportamento das populações para α=1.2 na presença do tumor e tratamento. ............ 58




Figura 43. Modelo PK/PD de um farmaco para diversas ordens fracionárias. ..................................... 60

Figura 44. Comportamento das populações para α=1.1 para o modelo PK/PD. .................................. 61



Figura 47. Dinâmica da população 1D para α=1.2. .............................................................................. 64




Figura 51. Dinâmica do modelo com tumor 1D para α=1.2. ................................................................. 67

Figura 52. Dinâmica do modelo com tumor 1D para α=1.1.. ................................................................ 67



Figura 55. Dinâmica do modelo com tumor e tratamento 1D para α=0.9. ............................................ 70



Figura 58. Dinâmica do modelo com PK/PD 1D para α=0.9. ................................................................ 72

Figura 59. Dinâmica do modelo com PK/PD 1D para α=0.7. ................................................................ 73

Figura 60. Modelo PK/PD em Simulink. ................................................................................................ 81

Figura 61. Blocos Drug Resistance2 e PD, respectivamente. .............................................................. 81

ix

Índice de tabelas

Tabela I. Parâmetros dos farmacos para o modelo PK/PD. ................................................................. 40

Tabela II. Parâmetros dos farmacos para o modelo PK/PD 1D............................................................ 46

Tabela III. Valores dos parâmetros utilizados para cada simulação. .................................................... 48

Tabela IV. Parâmetros das populações utilizados nas simulações de ordem fracionária. .................. 53

Tabela V. Parâmetros do tumor utilizados nas simulações de ordem fracionária. ............................... 53

Tabela VI. Parâmetros do tumor e do tratamento com inibidores de proteassoma utilizados na

simulação............................................................................................................................................... 57

Tabela VII. Parâmetros do tumor e do tratamento utilizados nas simulações com PK/PD para

derivadas fracionárias. .......................................................................................................................... 60

Tabela VIII. Parâmetros dos farmacos utilizadaos na simulação de derivadas fracionárias. ............... 61

Tabela IX. Parâmetros das populações utilizados nas simulações com 1D. ........................................ 63

Tabela X. Parâmetros das populações utilizados nas simulações com tumor em 1D. ........................ 66

Tabela XI. Parâmetros do tumor utilizados nas simulações em 1D. .................................................... 66

Tabela XII. Parâmetros do tumor e do tratamento com inibidores de proteassoma utilizados na

simulação em 1D com derivadas fracionárias....................................................................................... 69

Tabela XIII. Parâmetros do tumor e do tratamento utilizados nas simulações com PK/PD em 1D com

derivadas fracionárias. .......................................................................................................................... 72

Tabela XIV. Parâmetros das drogas utilizadaos na simulação em 1D com derivadas fracionárias. .... 72

x

Abreviaturas

1D Uma dimensão

BMP Bone Morphogenetic Protein

BMU Basic Multicelular Unit

c-fms colony factor macrophage stimulating

CSF-1 Colony Stimulating Factor 1

DKK-1 Dickkopf-related protein 1

EDP Equação Diferencial Parcial

IGF-I Insulin growth factor 1

IGF-II Insulin growth factor 2

IL-8 Interleukin 8

MCP-1 Monocyte Chemoattractant Protein 1

OPG Osteoprotegerin

PD Pharmacodynamic

PK Pharmacokinectic

PTH Parathyroid Hormone

RANK Receptor Activator of Nuclear Factor K

RANKL Receptor Activator of Nuclear Factor K Ligand

TGF-𝛽 Transforming Growth Factor 𝛽

1

1. Introdução

1.1. Considerações Iniciais

O esqueleto adulto é um órgão dinâmico e altamente especializado que sofre constante

remodelação ao longo do tempo. Apesar do propósito de tal processo ainda não ser totalmente

entendido, ele permite reparar danos existentes e evitar o envelhecimento [22].

Dois tipos de células altamente especializadas, geradas a partir da medula óssea, são

responsáveis pela remodelação do osso: os osteoclastos, que atuam na degradação do osso antigo,

e os osteoblastos, que formam o novo osso. Essas células interagem entre sí, de modo a afetar a sua

sintese e apoptose e contribuir para a homeostase do sistema. Algum problema metabólico do osso,

como a osteoporose, normalmente é causado por uma desordem no nascimento ou morte dessas

células [38].

O cancro é uma doença muito complexa que é caracterizada por uma população de células

que cresce e se divide de uma forma anormal. Tais células podem invadir e destruir tecidos

adjacentes e se espalhar para lugares distantes do corpo, processo conhecido como metástase [5].

O osso é um local comum de ocorrer metástase uma vez que apresenta características que

favorecem o crescimento e proliferação de células cancerígenas. Essas células acabam por interferir

no mecanismo de remodelamento ósseo, de modo a desequilibrar os processos de destruição e

reconstrução do osso. Quando a metástase acentua o fenômeno de destruição do osso ela é

classificada como osteolítica, como é o caso do cancro de mama e do mieloma múltiplo [6].

Figura 1. Interação das células na remodelação do osso (extraído de [3]).

O tratamento para a metástase óssea envolve dois aspetos: destruir o cancro através da

quimioterapia e hormônios e impedir a redução da densidade óssea através da regulação da

população de osteoclastos e osteoblastos. Essa última pode ser feita através do uso de certos tipos

de fármaco como bisfosfanatos, denosumab, inibidores de proteassoma, entre outros [20].

2

A terapia deve proporcionar o máximo de eficiência com o mínimo de toxidade para o

paciente, de modo a garantir o seu bem-estar.

Com isso em vista, o modelamento matemático desse complexo sistema pode vir a ser uma

ferramenta poderosa para um melhor compreendimento sobre o desenvolvimento da doença e

auxiliar na elaboração de novos mecanismos de tratamento que possam ser menos agressivos ao

organismo.

1.2. Proposta e objetivos

O estudo do cálculo fracionário tem ganho espaço nos últimos anos. Muitos fenômenos

físicos (mecânica, processamento de sinais e difusão) e biológicos (viscoelasticidade de membrana e

propagação de transmissões) apresentaram uma boa caracterização com o uso de ordens

fracionárias. A expansão de operações matemáticas através do calculo fracionário permite o

desenvolvimento de novas relações e modelamentos para sistemas físicos, químicos ou biológicos

complexos.

Assim, esta dissertação apresentará o estudo das derivadas fracionárias para a dinâmica de

remodelamento do osso descrita por equações diferenciais. Tal proposta é motivada por uma

curiosidade matemática e física em saber como que a dinâmica se comportará quando valores

fracionários para as derivadas forem utilizados no lugar dos tradicionais inteiros.

Será feita a aplicação de ordem fracionária em relação ao tempo para uma série de modelos

dinâmicos sem dimensões espaciais. O primeiro caso que será abordado é quando o sistema não

apresenta nenhum tumor, ou seja, uma pessoa saudável. Pretende-se estudar como que a população

de osteoclastos e osteoblastos interagem entre si e seu reflexo na densidade óssea e como que a

mudança do termo diferencial pode afetar tal comportamento.

O caso seguinte será incluir a doença mieloma ao sistema. Será possível analisar como que a

presença do tumor afeta a dinâmica de remodelamento do osso e suas respectivas consequências

para o organismo.

Dois modelos de tratamento serão estudados também. Um deles utiliza fármacos que visam

estimular a redução da apoptose dos osteoblastos na presença do tumor. Esses fármacos são

conhecidas como inibidores de proteassoma. O outro modelo aplicará fármacos que visam o aumento

da apoptose e redução da formação de osteoclastos, como os bisfosfanatos e denosumab. Além

disso, esse último modelo implementará também conceitos de farmacocinética e farmacodinâmica,

cada vez mais utilizados no desenvolvimento de terapias eficientes.

Por fim, todos esses modelos serão estendidos para o caso onde exista uma dimensão

espacial e ocorra difusão das células ao longo do osso. Esse processo busca uma caracterização

mais fiel do osso.

Uma conclusão e discussão dos resultados obtidos é feita no final dessa dissertação de modo

a destacar os pontos e resultados mais importantes.

3

1.3. Contributo da tese

A dissertação apresentada nesse documento contribui com a implementação e estudo de

derivadas fracionárias aos modelos dinâmicos de remodelação do osso na presença/ausência de

tumor/tratamento, bem como nos conceitos de farmacocinética e farmacodinâmica para os casos

discretos e com uma dimensão espacial.

1.4. Estrutura da dissertação

Esta dissertação foi elaborada em nove capítulos, de modo a que o leitor possa entender os

conceitos matemáticos básicos de derivadas fracionárias e biológicos, referentes à dinâmica de

remodelação óssea, para depois professeguir ao estudo dos modelos propostos. Tais capítulos

apresentam a seguinte divisão:

Capítulo 1: introduz a dissertação de uma maneira geral, abordando seu escopo, objetivos e

organização do relatório.

Capítulo 2: apresenta os conceitos e definições básicas do calculo fracionário. Desse modo o

leitor estará apto a entender o método aplicado para a resolução dos sistemas de equações

diferenciais de ordem fracionária.

Capítulo 3: explica os conceitos de estabilidade e comportamento dinâmico de sistemas para

diversas ordens fracionárias.

Capítulo 4: nessa parte é apresentado o método utilizado para a solução de equações

diferenciais fracionárias. O conceito é estendido também para equações diferenciais parciais

fracionárias.

Capítulo 5: expõe os conceitos biológicos que foram a base para a elaboração das equações

dinâmicas.

Capítulo 6: destina-se a explicar os conceitos de farmacocinética e farmacodinâmica e como

que eles são utilizados no planejamento de terapias.

Capítulo 7: apresenta os modelos dinâmicos e os resultados obtidos com o uso de derivadas

de ordem inteira. Além disso é feito uma análise entre esses resultados e seu significado, em

termos biológicos, que servirão de base para o estudo dos casos de ordem fracionária.

Capítulo 8: destina-se a implementação das ordens fracionárias aos modelos apresentados

no capítulo anterior, bem como a análise do resultado de suas simulações.

Capítulo 9: levantamento das conclusões obtidas nesta dissertação e apresentação de

possíveis próximos estudos.

Por fim, vale mencionar que todos os cálculos e simulações foram feitos com o uso do software

Matlab e Simulink.

4

2. Cálculo fracionário

Neste capítulo serão apresentados alguns conceitos e definições necessários para o estudo do

cálculo fracionário. Primeiramente é apresentada a definição da função gama e do operador D, bem

como algumas de suas propriedades, que são ferramentas fundamentais para a compreensão das

formulações de cálculo fracionário.

Em seguida são apresentadas as definições de Riemann-Liouville, Caputo e Grünwald-Letnikoff

de derivadas e integrais fracionárias, além do conceito de operador não-local relacionado. Tais

definições serão utilizadas para a elaboração do método de resolução de equações diferenciais

fracionárias.

Por fim, é feito uma apresentação das definições de Riemann-Liouville, Grünwald-Letnikoff e

Caputo para transformadas de Laplace de ordens fracionárias. Através dela será possível estudar a

dinâmica de sistemas de ordem fracionária através de funções de transferência, assunto discutido no

capítulo seguinte.

Esse capítulo tem como objetivo apenas expor as ferramentas necessárias para a compreensão

dos capítulos seguintes, bem como ser utilizado como instrumento de consulta. As demonstrações

das propriedades, bem como uma análise mais completa sobre os assuntos abordados podem ser

encontrados em [43], no qual esse capítulo foi baseado, e [33].

2.1. Função Gama

A função gama Γ está presente em diversos tópicos sobre o cálculo fracionário. A sua

definição, bem como algumas de suas propriedades importantes, é apresentada a seguir:

Definição 2.1

Define-se a função gama como:

{

Γ(𝑥) = ∫ 𝑒−𝑦

+∞

0

𝑦𝑥−1𝑑𝑦, 𝑥 ∈ ℝ+

Γ(𝑥) =Γ(𝑥 − ⌊𝑥⌋)

∏ (𝑥 + 𝑘)−⌊𝑥⌋−1𝑘=0

, 𝑥 ∈ ℝ− (2.1)

Propriedade 2.1

Γ(1) = ∫ 𝑒−𝑦𝑑𝑦

+∞

0

= [−𝑒−𝑦]0+∞ = 1 (2.2)

Propriedade 2.2

Γ(𝑥 + 1) = 𝑥Γ(𝑥) (2.3)

ou seja, a função gama é uma generalização do operador fatorial.

Propriedade 2.3

lim𝑥→0+

Γ(𝑥) = +∞ (2.4)

Definição 2.2

A combinação de um elemento a com um elemento b é definida como:

(𝑎𝑏) =

𝑎!

𝑏! (𝑎 − 𝑏)!, 𝑎, 𝑏 ∈ ℤ0

+ (2.5)

5

Através da função gama e suas propriedades, pode-se reescrever a operação de combinação

do seguinte modo:

(𝑎𝑏) =

{

Γ(𝑎 + 1)

Γ(𝑏 + 1)Γ(𝑎 − 𝑏 + 1), 𝑠𝑒 𝑎, 𝑏, 𝑎 − 𝑏 ∉ ℤ−

(−1)𝑏Γ(𝑏 − 𝑎)

Γ(𝑏 − 𝑎)Γ(−𝑎), 𝑠𝑒 𝑎 ∈ ℤ− ∧ 𝑏 ∈ ℤ0

+

0, 𝑠𝑒 [(𝑏 ∈ ℤ− ∨ 𝑏 − 𝑎 ∈ ℕ) ∧ 𝑎 ∉ ℤ−] ∨ (𝑎, 𝑏 ∈ ℤ− ∧ |𝑎| > |𝑏|)

(2.6)

2.2. Cálculo de ordens inteiras

Para o futuro estudo do cálculo fracionário é necessário utilizar algumas definições e

resultados do cálculo tradicional, de ordem inteira. Além disso, o operador D, associado a ordens 𝑛 ∈

ℤ, é introduzido:

Definição 2.3

𝐷𝑐 𝑡𝑛𝑓(𝑡) =

{

𝑑𝑛𝑓(𝑡)

𝑑𝑡𝑛, 𝑠𝑒 𝑛 ∈ ℕ

𝑓(𝑡), 𝑠𝑒 𝑛 = 0

∫ 𝐷𝑡𝑛+1𝑓(𝑡)𝑑𝑡𝑐

, 𝑠𝑒 𝑛 ∈ ℤ−𝑡

𝑐

(2.7)

𝐷𝑐𝑛𝑓(𝑡)𝑡

=

{

(−1)𝑛

𝑑𝑛𝑓(𝑡)

𝑑𝑡𝑛, 𝑠𝑒 𝑛 ∈ ℕ

𝑓(𝑡), 𝑠𝑒 𝑛 = 0

∫ 𝐷𝑐𝑛+1𝑓(𝑡)𝑑𝑡𝑡

, 𝑠𝑒 𝑛 ∈ ℤ−𝑐

𝑡

(2.8)

Propriedade 2.4

A derivada de ordem 𝑛 ∈ ℕ da função f(t) é dada por:

𝐷𝑛𝑓(𝑡) =𝑑𝑛𝑓(𝑡)

𝑑𝑡𝑛= lim

ℎ→0

∑ (−1)𝑘 (𝑛𝑘) 𝑓(𝑡 − 𝑘ℎ)𝑛

𝑘=0

ℎ𝑛 (2.9)

Propriedade 2.5

Dado uma integral de Riemann, descrita por:

∫𝑓(𝑡)𝑑𝑡

𝑡

𝑐

= limℎ→0+

∑ ℎ

⌊𝑡−𝑐ℎ⌋

𝑘=0

𝑓(𝑡 − 𝑘ℎ) (2.10)

onde ⌊𝑡−𝑐

ℎ⌋ é o maior inteiro igual ou inferior a

𝑡−𝑐

ℎ. E a fórmula de Cauchy para integrais indefinidas de

ordem 𝑛 ∈ ℕ de f(t):

𝐷𝑐 𝑡−𝑛𝑓(𝑡) = ∫

(𝑡 − 𝜏)𝑛−1

(𝑛 − 1)!

𝑡

𝑐

𝑓(𝜏)𝑑𝜏 (2.11)

𝐷𝑡 𝑐−𝑛𝑓(𝑡) = ∫

(𝜏 − 𝑡)𝑛−1

(𝑛 − 1)!

𝑐

𝑡

𝑓(𝜏)𝑑𝜏 (2.12)

pode-se obter:

𝐷𝑐 𝑡−𝑛𝑓(𝑡) = lim

ℎ→0+∑ ℎ

(𝑘ℎ)𝑛−1

(𝑛 − 1)!


𝑘=0

𝑓(𝑡 − 𝑘ℎ) (2.13)

6

Uma expressão semelhante é encontrada para o caso 𝐷𝑡 𝑐−𝑛𝑓(𝑡).

Propriedade 2.6

Se todas as derivadas existirem, a Lei dos Expoentes é dada por:

𝐷𝑐 𝑡𝑚 𝐷𝑡

𝑛𝑐 𝑓(𝑡) = 𝐷𝑡

𝑚+𝑛𝑐 𝑓(𝑡) (2.14)

para os seguintes casos:

𝑚, 𝑛 ∈ ℤ0+ (2.15)

𝑚, 𝑛 ∈ ℤ0− (2.16)

𝑚 ∈ ℤ + ∧ 𝑛 ∈ ℤ

− (2.17)

2.3. Cálculo de ordens fracionárias

2.3.1. Definições

No cálculo fracionário existem diversas definições alternativas para o conceito de derivada e

integral [33,43]. Neste tópico serão tratadas as três definições mais famosas, que são a de Riemann-

Liouville, Caputo e Grünwald-Letnikoff.

Definição 2.4 (Riemann-Liouville)

𝐷𝑡𝛼𝑓(𝑡)𝑐

=

{

∫(𝑡 − 𝜏)−𝛼−1

Γ(−𝛼)𝑓(𝜏)𝑑𝜏, 𝑠𝑒 𝛼 ∈ ℝ−

𝑡

𝑐

𝑓(𝑡), 𝑠𝑒 𝛼 = 0

𝑑⌈𝛼⌉

𝑑𝑡⌈𝛼⌉𝐷𝑡𝛼−⌈𝛼⌉𝑓(𝑡)𝑐

, 𝑠𝑒 𝛼 ∈ ℝ+

(2.18)

𝐷𝑐𝛼

𝑡 𝑓(𝑡) =

{

∫(𝜏 − 𝑡)−𝛼−1


𝑐

𝑡

𝑓(𝑡), 𝑠𝑒 𝛼 = 0

(−1)⌈𝛼⌉𝑑⌈𝛼⌉

𝑑𝑡⌈𝛼⌉𝐷𝑐𝛼−⌈𝛼⌉𝑓(𝑡)𝑡

, 𝑠𝑒 𝛼 ∈ ℝ+

(2.19)

Definição 2.5 (Caputo)

𝐷𝑡𝛼𝑓(𝑡)𝑐

=

{

∫(𝑡 − 𝜏)−𝛼−1


𝑡

𝑐

𝑓(𝑡), 𝑠𝑒 𝛼 = 0

𝐷𝑡𝛼−⌈𝛼⌉

𝑐

𝑑⌈𝛼⌉

𝑑𝑡⌈𝛼⌉𝑓(𝑡), 𝑠𝑒 𝛼 ∈ ℝ+

(2.20)

𝐷𝑐𝛼

𝑡 𝑓(𝑡) =

{

∫(𝜏 − 𝑡)−𝛼−1


𝑐

𝑡

𝑓(𝑡), 𝑠𝑒 𝛼 = 0

(−1)⌈𝛼⌉ 𝐷𝑐𝛼−⌈𝛼⌉

𝑡

𝑑⌈𝛼⌉

𝑑𝑡⌈𝛼⌉𝑓(𝑡), 𝑠𝑒 𝛼 ∈ ℝ+

(2.21)

Definição 2.6 (Grünwald-Letnikoff)

𝐷𝑡𝛼

𝑐 𝑓(𝑡) = lim

ℎ→0+

∑ (−1)𝑘 (𝛼𝑘) 𝑓(𝑡 − 𝑘ℎ)


𝑘=0

ℎ𝑎

(2.22)

7

𝐷𝑐𝛼

𝑡 𝑓(𝑡) = lim

ℎ→0+

∑ (−1)𝑘 (𝛼𝑘) 𝑓(𝑡 + 𝑘ℎ)


𝑘=0

ℎ𝑎

(2.23)

2.3.2. Operador não-local

As três definições anteriores compartilham de algo em comum: o operador D sempre

depende do limite de integração c e t, com exceção de quando 𝛼 ∈ ℝ0+. Ou seja, o operador D é um

operador não local.

Tal fato faz com que qualquer operação fracionária (excluindo as exceções citadas acima)

dependa de valores passados ou futuros da função. Nesse caso, as operações de derivada, por

exemplo, se assemelham mais às integrais do cálculo tradicional [43].

2.3.3. Operadores lineares

Pelas definições anteriores é possível notar, por inspeção, que D é um operador linear.

Propriedade 2.7

𝐷𝑡𝛼

𝑐 [𝑎𝑓(𝑡) + 𝑏𝑔(𝑡)] = 𝑎 𝐷𝑐

𝑡𝛼𝑓(𝑡) + 𝑏 𝐷𝑡

𝛼𝑐 𝑔(𝑡), 𝑜𝑛𝑑𝑒 𝑎, 𝑏 ∈ ℝ (2.24)

𝐷𝑐𝛼

𝑡 [𝑎𝑓(𝑡) + 𝑏𝑔(𝑡)] = 𝑎 𝐷𝑡

𝑐𝛼𝑓(𝑡) + 𝑏 𝐷𝑐

𝛼𝑡 𝑔(𝑡), 𝑜𝑛𝑑𝑒 𝑎, 𝑏 ∈ ℝ (2.25)

2.3.4. Relações entre as diferenes definições

Pelas definições apresentadas anteriormente pode-se notar claramente que as de Riemann-

Liouville e Caputo são diferentes para 𝛼 > 0, mas as definições de Riemann-Liouville e Grünwald-

Letnikoff apresentam o mesmo resultado (considerando que a função f satisfaça todas as condições

de cada) [43].

Propriedade 2.8

Se f(t) possui máximo de {0,⌊𝛼⌋} derivadas continuas e 𝐷 max{0,⌊𝛼⌋}𝑓(𝑡) é integrável, então

𝐷𝑡𝛼

𝑐 𝑓(𝑡) existe pelas definições de Riemann-Liouville e Grünwald-Letnikoff e ambas apresentam o

mesmo resultado.

Propriedade 2.9

Nas mesmas condições do propriedade acima, se 𝐷𝑡𝛼

𝑐 𝑓(𝑡) = 𝛾 pelas definições de Riemann-

Liouville e Grünwald-Letnikoff, então a de Caputo resulta em:

𝐷𝑡𝛼

𝑐 𝑓(𝑡) = 𝛾 − ∑

(𝑡 − 𝑐)−𝛼+𝑘𝑑𝑘𝑓(𝑐)𝑑𝑡𝑘

Γ(𝑘 − 𝛼 + 1)

⌈𝛼⌉−1

𝑘=0

(2.26)

2.3.5. Lei dos expoentes

A lei dos expoentes permanece válida para expoentes inteiros com as condições citadas

anteriormente.

Propriedade 2.10

𝐷𝑡𝛼

𝑐 𝐷𝑡

𝛽𝑐 𝑓(𝑡) = 𝐷𝑡

𝛼+𝛽𝑐 𝑓(𝑡) (2.27)

A igualdade acima é válida para as definições de Riemann-Liouville e Grünwald-Letnikoff no

caso da propriedade 2.6 e quando 𝛽 ≤ 0 ∧ 𝛼 + 𝛽 ≤ 0

8

2.3.6. Princípio da Memória Curta (“Short Memory”)

Propriedade 2.11

Quando a definição de Riemann-Liouville (ou Grünwald-Letnikoff) é utilizada, se |𝑓(𝑡)| <

𝑀, ∀𝑡 > 𝑐, então o erro 𝜀 causado pela aproximação:

𝐷𝑡𝛼

𝑐 𝑓(𝑡) ≈ 𝐷𝑡

𝛼𝑡−𝐿

𝑓(𝑡), 𝑡 > 𝑐 + 𝐿, 𝛼 > 0 (2.28)

é limitado a:

|𝜀| <

𝑀

𝐿𝛼|Γ(1 − 𝛼)| (2.29)

O mesmo acontece para a aproximação:

𝐷𝑐𝛼

𝑡 𝑓(𝑡) ≈ 𝐷𝑡+𝐿

𝛼𝑡 𝑓(𝑡), 𝑡 < 𝑐 − 𝐿, 𝛼 > 0 (2.30)

2.4. Transformada de Laplace no cálculo fracionário

2.4.1. Riemann-Liouville e Grünwald-Leitnikoff

Propriedade 2.12

A transformada de Laplace aplicada nas definições de Riemann-Liouville (2.4) é dada por:

ℒ[ 𝐷0 𝑡𝛼𝑓(𝑡)] =

{

𝑠𝛼𝐹(𝑠), 𝑠𝑒 𝛼 ∈ ℝ−

𝐹(𝑠), 𝑠𝑒 𝛼 = 0

𝑠𝛼 − ∑ 𝑠𝑘 𝐷𝑡𝛼−𝑘−1𝑓(0)0

⌈𝛼⌉−1

𝑘=0

, 𝑠𝑒 𝛼 ∈ ℝ+ (2.31)

Nas condições da propriedade 2.8 a transformada de Laplace de D quando aplicada a

definição de Grünwald-Letnikoff (2.22) também é dada por (2.31) [43].

2.4.2. Caputo

Propriedade 4.2.1

A transformada de Laplace aplicada à definição de Caputo é dada por:

ℒ[ 𝐷0 𝑡𝛼𝑓(𝑡)] =

{

𝑠𝛼𝐹(𝑠), 𝑠𝑒 𝛼 ∈ ℝ−

𝐹(𝑠), 𝑠𝑒 𝛼 = 0

𝑠𝛼 − ∑ 𝑠𝛼−𝑘−1 𝐷𝑡𝑘𝑓(0)0

⌈𝛼⌉−1

𝑘=0

, 𝑠𝑒 𝛼 ∈ ℝ+

(2.32)

9

3. Funções de transferência

Pelas características do projeto em questão é interessante estudar-se a dinâmica do sistema no

que diz respeito a sua estabilidade. Para ter-se uma melhor compreensão sobre o comportamento de

sua estabilidade para as diversas ordens de derivadas, esse capítulo realiza uma análise através da

representação de entrada e saída, ou seja, por funções de transferência.

Primeiramente é apresentado a definição da função de transferência para sistemas que

apresentam equações diferenciais fracionarias, tanto para modelos de uma entrada e uma saída

(SISO), quanto para de múltiplas entradas e múltiplas saídas (MIMO). Depois mostra-se o mapa de

estabilidade no plano complexo e como que a mudança na ordem do termo derivativo pode levar um

sistema à instabilidade.

Por fim, é realizado um estudo prático sobre uma função de transferência de um modelo dinâmico

e sua mudança de comportamento para diversas ordens de derivada.

Esse capítulo também visa apenas ilustrar os conceitos sem uma grande profundidade e as

demonstrações, bem como outras propriedades podem ser encontradas em [31,33,43].

3.1. Definições

3.1.1. Função de transferência SISO

Dado um sistema com apenas uma entrada e uma saída (SISO), onde ambos estão

relacionados por uma equação diferencial linear e invariante no tempo, a sua função de transferência

é dada pela razão da transformada de Laplace da saída e da entrada (caso as condições iniciais

sejam nulas). Essa definição também se aplica caso a equação diferencial possua ordens

fracionárias.

Definição 3.1

𝐺(𝑠) =

∑ 𝑏𝑘𝑠𝛽𝑘𝑚

𝑘=1

∑ 𝑎𝑘𝑠𝛼𝑘𝑛

𝑘=1

(3.1)

onde 𝛼𝑘, 𝛽𝑘 ≥ 0 são as ordens do denominador e numerador, respectivamente.

3.1.2. Função de transferência MIMO

Definição 3.2

Um sistema MIMO com m entradas e p saídas, linear e invariante no tempo, pode ser

interpretado como uma matriz 𝑝 × 𝑚 onde os seus elementos serão funções de transferência

fracionárias SISO. Cada elemento relaciona uma entrada com uma saída onde as outras entradas

são zero.

Essa condição só é possível caso as condições iniciais sejam nulas, uma vez que a definição

de função de transferência exige condições iniciais nulas.

3.2. Estabilidade

Propriedade 3.1

Dado um sistema G(s), pela definição (3.1) , ele é estável se:

10

∀𝑠:∑𝑎𝑘𝑠

𝛼𝑘

𝑛

𝑘=0

= 0, |∠𝑠| >𝜋

2 (3.2)

onde ∠𝑠 ∈ [−𝜋,+𝜋] [31].

Figura 2. (a) Mapa de estabilidade para uma função de transferência de α=q menor que 1 e (b) para α=q maior

que 1. (extraído de [31]).

Vale notar que para sistemas com ordem de derivada inteira tem-se o critério clássico de

estabilidade, onde todos os polos precisam ter a parte real negativa.

Para ordens onde 𝛼 > 2 o sistema é necessariamente instável. É por essa razão que o

projeto focará apenas os casos onde 0 < 𝛼 < 2.

3.3. Resposta no domínio do tempo e da frequência

3.3.1. Resposta em frequência de uma função de transferência genérica

Propriedade 3.2

Quando uma entrada 𝑢(𝑡) = 𝐴 sin(𝜔𝑡) é aplicada a um sistema estável G(s), dado por:

𝐺(𝑠) =

∑ 𝑏𝑘𝑠 𝑘𝛼𝑚

𝑘=0

∑ 𝑎𝑘𝑠𝑘𝛼 𝑛

𝑘=0

(3.3)

a saída, após ter passado pelo regime transiente, será 𝑦(𝑡) = |𝐺(𝑗𝜔)|𝐴 sin(𝜔𝑡 + ∠𝐺(𝑗𝜔)).

Esse resultado significa que a resposta em frequência de uma função de transferência

fracionária pode ser encontrada ao substituir s por jω, do mesmo modo que nos casos inteiros.

3.3.2. Estudo de 𝟏 (𝒔

𝒂)𝟐𝜶

+ 𝟐𝜻 (𝒔

𝒂)𝜶

+ 𝟏⁄

As raízes da função de transferência dada por:

𝐺(𝑠) =

1

(𝑠𝑎)2𝛼

+ 2𝜁 (𝑠𝑎)𝛼

+ 1, 𝛼 ∈]0,2[ (3.4)

são 𝑠𝛼 = 𝑎𝛼(−𝜁 ± √𝜁2 − 1). Logo, a estabilidade dessa função de transferência dependerá do valor 𝜉.

Quando 𝜁 ≥ 1, a função de transferência será sempre estável.

Quando 𝜁 ≤ −1, o sistema sempre será instável.

Quando 𝜁 = 0, as raízes são complexas e iguais a ±𝑗. Com isso, o sistema será estável se

𝛼 < 1.

Quando 0 < 𝜁 < 1 e 0 < 𝛼 ≤ 1, as raízes são:

11

𝑠𝛼 = 𝑎𝛼 (−𝜁 ± √𝜁2 − 1) = 𝑎𝛼𝑒

±𝑗 tan−1√1−𝜁2

−𝜁 (3.5)

ou seja, a função de transferência sempre será estável.

Quando −1 < 𝜁 < 0 e 0 < 𝛼 < 1 pode-se provar que o sistema será estável se −cos𝛼𝜋

2< 𝜁.

Quando 0 < 𝜁 < 1 e 1 < 𝛼 < 2 pode-se provar que o sistema será estável se 𝜁 > − cos𝛼𝜋

2.

Quando −1 < 𝜁 < 0 e 1 ≤ 𝛼 < 2 o sistema é instável.

Todos os casos apresentados acima podem ser resumidos no seguinte formato: o sistema será

estável se

𝜁 > − cos𝛼𝜋

2 (3.6)

Figura 3. Comportamento da frequência (extraído de [43]).

A figura acima mostra o mapa com as localizações das regiões onde ocorre nenhum, um ou

duas frequências de ressonância (o que corresponde a nenhum, um ou dois máximos locais no

ganho do diagrama de Bode) em função dos parâmetros 𝜁 e 𝛼. Vale destacar que para o caso de

uma função de transferência de ordem inteira ela terá uma ou nenhuma frequência de ressonância,

como era de se esperar [43].

A figura 4 mostra a resposta da função de transferência dada por:

𝐺(𝑠) =

1

𝑠2𝛼 + 3𝑠𝛼 + 1 (3.7)

a um degrau unitário para diversas ordens de α. Tal modelo possui 𝜁 > 1 e, portanto, ambas as raízes

reais, de modo que ele sempre será estável.

Desse modo fica bem nítida a mudança de comportamento do sistema. Em ordens menores

do que um percebe-se um maior amortecimento do sistema, onde a resposta é muito mais lenta e

sem oscilações. Ao subir o valor de α, o sistema passa a apresentar uma resposta mais rápida até o

caso em que ocorre oscilações. Tal simulação fora realizada com o uso do software Simulink e da

toolbox FOMCON [42].

12

Figura 4. Resposta do sistema 𝐺(𝑠) = 1 𝑠2𝛼 + 3𝑠𝛼 + 1⁄ a um degrau unitário para diversas ordens de derivada.

13

4. Método aplicado

Atualmente existem diversos métodos matemáticos que buscam calcular os valores de derivadas

e integrais para ordens não inteiras [33].

O método proposto por [34] une a técnica de diferenciação numérica de ordem inteira com as

definições de derivadas fracionárias através do uso de matrizes triangulares e da discretização da

variável independente. Com isso, é possível resolver numericamente equações diferenciais de ordens

inteiras e fracionárias.

Tal método apresenta alguns benefícios, como:

Uma abordagem uniforme à discretização de derivadas de ordem arbitrária reais, incluindo

ordens inteiras e derivadas esquerda e direita.

Uma abordagem uniforme à resolução numérica de equações diferenciais de ordem inteira e

fracionária.

Solução para problemas de equações diferenciais de ordem não inteira que possuem valores

iniciais e de fronteira diferentes de zero.

Possibilidade de solucionar equações diferenciais não lineares de ordem fracionária.

Pode-se resolver equações diferenciais parciais de ordem não inteira e equações com

atraso.

Como adversidade, pode-se destacar:

Simulações que visam um longo horizonte temporal resultam em matrizes de dimensões

muito elevadas e, consequentemente, um alto custo computacional.

Erros de aproximação numérica.

A solução de equações diferenciais parciais resulta em matrizes de dimensões muito

elevadas que podem ser maiores do que a capacidade de registro e armazenamento do

software utilizado.

4.1. Matrizes triangulares

O método de [34] utiliza matrizes que possuem uma estrutura especial. Tais matrizes são as

matrizes triangulares inferior e as matrizes triangulares superior e possuem o seguinte formato,

respectivamente:

𝐿𝑁 =

[ 𝜔0 0 0 0 ⋯ 0𝜔1 𝜔0 0 0 ⋯ 0𝜔2 𝜔1 𝜔0 0 ⋯ 0⋱ ⋱ ⋱ ⋱ ⋯ ⋯

𝜔𝑁−1 ⋱ 𝜔2 𝜔1 𝜔0 0𝜔𝑁 𝜔𝑁−1 ⋱ 𝜔2 𝜔1 𝜔0]

(4.1)

𝑈𝑁 =

[ 𝜔0 𝜔1 𝜔2 ⋱ 𝜔𝑁−1 𝜔𝑁0 𝜔0 𝜔1 ⋱ ⋱ 𝜔𝑁−10 0 𝜔0 ⋱ 𝜔2 ⋱0 0 0 ⋱ 𝜔1 𝜔2⋯ ⋯ ⋯ ⋯ 𝜔0 𝜔10 0 0 ⋯ 0 𝜔0 ]

(4.2)

Estas matrizes triangulares são completamente descritas por sua primeira coluna (linha).

Além disso, se as matrizes C e D são ambas triangulares superior (inferior) elas comutam entre si:

𝐶𝐷 = 𝐷𝐶 (4.3)

Denomina-se

Ω𝑁 = 𝐿𝑁 − 𝜔0𝐸, 𝜓𝑁 = 𝑈𝑁 −𝜔0𝐸 (4.4)

14

onde E é uma matriz identidade, pode-se escrever:

𝐿𝑁 = Ω𝑁 + 𝜔0𝐸, 𝑈𝑁 = 𝜓𝑁 +𝜔0𝐸 (4.4)

Define-se as matrizes (𝑁 + 1) × (𝑁 + 1) 𝐸𝑝+, 𝑝 = 1,…𝑁, com 1 na diagonal p acima da

diagonal principal e 0 em todo o resto, e 𝐸𝑝−, 𝑝 = 1,…𝑁, com 1 na diagonal p abaixo da diagonal

principal e 0 em todo o resto, pode-se provar que:

Ω𝑁 =∑𝜔𝑘𝐸𝑘

−

𝑁

𝑘=1

, 𝜓𝑁 =∑𝜔𝑘𝐸𝑘+

𝑁

𝑘=1

(4.5)

Ω𝑁𝑁+1 = Ο, 𝜓𝑁

𝑁+1 = Ο (4.6)

com isso, as matrizes inversas (𝐿𝑁)−1 e (𝑈𝑁)

−1 podem ser escritas na seguinte forma explicita:

(𝐿𝑁)−1 = 𝜔0

−1𝐸 − 𝜔0−2Ω𝑁 + 𝜔0

−3Ω𝑁−2 +⋯+ (−1)𝑁𝜔0

−𝑁−1Ω𝑁𝑁 (4.7)

(𝑈𝑁)−1 = 𝜔0

−1𝐸 − 𝜔0−2𝜓𝑁 + 𝜔0

−3𝜓𝑁−2 +⋯+ (−1)𝑁𝜔0

−𝑁−1𝜓𝑁𝑁 (4.8)

ao introduzir o seguinte polinômio,

𝜚𝑁(𝑧) = 𝜔0 + 𝜔1𝑧 + 𝜔2𝑧2 +⋯+ 𝜔𝑁𝑧

𝑁 (4.9)

e as equações (4.7) e (4.8), tem-se:

𝜚𝑁(𝐸1−) = 𝜔0𝐸 + 𝜔1𝐸1

− + 𝜔2(𝐸1−)2 +⋯+𝜔𝑁(𝐸1

−)𝑁 = 𝐿𝑁 (4.10)

𝜚𝑁(𝐸1+) = 𝜔0𝐸 + 𝜔1𝐸1

+ + 𝜔2(𝐸1+)2 +⋯+𝜔𝑁(𝐸1

+)𝑁 = 𝑈𝑁 (4.11)

Define-se a operação de truncamento, 𝑡𝑟𝑢𝑛𝑐𝑁(∙), que trunca a série de potências 𝜚(𝑧),

𝜚(𝑧) = ∑𝜔𝑘𝑧

𝑘

∞

𝑘=0

(4.12)

no polinômio 𝜚𝑁(𝑧),

𝑡𝑟𝑢𝑛𝑐𝑁(𝜚(𝑧)) ≝∑𝜔𝑘𝑧

𝑘

∞

𝑘=0

= 𝜚𝑁(𝑧) (4.13)

onde as seguintes propriedades são notadas:

𝑡𝑟𝑢𝑛𝑐𝑁(𝛾𝜆(𝑧)) = 𝛾𝑡𝑟𝑢𝑛𝑐𝑁(𝜆(𝑧)) (4.14)

𝑡𝑟𝑢𝑛𝑐𝑁(𝜆(𝑧) + 𝜇(𝑧)) = 𝑡𝑟𝑢𝑛𝑐𝑁(𝜆(𝑧)) + 𝑡𝑟𝑢𝑛𝑐𝑁(𝜇(𝑧)) (4.15)

𝑡𝑟𝑢𝑛𝑐𝑁(𝜆(𝑧)𝜇(𝑧)) = 𝑡𝑟𝑢𝑛𝑐𝑁(𝑡𝑟𝑢𝑛𝑐𝑁(𝜆(𝑧))𝑡𝑟𝑢𝑛𝑐𝑁(𝜇(𝑧))) (4.16)

Através desses resultados, é possível demonstrar que operações como adição, subtração,

multiplicação e inversão entre duas matrizes triangulares superiores (inferiores) podem ser

expressadas na forma de operações entre suas respectivas séries de potências.

4.2. Produto Kronecker de matrizes

O produto Kronecker 𝐴⨂𝐵 de duas matrizes 𝐴𝑛×𝑚 e 𝐵𝑝×𝑞,

𝐴 = [

𝑎11 𝑎12 ⋯ 𝑎1𝑚𝑎21 𝑎22 ⋯ 𝑎2𝑚⋮ ⋮ ⋱ ⋮𝑎𝑛1 𝑎𝑛2 ⋯ 𝑎𝑛𝑚

] , 𝐵 =

[ 𝑏11 𝑏12 ⋯ 𝑏1𝑞𝑏21 𝑏22 ⋯ 𝑏2𝑞⋮ ⋮ ⋱ ⋮𝑎𝑝1 𝑎𝑝2 ⋯ 𝑏𝑝𝑞]

(4.17)

é uma matriz 𝑛𝑝 × 𝑚𝑞 que possui a seguinte estrutura

𝐴⨂𝐵 = [

𝑎11𝐵 𝑎12𝐵 ⋯ 𝑎1𝑚𝐵𝑎21𝐵 𝑎22𝐵 ⋯ 𝑎2𝑚𝐵⋮ ⋮ ⋱ ⋮

𝑎𝑛1𝐵 𝑎𝑛2𝐵 ⋯ 𝑎𝑛𝑚𝐵

] (4.18)

15

Dentre as diversas propriedades do produto Kronecker, umas delas será muito útil: se as

matrizes A e B forem triangular inferior (superior), então 𝐴⨂𝐵 também é triangular inferior (superior).

Duas operações serão importantes para o método, que são 𝐸𝑛⨂𝐴𝑛×𝑚 e 𝐴𝑛×𝑚⨂𝐸𝑚 onde 𝐸𝑛 é

a matriz identidade de dimensões 𝑛 × 𝑛 e 𝐸𝑚 a de dimensões 𝑚 ×𝑚.

O resultado da primeira operação será uma matriz diagonal em blocos onde sua diagonal

repete a matriz A. Já o caso da segunda será uma matriz esparsa feita de 𝑛 × 𝑚 blocos diagonais.

Por exemplo:

𝐴 = [𝑎11 𝑎12 𝑎13𝑎21 𝑎22 𝑎23

] (4.19)

então,

𝐸2⨂𝐴 = [

𝑎11 𝑎12 𝑎13 0 0 0𝑎21 𝑎22 𝑎23 0 0 00 0 0 𝑎11 𝑎12 𝑎130 0 0 𝑎21 𝑎22 𝑎23

] (4.20)

𝐴⨂𝐸3 =

[ 𝑎11 0 0 𝑎12 0 0 𝑎13 0 00 𝑎11 0 0 𝑎12 0 0 𝑎13 00 0 𝑎11 0 0 𝑎12 0 0 𝑎13𝑎21 0 0 𝑎22 0 0 𝑎23 0 00 𝑎21 0 0 𝑎22 0 0 𝑎23 00 0 𝑎21 0 0 𝑎22 0 0 𝑎23]

(4.21)

4.3. Eliminadoras

Um outro tipo de matriz que será necessário é a matriz eliminadora. Tal matriz é obtida

através da eliminação de linhas da matriz identidade 𝐸𝑛×𝑛. A matriz 𝑆𝑟1,𝑟2,…,𝑟𝑘 é obtida omitindo as

linhas de número 𝑟1, 𝑟2, … , 𝑟𝑘 da matriz identidade 𝐸𝑛×𝑛.

Ao definir uma matriz 𝐴𝑛×𝑛, o produto de 𝑆𝑟1,𝑟2,…,𝑟𝑘𝐴 possuirá apenas as linhas de A de número

diferente de 𝑟1, 𝑟2, … , 𝑟𝑘. Por exemplo:

𝐴 = [

𝑎11 𝑎12 𝑎13𝑎21 𝑎22 𝑎23𝑎31 𝑎32 𝑎33

] , 𝑆1 = [0 1 00 0 1

] (4.22)

𝑆1𝐴 = [

𝑎21 𝑎22 𝑎23𝑎31 𝑎32 𝑎33

] , 𝐴𝑆1𝑇 = [

𝑎12 𝑎13𝑎22 𝑎23𝑎32 𝑎33

] , 𝑆1𝐴𝑆1𝑇 = [

𝑎22 𝑎23𝑎32 𝑎33

]

(4.23)

4.4. Derivadas fracionárias com matrizes triangulares superiores (inferiores)

4.4.1. Derivada fracionária esquerda

Considera-se uma função f(t), definida no intervalo [a,b], de modo que 𝑓(𝑡) ≡ 0 para 𝑡 < 𝑎,

tem-se que sua derivada fracionária de ordem real esquerda α (𝑛 − 1 ≤ 𝛼 < 𝑛) é:

𝐷𝑡𝛼 =

1

Γ(𝑛 − 𝛼)(𝑑

𝑑𝑡)𝑛

𝑎 ∫

𝑓(𝜏)𝑑𝜏

(𝑡 − 𝜏)𝛼−𝑛+1

𝑡

𝑎

, (𝑎 < 𝑡 < 𝑏) (4.24)

Através da metodologia de abordagem do método das diferenças finitas, define-se

primeiramente nós equidistantes de espaçamento ℎ: 𝑡𝑘 = 𝑘ℎ (𝑘 = 0,1, … , 𝑁), no intervalo [a,b], onde

16

𝑡0 = 𝑎 e 𝑡𝑁 = 𝑏. Aplica-se a aproximação recursiva da diferenciação fracionária de ordem α para os

pontos 𝑡𝑘, (tal aproximação baseia-se na definição de Grünwald-Letnikoff) tem-se:

𝐷𝑡𝑘𝛼 ≈

∇𝛼𝑓(𝑡𝑘)

ℎ𝛼= ℎ−𝛼∑(𝑗)−1 (

𝛼𝑗 ) 𝑓𝑘−𝑗, 𝑘 = 0,1, … , 𝑁

𝑘

𝑗=0

𝑎 (4.25)

Todos os N+1 pontos podem ser escritos na forma matricial:

[ ℎ−𝛼∇𝛼𝑓(𝑡0)

ℎ−𝛼∇𝛼𝑓(𝑡1)

ℎ−𝛼∇𝛼𝑓(𝑡2)⋮

ℎ−𝛼∇𝛼𝑓(𝑡𝑁−1)

ℎ−𝛼∇𝛼𝑓(𝑡𝑁) ]

=1

ℎ𝛼

[ 𝜔0

(𝛼) 0 0 0 ⋯ 0

𝜔1(𝛼) 𝜔0

(𝛼) 0 0 ⋯ 0

𝜔2(𝛼) 𝜔1

(𝛼) 𝜔0(𝛼) 0 ⋯ 0

⋱ ⋱ ⋱ ⋱ ⋯ ⋯𝜔𝑁−1

(𝛼) ⋱ 𝜔2(𝛼) 𝜔1

(𝛼) 𝜔0(𝛼) 0

𝜔𝑁(𝛼) 𝜔𝑁−1

(𝛼) ⋱ 𝜔2(𝛼) 𝜔1

(𝛼) 𝜔0(𝛼)]

[ 𝑓0𝑓1𝑓2⋮

𝑓𝑁−1𝑓𝑁 ]

(4.26)

𝜔𝑗(𝛼)

= (−1)𝑗 (𝛼𝑗 ) , 𝑗 = 0,1, … , 𝑁 (4.27)

Pode-se, portanto, definir a matriz 𝐵𝑁𝛼 como:

𝐵𝑁(𝛼)

=1

ℎ𝛼

[ 𝜔0

(𝛼) 0 0 0 ⋯ 0

𝜔1(𝛼) 𝜔0

(𝛼) 0 0 ⋯ 0

𝜔2(𝛼) 𝜔1

(𝛼) 𝜔0(𝛼) 0 ⋯ 0

⋱ ⋱ ⋱ ⋱ ⋯ ⋯𝜔𝑁−1

(𝛼) ⋱ 𝜔2(𝛼) 𝜔1

(𝛼) 𝜔0(𝛼) 0

𝜔𝑁(𝛼) 𝜔𝑁−1

(𝛼) ⋱ 𝜔2(𝛼) 𝜔1

(𝛼) 𝜔0(𝛼)]

(4.28)

Onde a sua multiplicação pela coluna de valores da função 𝑓𝑘 (𝑘 = 0,1, … , 𝑁) resulta na

aproximação da derivada fracionária 𝐷𝑎 𝑡𝑘𝛼 𝑓(𝑡), 𝑘 = 0,1, … , 𝑁. A função geradora de 𝐵𝑁

(𝛼) será:

𝛽𝛼(𝑧) = ℎ−𝛼(1 − 𝑧)𝛼 (4.29)

Como a matriz 𝐵𝑁𝛼 é triangular inferior e, considerando as correspondentes derivadas

fracionárias esquerda 𝐷𝑎 𝑡𝛼 e 𝐷𝑎

𝑡𝛽, pode-se determinar:

𝐷𝑎 𝑡𝛼 ( 𝐷𝑎

𝑡𝛽𝑓(𝑡)) = 𝐷𝑎

𝑡𝛼+𝛽

𝑓(𝑡) (4.30)

4.4.2. Derivada fracionária direita

Considera-se uma função f(t), definida no intervalo [a,b], de modo que 𝑓(𝑡) ≡ 0 para 𝑡 > 𝑏,

tem-se que sua derivada fracionária de ordem real direita α (𝑛 − 1 ≤ 𝛼 < 𝑛) é:

𝐷𝑏𝛼 =

(−1)𝑛

Γ(𝑛 − 𝛼)(𝑑

𝑑𝑡)𝑛

𝑡 ∫

𝑓(𝜏)𝑑𝜏

(𝜏 − 𝑡)𝛼−𝑛+1

𝑡

𝑎

, (𝑎 < 𝑡 < 𝑏) (4.31)

A abordagem é a mesma do tópico anterior, de modo a obter-se a matriz 𝐹𝑁(𝛼)

𝐹𝑁(𝛼)

=1

ℎ𝛼

[ 𝜔0

(𝛼) 𝜔1(𝛼) ⋱ ⋱ 𝜔𝑁−1

(𝛼) 𝜔𝑁(𝛼)

0 𝜔0(𝛼) 𝜔1

(𝛼) ⋱ ⋱ 𝜔𝑁−1(𝛼)

0 0 𝜔0(𝛼) 𝜔1

(𝛼) ⋱ ⋱⋯ ⋯ ⋯ ⋱ ⋱ ⋱0 ⋯ 0 0 𝜔0

(𝛼) 𝜔1(𝛼)

0 0 ⋯ 0 0 𝜔0(𝛼) ]

(4.32)

A função geradora da matriz 𝐹𝑁(𝛼)

é a mesma da matriz 𝐵𝑁(𝛼)

:

𝛽𝛼(𝑧) = ℎ−𝛼(1 − 𝑧)𝛼 (4.33)

17

4.5. Solução numérica de equações diferenciais

Nos tópicos anteriores foram apresentados as ferramentas e definições que serão utilizadas

para solucionar numericamente equações diferenciais fracionárias através do método proposto por

[34]. Vale ressaltar que a metodologia utilizada para definir o valor das derivadas em cada nó também

pode ser estendida para o cálculo de suas respectivas integrais de ordem fracionária. Esse tópico

abordará a metodologia para solucionar equações diferenciais de ordem fracionária.

Considere uma matriz (𝑁 + 1) × (𝑁 + 1) triangular inferior 𝐿𝑁 e uma triangular superior 𝑈𝑁 e,

ao numerar as linhas de 0 a N, pode-se adquirir as seguintes relações através do uso das matrizes

eliminadoras.

𝑆0 {

𝐿𝑁𝑈𝑁} 𝑆0

𝑇 = {𝐿𝑁−1𝑈𝑁−1

} (4.34)

𝑆𝑁 {

𝐿𝑁𝑈𝑁} 𝑆𝑁

𝑇 = {𝐿𝑁−1𝑈𝑁−1

} (4.35)

𝑆0,1,…,𝑘 {

𝐿𝑁𝑈𝑁} 𝑆0,1,…,𝑘

𝑇 = {𝐿𝑁−𝑘−1𝑈𝑁−𝑘−1

} (4.36)

𝑆𝑁−𝑘,𝑁−𝑘+1,…,𝑁 {

𝐿𝑁𝑈𝑁} 𝑆𝑁−𝑘,𝑁−𝑘+1,…,𝑁

𝑇 = {𝐿𝑁−𝑘−1𝑈𝑁−𝑘−1

} (4.37)

Ou seja, a multiplicação mantém a estrutura da matriz triangular e apenas reduz o seu

tamanho por k+1 linhas e k+1 colunas.

Para solucionar a equação diferencial é preciso primeiramente que as condições iniciais do

problema sejam utilizadas para transformar o problema em um com equações diferenciais de

condições iniciais nulas.

Depois, o sistema de equações algébricas é obtido substituindo todos os termos derivativos

(ordem inteira ou não) por suas respectivas matrizes discretas análogas (𝐵𝑁𝛼 para derivadas

esquerdas e 𝐹𝑁𝛼 para derivadas direitas).

Por exemplo, considere a equação diferencial fracionária de m termos, dado por:

∑𝑝𝑘(𝑡)𝐷

𝛼𝑘𝑦(𝑡) = 𝑓(𝑡),0 ≤ 𝛼1 < 𝛼2 < ⋯ < 𝛼𝑚

𝑛 − 1 < 𝛼𝑚 < 𝑛

𝑚

𝑘=1

(4.38)

onde 𝐷𝛼𝑘 é a derivada de ordem 𝛼𝑘 esquerda pela definição de Riemann-Liouville ou Caputo. Define-

se as seguintes variáveis:

𝑃𝑁(𝑘)= 𝑑𝑖𝑎𝑔(𝑝𝑘(𝑡0), 𝑝𝑘(𝑡1), … , 𝑝𝑘(𝑡𝑁)) (4.39)

𝑌𝑁 = (𝑦(𝑡0), 𝑦(𝑡1), … , 𝑦(𝑡𝑁))𝑇 (4.40)

𝐹𝑁 = (𝑓(𝑡0), 𝑓(𝑡1), … , 𝑓(𝑡𝑁))𝑇 (4.41)

com isso, pode-se escrever a equação (4.38) no seguinte formado:

∑𝑃𝑁

(𝑘)𝐵𝑁𝛼𝑘𝑌𝑁 = 𝐹𝑁

𝑚

𝑘=1

(4.42)

O sistema linear algébrico para determinar a resposta é obtido através do sistema (4.42)

omitindo as primeiras n linhas e substituindo o valor inicial nulo nas equações restantes. Isso é o

mesmo que realizar a seguinte operação com o uso dos eliminadoras:

{𝑆0,1,…,𝑛−1 {∑𝑃𝑁

(𝑘)𝐵𝑁𝛼𝑘

𝑚

𝑘=1

} 𝑆0,1,…,𝑛𝑇 } {𝑆0,1,…,𝑛−1𝑌𝑁} = 𝑆0,1,…,𝑛−1𝐹𝑁 (4.43)

18

A solução desse sistema, junto com os valores iniciais fornecidos, permite encontrar a

solução numérica da equação diferencial (4.38) sob condições iniciais nulas.

Caso 𝑝𝑘(𝑡) ≡ 𝑝𝑘, então o sistema (4.43) resulta na seguinte forma simplificada:

∑𝑝𝑘𝐵𝑁−𝑛

𝛼𝑘 {𝑆0,1,…,𝑛−1𝑌𝑁}

𝑚

𝑘=1

= 𝑆0,1,…,𝑛−1𝐹𝑁 (4.44)

Como é possível notar, um dos pontos positivos de se trabalhar com a matriz triangular é a

possibilidade de solucionar o problema através da resolução de um sistema de equações algébricas,

não havendo a necessidade de criar recursões para encontrar variáveis de valor desconhecido.

Para equações diferenciais não lineares o procedimento é o mesmo, basta apenas escrever uma

equação equivalente que possua condições iniciais nulas.

4.6. Solução numérica de equações diferenciais parciais

4.6.1. Introdução

Até o momento as equações diferenciais estudadas eram todas ordinárias, ou seja,

continham apenas funções de uma variável e suas respectivas derivadas. No entanto, parte desse

trabalho realiza o estudo do fenômeno da difusão ao longo de uma dimensão do osso e,

consequentemente, equações diferenciais parciais surgem.

A metodologia apresentada em [35,36] propõe uma solução numérica para equações

diferenciais parciais de ordem fracionárias. Trata-se de uma expansão da ideia apresentada em [34]

com o uso das matrizes triangulares para transformar a equação diferencial parcial em uma equação

algébrica de fácil solução.

4.6.2. Metodologia

A figura 5 mostra o modelo tradicional de discretização para a equação da difusão. Dois nós

na direção do tempo são utilizados para aproximar a derivada no tempo e três nós na direção do

espaço para a derivada no espaço.

Figura 5. Discretização para derivadas de ordem inteira. (extraído de [36]).

No entanto, como uma derivada de ordem fracionária não é um operador local, deve-se

considerar todos os nós, desde o passado, da sua respectiva variável. A figura 6 mostra a

representação da derivada parcial em relação ao tempo para uma rede de 5 nós para o tempo. A

figura 7 ilustra o caso onde exista derivada fracionária tanto em relação ao tempo, quanto ao espaço.

19

Figura 6. Discretização do tempo para derivadas de ordem fracionárias (extraído de [36]).

Figura 7. Discretização no tempo e no espaço para derivadas de ordem fracionárias bilaterais (extraído de [36]).

Considere os nós (𝑖ℎ, 𝑗𝜏), 𝑗 = 0,1,2, … , 𝑛 correspondendo a todos os nós em relação ao tempo

para a “i-ésima” discretização no espaço. Em [34] é demonstrado a seguinte relação para a derivada

fracionária de 𝑢(𝑥, 𝑡):

[𝑢𝑖,𝑛(𝛼)

𝑢𝑖,𝑛−1(𝛼)

… 𝑢𝑖,2(𝛼)

𝑢𝑖,1(𝛼)

𝑢𝑖,0(𝛼)] = 𝐵𝑛

(𝛼)[𝑢𝑖,𝑛 𝑢𝑖,𝑛−1 … 𝑢𝑖,2 𝑢𝑖,1 𝑢𝑖,0]𝑇 (4.45)

Para ser possível encontrar a derivada fracionária α em relação ao tempo de 𝑢(𝑥, 𝑡) para

todos os nós da malha, [36] propõe rearranjar todos os termos 𝑢𝑖,𝑗 no seguinte formato:

𝑢𝑛𝑚

= [𝑢𝑚,𝑛𝑢𝑚−1,𝑛⋯𝑢0,𝑛 𝑢𝑚,𝑛−1𝑢𝑚−1,𝑛−1⋯𝑢0,𝑛−1 … 𝑢𝑚,1𝑢𝑚−1,1⋯𝑢0,1 𝑢𝑚,0𝑢𝑚−1,0⋯𝑢0,0]𝑇 (4.46)

onde m é a discretização em relação ao espaço. Com isso, tem-se um vetor coluna onde cada linha

dessa coluna representa um nó em sua respectiva coordenada tempo-espaço.

A matriz que transforma o vetor 𝑈𝑛𝑚 no vetor 𝑈𝑡(𝛼)

de derivadas fracionárias de ordem α em

relação ao tempo pode ser obtida com o produto Kronecker da matriz 𝐵𝑛(𝛼)

com a matriz identidade

𝐸𝑚.

𝑇𝑚𝑛(𝛼)

= 𝐵𝑛(𝛼)⨂𝐸𝑚 (4.47)

Do mesmo modo a matriz que transforma o vetor 𝑈 no vetor 𝑈𝑥(𝛽)

de derivadas fracionárias de

ordem β em relação ao espaço pode ser obtida através do produto Kronecker entre a matriz

identidade 𝐸𝑛 com a matriz 𝑅𝑚(𝛽)

, que corresponde a uma matriz simétrica ordinária de Riesz [27].

𝑅𝑚(𝛽)

=1

ℎ𝛽

[ 𝜔0

(𝛽) 𝜔1(𝛽) 𝜔2

(𝛽) 𝜔3(𝛽) ⋯ 𝜔𝑚

(𝛽)

𝜔1(𝛽) 𝜔0

(𝛽) 𝜔1(𝛽) 𝜔2

(𝛽) ⋯ 𝜔𝑚−1(𝛽)

𝜔2(𝛽) 𝜔1

(𝛽) 𝜔0(𝛽) 𝜔1

(𝛽) ⋯ 𝜔𝑚−2(𝛽)

⋱ ⋱ ⋱ ⋱ ⋯ ⋯𝜔𝑚−1

(𝛽) ⋱ 𝜔2(𝛽) 𝜔1

(𝛽) 𝜔0(𝛽) 𝜔1

(𝛽)

𝜔𝑚(𝛽) 𝜔𝑚−1

(𝛽) ⋱ 𝜔2(𝛽) 𝜔1

(𝛽) 𝜔0(𝛽) ]

(4.48)

𝜔𝑘(𝛽)

=(−1)𝑘Γ(𝛽 + 1) cos(𝛽𝜋 2⁄ )

Γ(𝛽 2⁄ − 𝑘 + 1)Γ(𝛽 2⁄ + 𝑘 + 1), 𝑘 = 0,1, … ,𝑚 (4.49)

o que resulta em:

20

𝑆𝑚𝑛(𝛽)

= 𝐸𝑛⨂𝑅𝑚(𝛽)

(4.50)

Tendo essas aproximações para as derivadas parciais fracionárias para ambas as variáveis é

possível discretizar as equações de uma EDP e transformá-la em uma equação algébrica. Por

exemplo, a equação clássica da difusão, dada por:

𝐷0𝐶

𝑡𝛼𝑢 − 𝜆

𝜕𝛽𝑢

𝜕|𝑥|𝛽= 𝑓(𝑥, 𝑡) (4.51)

pode ser discretizada para a seguinte equação algébrica:

{𝐵𝑛(𝛼)⨂𝐸𝑚 − 𝜆𝐸𝑛⨂𝑅𝑚

(𝛽)}𝑢𝑛𝑚 = 𝑓𝑛𝑚 (4.52)

Vale ressaltar que esse método também exige que as condições iniciais e de fronteira sejam

nulos, caso contrário uma mudança de variável deve ser feita para que tal condição seja atendida.

Apesar de que o método discutido tenha abordado apenas duas dimensões (tempo e um

espaço), a ideia pode ser estendida para um número maior de dimensões, bastando apenas aplicar

repetidamente a representação das matrizes triangulares de derivadas fracionárias em combinação

com o produto Kronecker de matrizes.

4.7. Método do passo largo (“Large Step”)

A técnica apresentada anteriormente tem se mostrado eficiente, simples e de fácil aplicação.

No entanto, pode-se notar que ao buscar uma simulação em um intervalo de tempo elevado ou com

um nível de espaçamento muito baixo entre nós obter-se-ão matrizes de dimensões elevadas, fato

esse que pode vir a se tornar um problema.

Esse problema é causado porque as derivadas fracionárias não são operadores locais. Uma

das possíveis soluções para o problema é aplicar o princípio da memória curta (“Short Memory”)

[10,33]. No entanto, em [37] é apresentado um novo método, chamado de “Large Step”, onde o

problema de matrizes muito grandes pode ser solucionado.

O método consiste em solucionar a equação diferencial fracionária para um intervalo de

tempo pequeno (“Small Step”) de modo que as dimensões das matrizes dadas em [34,35] não sejam

um problema. Em seguida, soluciona-se novamente a equação diferencial fracionária para o intervalo

de tempo seguinte (não necessariamente da mesma dimensão que o “Small Step” anterior) tal que as

soluções obtidas do passo anterior são utilizadas e, com isso, tem-se o “Large Step”.

Por exemplo, dado o sistema abaixo:

𝐷0𝐶

𝑡𝛼𝑦(𝑡) = 𝑓(𝑦(𝑡), 𝑡), 0 < 𝛼 < 1 (4.53)

𝑦(0) = 0 (4.54)

pode-se obter facilmente a solução para um intervalo de tempo (0,a) através da metodologia de [33],

onde tem-se como valor final 𝑦𝑎 para 𝑡 = 𝑎. O problema (4.53) pode ser transformado em:

𝐷𝑎𝐶

𝑡𝛼𝑦(𝑡) = 𝑓(𝑦(𝑡), 𝑡) − 𝑃0

𝑎𝛼𝑦(𝑡), 𝑡 > 𝑎 (4.55)

𝑦(𝑎) = 𝑦𝑎 (4.56)

onde

𝑃0 𝑎𝛼𝑦(𝑡) =

1

Γ(1 − 𝛼)∫ (𝑡 − 𝜏)𝛼−1𝑦′(𝜏)𝑑𝜏𝑎

0

, 𝑡 > 𝑎 (4.57)

21

A expressão acima representa a contribuição do momento anterior de 𝑦(𝑡) no intervalo de

tempo (0,a) para a solução da equação diferencial em um novo intervalo (a,b).

Vale ressaltar que no caso da equação (4.55) as condições iniciais não são mais nulas e,

portanto, uma mudança de variável é necessária para que atenda as condições do método das

matrizes triangulares.

O método apresentado pode ser repetido quantas vezes for necessário até que o intervalo de

tempo requerido seja atingido. No caso de equações diferenciais fracionárias parciais o procedimento

é o mesmo.

Com isso é possível desviar de algumas limitações que o princípio da memória curta pode

causar e obter respostas com alta precisão para intervalos de tempo elevados.

22

5. Dinâmica do Osso

Esse capítulo possui o objetivo de apresentar como que funciona a dinâmica do osso, bem como

as células responsáveis. Primeiramente é mostrado o papel desempenhado por cada uma dessas

células, além de suas características, formação e morte.

O tópico seguinte explica como que funciona cada uma das cinco etapas do processo de

remodelamento do osso (ativação, reabsorção, reversão, formação e conclusão, respectivamente),

onde os principais agentes e seus efeitos são discutidos. Vale ressaltar que tais etapas ainda não são

totalmente compreendidas atualmente, de modo que algumas hipóteses são feitas [38].

Em seguida é apresentado a doença que será o foco dessa dissertação: a metástase óssea. Tal

doença provoca, dentre outros sintomas, uma queda na densidade óssea e em muitos casos pode

levar a morte. A compreensão sobre os mecanismos de reprodução de tais células cancerígenas,

bem como a sua interação com as células responsáveis pela dinâmica óssea é fundamental para o

desenvolvimento de novos tratamentos.

Por fim, é feito uma breve discussão sobre alguns métodos de tratamento que são atualmente

utilizados no combate a essa doença.

5.1. Regeneração do osso

O esqueleto é um órgão altamente especializado e dinâmico que sofre continua regeneração

[22]. Ele é constituído de células especializadas, tecidos de conexão, canais e cavidades. Durante o

processo de desenvolvimento, pequenas partes do osso de uma região são removidas e depositadas

em outro local, de modo a atingir o novo formato. Tal processo é chamado de modelamento.

Quando o osso atinge a maturidade, a regeneração continua de modo periódico no mesmo

local, trocando o osso antigo pelo novo [15]. Tal processo é chamado de remodelamento. O propósito

desse remodelamento no esqueleto adulto ainda não é totalmente conhecido mas sabe-se que ele

colabora para reparar danos por fadiga ou envelhecimento.

A remoção do osso é a tarefa desempenhada pelas células chamadas osteoclastos através

da acidificação e digestão proteolítica, enquanto que formação de um novo osso é feito pelos

osteoblastos através da secreção de osteóides, que serão convertidos, eventualmente, em um novo

osso. No entanto, tais tarefas não são feitas separadamente, em momentos distintos, de modo

independente, mas sim de modo conjunto.

5.2. Formação e morte dos osteoclastos e osteoblastos

Ambas as células são geradas através da medula óssea. Os precursores dos osteoblastos

são células mesenquimais (núcleo alongado e cromatina condensada) multipotentes, enquanto que

no caso dos osteoclastos, são células hematopoiéticas mononucleadas. Enquanto que os

precursores dos osteoclastos atingem o osso através da circulação, os dos osteoblastos comumente

chegam ao osso pela migração dos progenitores de tecidos conectivos vizinhos [4,40].

A medula óssea, através dos fatores de crescimento [1] e citocinas [7], controla o

desenvolvimento e diferenciação dos osteoclastos e osteoblastos. Além disso, alguns hormônios e

sinais [38] também podem estimular ou inibir o desenvolvimento dessas células.

23

Os precursores dos osteoclastos apresentam c-fms (fator estimulador de colônias de

macrófagos) e RANK (receptor ativador de NF-B). Na presença de CSF-1 (fator estimulador de

colônia 1) e RANKL (ligante de RANK), que são os respectivos ligantes de c-fms e RANK, podem se

diferenciar em osteoclastos. O fim da vida do osteoclasto é causado pelo fenômeno chamado de

apoptose (autodestruição da célula), e então é removido através de fagócitos [38].

No caso dos osteoblastos, suas células precursoras são estimuladas através proteínas

morfogenéticas do osso (BMP), vitamina D e alguns outros fatores. Os osteoblastos possuem o

receptor ativo para PTH (hormônio paratireoide) e RANKL, que é o ligante de RANK, na sua

superfície externa e, na presença de c-fms, promovem a fusão de diversos monócitos em

osteoclastos. Os osteoblastos ainda produzem osteoprotegerina (OPG), que inibe a formação de

osteoclastos. Osteoblastos também morrem de apoptose ou se diferenciam em osteócitos, que

possuem RANKL e, consequentemente, contribuem para o remodelamento do osso [38].

5.3. Remodelamento do osso

No esqueleto adulto todos os osteoclastos e osteoblastos pertencem a uma única estrutura

temporária, conhecida como uma unidade multicelular básica (BMU). Em outras palavras, BMU é um

conjunto de células que criam um compartimento de remodelamento do osso [15,38]. Tal processo de

remodelamento possui 5 fases, respectivamente: ativação, reabsorção, reversão, formação e

conclusão, que serão explicadas com base em [38].

Ativação: O primeiro estágio do remodelamento do osso consiste na identificação de um sinal

para o início do processo. Tal sinal pode ser de diversas formas, como dano estrutural no

osso ou hormônios (PTH e estrogênio). Em condições normais, osteócitos secretam TGF- 𝛽

(fator de transformação de crescimento 𝛽), o que inibe a formação de osteoclastos. Quando

ocorre uma danificação da matriz do osso, tais células sofrem apoptose, o que reduz os

níveis de TGF- 𝛽 e favorece a formação de osteoclastos. O hormônio PTH é um sinal

endócrino que visa manter os níveis de cálcio em homeostase. Os osteoblastos possuem

receptores para esse tipo de hormônio, o que acaba por gerar uma resposta que induz o

recrutamento de osteoclastos.

Reabsorção: Em resposta ao PTH, osteoblastos produzem a quimiocina MCP-1 (proteína

quimioatrativa de monócitos-1), o que é atrativa para os precursores dos osteoclastos. Além

disso, os níveis de CSF-1, RANKL e OPG também são regulados em resposta ao PTH

(aumento do primeiro e segundo e redução do terceiro) de modo a favorecer a proliferação de

osteoclastos.

Reversão: Nesta fase as lacunas no osso geradas pela reabsorção são preparadas para o

processo de formação do osso através da remoção de colágeno não digerido da fase

anterior. Tal tarefa é desempenhada por células mesenquimais e macrófagos.

Formação: Nesta etapa ocorre a formação de osso nos locais onde ele, anteriormente, sofreu

reabsorção. No entanto, os mecanismos e sinais que coordenam essa transição não são

totalmente compreendidos. Fatores liberados no osso durante a reabsorção como IGF-I,II

(fator de crescimento de insulina I e II, respectivamente) e TGF- 𝛽 podem estar envolvidos

nessa etapa, bem como a hipótese de que osteoclastos produzem fatores de acoplamento,

responsáveis por atrair osteoblastos até a região em questão. Tais fatores incluem

semaforinas e efrinas.

Conclusão: Quando uma quantidade igual ao do osso reabsorvido fora formada, o ciclo de

remodelamento está concluído e uma nova fase está pronta para ser iniciada. No entanto, os

sinais terminais que informam o fim do processo ainda não são totalmente conhecidos.

24

5.4. Tumor e seus efeitos

O surgimento do tumor é causado por um crescimento anormal de tecido ou divisão celular. O

crescimento do tumor pode fazer com que ele se espalhe em outros órgãos além do seu local de

nascimento e agravar o estado da doença.

No caso do osso, o surgimento do tumor normalmente é causado por uma falta de balanço

entre a absorção e formação óssea. A metástase óssea pode ser classificada em osteoblástica,

quando a formação de osso sobrepõe-se e osteolítica, quando a absorção óssea supera a sua

formação. Em muitos casos as células de cancro são diagnosticadas em pacientes antes de

diagnosticar a doença primária. Sendo assim, um melhor entendimento da especificidade e

patogênese da metástase irá contribuir para melhores terapias e qualidade de vida ao paciente [5].

Nesse trabalho será abordado a metástase osteolítica referente ao mieloma múltiplo.

A metástase de um tumor primário para outros órgãos requer uma sequência de

acontecimentos. A proliferação do tumor primário é causada por agentes autócrinos, que são

hormônios produzidos pela célula cancerígena e que agem nela mesma, ou por fatores locais de

crescimento, como TGF- 𝛽. Após as células se separarem do tumor primário elas invadem os vasos

sanguíneos e seguem para novos alvos.

Depois que as células cancerígenas alcançam o esqueleto, elas aderem ao osso e iniciam

uma nova colônia. Uma vez no osso, essas células passam a interagir com as outras células já

presentes (osteoclastos, osteoblastos e etc) promovendo a absorção do osso. Tal processo ocorre

através do mecanismo RANKL-dependente e RANKL-independente.

Como discutido anteriormente, durante a absorção do osso ocorre a produção de TGF- 𝛽, o

que acaba por estimular a proteína relacionada a PTH (PTHrP) nas células metásticas. Tal proteína

combina-se com os receptores PTH dos osteoblastos, reforçando a secreção de RANKL e,

consequentemente, a formação de osteoclastos e uma maior taxa de absorção do osso. Como

consequência, a atividade dos osteoclastos irá produzir mais TGF- 𝛽 e criar um ciclo vicioso [6].

Além disso, o TGF- 𝛽 podem induzir a IL-8 (interleucina 8) nas células metásticas, o que

promove um aumento direto na população de osteoclastos, e DKK-1 (Dickkopf-1), uma substância

que se relaciona diretamente ao número e gravidade das lesões ósseas, além de impedir a formação

osteoblástica [5].

5.5. Tratamento

O tratamento da metástase óssea depende de quais sintomas o paciente está sofrendo, além

de sua saúde geral. É comum a equipa de médicos trabalhar com o auxílio do paciente de modo a

determinar o melhor plano de tratamento. Os principais objetivos do tratamento são a eliminação das

células cancerígenas, controlar o crescimento do tumor, diminuir a dor e permitir que o paciente tenha

uma vida ativa. No entanto, atualmente ainda não existe uma cura total para a metástase óssea,

apesar de ser possível controlar o crescimento do tumor [20].

Os tratamentos mais comuns são apresentados a seguir. Tais tratamentos podem ser

divididos em dois blocos, um onde o paciente não apresenta nenhum sintoma e outro onde apresenta

25

sintomas. Vale destacar que o plano de tratamento feito para uma pessoa que acabou de ser

diagnosticada pela primeira vez será diferente de um caso de reincidência.

Pacientes que apresentam a doença em estágio primário e sem sintomas (chamado de

Mieloma latente) devem ser monitorados constantemente e, caso algum sintoma apareça, o

tratamento deve ser iniciado imediatamente.

No caso de o paciente apresentar algum sintoma, o tratamento deve buscar controlar a

doença, além de permitir uma terapia suportável para o paciente de modo a melhorar a sua qualidade

de vida. Tais tratamentos normalmente incluem terapia com fármacos, como a quimioterapia, com ou

sem o uso de esteroides, e são denominadas terapias direcionadas.

5.6. Terapias direcionadas

As terapias direcionadas buscam atacar especificamente genes do cancro, proteínas ou

tecidos que permitem a instalação e desenvolvimento da doença. Tais métodos apresentaram bons

resultados de modo a estimular a investigação de novas tipos de fármacos.

Fármacos como a lenalidomida e pomalidomide são conhecidas por impedir o crescimento do

tumor no osso. Tais fármacos agem de modo a estimular o organismo a atacar as células

cancerígenas através do bloqueio na formação de vasos sanguíneos, o que causa uma dificuldade

dessas células obterem nutrientes [25,26].

O bortezomib e o carfilzomib são classificados como inibidores de proteassoma e atuam de

modo a inibir certas enzimas, chamadas de proteassoma, que digerem proteínas nas células. Como

as células do Mieloma produzem muitas proteínas, acabam por ficarem vulneráveis a tal ação [16].

A utilização de fármacos que visam destruir células cancerígenas é conhecida como

quimioterapia. Tais fármacos são comumente administradas de forma intravenosa ou através de

cápsulas via oral por um determinado número de ciclos ao longo do tempo. O paciente pode receber

um fármaco por vez ou uma combinação de diferentes tipos. Alguns exemplos desses tipos de

fármacos incluem melfalano, etoposídeo e carmustina [2,14].

Outros tipos de fármacos como o bifosfanato, visa fortalecer o osso e evitar fraturas, de modo

a auxiliar no tratamento dos sintomas causados pelo Mieloma. Esse fármaco previne o aumento de

cálcio no sangue e diminui o risco de hipercalcificação [12].

Caso o tratamento utilizado não obtenha sucesso, a doença pode ser classificada como

avançada ou terminal e é importante apresentar apoio ao paciente. Além disso, deve-se fazer com

que esse esteja confortável fisicamente e livre de dores.

26

6. A farmacocinética e a farmacodinâmica

A farmacocinética e a farmacodinâmica são dois assuntos muito importantes no que diz respeito

a metodologia de um tratamento terapêutico. Cada medicamento prescrito a um paciente deve ser

adequado a sua situação clínica.

A dosagem de um regime de medicamentos depende de um básico entendimento da sua ação.

Quando um paciente apresenta determinados sinais e sintomas, deve-se avalia-lo e realizar um

diagnóstico. Com isso, e com os conhecimentos da ação do medicamento, é possível prescrever um

determinado regime de fármacos apropriado para a cura da doença e o bem-estar do paciente.

Tendo isso em vista, o farmacêutico precisa obter os seguintes conhecimentos sobre o

tratamento:

Necessidade de se utilizar um fármaco.

Escolha do fármaco apropriado.

Objetivos da terapia.

Design do regime (rota, dose, frequência, duração e outros).

Monitoramento e revisão.

Aconselhamento.

Após a decisão do uso de um medicamento, deve-se escolher, através do conhecimento da

farmacocinética e farmacodinâmica, a formulação apropriada dos itens listados anteriormente. Tais

itens dependem de algumas características do paciente, como a taxa de absorção, excreção,

distribuição e metabolismo.

O objetivo é garantir que o regime atinja o máximo de eficiência com o mínimo de toxicidade. De

modo que o monitoramento é feito para que o paciente receba informações apropriadas sobre o seu

estado.

6.1. Farmacocinética

6.1.1. Introdução

A farmacocinética provém uma base matemática para modelar o percurso do fármaco e seus

efeitos no corpo ao longo do tempo. Tal modelo leva em consideração os parâmetros de absorção,

distribuição, metabolismo e excreção [11].

A efetividade de uma dose é determinada pela concentração de fármaco no corpo do

paciente. Idealmente, a concentração de fármaco deve ser medida no local de sua ação (receptor),

no entanto, devido a inacessibilidade, a concentração do fármaco é normalmente medida através do

sangue onde plasma é gerado. Em alguns casos a urina, saliva ou outro fluido podem ser utilizado

para medição, onde assume-se que a concentração do fármaco nesses fluidos está em equilíbrio com

a concentração do receptor.

6.1.2. Modelos

6.1.2.1. Taxas de processo

Existem algumas técnicas para representar a farmacocinética de um farmaco. A mais comum

é assumir o corpo como uma série de compartimentos na qual o fármaco flui e é excretado. A

transferência de fármaco entre compartimentos é representada por uma taxa constante, que pode ser

definida como de zero ou primeira ordem.

27

A representação do modelo de ordem zero para a eliminação de um fármaco do corpo é dado

pela seguinte equação:

𝑑𝐴

𝑑𝑡= −𝑘 (6.1)

onde A é a concentração do fármaco e k é a taxa constante de ordem zero.

Nesse caso a taxa de eliminação é constante e independe da concentração A do fármaco no

corpo.

Já a representação do modelo de primeira ordem é dada por:

𝑑𝐴

𝑑𝑡= −𝑘∗𝐴 (6.2)

onde 𝑘∗ é a taxa constante de primeira ordem.

Nesse caso a taxa de eliminação do fármaco A depende de sua concentração no corpo. A

grande parte dos fármacos é eliminado dessa maneira.

6.1.2.2. Compartimentos

O modelo da farmacocinética possui uma estrutura hipotética que é utilizada para descrever o

destino do fármaco em um sistema biológico.

No modelo de apenas um compartimento o corpo é analisado como uma unidade

homogênea. Ou seja, o fármaco se distribui instantaneamente por todo o corpo e atinge o equilíbrio

entre tecidos. Tal modelo é esquematizado pela figura 8.

Figura 8. Modelo de um compartimento (extraído de [11]).

Vale ressaltar que a concentração do fármaco no plasma (𝐶𝑝) não é, necessariamente, igual a

concentração nos tecidos. No entanto, mudanças no plasma refletem, quantitativamente, nos tecidos.

No caso de dois compartimentos o corpo é “dividido” em um compartimento central e outro

periférico. Apesar de não possuir um significado fisiológico ou anatômico, assume-se que o

compartimento central é composto por tecidos altamente perfundidos (como coração, cérebro e

outros) e o periférico por tecidos menos perfundidos (pele, gordura e músculo).

Tal modelo assume que o fármaco entrará primeiramente no compartimento central e se

espalhará através deste e do periférico. No entanto o fármaco não atinge distribuição instantânea

entre os compartimentos. O modelo é representado pela figura 9.

A figura 10 mostra como a concentração de fármaco no plasma se altera ao longo do tempo.

Inicialmente a concentração diminui rapidamente até o momento em que os compartimentos atingem

o equilíbrio.

28

Figura 9. Modelo de dois compartimentos (extraído de [11]).

Figura 10. Alteração da concentração de um fármaco no plasma no modelo de dois compapartimentos (extraído

de [11]).

6.1.3. Aplicações

A farmacocinética para um modelo de apenas um compartimento de absorção e eliminação

de primeira ordem para uma administração oral será tratada neste tópico.

A concentração de fármaco que ainda não foi absorvida, 𝐶𝑔(𝑡), e a concentração do fármaco

no plasma, 𝐶𝑝(𝑡), são descritos pelos seguinte sistema de equações diferenciais.

𝑑𝐶𝑔(𝑡)

𝑑𝑡= −𝑘𝑎𝐶𝑔(𝑡) (6.3)

𝑑𝐶𝑝(𝑡)

𝑑𝑡= 𝑘𝑎𝐶𝑔(𝑡) − 𝑘𝑑𝐶𝑝(𝑡) (6.4)

onde 𝑘𝑎 e 𝑘𝑑 são as taxas de absorção e eliminação, respectivamente.

Para uma administração oral de apenas uma dose de concentração inicial 𝐶𝑔(0) = 𝐶0 pode-

se obter a concentração do plasma, 𝐶𝑝(𝑡), através da transformada de Laplace.

𝐶𝑝(𝑠) = 𝐶0

𝑘𝑎(𝑠 + 𝑘𝑎)(𝑠 + 𝑘𝑑)

(6.5)

Pela equação acima, a resposta no tempo para uma administração de dose única é dada por:

𝐶𝑝(𝑡) = 𝐶0

𝑘𝑎𝑘𝑎 − 𝑘𝑑

(𝑒−𝑘𝑑𝑡 − 𝑒−𝑘𝑎𝑡) (6.6)

29

A figura 11 mostra como varia a concentração de um fármaco no plasma após uma única

administração oral. Os parâmetros foram 𝐶0 = 0.008, 𝑘𝑎 = 3 e 𝑘𝑑 = 0.5.

Figura 11. Modelo farmacocinético para a administração de uma dose oral (extraído de [9]).

Para uma administração de diversas doses, assumindo todas com concentrações iniciais 𝐶0 e

administradas em intervalos de tempo constante, 𝜏, a concentração do fármaco no plasma para a “n-

ésima” dose será:

𝐶𝑝(𝑛, 𝑡

′) = 𝐶0𝑘𝑎

𝑘𝑎 − 𝑘𝑑(1 − 𝑒−𝑛𝑘𝑑𝜏

1 − 𝑒−𝑘𝑑𝜏𝑒−𝑘𝑑𝑡

′−1 − 𝑒−𝑛𝑘𝑎𝜏

1 − 𝑒−𝑘𝑎𝜏𝑒−𝑘𝑎𝑡

′) (6.7)

onde 𝑡′ = 𝑡 − (𝑛 − 1)𝜏 representa o tempo corrido após a “n-ésima” dose.

Após um elevado número de doses a concentração atinge o seu estado estacionário, dado

por 𝐶𝑝𝑠𝑠 =𝐶0

𝜏𝑘𝑑.

A figura 12 mostra o modelo farmacocinético para múltiplas doses de administração diária

(𝜏 = 1, 𝐶0 = 0.008, 𝑘𝑎 = 3, 𝑘𝑑 = 0.5).

Figura 12. Modelo farmacocinético para múltiplas doses (extraído de [9]).

6.2. Farmacodinâmica

6.2.1. Introdução

A farmacodinâmica estuda a inter-relação da concentração de um fármaco e a estrutura alvo,

bem como o respectivo mecanismo de ação. Isso permite estabelecer as faixas de doses apropriadas

30

para os pacientes, bem como para comparar a potência, a eficácia e a segurança de um fármaco com

outro.

O estudo da farmacodinâmica baseia-se no conceito da ligação fármaco-receptor [17].

Quando um fármaco se liga a seu receptor, pode ocorrer uma resposta como consequência. Em

algum momento todos os receptores podem estar ligados ao fármaco, o que gera uma resposta

máxima.

Tal resposta quando é desencadeado por diversas células pode criar um efeito a nível de

órgão ou até mesmo do paciente. Consequentemente, é útil dispor de um modelo que possa

descrever a ligação do fármaco com o seu receptor para prever o seu efeito a níveis moleculares,

celulares e do organismo como um todo.

6.2.2. Modelamento

Do mesmo modo que em [32], o efeito do fármaco, 𝑑(𝑡), é descrito por uma função de Hill

dada por:

𝑑(𝑡) =

𝐶𝑝(𝑡)

𝐶50(𝑡) + 𝐶𝑝(𝑡) (6.8)

onde 𝐶50(𝑡) representa a concentração para 50% do efeito máximo. O valor de 𝑑 está entre 0 e 1, tal

que zero representa nenhum efeito do fármaco e um representa o efeito máximo.

A resistência ao fármaco pode ser descrita por:

𝐶50(𝑡) = 𝑓(𝑡)𝐶50𝑏𝑎𝑠𝑒 (6.9)

𝑓(𝑡) = 1 + 𝐾𝑟 ∫ max {0, 𝐿𝑟 − 𝐶𝑝(𝜏)}𝑑𝜏

𝑡

0

(6.10)

onde a concentração baixa de fármaco no plasma induz uma maior resistência ao fármaco. Com isso,

o mesmo efeito do fármaco só pode ser atingido com uma maior concentração no plasma.

A variável 𝐶50𝑏𝑎𝑠𝑒 representa o valor inicial de 𝐶50, enquanto que 𝐾𝑟 é a capacidade de aumento

da resistência ao fármaco no caso de uma queda na concentração do plasma. 𝐿𝑟 é o limiar de

concentração do plasma que permite um aumento na resistência ao fármaco.

A figura 13 mostra o efeito do fármaco ao longo do tempo. Pode-se notar que 𝐶50(𝑡) = 0.1.

Figura 13. Modelo farmacodinâmico de um fármaco (extraído de [9]).

31

6.2.3. Efeito de múltiplos fármacos

O uso de múltiplos fármacos pode afetar múltiplos alvos ou doenças simultaneamente. O uso

de múltiplos fármacos com mecanismos ou modos de ação distintos pode ser direcionado para um

único alvo, de modo a aumentar a eficiência do tratamento (ação sinergista).

Tal ação sinergista pode trazer boas consequências como:

O aumento do efeito terapêutico

Diminuição da dosagem e, consequentemente, da toxicidade

Diminuição da resistência ao fármaco

Tais benefícios estimularam o estudo a respeito dessa metodologia e sua aplicação para o tratamento

de doenças.

De modo a quantificar o efeito sinergista (ou antagonista) da combinação de fármacos, o

índice de combinação 𝐶𝐼, definido em [8], pode ser descrito pela seguinte equação:

𝐶𝐼 =

(𝐶𝑝∗)1

(𝐶𝑝)1+(𝐶𝑝

∗)2

(𝐶𝑝)2=

(𝐶𝑝∗)1

(𝐶50)1 𝑑 (1 − 𝑑)⁄+

(𝐶𝑝∗)2

(𝐶50)2 𝑑 (1 − 𝑑)⁄ (6.11)

onde (𝐶𝑝∗)𝑖 é a concentração necessária do fármaco 𝑖 quando combinada com outro fármaco para

produzir o mesmo efeito 𝑑 da concentração 𝐶𝑝𝑖 do fármaco 𝑖 quando tomada sozinha.

Se 𝐶𝐼 = 1, a combinação de fármacos apresenta efeito aditivo, o que significa que os efeitos

dos fármacos são mutualmente excludentes. Quando 0 < 𝐶𝐼 < 1, a combinação apresenta caráter

sinergista. Caso 𝐶𝐼 > 1, o efeito é antagonista.

Pode-se rearranjar a equação (6.11) de modo a isolar o termo referente ao efeito do fármaco

combinado, obtendo:

𝑑12 =[(𝐶𝑝

)1(𝐶50)1

+(𝐶𝑝

)2(𝐶50)2

]

𝐶𝐼 + [(𝐶𝑝

)1(𝐶50)1

+(𝐶𝑝

)2(𝐶50)2

]

(6.12)

Essa combinação de modelos farmacocinético e farmacodinâmico dá origem ao chamado

sistema PK/PD para um fármaco.

32

7. Estudo dos modelos dinâmicos

Neste capítulo é apresentado o estudo do modelo dinâmico para diversas situações clínicas, bem

como a ligação entre os diversos parâmetros das equações e os seus respectivos significados

biológicos, que também servirão de base para quando forem aplicadas as derivadas fracionárias.

Primeiramente o estudo será focado em um modelo com nenhuma dimensão espacial. Será

apresentado um sistema sem a presença do tumor, o que corresponde a uma pessoa sadia,

desenvolvido em [19]. Em seguida é feita a implementação da doença Mieloma e analise de suas

consequências para a dinâmica do processo, com base nos estudos feitos em [3] e [9]. Por fim, o

estudo de dois modelos de tratamento e suas respectivas consequências na dinâmica das células. O

primeiro modelo representará a aplicação dos inibidores de proteassoma, que visam inibir a apoptose

dos osteoblastos e eliminar as células cancerígenas [3]. O segundo modelo consiste na aplicação de

conceitos de farmacocinética e farmacodinâmica onde busca-se diminuir a atividade absorvedora dos

osteoclastos bem como a eliminação do tumor [9].

Em seguida, estende-se os modelos dinâmicos apresentados para o caso onde exista uma

dimensão espacial, bem como o fenômeno da difusão.

Por fim, os valores dos parâmetros que serão utilizados nas equações serão os mesmos, ou

próximos dos propostos por [19], [3] e [9]. Tais valores foram escolhidos por reproduzirem um

comportamento qualitativo do que se sabe da dinâmica das células. No entanto, vale ressaltar que

não há medições feitas para alguns desses parâmetros, uma vez que o estudo in vivo pode ser

invasivo, além de necessitar de um comprometimento por parte do paciente em ter o seu quadro

clínico acompanhado por um extenso período, o que dificulta bastante os estudos e medições de tais

parâmetros.

7.1. Caso sem tumor e sem tratamento 0D

O primeiro caso a ser analisado será a dinâmica da população de osteoclastos e osteoblastos

e seu resultado na densidade óssea de uma pessoa na ausência de tumor e tratamento. Para tal, fora

utilizado o modelo dinâmico proposto por [19] para um ponto discreto. Nele, tem-se uma equação

para a dinâmica dos osteoclastos, responsáveis pela destruição do osso, outra para os osteoblastos,

responsáveis pela reconstrução do osso, e, por fim, a equação que descreve a dinâmica da

densidade óssea. As equações são dadas por:

𝑑𝐶(𝑡)

𝑑𝑡= 𝛼1𝐶(𝑡)

𝑔11𝐵(𝑡)𝑔21 − 𝛽1𝐶(𝑡) (7.1)

𝑑𝐵(𝑡)


𝑔12𝐵(𝑡)𝑔22 − 𝛽2𝐵(𝑡) (7.2)

𝑑𝑧(𝑡)

𝑑𝑡= −𝑘1max{0, 𝐶(𝑡) − 𝐶𝑠𝑠} + 𝑘2max{0, 𝐵(𝑡) − 𝐵𝑠𝑠} (7.3)

Os números de células de osteoclastos 𝐶(𝑡) e de osteoblastos 𝐵(𝑡) são as variáveis do

sistema para um determinado instante t, e possuem 𝐶(0) = 𝐶0 e 𝐵(0) = 𝐵0 como valores iniciais. A

produção de cada uma dessas células afeta o desenvolvimento de si mesmas e/ou do outro tipo de

célula. Esses efeitos são denominados autócrinos e parácrinos, representadas por 𝑔11, 𝑔22, 𝑔12 e

𝑔21respectivamente. Tais parâmetros são adimensionais.

33

As características de RANK/RANKL/OPG estão relacionadas a 𝑔21 como a razão entre o

RANKL e OPG produzidos por osteoblastos, que tem efeito parácrino na população de osteoclastos.

Seu valor é assumido negativo para representar o efeito predominante do OPG sobre RANKL no

estado estacionário.

O efeito autócrino dos osteoblastos, 𝑔22, é considerado pequeno (pouca influência na

maturação de novos osteoblastos), de modo a possuir o valor 0 ou próximo disso.

Como o efeito parácrino dos osteoclastos regula a formação de osteoblastos, 𝑔12, seu valor é

assumido positivo. Por fim, o autócrino dos osteoclastos, 𝑔11, apresenta uma grande influência no

comportamento dinâmico do sistema.

As variáveis 𝛼1, 𝛼2, 𝛽1 e 𝛽2 representam a taxa de produção e morte de cada tipo de célula,

respectivamente e possuem como unidade 𝑑𝑖𝑎−1.

A densidade do osso (expresso em porcentagem) para um instante t é representado pela

variável 𝑧(𝑡) e é resultado da interação entre as quantidades de cada tipo de célula. Pelo modelo, a

taxa de destruição e formação do osso é proporcional ao número de osteoclastos e osteoblastos que

excedem o valor de seu estado estacionário 𝐶𝑠𝑠 e 𝐵𝑠𝑠, respectivamente. Os valores 𝑘1 e 𝑘2

representam a atividade de destruição e formação do osso e apresentam as unidades

% 𝑐é𝑙𝑢𝑙𝑎𝑠−1𝑑𝑖𝑎−1 onde a porcentagem indica a mudança em relação ao estado inicial.

O estado estacionário é obtido quando 𝑑𝐶(𝑡) 𝑑𝑡⁄ = 𝑑𝐵(𝑡) 𝑑𝑡⁄ = 0. Esse sistema possui o

estado estacionário não trivial dado pela seguinte expressão:

𝐶𝑠𝑠 = (𝛽1𝛼1)

(1−𝑔22)Ω

(𝛽2𝛼2)

𝑔21Ω

(7.4)

𝐵𝑠𝑠 = (𝛽1𝛼1)

𝑔12Ω

(𝛽2𝛼2)

(1−𝑔11)Ω

(7.5)

onde Ω = 𝑔12𝑔21 − (1 − 𝑔11)(1 − 𝑔22). A estabilidade do estado estacionário de um sistema não linear

para pequenas perturbações pode ser calculada analiticamente através de sua linearização por série

de Taylor. O cálculo do Jacobiano para os pontos de estado estacionário será:

𝐽(𝐶, 𝐵) = [

𝛼1𝑔11𝐶𝑔11−1𝐵𝑔21 − 𝛽1 𝛼1𝑔21𝐶

𝑔11𝐵𝑔21−1

𝛼2𝑔12𝐶𝑔12−1𝐵𝑔22 𝛼2𝑔22𝐶

𝑔12𝐵𝑔22−1 − 𝛽2] (7.6)

para a solução não trivial,

𝐽(𝐶𝑠𝑠, 𝐵𝑠𝑠) = [

𝛽1(𝑔11 − 1) 𝛽1𝑔12 𝐶𝑠𝑠 𝐵𝑠𝑠⁄

𝛽2𝑔21 𝐵𝑠𝑠 𝐶𝑠𝑠⁄ 𝛽2(𝑔22 − 1)] (7.7)

O sistema possuirá soluções periódicas se 𝑡𝑟(𝐽) = 0, portanto:

𝑡𝑟(𝐽) = Ψ = 𝛽1(𝑔11 − 1) + 𝛽2(𝑔22 − 1) = 0 (7.8)

Para Ψ < 0 o sistema terá oscilações amortecidas que convergirão para o estado estacionário

não trivial 𝐶𝑠𝑠 e 𝐵𝑠𝑠. Se Ψ > 0 (onde Ψ é próximo de 0) o sistema ficará instável e divergirá do estado

estacionário não trivial.

A figura 14 mostra a dinâmica dos osteoclastos, osteoblastos e da densidade óssea para uma

condição inicial de 𝐶0 = 𝐶𝑠𝑠 + 10, 𝐵0 = 𝐵𝑠𝑠 e 𝑧0 = 100%. Os valores dos parâmetros das equações

foram selecionados com base em [19] e são dados por 𝛼1 = 3 dia−1, 𝛼2 = 4 dia

−1, 𝛽1 = 0.2 dia−1, 𝛽2 =

0.02 dia−1, 𝑔11 = 0.5, 𝑔12 = 1, 𝑔21 = −0.5, 𝑔22 = 0. Para as variáveis da equação óssea tem-se 𝑘1 =

34

0.24% 𝑐é𝑙𝑢𝑙𝑎𝑠−1𝑑𝑖𝑎−1 e 𝑘2 = 0.0017% 𝑐é𝑙𝑢𝑙𝑎𝑠−1𝑑𝑖𝑎−1. Nesse caso Ψ = −0.12 mas não é

suficientemente perto de zero para exibir oscilações amortecidas.

Figura 14. Dinâmica das populações provocadas por um distúrbio na população de osteoclastos.

Pode-se notar que o distúrbio causado pelo maior número de osteoclastos provoca um

aumento na produção de osteoblastos e, consequentemente, uma diminuição na população de

osteoclastos. A densidade óssea sofre uma queda inicial causada pela ação destrutiva dos

osteoclastos mas retoma o seu valor inicial no momento que o sistema começa a convergir para o

seu estado estacionário. As linhas escuras mostram o valor do estados estacionários 𝐶𝑠𝑠 = 1.06 e

𝐵𝑠𝑠 = 212.1.

A figura 15 ilustra o resultado obtido para uma simulação com outros valores de parâmetros.

Os valores utilizados dessa vez também foram extraídos de [19] e são 𝛼1 = 3 dia−1, 𝛼2 = 4 dia−1, 𝛽1 =

0.2 dia−1, 𝛽2 = 0.02 dia−1, 𝑔11 = 1.1, 𝑔12 = 1, 𝑔21 = −0.5, 𝑔22 = 0, 𝑘1 = 0.0748% 𝑐é𝑙𝑢𝑙𝑎𝑠−1𝑑𝑖𝑎−1, 𝑘2 =

0.0006395% 𝑐é𝑙𝑢𝑙𝑎𝑠−1𝑑𝑖𝑎−1. Nesse caso ocorre oscilações estáveis periódicas entre as populações

de osteoclastos e osteoblastos (Ψ = 0). A oscilação é estimulada novamente por um aumento de dez

unidades na quantidade de osteoclastos acima do estado estacionário. O resultado é uma constante

remodelação na densidade óssea sobre o seu valor normalizado de 100%.

Figura 15. Dinâmica das populações com remodelação constante.

35

7.2. Caso com tumor e sem tratamento 0D

O modelo anteriormente estudado referia-se ao sistema dinâmico entre os osteoclastos e

osteoblastos na ausência de um tumor. De um modo resumido, pode-se dizer que um distúrbio

causado por um aumento da população de osteoclastos provoca o respectivo aumento da população

de osteoblastos. Esse aumento provoca uma queda no número de osteoclastos e o sistema busca

convergir para o estado estacionário. A densidade óssea segue essa dinâmica, tendo o seu valor

reduzido quando a proporção de osteoclastos, responsáveis pela destruição do osso, acima do

estado estacionário é superior ao de osteoblastos, responsáveis pela construção do osso, ou

aumentada em caso contrário.

Nesse tópico é feito o estudo do sistema dinâmico das células na presença da doença óssea

mieloma. Tal doença é discutida em mais detalhes no capítulo 5 mas, de uma maneira simplificada,

pode-se dizer que a presença do tumor altera o equilíbrio entre as células através de um aumento na

amplitude de oscilação dos osteoclastos e na diminuição da amplitude de osteoblastos. Como

resultado, tem-se uma queda da formação óssea e uma consequente queda de sua densidade.

[3] adiciona o tumor à dinâmica proposta em [19] de modo a afetar os parâmetros autócrinos

e parácrinos. A equação do tumor possui o formato de Gompertz, dado por:

𝑑𝑇(𝑡)

𝑑𝑡= 𝛾𝑇𝑇(𝑡) log (

𝐿𝑇𝑇(𝑡)

) (7.9)

onde 𝑇(𝑡) é a quantidade de células cancerígenas para um determinado tempo t. No entanto, tal

modelo não é definido em zero e pode ser substituído pelo seguinte modelo analítico equivalente [9]:

𝑑𝑇(𝑡)

𝑑𝑡= 𝑟𝑇𝑇(𝑡)

𝑎 (1 −𝑇(𝑡)

𝐿𝑇)

𝑏

(7.10)

onde 𝑎, 𝑏 𝑒 𝑟𝑇 são parâmetros adimensionais relacionados ao crescimento do tumor.

O novo sistema dinâmico adquire o seguinte formato, sendo 𝐶(𝑡) e 𝐵(𝑡) as densidades de

osteoclastos e osteoblastos, como anteriormente e 𝑧(𝑡) a densidade óssea, expressa em

porcentagem:

𝑑𝐶(𝑡)


𝑔11(1+𝑟11𝑇(𝑡)𝐿𝑇

)𝐵(𝑡)


)− 𝛽1𝐶(𝑡) (7.11)

𝑑𝐵(𝑡)


𝑔12 (1+𝑟12𝑇(𝑡)𝐿𝑇

)⁄𝐵(𝑡)

(𝑔22−𝑟22𝑇(𝑡)𝐿𝑇

)− 𝛽2𝐵(𝑡) (7.12)

𝑑𝑇(𝑡)

𝑑𝑡= 𝑟𝑇𝑇(𝑡)

𝑎 (1 −𝑇(𝑡)

𝐿𝑇)

𝑏

(7.13)

𝑑𝑧(𝑡)


Os parâmetros 𝛼1, 𝛼2, 𝛽1, 𝛽2, 𝑔11, 𝑔12, 𝑔21, 𝑔22, 𝑘1 e 𝑘2 são os mesmo de anteriormente. Os

novos parâmetros são 𝐿𝑇, que representa a capacidade máxima que o tumor pode atingir, 𝑟𝑇 que é a

velocidade de crescimento do tumor (assumiu-se que esse valor independe da perda de osso), 𝑎 e 𝑏

também são parâmetros que influenciam o crescimento e 𝑟11, 𝑟12, 𝑟21, 𝑟22 que são constantes positivas

responsáveis pela intensidade com que o tumor afeta os efeitos autócrinos e parácrinos. Tanto em [3]

quanto em [9] não são apresentadas as unidades para esses parâmetros, de modo que nessa

dissertação elas também estarão ausentes.

36

Nesse novo modelo pode-se notar que a presença de tumor promove um aumento no efeito

autócrino dos osteoclastos (𝑔11 (1 + 𝑟11𝑇(𝑡)

𝐿𝑇) > 𝑔11, se 𝑔11 > 0), devido ao mecanismo RANKL-

independente. Uma redução no efeito parácrino dos osteoclastos (−𝑔21 > −𝑔21 (1 + 𝑟21𝑇(𝑡)

𝐿𝑇) ,

assumindo 𝑔21 < 0 ), através do aumento de RANKL e redução de OPG (causada pela liberação de

PTHrP). Uma redução no parácrino dos osteoblastos (𝑔12 (1 + 𝑟12𝑇(𝑡)

𝐿𝑇)⁄ > 𝑔12, se 𝑔12 > 0) e uma

redução no efeito autócrino dos osteoblastos (𝑔22 − 𝑟22𝑇(𝑡)

𝐿𝑇< 𝑔22, se 𝑔22 > 0), através do inibidor da

diferenciação e ativação dos osteoblastos.

A presença do tumor faz com que o sistema tenha oscilações em torno de um novo estado

estacionário, dado por:


(1−𝑔22−𝑟22)Λ

(𝛽2𝛼2)

𝑔21(1+𝑟21)Λ

(7.15)


(𝑔12 (1+𝑟12)⁄ )Λ

(𝛽2𝛼2)

1−𝑔11(1+𝑟11)Λ

(7.16)

𝑇𝑠𝑠 = 𝐿𝑇 (7.17)

onde Λ = (𝑔12 (1 + 𝑟12)⁄ )(𝑔21(1 + 𝑔21)) − (1 − 𝑔11(1 + 𝑟11))(1 − 𝑔22 + 𝑟22).

O estudo da estabilidade desse novo sistema é novamente realizado através da análise do

Jacobiano para os novos pontos de equilíbrio. Seguindo a mesma metodologia utilizada em (colocar a

referência do jacobiano anterior), define-se a seguinte equação:

Φ = 𝛽1(𝑔11(1 + 𝑟11) − 1) + 𝛽2(𝑔22 − 𝑟22 − 1) (7.18)

O modelo terá ciclos estáveis se Φ = 0. Se Φ < 0 o sistema apresentará oscilações

amortecidas que convergirá para o estado estacionário e oscilações instáveis divergindo do estado

estacionário para o caso de Φ > 0 (onde Φ é próximo de 0).

A figura 16 mostra o resultado de uma simulação com Φ < 0. Os parâmetros utilizados na

equação foram 𝛼1 = 3 células/dia, 𝛼2 = 4 dia−1, 𝛽1 = 0.2 dia−1, 𝛽2 = 0.02 dia

−1, 𝑔11 = 1.1, 𝑔12 = 1,

𝑔21 = −0.5, 𝑔22 = 0, 𝑘1 = 0.0748% 𝑐é𝑙𝑢𝑙𝑎𝑠−1𝑑𝑖𝑎−1, 𝑘2 = 0.0006395% 𝑐é𝑙𝑢𝑙𝑎𝑠−1𝑑𝑖𝑎−1, 𝑎 = 2, 𝑏 = 3.1

𝑟11 = 0.005, 𝑟12 = 0.0, 𝑟21 = 0.0, 𝑟22 = 0.2, 𝑟𝑇 = 0.005 e 𝐿𝑇 = 100%. Os valores iniciais para cada

população foram 𝐶(0) = 15.0, 𝐵(0) = 316.0 e 𝑇(0) = 1. Tais valores foram implementados em [3] de

modo que o sistema apresentasse um comportamento qualitativo semelhante ao que se sabe sobre a

doença.

Nesse caso, o tumor se desenvolve ao longo do tempo e converge para o seu tamanho

máximo, 𝐿𝑇. A presença do tumor causa uma alteração no equilíbrio da população das células,

resultando em um aumento da amplitude das oscilações dos osteoclastos e uma diminuição na de

osteoblastos (se comparado ao modelo sem a presença da doença), o que é uma das características

da Mieloma.

No modelo com a doença, o número de osteoclastos e osteoblastos sofre um aumento e

converge para um novo estado estacionário não trivial 𝐶𝑠𝑠 = 5.0 e 𝐵𝑠𝑠 = 316.0. Esse aumento pode

ser entendido como uma tentativa do corpo humano de tentar buscar o equilíbrio novamente entre a

porcentagem de destruição e reconstrução da densidade óssea causado pelas células.

37

Figura 16. Dinâmica das populações na presença de um tumor.

A densidade óssea apresenta um ligeiro aumento inicial por conta do aumento populacional

citado anteriormente. No entanto, a medida que o tumor se desenvolve e afeta a amplitude das

oscilações das células, a densidade óssea diminui até atingir 0. Como reflexo dessa diminuição, a

amplitude das oscilações também reduz.

7.3. Caso com tumor e com tratamento de inibidores de proteassoma 0D

No tópico anterior fez-se o estudo do modelo dinâmico proposto por [3]. Através dele, é

possível analisar o comportamento dos osteoclastos e osteoblastos na presença do tumor para

derivadas de primeira ordem.

De um modo geral, pode-se dizer que o tumor provoca um aumento na amplitude de

oscilação dos osteoclastos e uma diminuição na de osteoblastos (se levado em comparação ao

modelo sem tumor). Esse desequilíbrio acentua-se cada vez mais a medida que o tumor se

desenvolve e, uma vez que as atividades de destruição do osso superam a de reconstrução, a

densidade óssea de um modo geral começa a sofrer uma redução.

A aplicação de um tratamento deve ser feita o quanto antes para evitar que o tumor atinja um

estado muito evoluído e a densidade óssea não seja muito comprometida. Tendo isso em foco, o

tratamento deve possuir duas características que são: eliminar a presença de tumor e permitir uma

restauração na densidade óssea, bem como no seu processo de remodelamento.

Existem muito tratamentos terapêuticos para a doença mieloma, tanto de modo a afetar as

células cancerígenas diretamente quanto indiretamente, através dos osteoclastos e osteoblastos. No

entanto, uma resposta geral para cada tipo de tratamento é algo difícil de se determinar através de

resultados de estudos in vivo.

Nesse tópico é feito o estudo do modelo de tratamento proposto em [3]. Nele é realizado o

modelamento dos efeitos dos inibidores proteassoma para a doença óssea Mieloma.

Os inibidores de proteassoma são caracterizados por apresentarem atuação direta nas

células cancerígenas, através de sua eliminação, pois são capazes de inibir a ação de certas enzimas

que deterioram proteínas responsáveis por manter o processo de divisão celular sob controle. Além

38

disso, afetam a formação de osteoblastos, através da diminuição de sua apoptose e,

consequentemente, aumento da formação óssea. Tal diminuição é causada por uma degradação dos

receptores de vitamina D dos proteassoma (VDR), o que implica em um aumento na diferenciação de

osteoblastos, que são estimulados por essa vitamina [13,16,21,29,30]. Com esse tratamento, as

equações apresentam o seguinte formato:

𝑑𝐶(𝑡)



)𝐵(𝑡)


)− 𝛽1𝐶(𝑡) (7.19)

𝑑𝐵(𝑡)


𝑔12 (1+𝑟12𝑇(𝑡)𝐿𝑇

)⁄𝐵(𝑡)


)− (𝛽2 − 𝑉1(𝑡))𝐵(𝑡) (7.20)

𝑑𝑇(𝑡)

𝑑𝑡= (𝑟𝑇 − 𝑉2(𝑡))𝑇(𝑡)

𝑎 (1 −𝑇(𝑡)

𝐿𝑇)

𝑏

(7.21)

𝑑𝑧(𝑡)


A diferença desse modelo para o estudado anteriormente é a presença dos parâmetros 𝑉1(𝑡)

e 𝑉2(𝑡), responsáveis por diminuir a apoptose dos osteoblastos e eliminar células cancerígenas,

respectivamente. Eles são definidos da seguinte maneira:

{

𝑉1(𝑡) = 0 se 𝑡 < 𝑡𝑠𝑡𝑎𝑟𝑡𝑉1(𝑡) = 𝑣1 se 𝑡 ≥ 𝑡𝑠𝑡𝑎𝑟𝑡𝑉2(𝑡) = 0 se 𝑡 < 𝑡𝑠𝑡𝑎𝑟𝑡𝑉2(𝑡) = 𝑣2 se 𝑡 ≥ 𝑡𝑠𝑡𝑎𝑟𝑡

(7.23)

onde 𝑡𝑠𝑡𝑎𝑟𝑡 indica o momento de início do tratamento e 𝑣1 e 𝑣2 são variáveis de intensidade do

tratamento.

A figura 17 mostra o resultado de uma simulação realizada com os mesmos parâmetros e

condições iniciais do tópico 7.2, apenas com a adição de 𝑣1 = 0.001 e 𝑣2 = 0.008, extraídos de [3].

Figura 17. Dinâmica das populações na presença de um tumor e tratamento.

Nesse caso, o tumor inicia-se em 𝑡 = 0 e desenvolve-se normalmente, através do modelo

estudado anteriormente, onde as células com tumor aumentam, junto com um aumento na amplitude

de oscilação dos osteoclastos e uma diminuição no de osteoblastos. A partir de 𝑡 = 𝑡𝑠𝑡𝑎𝑟𝑡 = 600 o

tratamento com os inibidores de proteassoma inicia-se e o tumor é reduzido gradativamente.

39

Com isso, a população de osteoclastos sofre uma diminuição, porém, pelos resultados da

simulação, não é suficiente para garantir um aumento na densidade óssea. Esse aumento só ocorre

após o número de osteoblastos começar a subir.

7.4. Caso com tumor e com tratamento utilizando PK/PD 0D

No tópico anterior foi apresentado um modelo de tratamento terapêutico para a doença óssea

Mieloma. Tal tratamento consiste na aplicação de inibidores de proteassoma, de modo a estimular

uma maior diferenciação da população de osteoblastos e formação do osso, ao mesmo tempo em

que inibe a diferenciação de células cancerígenas [13]. Os resultados foram condizentes com os

estudos feitos em [13,16,21,29,30], de modo que o tumor fora eliminado e o indivíduo recuperou a

dinâmica original das células.

No entanto, em [9] é proposto um outro modelo dinâmico de tratamento terapêutico, no qual,

além de incluir os modelos farmacocinéticos e farmacodinâmicos, utiliza outros tipos de fármacos.

Tendo em vista que o tratamento para uma doença tão severa como o Mieloma possui diversos

métodos de terapia [20], é interessante realizar o estudo de novas dinâmicas e os respectivos efeitos

no organismo de modo a melhorar a compreensão e escolha da terapia adequada a um determinado

paciente.

Diferentemente do modelo de terapia apresentado em [3], que busca um maior estímulo na

produção de osteoblastos, o modelo de [9] apresenta um tratamento através da inibição da atividade

de destruição do osso causada pelos osteoclastos, além da eliminação do tumor. Outro ponto

interessante é a inclusão dos conceitos de farmacocinética e farmacodinâmica, utilizados atualmente

de forma a melhorar a resposta do indivíduo ao tratamento.

A terapia que busca inibir a atividade dos osteoclastos costuma utilizar bisfosfanatos e

denosumab, que atuam de modos distintos na reabsorção do osso. Os bisfosfanatos são fármacos

desenvolvidas para o combate a osteoporose através da mineralização da superfície óssea, de modo

a impedir a absorção do osso por parte dos osteoclastos, além da liberação de algumas enzimas

específicas que causam um aumento na apoptose dessas células [27]. No caso da denosumab a sua

ação inibe RANK-ligante e, consequentemente, a maturação dos osteoclastos [6]. Com isso, as ações

de ambas os fármacos são complementares. Os bisfosfanatos inibem as atividades e estimulam a

morte de osteoclastos, enquanto que denosumab impede a diferenciação de novas células.

O tratamento do tumor é realizado através de quimioterapia e terapia de hormônios que visam

eliminar as células cancerígenas, de modo a eliminar a sua capacidade de crescimento e divisão [14].

Dependendo do modo de atuação do fármaco, ele pode ter um caminho inibitório ou

estimulador. Tais modelos podem ser descritos no seguinte formato, respectivamente:

𝐼𝑒(𝑡) = 1 − 𝐾𝑖𝑑(𝑡) (7.24)

𝑆𝑒 = 1 + 𝐾𝑠𝑑(𝑡) (7.25)

onde 𝐾𝑖 e 𝐾𝑠 são constantes adimensionais positivas, que representam o efeito máximo do fármaco

em determinado mecanismo e 𝑑(𝑡) é o efeito do modelo PK/PD para uma ou múltiplas doses do

fármaco apresentado no capítulo 6.

40

O sistema dinâmico do osso com o tumor e com o tratamento descrito acima apresenta o

seguinte formato:

𝑑𝐶(𝑡)



)𝐵(𝑡)


)(1+𝐾𝑠1𝑑1(𝑡)) − (1 + 𝐾𝑠2𝑑2(𝑡))𝛽1𝐶(𝑡) (7.26)

𝑑𝐵(𝑡)


𝑔12 (1+𝑟12𝑇(𝑡)𝐿𝑇

)⁄𝐵(𝑡)


)− 𝛽2𝐵(𝑡) (7.27)

𝑑𝑇(𝑡)

𝑑𝑡= (1 − 𝐾𝑖34𝑑𝑐34(𝑡))𝑟𝑇𝑇(𝑡)

𝑎 (1 −𝑇(𝑡)

𝐿𝑇)

𝑏

(7.28)

𝑑𝑧(𝑡)


O efeito farmacodinâmico ao longo do tempo do denosumab, bisfosfanatos e quimioterapia

(que é uma combinação de fármacos) estão representados por 𝑑1(𝑡), 𝑑2(𝑡) e 𝑑𝑐34(𝑡),

respectivamente.

O efeito da quimioterapia atua diretamente na eliminação das células cancerígenas e

impedindo que o tumor se desenvolva. Nesse modelo ela é interpretada como a combinação de dois

fármacos.

A figura 18 mostra o resultado da simulação para os mesmos parâmetros e condições iniciais

do tópico 7.2. O parâmetro de cada fármaco é dado pela tabela I.

Tabela I. Parâmetros dos fármacos para o modelo PK/PD.

Combinação 𝑑34

Fármaco 𝑑3 𝑑4 𝑑1 𝑑2

𝜏 1 1 1 1

𝑘𝑑 1 1 1 0.5

𝑘𝑎 0.01 0.01 0.01 0.01

𝐶0 0.008 0.008 0.008 0

𝐶0𝑏𝑎𝑠𝑒 0.002 0.002 0.002 0.002

𝐾𝑟 0 0 0 0

𝐿𝑟 0 0 0 0

𝐾𝑠,𝑖 2 0.1 0.1

𝐶𝐼 1 - -

Tais valores foram propostos em [9] de modo a eliminar o tumor e restaurar o ciclo de

remodelamento do osso.

41

Figura 18. Modelo de tratamento com o uso de PK/PD.

Pode-se notar que o tumor começa a se desenvolver e, no momento que atinge um certo

patamar, o tratamento com os fármacos inicia-se e, aos poucos, as células cancerígenas são

eliminadas, chegando a ficar apenas resíduos que devem acabar sendo destruidos no futuro. As

populações de osteoclastos e osteoblastos recuperam o seu ciclo normal, através do aumento da

apoptose dos osteoclastos causado pelos bisfosfanatos e do bloqueio de sua maturação, motivado

pelo denosumab. A densidade óssea começa a ser vagarosamente aumentada ao passar do tempo.

Tais resultados condizem com o objetivo da terapia e com estudos feitos com esses

fármacos, tanto no que diz respeito a eliminação da doença através da quimioterapia, quanto na

recuperação do ciclo dinâmico das células e da densidade [2,12,14,18,41].

Diferentemente do modelo de tratamento anterior, pode-se notar uma redução mais lenta do

tumor nos momentos iniciais da terapia. Isso deve-se ao fato de que as concentrações dos fármacos

no organismo ainda estarem baixas. Esse comportamento reflete os conceitos de farmacocinética e

farmacodinâmica e como eles estão relacionados aos processos terapêuticos.


O modelo apresentado em [19] propõe um modelamento matemático para a dinâmica das

células do osso. As equações em questão estudam a interação das células para um ponto discreto do

osso.

Em [3] estende-se o conceito apresentado para o caso de uma dimensão espacial. Essa nova

dimensão espacial é introduzida através de um modelo de difusão na região Ω = [0 1]. Nesse novo

modelo as varáveis são 𝐶(𝑡, 𝑥) para a população de osteoclastos, 𝐵(𝑡, 𝑥) para a de osteoblastos e

𝑧(𝑥, 𝑡) para a densidade do osso, em referência ao tempo t e distância x. Nesse caso as equações

ficam:

𝜕𝐶(𝑡, 𝑥)

𝜕𝑡= 𝜎1

𝜕2

𝜕𝑡2𝐶(𝑡, 𝑥) + 𝛼1𝐶(𝑡, 𝑥)

𝑔11𝐵(𝑡, 𝑥)𝑔21 − 𝛽1𝐶(𝑡, 𝑥) (7.30)

𝜕𝐵(𝑡, 𝑥)

𝜕𝑡= 𝜎2

𝜕2

𝜕𝑡2𝐵(𝑡, 𝑥) + 𝛼2𝐶(𝑡, 𝑥)

𝑔12𝐵(𝑡, 𝑥)𝑔22 − 𝛽2𝐵(𝑡, 𝑥) (7.31)

𝜕𝑧(𝑡, 𝑥)

𝜕𝑡= 𝜎3

𝜕2

𝜕𝑡2𝑧(𝑡, 𝑥) − 𝑘1(𝑥)max{0, 𝐶(𝑡, 𝑥) − 𝐶𝑠𝑠} + 𝑘2(𝑥)max{0, 𝐵(𝑡, 𝑥) − 𝐵𝑠𝑠} (7.32)

onde as condições iniciais são 𝐶(0, 𝑥) = 𝐶0(𝑥), 𝐵(0, 𝑥) = 𝐵0(𝑥) e 𝑧(0, 𝑥) = 𝑧0. As condições de

fronteira são dadas por:

42

𝜕

𝜕𝑡𝐶(𝑡, 0) =

𝜕

𝜕𝑡𝐶(𝑡, 1) = 0 (7.33)

𝜕

𝜕𝑡𝐵(𝑡, 0) =

𝜕

𝜕𝑡𝐵(𝑡, 1) = 0 (7.34)

𝜕

𝜕𝑡𝑧(𝑡, 0) =

𝜕

𝜕𝑡𝑧(𝑡, 1) = 0 (7.35)

Os estados estacionários não triviais no caso adimensional 𝐶𝑠𝑠 e 𝐵𝑠𝑠 são os mesmos para o

caso de uma dimensão, sendo considerados funções constantes 𝐶𝑠𝑠(𝑥) = 𝐶𝑠𝑠 e 𝐵𝑠𝑠(𝑥) = 𝐵𝑠𝑠 em Ω.

Figura 19. Distribuição inicial da população de osteoclastos (extraído de [3]).

Para a simulação do modelo apresentado, fora escolhido uma distribuição inicial de

osteoclastos 𝐶(0, 𝑥) dado pela figura 19. A distribuição inicial de osteoblastos assumiu-se constante e

dada por 𝐵0(0, 𝑥) = 𝐵𝑠𝑠 = 231.7. Os parâmetros das equações foram 𝛼1 = 3 𝑑𝑖𝑎−1, 𝛼2 = 4 𝑑𝑖𝑎−1, 𝛽1 =

0.2 𝑑𝑖𝑎−1, 𝛽2 = 0.02 𝑑𝑖𝑎−1, 𝑔11 = 1.1, 𝑔12 = 1, 𝑔21 = −0.5, 𝑔22 = 0, 𝜎1 = 0.000001, 𝜎2 = 0.000001.

Para as variáveis da equação óssea tem-se 𝑧0(0, 𝑥) = 100%, 𝜎3 = 0.000001, 𝑘1(𝑥) =

0.45 % 𝑐é𝑙𝑢𝑙𝑎𝑠−1𝑑𝑖𝑎−1e 𝑘2(𝑥) = 0.0048 % 𝑐é𝑙𝑢𝑙𝑎𝑠−1𝑑𝑖𝑎−1. Tais valores também foram extraídos de

[3].

Figura 20. Comportamento dinâmico para o caso com uma dimensão espacial.

43

Pode-se notar que o sistema apresenta oscilações regulares no espaço e no tempo,

característico do remodelamento do osso, e aos poucos converge para o seu estado estacionário.


Para o modelo com uma dimensão assumiu-se que as células do tumor também estão

difundindo em Ω. As variáveis 𝐶(𝑡, 𝑥) e 𝐵(𝑡, 𝑥) representam as mesmas populações do modelo

anterior e 𝑇(𝑡, 𝑥) é a densidade de células do tumor para um determinado momento t e espaço x. As

equações para o modelo são dadas por:

𝜕𝐶(𝑡, 𝑥)

𝜕𝑡= 𝜎1

𝜕2

𝜕𝑡2𝐶(𝑡, 𝑥) + 𝛼1𝐶(𝑡, 𝑥)

𝑔11(1+𝑟11𝑇(𝑡,𝑥)𝐿𝑇

)𝐵(𝑡, 𝑥)


)− 𝛽1𝐶(𝑡, 𝑥) (7.36)

𝜕𝐵(𝑡, 𝑥)

𝜕𝑡= 𝜎1

𝜕2

𝜕𝑡2𝐵(𝑡, 𝑥) + 𝛼2𝐶(𝑡, 𝑥)

𝑔12 (1+𝑟12𝑇(𝑡,𝑥)𝐿𝑇

)⁄𝐵(𝑡, 𝑥)

(𝑔22−𝑟22𝑇(𝑡,𝑥)𝐿𝑇

)− 𝛽2𝐵(𝑡, 𝑥) (7.37)

𝜕𝑇(𝑡, 𝑥)

𝜕𝑡= 𝜎4

𝜕2

𝜕𝑡2𝑇(𝑡, 𝑥) + 𝑟𝑇𝑇(𝑡, 𝑥)

𝑎 (1 −𝑇(𝑡, 𝑥)

𝐿𝑇)

𝑏

(7.38)

𝜕𝑧(𝑡, 𝑥)

𝜕𝑡= 𝜎3

𝜕2

𝜕𝑡2𝑧(𝑡, 𝑥) − 𝑘1max{0, 𝐶(𝑡, 𝑥) − 𝐶𝑠𝑠} + 𝑘2max{0, 𝐵(𝑡, 𝑥) − 𝐵𝑠𝑠} (7.39)

com as seguintes condições de fronteira:

𝜕

𝜕𝑡𝐶(𝑡, 0) =

𝜕

𝜕𝑡𝐶(𝑡, 1) = 0 (7.40)

𝜕

𝜕𝑡𝐵(𝑡, 0) =

𝜕

𝜕𝑡𝐵(𝑡, 1) = 0 (7.41)

𝜕

𝜕𝑡𝑇(𝑡, 0) =

𝜕

𝜕𝑡𝑇(𝑡, 1) = 0 (7.42)

𝜕

𝜕𝑡𝑧(𝑡, 0) =

𝜕

𝜕𝑡𝑧(𝑡, 1) = 0 (7.43)

e condições iniciais dadas por 𝐶(0, 𝑥) = 𝐶0(𝑥), 𝐵(0, 𝑥) = 𝐵0(𝑥), 𝑧(0, 𝑥) = 𝑧0(𝑥) e 𝑇(0, 𝑥) = 𝑇0(𝑥). Vale

ressaltar que tais condições de fronteira assumem que não ocorra mudança nos valores das variáveis

nos extremos, no entanto tal facto não condiz com a realidade, pois o tumor pode propagar-se por

todo o osso.

Para a simulação os seguintes parâmetros foram utilizados: 𝛼1 = 3 𝑑𝑖𝑎−1, 𝛼2 = 4 𝑑𝑖𝑎−1, 𝛽1 =

0.2 𝑑𝑖𝑎−1, 𝛽2 = 0.02 𝑑𝑖𝑎−1, 𝑔11 = 1.1, 𝑔12 = 1, 𝑔21 = −0.5, 𝑔22 = 0, 𝑎 = 2, 𝑏 = 3, 𝑟11 = 0.005, 𝑟12 = 0.0,

𝑟21 = 0.0, 𝑟22 = 0.2, 𝑟𝑇 = 0.004, 𝐿𝑇 = 100, 𝜎1 = 0.000001, 𝜎2 = 0.000001, 𝜎3 = 0.000001, 𝜎4 =

0.000001, 𝑘1 = 0.45% 𝑐é𝑙𝑢𝑙𝑎𝑠−1𝑑𝑖𝑎−1 e 𝑘2 = 0.0048% 𝑐é𝑙𝑢𝑙𝑎𝑠

−1𝑑𝑖𝑎−1.

A distribuição inicial das populações de osteoclastos e osteoblastos são os mesmos do tópico

7.5 e o tumor é inicialmente pequeno e localizado no lado direito de Ω = [0,1].

Ao longo do tempo o tumor propaga-se para o resto do osso e, com isso, afeta a dinâmica

das células presentes. Ocorre um aumento na produção de osteoclastos e uma diminuição na de

osteoblastos, onde ambas buscam atingir o novo estado estacionário. Como resultado, a ação de

destruição do osso supera a de reconstrução e a densidade óssea diminui até atingir o valor nulo.

Apesar de não apresentado na simulação, espera-se que o tumor continue propagando-se até a

eliminação total do osso.

44

Figura 21. Comportamento dinâmico na presença de tumor para uma dimensão espacial.


O modelo com o tratamento é obtido através da adição das funções dependente do tempo

𝑉1(𝑡) e 𝑉2(𝑡) como no caso adimensional. As equações ficam:

𝜕𝐶(𝑡, 𝑥)

𝜕𝑡= 𝜎1

𝜕2

𝜕𝑡2𝐶(𝑡, 𝑥) + 𝛼1𝐶(𝑡, 𝑥)


)𝐵(𝑡, 𝑥)


)− 𝛽1𝐶(𝑡, 𝑥) (7.44)

𝜕𝐵(𝑡, 𝑥)

𝜕𝑡= 𝜎2

𝜕2

𝜕𝑡2𝐵(𝑡, 𝑥) + 𝛼2𝐶(𝑡, 𝑥)

𝑔12 (1+𝑟12𝑇(𝑡,𝑥)𝐿𝑇

)⁄𝐵(𝑡, 𝑥)


)− (𝛽2 − 𝑉1(𝑡))𝐵(𝑡, 𝑥) (7.45)

𝜕𝑇(𝑡, 𝑥)

𝜕𝑡= 𝜎4

𝜕2

𝜕𝑡2𝑇(𝑡, 𝑥)(𝑟𝑇 − 𝑉2(𝑡))𝑇(𝑡)

𝑎 (1 −𝑇(𝑡)

𝐿𝑇)

𝑏

(7.46)

𝜕𝑧(𝑡, 𝑥)

𝜕𝑡= 𝜎3

𝜕2


com as seguintes condições de fronteira:

𝜕

𝜕𝑡𝐶(𝑡, 0) =

𝜕

𝜕𝑡𝐶(𝑡, 1) = 0 (7.48)

𝜕

𝜕𝑡𝐵(𝑡, 0) =

𝜕

𝜕𝑡𝐵(𝑡, 1) = 0 (7.49)

𝜕

𝜕𝑡𝑇(𝑡, 0) =

𝜕

𝜕𝑡𝑇(𝑡, 1) = 0 (7.50)

𝜕

𝜕𝑡𝑧(𝑡, 0) =

𝜕

𝜕𝑡𝑧(𝑡, 1) = 0 (7.51)

e condições iniciais dadas por 𝐶(0, 𝑥) = 𝐶0(𝑥), 𝐵(0, 𝑥) = 𝐵0(𝑥), 𝑧(0, 𝑥) = 100% e 𝑇(0, 𝑥) = 𝑇0(𝑥).

45

A figura 22 mostra o resultado da simulação para esse novo caso. Todos os parâmetros e

condições iniciais são os mesmos do caso com tumor e sem tratamento e 𝑣1 = 0.0001, 𝑣2 = 0.006. O

início do tratamento é dado por 𝑡𝑠𝑡𝑎𝑟𝑡 = 600 dias.

Figura 22. Comportamento dinâmico na presença de tumor e tratamento com uma dimensão espacial.

Essa forma de tratamento corresponde a fármacos como os inibidores de proteassoma, os

quais promovem a produção de osteoblastos e inibem o crescimento do tumor.

Os gráficos mostram as populações de osteoclastos e osteoblastos recuperando o seu ciclo

regular após a eliminação da doença. O tumor acaba sendo quase totalmente extinguido pelo

tratamento, sobrando apenas alguns pequenos pontos residuais, antes que todo o osso seja

destruído (lado esquerdo é recuperado, mas não o lado direito). Tais resultados são coerentes com

os estudos apresentados em [13,16].

No entanto, vale ressaltar que não houve recuperação da densidade óssea perdida,

provavelmente como reflexo dos parâmetros escolhidos para a simulação. Espera-se que tal

recuperação ocorra ao longo de um grande período. Além disso, aparentemente o tempo para os

inibidores de proteassoma recuperarem a densidade óssea perdida é muito maior do que o tempo

necessário para a eliminação do tumor, o que pode ter implicações no período de tratamento.


O sistema dinâmico para o modelo PK/PD também pode ser estendido para uma dimensão

espacial. Nesse caso, considera-se que o efeito farmacodinâmico, 𝑑(𝑡), de um determinado fármaco

será o mesmo para todo o comprimento do osso em um determinado instante. Esse novo modelo é

descrito como:

𝜕𝐶(𝑡, 𝑥)

𝜕𝑡= 𝜎1

𝜕2

𝜕𝑡2𝐶(𝑡, 𝑥) + 𝛼1𝐶(𝑡, 𝑥)


)𝐵(𝑡, 𝑥)


)(1+𝐾𝑠1𝑑1(𝑡)) − (1 + 𝐾𝑠2𝑑2(𝑡))𝛽1𝐶(𝑡, 𝑥) (7.52)

𝜕𝐵(𝑡, 𝑥)

𝜕𝑡= 𝜎2

𝜕2

𝜕𝑡2𝐵(𝑡, 𝑥) + 𝛼2𝐶(𝑡, 𝑥)

𝑔12 (1+𝑟12𝑇(𝑡,𝑥)𝐿𝑇

)⁄𝐵(𝑡, 𝑥)


)− 𝛽2𝐵(𝑡, 𝑥) (7.53)

46

𝜕𝑇(𝑡)

𝜕𝑡= 𝜎4

𝜕2

𝜕𝑡2𝑇(𝑡, 𝑥) + (1 − 𝐾𝑖34𝑑𝑐34(𝑡))𝑟𝑇𝑇(𝑡, 𝑥)

𝑎 (1 −𝑇(𝑡, 𝑥)

𝐿𝑇)

𝑏

(7.54)

𝜕𝑧(𝑡)

𝜕𝑡= 𝜎3

𝜕2


Para a simulação desse caso foram utilizados os mesmos parâmetros e condições inicias do

tópico 7.5. O parâmetro de cada fármaco é dado pela tabela II.

Tabela II. Parâmetros dos fármacos para o modelo PK/PD 1D.

Combinação 𝑑34


𝜏 1 1 1 1

𝑘𝑑 1 1 1 0.5

𝑘𝑎 0.01 0.01 0.01 0.01

𝐶0 0.008 0.008 0.008 0

𝐶0𝑏𝑎𝑠𝑒 0.002 0.002 0.002 0.002

𝐾𝑟 0 0 0 0

𝐿𝑟 0 0 0 0

𝐾𝑠,𝑖 2.3 0.1 0.1

𝐶𝐼 1 - -

Figura 23. Modelo de tratamento com o uso de PK/PD com uma dimensão espacial.

Pode-se notar que, de modo semelhante ao tratamento com inibidores de proteassoma,

houve a redução do tumor a tempo que ele não destruísse todo o osso, sendo que apenas uma parte

fora perdida. Após o início do tratamento, ocorre um pequeno atraso até que o tumor comece a ser

reduzido. Tal como no modelo sem dimensão espacial, esse atraso é devido a farmacocinética e a

farmacodinâmica que reproduzem um efeito inibidor pequeno enquanto as concentrações dos

fármacos no corpo do paciente ainda são pequenas.

O tratamento acaba por afetar também as células que não tinham sido atingidas ainda pela

doença, de modo que passam a convergir para um novo estado estacionário. Desse modo, pode-se

notar uma pequena queda em um dos extremos do osso. Tal acontecimento sugere que a dosagem

47

dos fármacos pode ter sido um pouco elevada e acabou reduzindo por demais a quantidade de

osteoclastos e osteoblastos. No entanto, de um modo geral, pode-se dizer que o tratamento foi

efetivo e coerente com os resultados apresentados em [2,12,14,18,41].

8. Aplicação de derivadas fracionárias

Nesse capítulo será apresentado o estudo das derivadas fracionárias nos modelos dinâmicos do

osso presentes no capítulo 7. Tal processo será feito através da substituição dos termos derivativos

em relação ao tempo por um de ordem fracionária, tanto para o modelo de dimensão espacial nula,

quanto para o de uma dimensão espacial.

Os valores das variáveis das equações presentes serão escolhidos com bases nos trabalhos de

[19], [3] e [9] para as simulações, conforme feito anteriormente. Com isso, busca-se um melhor grau

de comparabilidade entre os diversos modelos de ordem fracionária e o de ordem inteira. Apenas em

alguns casos alguma variável precisará ter o seu valor alterado para que o modelo seja, do ponto de

vista biológico, relevante. Além disso, um gráfico comparativo entre diversas ordens de derivada para

um mesmo conjunto de parâmetros será feito em cada caso, com o objetivo de ilustrar melhor as

diferenças de comportamento do sistema. Recorda-se que muito dos valores de parâmetros utilizados

são estimações feitas de modo que o sistema apresente um comportamento qualitativo próximo ao

que se conhece sobre a dinâmica.

A solução das equações diferenciais de ordem fracionária será feita através do método numérico

das matrizes triangulares em combinação com o método do passo largo (“Large Step”). Esses

métodos estão presentes em [34,36,37] e explicados no capítulo 4. É importante salientar que tais

modelos dinâmicos não apresentam condições iniciais nulas e, portanto, uma mudança de variável

será necessária para a aplicação do método das matrizes triangulares.

O modelo PK/PD será feito em Simulink com o auxílio da toolbox FOMCON [42] e os detalhes de

sua implementação são dados em anexo.


O modelo sem tumor e sem tratamento, apresentado em [19], é dado por (7.1), (7.2) e (7.3). A

aplicação das derivadas fracionárias consiste na substituição do termo derivativo de (7.1), (7.2) e

(7.3) pela definição de derivada de Caputo, dado por:

𝐷𝑐 𝑡𝛼𝑓(𝑡) =

{

∫

(𝑡 − 𝜏)−𝛼−1

Γ(−𝛼)𝑓(𝜏)𝑑𝜏, 𝑠𝑒 𝛼 𝜖 ℝ−

𝑡

𝑐

𝑓(𝑡), 𝑠𝑒 𝛼 = 0

𝐷𝑐 𝑡𝛼−⌈𝛼⌉ 𝑑

⌈𝛼⌉

𝑑𝑡⌈𝛼⌉𝑓(𝑡), 𝑠𝑒 𝛼 𝜖 ℝ+

(8.1)

onde o termo 𝛼 ∈ ℝ é a ordem da derivada. Caso 𝛼 = 1, tem-se a definição tradicional de derivada.

Com essa substituição, o sistema de equações terá o seguinte formato:

𝐷0 𝑡𝛼𝐶(𝑡) = 𝛼1𝐶(𝑡)

𝑔11𝐵(𝑡)𝑔21 − 𝛽1𝐶(𝑡) (8.2)

𝐷0 𝑡𝛼𝐵(𝑡) = 𝛼2𝐶(𝑡)

𝑔12𝐵(𝑡)𝑔22 − 𝛽2𝐵(𝑡) (8.3)

𝐷0 𝑡𝛼𝑧(𝑡) = −𝑘1max{0, 𝐶(𝑡) − 𝐶𝑠𝑠} + 𝑘2max{0, 𝐵(𝑡) − 𝐵𝑠𝑠} (8.4)

48

Sabe-se que os pontos de equilíbrio não triviais para o sistema (8.2), (8.3) e (8.4) serão os

mesmos que o do sistema (7.1), (7.2) e (7.3) para todo 𝛼 ∈ (0,1) [22]. Os pontos de equilíbrio serão:


(1−𝑔22)Ω

(𝛽2𝛼2)

𝑔21Ω

(8.5)


𝑔12Ω

(𝛽2𝛼2)

(1−𝑔11)Ω

(8.6)

onde Ω = 𝑔12𝑔21 − (1 − 𝑔11)(1 − 𝑔22). Para o caso de 𝛼 maior que 1 tal condição não é

necessariamente verdadeira e, portanto, os pontos serão determinados com base nos resultados das

simulações.

A solução desse sistema não-linear de ordem fracionária é encontrada através do método

numérico de [34]. Como as equações (8.2), (8.3) e (8.4) não possuem condições iniciais nulas, é

necessário realizar uma mudança de variáveis de modo que o novo sistema apresente tais condições

iniciais nulas. Em anexo é descrito o modo de se fazer isso.

A tabela III mostra os valores dos parâmetros, bem como a ordem da derivada utilizada em

cada simulação. Tais parâmetros são os mesmos utilizados em [19] de modo a facilitar a comparação

dos modelos. As figuras 24-29 mostram o comportamento da população de osteoclastos e

osteoblastos, bem como a variação na densidade óssea.

Como as populações de osteoblastos e osteoclastos passam a demorar tempos distintos para

atingir o seu estado estacionário, permanecendo acima ou abaixo desse valor por mais tempo, para

cada valor de 𝛼, o seu efeito de atividade na reabsorção e reconstrução do osso, 𝑘1e 𝑘2, devem ser

ajustados para que a densidade óssea atinja o seu valor inicial.

Tabela III. Valores dos parâmetros utilizados para cada simulação.

𝛼 𝐶0 𝐵0 𝛼1 𝛼2 𝛽1 𝛽2 𝑔11 𝑔12 𝑔21 𝑔22 𝑘1 𝑘2

1.3 5.0 210 3.0 4.0 0.2 0.02 0.3 1.0 -0.5 0.0 0.24 0.00249

1.2 6.0 211 3.0 4.0 0.2 0.02 0.4 1.0 -0.5 0.0 0.24 0.00215

1.1 11 211 3.0 4.0 0.2 0.02 0.4 1.0 -0.5 0.0 0.24 0.00173

0.9 11.06 212.1 3.0 4.0 0.2 0.02 0.5 1.0 -0.5 0.0 0.24 0.0012

0.8 11.06 212.1 3.0 4.0 0.2 0.02 0.5 1.0 -0.5 0.0 0.24 0.0012

0.7 11.05 210.0 3.0 4.0 0.2 0.02 0.3 1.0 -0.5 0.0 0.24 0.0012

A figura 24 mostra o comportamento do sistema para 𝛼 = 1.3. Nota-se uma amplitude de

oscilação maior nas populações se comparado ao modelo de 𝛼 = 1.0, além de um tempo menor para

o sistema atingir o seu estado estacionário, característico de um sistema com menor amortecimento.

Nesse caso a população de osteoclastos atinge o seu estado estacionário em aproximadamente 60

dias, enquanto que os osteoblastos e a densidade do osso demoram cerca de 190 dias. A redução no

efeito autócrino dos osteoclastos, 𝑔11,bem como o valor do distúrbio da população inicial de

osteoclastos, 𝐶0, foram necessários uma vez que valores mais elevados provocam como resposta um

aumento muito grande da população de osteoblastos. Esse aumento leva a uma queda muito brusca

na população de osteoclastos e faz com que esta oscile e atinja valores negativos, o que é, em

termos biológicos, errado.

49

A simulação apresentada pela figura 25 corresponde a 𝛼 = 1.2. O sistema ainda apresenta

oscilações da população de osteoclastos, porém já é possível aumentar um pouco o valor de 𝑔11 e 𝐶0

em comparação ao modelo de 𝛼 = 1.3 sem que a população de osteoclastos atinja valores negativos.

Nesse caso a população de osteoclastos demora cerca de 90 dias para atingir o estado estacionário,

enquanto que a de osteoblastos leva 190 dias e apresenta um pico populacional de aproximadamente

290 unidades. Tal aumento em relação a 𝛼 = 1.3 deve-se aos novos valores de 𝑔11 e 𝐶0, que são

resultam em uma maior quantidade de osteoclastos e, como resposta, um maior aumento na

população de osteoblastos.

Figura 24. Dinâmica das populações para α=1.3.


O modelo de 𝛼 = 1.1 é ilustrado pela figura 26. Ele possui um nível de amortecimento maior

do que o modelo apresentado anteriormente, de modo que é possível aumentar o valor inicial do

distúrbio da população de osteoclastos, 𝐶0, sem afetar a coerência do sistema.


50

No caso da figura 27 tem-se o resultado da simulação para 𝛼 = 0.9. Pode-se notar um maior

amortecimento do sistema em relação aos modelos anteriores, onde a população de osteoclastos

sofre uma pequena oscilação abaixo do seu valor estacionário e estabiliza-se por volta de 130 dias. A

população de osteoblastos apresenta um pico de aproximadamente 390 unidades e demora cerca de

550 dias para tingir o seu estado estacionário, devido a característica um pouco mais lenta do

sistema. Já a densidade óssea recupera o seu valor inicial em aproximadamente 510 dias, uma vez

que a população de osteoblastos acima do estado estacionário após esse período é extremamente

pequena e acaba por não exercer um efeito construtivo significativo no osso.


Na figura 28 o modelo de 𝛼 = 0.8 é apresentado. Diferentemente dos modelos anteriores,

agora a população de osteoclastos não apresenta mais uma oscilação em relação ao seu valor

estacionário, mas sim uma lenta convergência de cerca de 900 dias até atingi-lo. A população de

osteoblastos mantém um comportamento semelhante ao modelo anterior, mas apresenta um tempo

muito maior até atingir o seu ponto de equilíbrio, que acontece por volta de 1750 dias.


No caso de 𝛼 = 0.7, ilustrado na figura 29, o sistema passa a ficar muito amortecido. Em

razão disso, o valor autócrino dos osteoclastos precisou ser reduzido, de modo que o sistema não

demorasse um tempo muito longo para atingir o estado estacionário e ficasse biologicamente

irrelevante. Com essa redução a capacidade dos osteoclastos de estimularem a produção de mais

osteoclastos é afetada, o que faz com que a população atinja o estado estacionário mais

rapidamente.

51

Essa redução do parâmetro acaba por também afetar a população de osteoblastos. Como

ambas as populações estão diretamente relacionadas entre si, o menor número de osteoclastos exige

uma quantidade menor de osteoblastos para controlar a reabsorção do osso. Tal efeito fica evidente

ao observar o pico da população de osteoblastos, que é menor em relação aos modelos anteriores.


A figura 30 ilustra uma simulação comparativa entre diversos valores dos expoentes das

derivadas para a população de osteoclastos e osteoblastos. Com isso, fica mais fácil a visualização

da mudança de comportamento do sistema. Os parâmetros utilizados para essa simulação foram

𝛼1 = 3, 𝛼2 = 4, 𝛽1 = 0.2, 𝛽2 = 0.02, 𝑔11 = 0.5, 𝑔12 = 1, 𝑔21 = −0.5, 𝑔22 = 0, 𝐶(0) = 4 e 𝐵(0) = 212.

Figura 30. Comparação do comportamento de diversas ordens fracionárias para um mesmo parâmetro de

equações

Pode-se notar que, de uma forma geral, o comportamento qualitativo do sistema não mudou.

O aumento inicial na população de osteoclastos provoca um consequente aumento da população de

osteoclastos. Esse aumento faz com que o número de osteoclastos se reduza e o sistema convirja

para o estado estacionário.

No entanto, cada vez que o valor α da derivada aproxima-se de zero, o sistema torna-se mais

amortecido (comportamento já esperado conforme [31,33,43]) e demora mais para atingir o estado

estacionário. Quando o expoente da derivada passa a ser maior que um ocorre o oposto, o sistema

começa a convergir mais rapidamente para o estado estacionário, no entanto acaba apresentando

maiores oscilações

Por fim, para ilustrar melhor o efeito do amortecimento, a figura 31 apresenta o resultado de

uma simulação para 𝛼 = 0.8 com os mesmos parâmetros utilizados na simulação da figura 15. Neste

caso o sistema, ao invés de entrar em um estado oscilatório constante conforme o modelo com 𝛼 =

52

1.0, acaba por convergir para o estado estacionário 𝐶𝑠𝑠 = 1.16 e 𝐵𝑠𝑠 = 231.72, após uma pequena

oscilação por parte da população de osteoclastos.

Figura 31. Sistema não apresenta comportamento cíclico quando α=0.8.


O modelo com o tumor e sem tratamento é dado por (7.11), (7.12), (7.13) e (7.14). A

aplicação das derivadas fracionárias segue a mesma metodologia apresentada anteriormente, com a

substituição do termo derivativo pela definição de Caputo. Com isso, o sistema anterior apresenta o

seguinte formato com as derivadas fracionárias:



)𝐵(𝑡)


)− 𝛽1𝐶(𝑡)

(8.7)


𝑔12 (1+𝑟12𝑇(𝑡)𝐿𝑇

)⁄𝐵(𝑡)


)− 𝛽2𝐵(𝑡)

(8.8)

𝐷0 𝑡𝛼𝑇(𝑡) = 𝑟𝑇𝑇(𝑡)

𝑎 (1 −𝑇(𝑡)

𝐿𝑇)

𝑏

(8.9)


As condições iniciais para esse caso são 𝐶(0) = 𝐶0, 𝐵(0) = 𝐵0, 𝑇(0) = 𝑇0 e 𝑧(0) = 𝑧0. Para

𝛼 ∈ (0,1) o sistema terá os pontos de equilíbrio expressos por:


(1−𝑔22−𝑟22)Λ

(𝛽2𝛼2)

𝑔21(1+𝑟21)Λ

(8.11)


(𝑔12 (1+𝑟12)⁄ )Λ

(𝛽2𝛼2)

1−𝑔11(1+𝑟11)Λ

(8.12)

𝑇𝑠𝑠 = 𝐿𝑇 (8.13)

onde Λ = (𝑔12 (1 + 𝑟12)⁄ )(𝑔21(1 + 𝑔21)) − (1 − 𝑔11(1 + 𝑟11))(1 − 𝑔22 + 𝑟22). Em anexo é apresentado a

mudança de variável necessária para a aplicação do método das matrizes triangulares.

As figuras 32-36 mostram o resultado das simulações para os parâmetros e condições inicias

apresentados na tabela IV e V. Tais valores foram baseados em [3] com exceção dos parâmetros da

equação dinâmica do tumor, que foram os de [9], de modo a facilitar a comparação entre modelos.

Alguns desses parâmetros precisaram ser alterados para determinadas ordens de derivadas, como 𝑘1

e 𝑘2, que estão diretamente relacionados com as amplitudes das oscilações de cada população, 𝑟𝑇,

responsável pela velocidade com que o tumor se desenvolve e 𝑔11, que possui uma forte influência

na dinâmica do sistema.

53

Tabela IV. Parâmetros das populações utilizados nas simulações de ordem fracionária.

𝛼 𝐶0 𝐵0 𝛼1 𝛼2 𝛽1 𝛽2 𝑔11 𝑔12 𝑔21 𝑔22 𝑘1 𝑘2

1.3 8.0 302.0 3.0 4.0 0.2 0.02 1.08 1.0 -0.5 0.0 0.0748 0.0003395

1.2 10.0 302.0 3.0 4.0 0.2 0.02 1.08 1.0 -0.5 0.0 0.0748 0.0003395

1.1 15.0 316.0 3.0 4.0 0.2 0.02 1.1 1.0 -0.5 0.0 0.0748 0.0003395

0.9 17.8 460.6 3.0 4.0 0.2 0.02 1.163 1.0 -0.5 0.0 0.1248 0.0005095

0.8 22.1 661.1 3.0 4.0 0.2 0.02 1.21 1.0 -0.5 0.0 0.1248 0.0005095

0.7 30.6 1028.5 3.0 4.0 0.2 0.02 1.244 1.0 -0.5 0.0 0.1248 0.0005095

Tabela V. Parâmetros do tumor utilizados nas simulações de ordem fracionária.

𝛼 𝑇0 a b 𝑟𝑇 𝐿𝑇 𝑟11 𝑟12 𝑟21 𝑟22

1.3 1.0 2.0 3.1 0.0025 100 0.005 0.0 0.0 0.2

1.2 1.0 2.0 3.1 0.0033 100 0.005 0.0 0.0 0.2

1.1 1.0 2.0 3.1 0.005 100 0.005 0.0 0.0 0.2

0.9 1.0 2.0 3.1 0.01 100 0.005 0.0 0.0 0.2

0.8 1.0 2.0 3.1 0.02 100 0.005 0.0 0.0 0.2

0.7 1.0 2.0 3.1 0.026 100 0.005 0.0 0.0 0.2

O primeiro modelo analisado, correspondente a 𝛼 = 1.3, é apresentado pela figura 32. Pode-

se notar que, mesmo com um distúrbio inicial pequeno na população de osteoclastos e um menor

valor no efeito autócrino do mesmo, o sistema não apresenta convergência para o seu estado

estacionário, mas sim oscilações cada vez maiores que fazem o modelo divergir. O tumor apresenta

um desenvolvimento mais acelerado, mesmo com o parâmetro 𝑟𝑇 menor e leva aproximadamente

800 dias para atingir o seu valor máximo. Além disso, sua presença colabora para o aumento da

população de osteoclastos e inibição da de osteoblastos e, como consequência, uma maior

reabsorção do osso e queda na densidade óssea.

Figura 32. Comportamento das populações para α=1.3 na presença do tumor.

A simulação com 𝛼 = 1.2 é ilustrado através da figura 33. Agora o modelo apresenta uma

convergência para o seu estado estacionário, mesmo com um distúrbio em 𝐶0 maior. O tempo de

desenvolvimento completo do tumor fica menor, se comparado a 𝛼 = 1.3, uma vez que também leva

aproximadamente 800 dias porém com um 𝑟𝑇 maior.

54

O caso de 𝛼 = 1.1 é apresentado pela figura 34. Novamente os parâmetros 𝑔11 e 𝐶0 puderem

ser aumentados, se assemelhando mais ao modelo de 𝛼 = 1.0, sem que o sistema adquirisse um

comportamento divergente. O tempo de convergência ainda continua menor, característico de

sistemas com ordem fracionária maiores que um.



A figura 35 mostra o modelo de 𝛼 = 0.9. Nesse caso o nível de amortecimento do sistema

passa a ser bem mais elevado. Com isso, é necessário aumentar o valor de 𝑔11 de modo

considerável, bem como o distúrbio inicial 𝐶0, para estimular uma maior oscilação do sistema e a

relação entre as populações fique mais evidente.

Devido ao maior amortecimento do sistema, as amplitudes das oscilações ficam menores e a

convergência mais lenta.

A simulação para 𝛼 = 0.8, representado pela figura 36, apresenta um comportamento

semelhante ao de 𝛼 = 0.9, onde o grau de amortecimento elevado exige um aumento maior dos

parâmetros 𝑔11, 𝐶0 e 𝑟𝑇 para induzir oscilações no sistema. O tumor apresenta um desenvolvimento

muito mais lento, no qual leva mais de 3000 dias para atingir o seu tamanho máximo.

O último modelo é ilustrado pela figura 37. Para esse caso, onde 𝛼 = 0.7, o nível de

amortecimento é extremamente elevado e resulta em uma menor quantidade de oscilações, bem

como uma convergência mais devagar, mesmo com alguns parâmetros aumentados.

55




Na figura 38 é feito um gráfico comparativo entre os diversos comportamentos do sistema

para certas ordens de derivadas fracionárias.

Como ordens de derivada mais elevadas (maiores que um) divergem mais facilmente, mesmo

para valores de parâmetros menores, optou-se por subdividir a comparação em dois blocos. Sendo

que no caso onde 0.5 ≤ 𝛼 ≤ 0.95 os parâmetros utilizados foram 𝛼1 = 3.0, 𝛼2 = 4.0, 𝛽1 = 0.2, 𝛽2 =

0.02, 𝑔11 = 1.1, 𝑔12 = 1.0, 𝑔21 = −0.5, 𝑔22 = 0, 𝑟11 = 0.005, 𝑟12 = 0, 𝑟21 = 0, 𝑟22 = 0.2, 𝑎 = 2.0, 𝑏 = 3.1,

𝐿𝑇 = 100, 𝑟𝑇 = 0.006, 𝐶(0) = 15, 𝐵(0) = 360 e 𝑇(0) = 1.0. No caso 1.0 ≤ 𝛼 ≤ 1.3, 𝑟𝑇 = 0.002.

Novamente fica evidente a diferença de comportamento para as diversas ordens fracionárias.

Enquanto que em alguns casos o sistema não oscila, em outros é possível ver pequenas oscilações

que convergem para o estado estacionário e até mesmo divergência, como em 𝛼 ∈ (1.2,1.3).

56

Figura 38. Gráficos comparativos entre diversas ordens fracionárias. Os três primeiros correspondem a α≤0.95 e

os três últimos para 1≤α≤1.3.

Pode-se notar que a velocidade de crescimento do tumor é fortemente afetada pela ordem de

sua derivada. Quanto menor o seu valor, mais tempo o tumor demora para atingir o seu estado

estacionário 𝑇𝑠𝑠 = 𝐿𝑇. Por isso, o parâmetro 𝑟𝑇, responsável pela taxa de crescimento do tumor, teve

de ser aumentado conforme a ordem da derivada diminuía de modo a garantir que o valor 𝐿𝑇 fosse

atingido em um tempo aceitável em termos biológicos. Vale ressaltar que essa mudança no valor do

parâmetro não trouxe profundas modificações na característica geral da curva de crescimento do

tumor.

As populações apresentam o mesmo comportamento qualitativo que no modelo original. A

presença do tumor provoca um maior estimulo da produção de osteoclastos através do aumento de

RANKL e redução de OPG, de modo que essa passe a convergir para um novo estado estacionário.

Já a população de osteoblastos é afetada negativamente pela presença do tumor, devido a sua ação

inibidora, uma vez que a diferenciação e maturação dessas células não é mais proporcional a dos

osteoclastos.

Essa alteração no equilíbrio entre as duas populações faz com que a atividade de destruição

do osso causada pelos osteoclastos seja superior ao de formação do osso, por parte dos

osteoblastos. Como efeito, a densidade óssea acaba por reduzir-se ao longo do tempo, tendendo ao

57

valor nulo. Naturalmente, o paciente acaba falecendo muito antes disso devido a diversas

complicações do organismo.

A frequência das oscilações também diminui à medida que a ordem da derivada diminui e o

motivo é a característica de maior amortecimento para sistemas de menor grau. Uma das maneiras

de modificar essa frequência é através de 𝑔11 O aumento no valor do parâmetro serve para provocar

as oscilações no sistema, uma vez que valores baixos fazem com que o sistema não oscile e a

relação entre as duas populações de células não fique muito evidente.


O modelo de equações diferenciais com a presença do tumor, expresso pela equação

analítica de [9], e do tratamento de proteassoma [3], é dado por (7.19), (7.20), (7.21) e (7.22). A

aplicação de derivadas fracionárias para o modelo, através da definição de Caputo, resulta em:



)𝐵(𝑡)


)− 𝛽1𝐶(𝑡)

(8.14)


𝑔12 (1+𝑟12𝑇(𝑡)𝐿𝑇

)⁄𝐵(𝑡)


)− (𝛽2 − 𝑉1(𝑡))𝐵(𝑡)

(8.15)

𝐷0 𝑡𝛼𝑇(𝑡) = (𝑟𝑇 − 𝑉2(𝑡))𝑇(𝑡)

𝑎 (1 −𝑇(𝑡)

𝐿𝑇)

𝑏

(8.16)


O procedimento para tornar as equações acima com condições iniciais nulas é o mesmo

mostrado no tópico anterior, de modo que será omitido dessa vez.

A tabela VI mostra os parâmetros do tumor e do tratamento utilizados para a simulação de

cada ordem de derivada enquanto que os outros são fornecidos pela tabela IV. Os resultados das

simulações são dados pelas figuras 38-41.

Tabela VI. Parâmetros do tumor e do tratamento com inibidores de proteassoma utilizados na simulação.

𝛼 𝑇0 a b 𝑟𝑇 𝐿𝑇 𝑟11 𝑟12 𝑟21 𝑟22 𝑣1 𝑣2 𝑡𝑠𝑡𝑎𝑟𝑡

1.2 1.0 2.0 3.1 0.0033 100 0.005 0.0 0.0 0.2 0.0015 0.006 400

1.1 1.0 2.0 3.1 0.005 100 0.005 0.0 0.0 0.2 0.001 0.008 500

0.9 1.0 2.0 3.1 0.01 100 0.005 0.0 0.0 0.2 0.0009 0.017 900

0.8 1.0 2.0 3.1 0.02 100 0.005 0.0 0.0 0.2 0.00005 0.024 1500

O primeiro modelo apresentado, correspondendo a 𝛼 = 1.2, é ilustrado pela figura 39. Nele o

tumor se desenvolve livremente durante um período de 400 dias, momento no qual o tratamento com

inibidores de proteassoma tem início. A terapia consegue inibir o desenvolvimento do tumor de modo

que apenas permanece uma quantidade residual após 1500 dias. A dinâmica do osso é recuperada e

a densidade óssea começa a apresentar sinais de recuperação após 1000 dias de simulação.

Na figura 40 tem-se o resultado da simulação para 𝛼 = 1.1. Nesse caso foi implementado um

aumento no parâmetro 𝑔11 e nas condições inicias, uma vez que não causava instabilidade do

sistema. O modelo apresenta um comportamento bem semelhante ao anterior, apesar de ser um

pouco mais lento. Mesmo com um aumento no efeito 𝑣2, responsável pela redução do tumor, as

células cancerígenas demoraram um pouco mais de tempo para serem eliminadas.

58

Figura 39. Comportamento das populações para α=1.2 na presença do tumor e tratamento.


A figura 41 ilustra o caso de 𝛼 = 0.9. Nesse caso o nível de amortecimento passa a ficar

muito mais elevado e o sistema passa a ter um desenvolvimento mais lento. O tratamento inicia-se

depois de 1000 dias, quando o tumor já está em um estágio mais evoluído, e leva cerca de 1500 dias

para extinguir a doença. O sistema apresenta uma recuperação muito mais lenta do seu ciclo se

comparado aos modelos anteriores.


59

O modelo de 𝛼 = 0.8 é dado pela figura 42. Nota-se uma grande redução do número de

oscilações entre os osteoclastos e osteoblastos. O tratamento consegue reduzir o tumor e evitar a

perca de mais osso, no entanto o sistema acaba por atingir seu novo estado estacionário e a

recuperação óssea não acontece.

Com esses resultados pode-se concluir que os inibidores de proteassoma possuem um

grande efeito na redução do tumor, bem como na prevenção da diminuição da densidade óssea. Os

resultados são consistentes com os estudos clínicos realizados em [13,16,20,27,28] e com os

resultados para o caso de 𝛼 = 1.0 feito no capítulo 7, onde o tumor fora reduzido e a perda óssea

evitada através da redução da apoptose dos osteoblastos.


Os modelos de ordem mais elevada permitiram uma quantidade maior de oscilações e

velocidade de convergência, o que resultou em uma recuperação óssea no período de simulação

apresentado. Nos casos de ordens menores, devido a seu caráter mais lento e amortecido, o número

de oscilações diminuiu e o tratamento acaba por não conseguir restaurar a massa óssea perdida.

Novamente é possível notar o papel decisivo de 𝑔11 para a dinâmica do sistema. Para cada

ordem de derivada é necessário o seu ajuste de modo que o comportamento dinâmico apresente um

caráter qualitativo adequado. Os parâmetros 𝑟𝑇 e 𝑣2, estão relacionados com a velocidade de

crescimento e destruição do tumor e seus valores precisam ser maiores para os casos onde 𝛼 é

menor. Por fim, 𝑣1 diz respeito a redução da apoptose dos osteoblastos e não pode ser muito elevado

pois acaba influenciando a produção de muitos osteoblastos e, consequentemente, a redução

abrupta de osteoclastos. Como resultado, a densidade óssea pode sofrer um aumento muito grande e

levar o paciente a um novo quadro clínico.


A aplicação das derivadas fracionárias inicia-se no modelo farmacocinético apresentado no

capitulo 6. As concentrações de fármacos que ainda faltam ser absorvidos e presentes no plasma

sanguíneo, 𝐶𝑔 e 𝐶𝑝, são escritas no seguinte formato:

𝐷0 𝑡𝛼𝐶𝑔(𝑡) = −𝑘𝑎𝐶𝑔(𝑡) (8.18)

𝐷0 𝑡𝛼𝐶𝑝(𝑡) = 𝑘𝑎𝐶𝑔(𝑡) − 𝑘𝑑𝐶𝑝(𝑡) (8.19)

60

e a transformada de Laplace, para a ingestão de um fármaco oral, de 𝐶𝑝 é dada por:


𝑘𝑎(𝑠𝛼 + 𝑘𝑎)(𝑠

𝛼 + 𝑘𝑑) (8.20)

A formulação da farmacodinâmica é a mesma apresentada no capítulo 7 e a metodologia

para encontrar o modelo PK/PD é dada em Anexo nessa dissertação. A figura 42 ilustra a

concentração de um fármaco 𝑑(𝑡) de uma série de administrações cíclicas para diversas ordens

fracionárias (𝐶0 = 0.008, 𝑘𝑎 = 0.01, 𝑘𝑑 = 1, 𝐶50 = 0.002 e 𝜏 = 1).

Figura 43. Modelo PK/PD de um fármaco para diversas ordens fracionárias.

O modelo dinâmico do osso, baseado em [9], com equações fracionárias é dado por:



)𝐵(𝑡)


)(1+𝐾𝑠1𝑑1(𝑡)) − (1 + 𝐾𝑠2𝑑2(𝑡))𝛽1𝐶(𝑡) (8.21)


𝑔12 (1+𝑟12𝑇(𝑡)𝐿𝑇

)⁄𝐵(𝑡)


)− 𝛽2𝐵(𝑡)

(8.22)

𝐷0 𝑡𝛼𝑇(𝑡) = (1 − 𝐾𝑖34𝑑𝑐34(𝑡))𝑟𝑇𝑇(𝑡)

𝑎 (1 −𝑇(𝑡)

𝐿𝑇)

𝑏

(8.23)


Os parâmetros e condições iniciais das equações são fornecidos pela tabela IV e VII. O

parâmetro de cada fármaco é dado pela tabela VIII. Tais valores são baseados em [3] e [9], sendo

que apenas alguns foram alterados conforme a necessidade de cada ordem fracionária.

Tabela VII. Parâmetros do tumor e do tratamento utilizados nas simulações com PK/PD para derivadas

fracionárias.

𝛼 𝑇0 a b 𝑟𝑇 𝐿𝑇 𝑟11 𝑟12 𝑟21 𝑟22 𝐾𝑠1 𝐾𝑠2 𝐾𝑖34 𝑡𝑠𝑡𝑎𝑟𝑡

1.1 1.0 2.0 3.1 0.005 100 0.005 0.0 0.0 0.2 0.001 0.001 1.8 500

0.9 1.0 2.0 3.1 0.01 100 0.005 0.0 0.0 0.2 0.1 0.01 2 700

0.8 1.0 2.0 3.1 0.02 100 0.005 0.0 0.0 0.2 0.001 0.001 1.3 1500

61

Tabela VIII. Parâmetros das fármacos utilizadaos na simulação de derivadas fracionárias.

Combinação 𝑑34


𝜏 1 1 1 1

𝑘𝑑 1 1 1 0.5

𝑘𝑎 0.01 0.01 0.01 0.01

𝐶0 0.008 0.008 0.008 0

𝐶0𝑏𝑎𝑠𝑒 0.002 0.002 0.002 0.002

𝐾𝑟 0 0 0 0

𝐿𝑟 0 0 0 0

𝐶𝐼 1 - -

O primeiro modelo analisado é para o caso onde 𝛼 = 1.1, ilustrado pela figura 44. Nota-se um

comportamento mais oscilatório se comparado aos outros modelos, mesmo utilizando valores de

parâmetros menores. A quimioterapia é eficiente em eliminar o tumor, no entanto acaba por reduzir

bastante a população de osteoclastos e osteoblastos, de modo que não ocorre uma recuperação da

densidade óssea perdida a curto prazo.

Figura 44. Comportamento das populações para α=1.1 para o modelo PK/PD.

A figura 45 mostra o resultado quando 𝛼 = 0.9. O tratamento proporciona uma brusca queda

no tumor, o que sugere uma forte ação da quimioterapia na eliminação das células metasticas. Esse

facto é importante na hora de decidir as dosagens dos fármacos utilizadas para evitar uma grande

toxicidade ao paciente. As populações de osteoclastos e osteoblastos recuperam o seu ciclo natural e

a densidade óssea começa a aumentar ligeiramente.

O ultimo modelo é dado pela figura 46 onde 𝛼 = 0.8. Nesse caso o grande amortecimento do

sistema faz com que as células convirjam para o seu novo estado estacionário com poucas

oscilações. O cancro é reduzido para um estado residual, porém não é totalmente eliminado.

62



Os modelos com derivadas fracionárias apresentaram um comportamento qualitativo

semelhante ao caso de ordem inteira. Nota-se que o os parâmetros de efeito máximo dos fármacos

apresentam uma forte influencia no comportamento oscilatório das populações uma vez que afetam a

taxa de apoptose e diferenciação dos osteoclastos.

Em todos os modelos ocorreu a redução do tumor e a recuperação da dinâmica natural das

células, no entanto não fora uninimidade a recuperação da densidade óssea. Como trabalho futuro

pode ser feito uma análise mais profunda sobre os efeitos da farmacocinética e farmacodinâmica,

como alterações na frequência de dosagem, concentrações iniciais de fármacos e a implementação

de outras combinações de fármacos.

63


A aplicação das derivadas fracionárias segue a mesma metodologia apresentada nos outros

tópicos. Utilizando o modelo dinâmico do osso para o caso sem tumor e sem tratamento apresentado

no capítulo 7 e, considerando o osso com uma dimensão espacial normalizada Ω = [0 1] que sofre

difusão, pode-se escrever:

𝐷0 𝑡𝛼𝐶(𝑡, 𝑥) = 𝜎1

𝜕2

𝜕𝑡2𝐶(𝑡, 𝑥) + 𝛼1𝐶(𝑡, 𝑥)

𝑔11𝐵(𝑡, 𝑥)𝑔21 − 𝛽1𝐶(𝑡, 𝑥) (8.25)

𝐷0 𝑡𝛼𝐵(𝑡, 𝑥) = 𝜎2

𝜕2

𝜕𝑡2𝐵(𝑡, 𝑥) + 𝛼2𝐶(𝑡, 𝑥)

𝑔12𝐵(𝑡, 𝑥)𝑔22 − 𝛽2𝐵(𝑡, 𝑥) (8.26)

𝐷0 𝑡𝛼𝑧(𝑡, 𝑥) = 𝜎3

𝜕2

𝜕𝑡2𝑧(𝑡, 𝑥) − 𝑘1(𝑥)max{0, 𝐶(𝑡, 𝑥) − 𝐶𝑠𝑠} + 𝑘2(𝑥)max{0, 𝐵(𝑡, 𝑥) − 𝐵𝑠𝑠} (8.27)

com as seguintes condições iniciais 𝐶(0, 𝑥) = 𝐶0(𝑥), 𝐵(0, 𝑥) = 𝐵0(𝑥) e 𝑧(0, 𝑥) = 𝑧0. As condições de

fronteira são dadas por:

𝜕

𝜕𝑡𝐶(𝑡, 0) =

𝜕

𝜕𝑡𝐶(𝑡, 1) = 0 (8.28)

𝜕

𝜕𝑡𝐵(𝑡, 0) =

𝜕

𝜕𝑡𝐵(𝑡, 1) = 0 (8.29)

𝜕

𝜕𝑡𝑧(𝑡, 0) =

𝜕

𝜕𝑡𝑧(𝑡, 1) = 0 (8.30)

No entanto, as condições iniciais nulas não são atendidas. A mudança de variável necessária

é dada em anexo.

A tabela IX apresenta os valores dos parâmetros utilizados nas simulações e seus resultados

são mostrados nas figuras a seguir. A distribuição inicial de osteoblastos foi assumida constante e

igual ao seu estado estacionário, já os osteoclastos apresentaram uma distribuição inicial ilustrada

pela figura 19.

Tabela IX. Parâmetros das populações utilizados nas simulações com 1D.

𝛼 𝜎1 𝜎2 𝜎3 𝛼1 𝛼2 𝛽1 𝛽2 𝑔11 𝑔12 𝑔21 𝑔22 𝑘1 𝑘2

1.2 10−5 10−5 10−5 3.0 4.0 0.2 0.02 1.08 1.0 -0.5 0.0 0.45 0.0053

1.1 10−5 10−5 10−5 3.0 4.0 0.2 0.02 1.0 1.0 -0.5 0.0 0.45 0.005

0.9 10−5 10−5 10−5 3.0 4.0 0.2 0.02 1.163 1.0 -0.5 0.0 0.45 0.0048

0.8 10−5 10−5 10−5 3.0 4.0 0.2 0.02 1.21 1.0 -0.5 0.0 0.45 0.003

64

Figura 47. Dinâmica da população 1D para α=1.2.



65


Pode-se notar que nos casos onde o valor da ordem de derivada é maior o sistema apresenta

um grande número de oscilações, mesmo para valores menores de 𝑔11. Tal comportamento era

esperado, uma vez que o nível de amortecimento nesses casos é menor. Já para os modelos onde a

ordem da derivada é menor, optou-se por manter o mesmo valor de 𝑔11 e ilustrar a diminuição cada

vez maior das oscilações até atingir o estado estacionário. Os valores de 𝑘1 e 𝑘2 precisaram ser

ajustados com base na amplitude das oscilações.

O comportamento qualitativo do sistema não sofreu alteração. O distúrbio inicial na população

de osteoclastos provoca um aumento na população de osteoblastos que, como consequência, gera a

redução do número de osteoclastos. O sistema segue essa dinâmica até atingir o estado estacionário

e a densidade óssea recupera o seu valor inicial.


O modelo com o tumor e sem tratamento para a derivada fracionária pode ser descrito, com

base no sistema do capítulo 7, por:


𝜕2

𝜕𝑡2𝐶(𝑡, 𝑥) + 𝛼1𝐶(𝑡, 𝑥)


)𝐵(𝑡, 𝑥)


)− 𝛽1𝐶(𝑡, 𝑥) (8.31)


𝜕2

𝜕𝑡2𝐵(𝑡, 𝑥) + 𝛼2𝐶(𝑡, 𝑥)

𝑔12 (1+𝑟12𝑇(𝑡,𝑥)𝐿𝑇

)⁄𝐵(𝑡, 𝑥)


)− 𝛽2𝐵(𝑡, 𝑥) (8.32)

𝐷0 𝑡𝛼𝑇(𝑡, 𝑥) = 𝜎4

𝜕2

𝜕𝑡2𝑇(𝑡, 𝑥) + 𝑟𝑇𝑇(𝑡, 𝑥)

𝑎 (1 −𝑇(𝑡, 𝑥)

𝐿𝑇)

𝑏

(8.33)


𝜕2


66

com condições iniciais (0, 𝑥) = 𝐶0(𝑥), 𝐵(0, 𝑥) = 𝐵0(𝑥), 𝑧(0, 𝑥) = 100 e 𝑇(0, 𝑥) = 𝑇0(𝑥). As condições

de fronteira são dadas por:

𝜕

𝜕𝑡𝐶(𝑡, 0) =

𝜕

𝜕𝑡𝐶(𝑡, 1) = 0 (8.35)

𝜕

𝜕𝑡𝐵(𝑡, 0) =

𝜕

𝜕𝑡𝐵(𝑡, 1) = 0 (8.36)

𝜕

𝜕𝑡𝑇(𝑡, 0) =

𝜕

𝜕𝑡𝑇(𝑡, 1) = 0 (8.37)

𝜕

𝜕𝑡𝑧(𝑡, 0) =

𝜕

𝜕𝑡𝑧(𝑡, 1) = 0 (8.38)

no entanto, vale reforçar que tais condições não representam fielmente a realidade.

Novamente é necessário realizar as mudanças de variáveis para que o sistema apresente

condições iniciais nulas. Os passos são os mesmos descritos anteriormente e estão dados em anexo.

Tabela X. Parâmetros das populações utilizados nas simulações com tumor em 1D.

𝛼 𝜎1 𝜎2 𝜎3 𝜎4 𝛼1 𝛼2 𝛽1 𝛽2 𝑔11 𝑔12 𝑔21 𝑔22 𝑘1 𝑘2

1.2 10−5 10−5 10−5 10−5 3.0 4.0 0.2 0.02 1.05 1.0 -0.5 0.0 0.45 0.0043

1.1 10−5 10−5 10−5 10−5 3.0 4.0 0.2 0.02 1.0 1.0 -0.5 0.0 0.45 0.0043

0.9 10−5 10−5 10−5 10−5 3.0 4.0 0.2 0.02 1.1 1.0 -0.5 0.0 0.45 0.0043

0.8 10−5 10−5 10−5 10−5 3.0 4.0 0.2 0.02 1.1 1.0 -0.5 0.0 0.45 0.0043

Tabela XI. Parâmetros do tumor utilizados nas simulações em 1D.

𝛼 𝑇0 a b 𝑟𝑇 𝐿𝑇 𝑟11 𝑟12 𝑟21 𝑟22

1.2 1.0 2.0 3.1 0.0033 100 0.005 0.0 0.0 0.2

1.1 1.0 2.0 3.1 0.005 100 0.005 0.0 0.0 0.2

0.9 1.0 2.0 3.1 0.01 100 0.005 0.0 0.0 0.2

0.8 1.0 2.0 3.1 0.02 100 0.005 0.0 0.0 0.2

A tabela X e XI mostra os parâmetros utilizados para as simulações. A distribuição inicial de

osteoclastos e osteoblastos são os mesmos utilizados no capítulo anterior e o tumor encontra-se

inicialmente localizado no lado direito de Ω. Vale relembrar que tais valores foram baseados em [3].

67

Figura 51. Dinâmica do modelo com tumor 1D para α=1.2.

Figura 52. Dinâmica do modelo com tumor 1D para α=1.1..

68



A medida que o tumor se desenvolve e propaga pelo resto do osso acaba por afetar a

dinâmica das células e faz com que elas estabilizem em um novo estado estacionário. Para ordens

fracionárias elevadas o sistema tende a oscilar mais e por isso reduziu-se o valor de 𝑔11 de modo que

houvesse uma convergência da dinâmica. Para o caso de ordens fracionárias menores que um o

número de oscilações diminui muito devido a característica mais amortecida.

Nota-se também uma velocidade de desenvolvimento do tumor proporcional a ordem da

derivada. Sendo muito mais rápido nos casos onde 𝛼 é maior.


O primeiro modelo de tratamento com derivadas fracionárias é o de uso dos inibidores de

proteassoma. Tal modelo é obtido através da adição dos termos 𝑉1(𝑡) e 𝑉2(𝑡) nas equações do tópico

anterior. Esses parâmetros representam a ação da terapia, iniciada em 𝑡𝑠𝑡𝑎𝑟𝑡, no organismo.

69


𝜕2

𝜕𝑡2𝐶(𝑡, 𝑥) + 𝛼1𝐶(𝑡, 𝑥)


)𝐵(𝑡, 𝑥)


)− 𝛽1𝐶(𝑡, 𝑥) (8.39)


𝜕2

𝜕𝑡2𝐵(𝑡, 𝑥) + 𝛼2𝐶(𝑡, 𝑥)

𝑔12 (1+𝑟12𝑇(𝑡,𝑥)𝐿𝑇

)⁄𝐵(𝑡, 𝑥)


)− (𝛽2 − 𝑉1(𝑡))𝐵(𝑡, 𝑥) (8.40)


𝜕2

𝜕𝑡2𝑇(𝑡, 𝑥) + (𝑟𝑇 − 𝑉2(𝑡))𝑇(𝑡, 𝑥)

𝑎 (1 −𝑇(𝑡, 𝑥)

𝐿𝑇)

𝑏

(8.41)


𝜕2


As condições de fronteira são dadas por:

𝜕

𝜕𝑡𝐶(𝑡, 0) =

𝜕

𝜕𝑡𝐶(𝑡, 1) = 0 (8.43)

𝜕

𝜕𝑡𝐵(𝑡, 0) =

𝜕

𝜕𝑡𝐵(𝑡, 1) = 0 (8.44)

𝜕

𝜕𝑡𝑇(𝑡, 0) =

𝜕

𝜕𝑡𝑇(𝑡, 1) = 0 (8.45)

𝜕

𝜕𝑡𝑧(𝑡, 0) =

𝜕

𝜕𝑡𝑧(𝑡, 1) = 0 (8.46)

e as condições iniciais dadas por 𝐶(0, 𝑥) = 𝐶0(𝑥), 𝐵(0, 𝑥) = 𝐵0(𝑥), 𝑧(0, 𝑥) = 100% e 𝑇(0, 𝑥) = 𝑇0(𝑥).

.

Tabela XII. Parâmetros do tumor e do tratamento com inibidores de proteassoma utilizados na simulação em 1D

com derivadas fracionárias.

𝛼 𝑇0 a b 𝑟𝑇 𝐿𝑇 𝑟11 𝑟12 𝑟21 𝑟22 𝑣1 𝑣2 𝑡𝑠𝑡𝑎𝑟𝑡

0.9 1.0 2.0 3.1 0.01 100 0.005 0.0 0.0 0.2 0.0009 0.017 600

0.8 1.0 2.0 3.1 0.02 100 0.005 0.0 0.0 0.2 0.00005 0.024 600

0.7 1.0 2.0 3.1 0.02 100 0.005 0.0 0.0 0.2 0.00001 0.031 600

A tabela XII apresenta os parâmetros utilizados nas simulações para o tratamento. Os outros

parâmetros são 𝛼1 = 3 , 𝛼2 = 4, 𝛽1 = 0.2, 𝛽2 = 0.02, 𝑔11 = 1.1, 𝑔12 = 1, 𝑔21 = −0.5, 𝑔22 = 0, 𝜎1 =

0.000001, 𝜎2 = 0.000001. Para as variáveis da equação óssea tem-se 𝜎3 = 0.000001, 𝑘1(𝑥) = 0.45,

𝑘2(𝑥) = 0.0048 e as condições iniciais são as mesmas dos tópicos anteriores. As figuras 55-57

mostram os respectivos resultados.

70

Figura 55. Dinâmica do modelo com tumor e tratamento 1D para α=0.9.


71


Os gráficos apresentam o tumor se desenvolvendo até o momento em que se inicia o

tratamento. A terapia mostra-se eficaz ao reduzir o tumor, levando-o a um valor residual e evitando

uma maior perda do osso. A apoptose dos osteoblastos é reduzida e as populações de osteoclastos e

osteoblastos retornam aos seus respectivos estados estacionários na ausência de tumor, no entanto

não ocorre aumento da densidade óssea.

Sendo assim, os resultados apresentam um comportamento qualitativo semelhante ao de ordem

𝛼 = 1, onde a grande diferença está em relação ao número de oscilações e o grau de amortecimento,

como discutido em outros tópicos.


O último modelo a ser analisado será o tratamento com os fármacos bisfosfanato e

denosumab, junto com a quimioterapia. Nesse caso considerou-se que a concentração de cada

fármaco exerce um mesmo efeito por todo o osso para um determinado instante t. O modelo

farmacocinético e farmacodinâmico é o mesmo apresentado no tópico 8.4 e o cálculo do efeito de

cada fármaco é feito através de um modelo Simulink apresentado em anexo. Por fim, as equações

dinâmicas do osso são:


𝜕2

𝜕𝑡2𝐶(𝑡, 𝑥) + 𝛼1𝐶(𝑡)


)𝐵(𝑡)


)(1+𝐾𝑠1𝑑1(𝑡)) − (1 + 𝐾𝑠2𝑑2(𝑡))𝛽1𝐶(𝑡) (8.47)


𝜕2

𝜕𝑡2𝐵(𝑡, 𝑥) + 𝛼2𝐶(𝑡)

𝑔12 (1+𝑟12𝑇(𝑡)𝐿𝑇

)⁄𝐵(𝑡)


)− 𝛽2𝐵(𝑡) (8.48)


𝜕2

𝜕𝑡2𝑇(𝑡, 𝑥) + (1 − 𝐾𝑖34𝑑𝑐34(𝑡))𝑟𝑇𝑇(𝑡)

𝑎 (1 −𝑇(𝑡)

𝐿𝑇)

𝑏

(8.49)


𝜕2

𝜕𝑡2𝑧(𝑡, 𝑥) − 𝑘1max{0, 𝐶(𝑡) − 𝐶𝑠𝑠} + 𝑘2max{0, 𝐵(𝑡) − 𝐵𝑠𝑠} (8.50)

onde 𝑑1(𝑡), 𝑑2(𝑡) e 𝑑34(𝑡) são os efeitos farmacodinâmicos de denosumab, bisfosfanato e da

quimioterapia (interpretada como a combinação de dois fármacos) e 𝐾𝑠1, 𝐾𝑠2, 𝐾𝑖34 os efeitos máximos

correspondentes.

As condições iniciais são dadas por (0, 𝑥) = 𝐶0(𝑥), 𝐵(0, 𝑥) = 𝐵0(𝑥), 𝑧(0, 𝑥) = 100% e 𝑇(0, 𝑥) =

𝑇0(𝑥). As condições de fronteira são:

𝜕

𝜕𝑡𝐶(𝑡, 0) =

𝜕

𝜕𝑡𝐶(𝑡, 1) = 0 (8.51)

72

𝜕

𝜕𝑡𝐵(𝑡, 0) =

𝜕

𝜕𝑡𝐵(𝑡, 1) = 0 (8.52)

𝜕

𝜕𝑡𝑇(𝑡, 0) =

𝜕

𝜕𝑡𝑇(𝑡, 1) = 0 (8.53)

𝜕

𝜕𝑡𝑧(𝑡, 0) =

𝜕

𝜕𝑡𝑧(𝑡, 1) = 0 (8.54)

Os parâmetros e condições iniciais utilizados para a simulação são dados pela tabela X, XIII e

XIV e foram baseados nos valores apresentados em [9]. As figuras 58 e 59 mostram os resultados.

Tabela XIII. Parâmetros do tumor e do tratamento utilizados nas simulações com PK/PD em 1D com derivadas

fracionárias.

𝛼 𝑇0 a b 𝑟𝑇 𝐿𝑇 𝑟11 𝑟12 𝑟21 𝑟22 𝐾𝑠1 𝐾𝑠2 𝐾𝑖34 𝑡𝑠𝑡𝑎𝑟𝑡

0.9 1.0 2.0 3.1 0.01 100 0.005 0.0 0.0 0.2 0.001 0.001 2 700

0.8 1.0 2.0 3.1 0.02 100 0.005 0.0 0.0 0.2 0.001 0.001 1.8 600

Tabela XIV. Parâmetros dos fármacos utilizadaos na simulação em 1D com derivadas fracionárias.

Combinação 𝑑34


𝜏 1 1 1 1

𝑘𝑑 1 1 1 0.5

𝑘𝑎 0.01 0.01 0.01 0.01

𝐶0 0.008 0.008 0.008 0

𝐶0𝑏𝑎𝑠𝑒 0.002 0.002 0.002 0.002

𝐾𝑟 0 0 0 0

𝐿𝑟 0 0 0 0

𝐶𝐼 1 - -

Figura 58. Dinâmica do modelo com PK/PD 1D para α=0.9.

73

Figura 59. Dinâmica do modelo com PK/PD 1D para α=0.7.

Os fármacos bisfosfanato e denosumab conseguiram reduzir o constante aumento da

população de osteoclastos e, com isso, recuperar o ciclo natural das células. No entanto, essas

populações passam a convergir para um novo estado estacionário, abaixo dos níveis iniciais, o que

sugere que as dosagens e/ou frequência de toma fora elevado. No entanto, vale destacar que essa

redução não trouxe graves consequências para a densidade óssea. A quimioterapia conseguiu parar

o desenvolvimento do tumor e reduzi-lo a um estado residual de modo bem rápido, o que sugere uma

alta eficiência desse tratamento. O comportamento apresentado em cada um desses casos está

coerente com os estudos apresentados em [2,12,14,18,41] no que diz respeito a redução do tumor e

recuperação do ciclo natural das células.

Em relação aos termos derivativos, assim como nos outros exemplos, ocorre maiores

oscilações em ordens elevadas e maior amortecimento para ordens menores. No caso da simulação

com 𝛼 = 0.7 o tratamento iniciasse antes do tumor alcançar um nível muito elevado e, com isso, não

ocorre uma grande perda de osso, além de ser possível notar sua ligeria recuperação. Isso reforça a

importância e benefício de se diagnosticar a doença e iniciar o seu tratamento o quanto antes.

74

9. Conclusão

A dissertação apresentou o estudo da dinâmica de remodelamento do osso através de uma série

de equações diferenciais. Com isso, fora possível analisar o comportamento da população de

osteoclastos e osteoblastos para diferentes ambientes biológicos, como na presença/ausência da

doença Mieloma e de tratamento.

A aplicação de derivadas fracionárias não alterou o caráter qualitativo dos modelos, no qual as

principais mudanças observadas foram o aumento ou diminuição do seu nível de amortecimento, que

acaba refletindo no número de oscilações, velocidade de convergência e estabilidade. Derivadas

fracionárias de ordem mais elevada apresentaram uma maior quantidade de oscilações, velocidade

de convergência mais rápida e uma menor estabilidade. Já as ordens menores apresentam o

comportamento oposto, com uma redução na amplitude de oscilações, diminuição na velocidade de

convergência e um aumento na estabilidade. Esses comportamentos estão de acordo com a literatura

sobre derivadas fracionárias [31,33,43].

O primeiro modelo examinou os papéis cooperativos entre os reguladores autócrinos e parácrinos

no controle do remodelamento ósseo. O modelo é mais sensível ao efeito autócrino dos osteoclastos,

o que reflete o fato de que sua atividade de destruição óssea é muito ativa e essas células são

rapidamente recrutadas e removidas. Por outro lado, a população de osteoblastos é muito menos

ativa e uma mudança em seu número ocorre de forma mais lenta. Como consequência, muito mais

osteoblastos precisam ser recrutados para remodelar determinado ponto do osso.

A introdução da doença Mieloma múltiplo provocou o aumento da população de osteoclastos e

osteoblastos, que reflete uma tentativa do organismo em manter a atividade de destruição e

remodelação do osso em equilíbrio. No entanto, tal processo é seguido por uma diminuição na

oscilação de osteoblastos (se comparado ao caso onde não exista tumor) e resulta numa queda

significativa de sua atividade. Como consequência, ocorre uma queda da densidade óssea, sintoma

característico da doença.

Dois modelos de tratamento foram analisados, um utilizando fármacos conhecidos como

inibidores de proteassoma e outro com bisfosfanatos e denosumab. O primeiro apresentou um efeito

direto no crescimento nas células cancerígenas e na formação de osteoblastos, de modo a prevenir a

doença Mieloma. Tal comportamento está de acordo com os estudos realizados in vivo. Os

osteoclastos e o tumor foram reduzidos ao longo do tempo de tratamento, no entanto um aumento na

população de osteoblastos e formação óssea apresentam requerem um tempo maior para

acontecerem, o que sugere que esse tipo de tratamento leva muito mais tempo para substituir o osso

perdido do que para eliminar a doença.

O segundo método de tratamento inclui também os conceitos de farmacocinética e

farmacodinâmica, como a administração de doses, análise da concentração de fármacos no

organismo, combinação de fármacos e análise de seus efeitos. Os fármacos utilizados nesse caso

estimularam a apoptose de osteoclastos e uma redução na sua diferenciação (caso do denosumab e

bisfosfanato) e destruíram as células cancerígenas (quimioterapia). Nessa caso houve um atraso nos

efeitos descritos ao iniciar o tratamento, devido ao fato de que as concentrações de fármacos no

organismo ainda estarem baixas e gerarem pouco efeito no corpo. A medida que a concentração de

75

tais fármacos no plasma sanguíneo aumenta, pode-se notar uma potencialização do seu efeito na

dinâmica das células. O tumor acaba sendo extinto, a população de osteoclastos é reduzida e a

destruição do osso é controlada, o que está de acordo com os estudos realizados in vivo.

Todos os modelos também foram estendidos para o caso de uma dimensão espacial, com o

objetivo de trazer uma maior semelhança ao ambiente real do osso. Todos os casos, incluindo os

modelos com derivadas fracionárias, apresentaram um comportamento semelhante aos sistemas

sem dimensão espacial.

Pode-se concluir que, mesmo de uma forma simples, modelar os processos de regulação entre

osteoclastos e osteoblastos resulta em comportamentos complexos e não lineares. Com isso, é

praticamente impossível recriá-los in vitro, o que gera uma grande dificuldade no desenvolvimento de

novas terapias para o combate a essa doença.

Vale salientar algumas limitações que esses modelos apresentam, como:

O modelo apenas considera dois tipos de células na dinâmica. Sabe-se que outras células,

como os osteócitos também são importantes na dinâmica de remodelação [19].

Somente a formação de osteoclastos e osteoblastos é regulado por fatores parácrinos e

autócrinos, enquanto que a atividade das células e morte são assumidos proporcionais ao

número de células [19].

Os parâmetros que indicam o efeito autócrino e parácrino incluem a ação de diversos fatores

biológicos [3].

O processo de remodelação é ativado através de um aumento no número inicial de

osteoclastos acima do seu estado estacionário, no entanto não descreve como que ocorre

esse aumento em determinado local do osso.

Por último, assume-se que o crescimento do tumor é independete da perca óssea, no

entanto, conforme explicado, eles estão acoplados [15].

9.1. Trabalho futuro

Muitas adições podem ser feitas aos modelos propostos nessa dissertação. No âmbito do

cálculo fracionário pode-se realizar o estudo do comportamento do sistema quando se aplicam

derivadas fracionárias em relação ao espaço nos modelos que involvem uma dimensão, pois sabe-se

que há fenômenos de difusão anômala em sistemas biológicos e estes são modelados com derivadas

fracionárias [24]. Outra oportunidade seria o uso de ordens fracionárias variantes no tempo [42].

Outras possibilidades seriam incluir o efeito de hormônios, como o estrogênio, aos modelos

uma vez que se sabe a sua importância na reconstrução óssea e no desenvolvimento da metástase.

Muitos tratamentos utilizam hormônios esteroides para combater as doenças, de modo que sua

implementação pode trazer grandes contribuições aos estudos.

A inclusão dos conceitos de farmacocinética e farmacodinâmica ao tratamento com inibidores

de proteassoma também permitira uma melhor comparação com a terapia de bisfosfanatos e

denosumab.

Além disso, tais modelos farmacocinéticos e farmacodinâmicos possibilitam prever a

resistência do organismo ao fármaco, no entanto esse fato não fora deveras explorado e pode vir a

76

ser de grande ajuda na determinação da quantidade de fármacos e intervalos de tempo que devem

ser administrados de modo a maximizar os efeitos da terapia com o mínimo de toxicidade para o

paciente.

Implementar mecanismos que ativam o processo de remodelamento ósseo através de sinais

como hormônios e estímulo mecânico também contribuirá para uma maior fidelidade do modelo ao

caso real.

Uma outra abordagem possível seria incluir a ligação que existe entre a queda da densidade

óssea e o seu efeito no crescimento do tumor. Além de outros tipos de células que também

participam diretamento do processo de remodelamento ou indiretamente, através do estimulo a

produção de osteoclastos e osteoblastos.

Por fim, a expansão do modelo de uma dimensão espacial para o caso de duas e três

dimensões espaciais. Com isso, será possível uma simulação mais fiel ao ambiente em que o osso

está inserido.

77

Referências

[1] Abe, E., Yamamoto, M., Taguchi, Y., Czernick, B. L., O’Brien, C., A., Economides, A. N., Stahl, N.,

Jilka, R. L., and Manolagas, S. C. (2000). Essential requirement of BMPs-2/4 for both osteoblast and

osteoclast formation in murine bone marrow cultures from adult mice: antagonism by Noggin. J. Bone

Miner. Res., 15:663-673.

[2] Alexanian, R., Bergsargel, D. E., Migliore, P. J., Vaughn, W. K., and Howe C. D. (1968).

Melphalam therapy for plasma cell myeloma. Journal of Hematology, 31:1-10.

[3] Ayati, B. P., Edwards, C. M., Webb, G. F., and Wikswo, J. P. (2010). A mathematical model of

bone remodeling dynamics for normal bone cell populations and myeloma bone disease. Biology

direct, 5:28.

[4] Bilezikian, J. P., Raisz, L. G., and Rodan, G. A. (2002). Bone formation. In Principles of Bone

Biology, chapter 3. Academic Press, 2 edition.

[5] Bussard, K. M., Gay, C. V., and Mastro, A. M. (2008). The bone microenvironment in mestastase;

what is special about bone?. Cancer Metastasis rev., 27:41-55.

[6] Casimiro, S., Guise, T. A., and Chirgwin, J. (2009). The critical role of the bone microenvironmente

in cancer metastases. Mollecular and Cellular Endocrinology, 310:71-81.

[7] Cho, T.-J., Kim, J. A., Chung, C. Y., Yoo, W. J., Gerstendfeld, L. C., Einhorn, T. A., and Choi, I. H.

(2007). Expression and role of interleukin-6 in distraction osteogenesis. Calcif. Tissue Int., 80:192-200.

[8] Chou, T. -C. (2006). Theoretical basis, experimental desing and computerized simulation of

synergism and antagonism in drug combination studies. Pharmacological reviews, 58(3):621-681.

[9] Coelho, R. M. M. A. V. (2015). CancerSys – multiscale modeling for personalized therapy of bone

metastasis. EXPL/EMS-SIS/1954/2013.

[10] Deng., W. (2007). Short memory principle and a predictor-corrector approach for fractional

differential equations. Journal of Computational and Applied Mathematics, 206:174-188.

[11] Dhillon, S. and Kostrzewski, A. (2006). Basic pharmacokinetics. In Clinical Pharmacokinetics,

chapter 1. Pharmaceutical Press, 1 edition.

[12] Drake, M. T., Clarcke, B. L., and Khosla, S. (2008). Bisphosphonates: mechanism of action and

role in clinical practice. National Institute of Health, 83(9):1032-1045.

[13] Edwards, C. M., Lwin, S. T., Fowler, J. A., Oyajobi B. O., Zhuang J., Bates, A. L., and Mundy, G.

R. (2009). Myeloma cells exhibit an increase in proteasome activity and an enhanced response to

proteasome inhibition in the bone marrow microenvironment in vivo. Am. J. Hematol, 84(5):268-272.

[14] Fermand, J.-P., Levy, Y., Gerota, J., Benbunan, M., Cosset, J. M., Caastaigne, S., Seligmann, M.,

and Brouet, J.-C. (1989). Treatment of aggressive multiple myeloma by high-dose chemotherapy and

total body irradiation followed by blood stem cells autologous graft. Journal of Hematology, 73(1):20-

23.

78

[15] Frost., H., M. (1973). Bone remodeling and its relationship to metabolic bone disease. Charles C.

Thomas, Springfield, MA.

[16] Garret, I. R., Chen, D., Gutierrez, G., Zhao, M., Escobedo, A., Rossini, G., Harris, S. E., Gallwitz,

W., Kim, K. B., Hu, S., Crew, C. M., and Mundy, G. R. (2003). Selective inhibitors of the osteoblasts

proteasome stimulate bone formation in vivo and in vitro. J. Clin. Invest. 111(11):1771-1782.

[17] Golan, D. E., Tashjian, A. H., Armstrong, E. J., and Armstrong, A., W. (2012).

Pharmacodynamics. In Principles of Pharmacology, chapter 2. Wolters Kluwer, 3 edition.

[18] Henry, D. H., Costa, L., Goldwaser, F., and Hirsh, V. (2011). Randomized, Double-blind study of

denosumab versus zoledronic acid in the treatment of bone metastases in patients with advanced

cancer (excluding breast and prostate cancer) or multiple myeloma. J. Clin. Oncol., 29:1125-1132.

[19] Komarova, S. V., Smith, R. J., Dixon, S., Sims, S. M., and Wahl, L. M. (2003). Mathematical

model predicts a critical role for osteoclast autocrine regulation in the control of bone remodeling.

Bone, 33(2):206-215

[20] Kyle, R. A., and Rajkumar, S. V. (2009). Treatment of multiple myeloma: a comprehensive review.

National Institute of Health, 9(4): 278-288.

[21] LeBlanc, R., Catley, L. P., Hideshima, T., Lentzsch, S., Mitsiades, C. S., Mitsiades, N., Neuberg,

D., Goloubeva, O., Pien, C. S., Adams, J., Gupta, D., Richardson, P. G., Munshi, N. C., and Anderson,

K. C. (2002). Proteasome inhibitor PS-341 inhibits human myeloma cell growth in vivo and prolongs

survival in a murine model. Cancer research, 62:4996-5000.

[22] Li, Y., Chen, Y., and Podlubny, I. (2010). Stability of fractional order nonlinear dynamic systems>

Lyapunov direct method and generalized Mittag-Leffler stability. Computers and Mathematics with

applications, 59:1810-1821.

[23] Manolagas, S., C. (2000). Birth and death of bone cells: basic regulatory mechanisms and

implications for the pathogenesis and treatment of osteoporosis. Endocrine reviews, 21:115-137.

[24] Margin, R. L. (2004). Fractional calculus in bioengineering. Biomedical Engineering, 32(1):1-104.

[25] McCarthy, P. L., Owzar, K., and Hofmeister, C. G. (2012). Lenalidomide after stem-cell

transplantation for multiple myeloma. New England Journal of Medicine, 366:1770-1781.

[26] McCurdy, A. R., and Lacy, M. Q. (2013). Pomalidomide and its clinical potential for relapsed or

refractory multiple myeloma: an update for the hematologist. Ther. Adv. Hematol., 4(3):211-216.

[27] Mundy, G. R. (2002). Metastasis to bone: causes, consequences and therapeutic opportunities.

Nat. Rev. Cancer, 2(8): 584-593.

[28] Ortigueira, M., D. (2006). Riesz potential operators and inverses via fractional centred derivatives.

International Journal of Mathematics and Mathematical Sciences, 1-12 Article ID 48391.

79

[29] Oyajobi, B. O., Garrett, I. R., Gupta, A., Flores, A., Esparza, J., Muñoz, S., Zhao, M., and Mundy,

G. R. (2007). Stimulation of new bone formation by the proteasome inhibitor, bortezomib: implications

for myeloma bone disease. British J. Haematology, 139(3):434-438.

[30] Pennisi, A., Li, X., Ling, W., Khan, S., Zangari, M., and Yaccoby, S. (2008). The proteasome

inhibitor, bortezomib suppresses primary myeloma and stimulates bone formation in myelomatous and

nonmyelomatous bones in vivo. Am. J. Hematol, 84:6-14.

[31] Petrás, I. (2009). Stability of fractional-order systems with rational orders: a survey. Fractional

Calculus and applied analysis, 12(3):269-298.

[32] Pinheiro, J. V., Lemos, J. M., and Vinga, S. (2011). Nonlinear MPC of HIV-1 infection with

periodics inputs. IEEE Conference on Decision and Control and European Control Conference, pages

65-70.

[33] Podlubny, I. (1999). Fractional derivatives and integrals. In Fractional Differential Equations,

chapter 2. Academic Press, 1 edition.

[34] Podlubny, I. (2000) Matrix approach to discrete fractional calculus. Fractional Calculus and

applied analysis, 3(4):359-386.

[35] Podlubny, I., Chechkin, A., Skovranek, T., Chen, Y., and Jara, B. V. (2008). Matrix approach to

discretization of ODEs and PDEs of arbitrary real order. Accessed in March/April of 2015 in

http://www.mathworks.com/matlabcentral/fileexchange/22071

[36] Podlubny, I., Chechkin, A., Skovranek, T., Chen, Y., and Jara, B., M., V. (2009). Matrix approach

to discrete fractional calculus II: partial fractional differential equations. Journal of Computational

Physics, 228:3137-3153.

[37] Podlubny, I., Skovranek, T., Jara, B., M., V., Petras, I., Verbitsky, V., and Chen, Y. (2013). Matrix

approach to discrete fractional calculus III: non-equidistant grids, variable step length and distributed

orders. Royal Society, 3(4):359-386.

[38] Raggatt, L. J., and Partridge, N. C. (2010). Cellular and molecular mechanisms of bone

remodeling. Journal of Biological Chemistry, 285(33):25103-25108.

[39] Solari, F., Flamant, F., Cherel, Y., Wyers, M. and Jurdic, P. (1996). The osteoclasts generation:

an in vitro and in vivo study with a genetically labelled avian monocytic cell line. Jornal of Cell Science,

109:1203-1213.

[40] Suda, T., Takahashi, N., Udagawa, N., Jimi, E., Gillespie, M. T., and Martin, T. J. (1999).

Modulation of osteoclasts differentiation and function by the new members of the tumor necrosis factor

receptor and ligand families. Endocrine reviews, 20(3):345-357.

[41] Steger, G. G., and Bartsch, R. (2011). Denosumab for the treatment of bone metastases in breast

cancer: evidence and opinion. Ther. Adv. Med. Oncol., 3(5):233-243.

80

[42] Tepljakov, A., Petlenkov, E., and Belikov, J. (2011). FOMCON toolbox. Accessed in May of 2015

in http://www.fomcon.net/

[43] Valério, D., and Costa, J. S. (2013). Fractional calculus: real orderns. In An Introduction to

Fractional Control, chapter 2. Institution of Engineering and Technology, 1 edition.

http://www.fomcon.net/

81

Anexo

Implementação do modelo PK/PD em Simulink

A figura a seguir ilustra o modelo em Simulink utilizado, junto com o toolbox FOMCON, para

calcular o efeito de uma determinado fármaco no organismo através de conceitos de farmacocinética

e farmacodinâmica.

Figura 60. Modelo PK/PD em Simulink.

O Pulse Generator1 representa a administração oral cíclica do fármaco, enquanto que LTI PK

Oral2 calcula a concentração do fármaco no plasma, 𝐶𝑝(𝑠), dada pela seguinte função de

transferência fracionária:


𝑘𝑎(𝑠𝛼 + 𝑘𝑎)(𝑠

𝛼 + 𝑘𝑑) (A.1)

onde 𝐶0 é a concentração inicial de uma única dose e 𝑘𝑎 e 𝑘𝑑 são taxas de absorção e secreção do

corpo.

A figura 61 mostra o interior do bloco Drug Resistance2 e PD, respectivamente.

Figura 61. Blocos Drug Resistance2 e PD, respectivamente.

O lado esquerdo implementa as equações dadas por:

𝐶50(𝑡) = 𝑓(𝑡)𝐶50𝑏𝑎𝑠𝑒 (A.2)

𝑓(𝑡) = 1 + 𝐾𝑟 ∫ max {0, 𝐿𝑟 − 𝐶𝑝(𝜏)}𝑑𝜏

𝑡

0

(A.3)

82

que representam a resistência do organismo ao efeito do fármaco. O lado direito implementa o

modelo farmacodinâmico dado pelas equações:

𝑑(𝑡) =

𝐶𝑝(𝑡)

𝐶50(𝑡) + 𝐶𝑝(𝑡) (A.4)

Através desse modelo é possível encontrar os respectivos efeitos de cada fármaco no organismo,

𝑑(𝑡), que são implementados nas equações dinâmicas do tratamento com bisfosfanatos, denosumab

e quimioterapia.

Mudança de variável para o caso sem tumor e sem tratamento 0D

Define-se as seguintes variáveis:

Θ(t) = 𝐶(𝑡) − 𝐶0 (A.5)

Ψ(t) = 𝐵(𝑡) − 𝐵0 (A.6)

Φ(t) = 𝑧(𝑡) − 𝑧0 (A.7)

No caso de 𝛼 maior que um é necessário a implementação de mais uma condição inicial,

referente a primeira derivada de cada variável, de modo que tem-se:

𝑑𝐶0𝑑𝑡

=𝑑𝐵0𝑑𝑡

=𝑑𝑧0𝑑𝑡

= 𝜌 = 0 (A.8)

O modelo apresentado em (8.2), (8.3) e (8.4) pode ser então reescrito do seguinte formato, no

caso de 𝛼 ≤ 1:

𝐷0 𝑡𝛼Θ(𝑡) = 𝛼1(Θ(𝑡) + 𝐶0)

𝑔11(Ψ(t) + 𝐵0)𝑔21 − 𝛽1(Θ(𝑡) + 𝐶0) (A.9)

𝐷0 𝑡𝛼Ψ(𝑡) = 𝛼2(Θ(𝑡) + 𝐶0)

𝑔12(Ψ(t) + 𝐵0)𝑔22 − 𝛽2(Ψ(𝑡) + 𝐵0) (A.10)

𝐷0 𝑡𝛼Φ(t) = −𝑘1max{0, (Θ(𝑡) + 𝐶0) − 𝐶𝑠𝑠} + 𝑘2max{0, (Ψ(t) + 𝐵0) − 𝐵𝑠𝑠} (A.11)

Já para 1 < 𝛼 ≤ 2 deve-se introduzir a segunda condição inicial. Para tanto, considere 𝛾 =

𝛼 − 1. Pode-se reescrever a equação do seguinte modo:

𝐷0 𝑡𝛾(𝑑

𝑑𝑡Θ(𝑡) + 𝜌) = 𝛼1(Θ(𝑡) + 𝐶0)

𝑔11(Ψ(t) + 𝐵0)𝑔21 − 𝛽1(Θ(𝑡) + 𝐶0) (A.12)

𝐷0 𝑡𝛾(𝑑

𝑑𝑡Ψ(𝑡) + 𝜌) = 𝛼2(Θ(𝑡) + 𝐶0)

𝑔12(Ψ(t) + 𝐵0)𝑔22 − 𝛽2(Ψ(𝑡) + 𝐵0) (A.13)

𝐷0 𝑡𝛾(𝑑

𝑑𝑡Φ(𝑡) + 𝜌) = −𝑘1max{0, (Θ(𝑡) + 𝐶0) − 𝐶𝑠𝑠} + 𝑘2max{0, (Ψ(t) + 𝐵0) − 𝐵𝑠𝑠} (A.14)

onde agora apresenta condições iniciais nulas e pode ser resolvido pelo método das matrizes

triangulares. Para obter o valor correto das populações das células, bem como a densidade do osso,

basta utilizar a relação em (A.5), (A.6) e (A.7) entre as respectivas variáveis.

Mudança de variável para o caso com tumor e sem tratamento 0D

Como o modelo não apresenta condições iniciais nulas, uma mudança de variável é

necessária para a aplicação do método das matrizes triangulares. Definindo as seguintes variáveis

como:

Θ(t) = 𝐶(𝑡) − 𝐶0 (A.15)

Ψ(t) = 𝐵(𝑡) − 𝐵0 (A.16)

83

Π(𝑡) = 𝑇(𝑡) − 𝑇0 (A.17)

Φ(t) = 𝑧(𝑡) − 𝑧0 (A.18)

e para o caso de 𝛼 > 1 :

𝑑𝐶0𝑑𝑡

=𝑑𝐵0𝑑𝑡

=𝑑𝑇0𝑑𝑡

=𝑑𝑧0𝑑𝑡

= 𝜌 = 0 (A.19)

o sistema (8.7), (8.8), (8.9) e (8.10) pode ser reescrito no seguinte formato para 𝛼 ≤ 1:

𝐷0 𝑡𝛼Θ(𝑡) = 𝛼1(Θ(t) + 𝐶0)

𝑔11(1+𝑟11(Π(𝑡)+𝑇0)

𝐿𝑇)(Ψ(t) + 𝐵0)

𝑔21(1+𝑟21(Π(𝑡)+𝑇0)

𝐿𝑇)− 𝛽1(Θ(t) + 𝐶0)

(A.20)

𝐷0 𝑡𝛼Ψ(𝑡) = 𝛼2(Θ(t) + 𝐶0)

𝑔12 (1+𝑟12(Π(𝑡)+𝑇0)

𝐿𝑇)⁄(Ψ(t) + 𝐵0)

(𝑔22−𝑟22(Π(𝑡)+𝑇0)

𝐿𝑇)− 𝛽2(Ψ(t) + 𝐵0)

(A.21)

𝐷0 𝑡𝛼Π(𝑡) = 𝑟𝑇(Π(𝑡) + 𝑇0)

𝑎 (1 −(Π(𝑡) + 𝑇0)

𝐿𝑇)

𝑏

(A.22)

𝐷0 𝑡𝛼Φ(𝑡) = −𝑘1max{0, (Θ(t) + 𝐶0) − 𝐶𝑠𝑠} + 𝑘2max{0, (Ψ(t) + 𝐵0) − 𝐵𝑠𝑠} (A.23)

e para o caso 1 < 𝛼 ≤ 2

𝐷0 𝑡𝛾(𝑑

𝑑𝑡Θ(𝑡) + 𝜌) = 𝛼1(Θ(t) + 𝐶0)

𝑔11(1+𝑟11(Π(𝑡)+𝑇0)

𝐿𝑇)(Ψ(t) + 𝐵0)

𝑔21(1+𝑟21(Π(𝑡)+𝑇0)

𝐿𝑇)− 𝛽1(Θ(t) + 𝐶0) (A.24)

𝐷0 𝑡𝛾(𝑑

𝑑𝑡Ψ(𝑡) + 𝜌) = 𝛼2(Θ(t) + 𝐶0)

𝑔12

(1+𝑟12(Π(𝑡)+𝑇0)

𝐿𝑇)⁄(Ψ(t) + 𝐵0)

(𝑔22−𝑟22

(Π(𝑡)+𝑇0)𝐿𝑇

)− 𝛽

2(Ψ(t) + 𝐵0) (A.25)

𝐷0 𝑡𝛾(𝑑

𝑑𝑡Π(𝑡) + 𝜌) = 𝑟𝑇(Π(𝑡)+ 𝑇0)

𝑎 (1 −(Π(𝑡)+ 𝑇0)

𝐿𝑇)

𝑏

(A.26)

𝐷0 𝑡𝛾(𝑑

𝑑𝑡Φ(𝑡) + 𝜌) = −𝑘1max{0, (Θ(t) + 𝐶0) − 𝐶𝑠𝑠} + 𝑘2max{0, (Ψ(t) + 𝐵0) − 𝐵𝑠𝑠} (A.27)

onde 𝛾 = 𝛼 − 1. Agora é possível aplicar o método das matrizes triangulares. O valor correto das

populações é obtido com a substituição dos valores encontrados em na solução pelos da (A.15-18).

Mudança de variável para o caso sem tumor e sem tratamento 1D

Realizando a seguinte mudança de variável:

Θ(t, x) = 𝐶(𝑡, 𝑥) − 𝐶(0, 𝑥) (A.28)

Ψ(t, x) = 𝐵(𝑡, 𝑥) − 𝐵(0, 𝑥) (A.29)

Φ(t, x) = 𝑧(𝑡, 𝑥) − 𝑧(0, 𝑥) (A.30)

e no caso de 𝛼 > 1:

𝑑𝐶(0, 𝑥)

𝑑𝑡=𝑑𝐵(0, 𝑥)

𝑑𝑡=𝑑𝑧(0, 𝑥)

𝑑𝑡= 𝜌(𝑥) = 0 (A.31)

o sistema (8.25), (8.26) e (8.27) pode ser reescrito do seguinte formato para 𝛼 ≤ 1:

𝐷0 𝑡𝛼Θ(t, x) = 𝜎1

𝜕2

𝜕𝑡2(Θ(𝑡, 𝑥) + 𝐶(0, 𝑥)) + 𝛼1(Θ(𝑡, 𝑥) + 𝐶(0, 𝑥))

𝑔11(Ψ(t, x) + 𝐵(0, 𝑥))𝑔21

− 𝛽1(Θ(𝑡, 𝑥) + 𝐶(0, 𝑥))

(A.32)

𝐷0 𝑡𝛼Ψ(t, x) = 𝜎2

𝜕2

𝜕𝑡2(Ψ(t) + 𝐵(0, 𝑥)) + 𝛼2(Θ(𝑡, 𝑥) + 𝐶(0, 𝑥))

𝑔12(Ψ(t) + 𝐵(0, 𝑥))𝑔22

− 𝛽2(Ψ(t) + 𝐵(0, 𝑥))

(A.33)

𝐷0 𝑡𝛼Φ(t, x) = 𝜎3

𝜕2

𝜕𝑡2(Φ(t, x) + 𝑧(0, 𝑥)) − 𝑘1(𝑥)max{0, (Θ(𝑡, 𝑥) + 𝐶(0, 𝑥)) − 𝐶𝑠𝑠}

+ 𝑘2(𝑥)max{0, (Ψ(t) + 𝐵(0, 𝑥)) − 𝐵𝑠𝑠}

(A.44)

84


𝐷0 𝑡𝛾(𝑑

𝑑𝑡Θ(t, x) + 𝜌(𝑥)) = 𝜎1

𝜕2

𝜕𝑡2(Θ(𝑡, 𝑥) + 𝐶(0, 𝑥)) + 𝛼1(Θ(𝑡, 𝑥) + 𝐶(0, 𝑥))

𝑔11(Ψ(t, x) + 𝐵(0, 𝑥))𝑔21 − 𝛽1(Θ(𝑡, 𝑥) + 𝐶(0, 𝑥)) (A.45)

𝐷0 𝑡𝛾(𝑑

𝑑𝑡Ψ(t, x) + 𝜌(𝑥)) = 𝜎2

𝜕2

𝜕𝑡2(Ψ(t) + 𝐵(0, 𝑥)) + 𝛼2(Θ(𝑡, 𝑥) + 𝐶(0, 𝑥))

𝑔12(Ψ(t) + 𝐵(0, 𝑥))𝑔22 − 𝛽2(Ψ(t) + 𝐵(0, 𝑥)) (A.46)

𝐷0 𝑡𝛾(𝑑

𝑑𝑡Φ(t, x) + 𝜌(𝑥)) = 𝜎3

𝜕2

𝜕𝑡2(Φ(t, x) + 𝑧(0, 𝑥)) − 𝑘1(𝑥)max{0, (Θ(𝑡, 𝑥) + 𝐶(0, 𝑥)) − 𝐶𝑠𝑠}

+ 𝑘2(𝑥)max{0, (Ψ(t) + 𝐵(0, 𝑥)) − 𝐵𝑠𝑠} (A.47)

onde 𝛾 = 𝛼 − 1.

Mudança de variável para o caso com tumor e sem tratamento 1D

Realizando a seguinte mudança de variável:

Θ(t, x) = 𝐶(𝑡, 𝑥) − 𝐶(0, 𝑥) (A.48)

Ψ(t, x) = 𝐵(𝑡, 𝑥) − 𝐵(0, 𝑥) (A.49)

Π(𝑡, 𝑥) = 𝑇(𝑡, 𝑥) − 𝑇(0, 𝑥) (A.50)

Φ(t, x) = 𝑧(𝑡, 𝑥) − 𝑧(0, 𝑥) (A.51)

e para 𝛼 > 1

𝑑𝐶(0, 𝑥)

𝑑𝑡=𝑑𝐵(0, 𝑥)

𝑑𝑡=𝑑𝑇(0, 𝑥)

𝑑𝑡=𝑑𝑧(0, 𝑥)

𝑑𝑡= 𝜌(𝑥) = 0 (A.52)

o sistema (8.31), (8.32), (8.33) e (8.34) pode ser reescrito do seguinte formato para 𝛼 ≤ 1:

𝐷0 𝑡𝛼Θ(𝑡) = 𝜎1

𝑑2

𝑑𝑡2(Θ(t) + 𝐶(0, 𝑥)) + 𝛼1(Θ(t) + 𝐶(0,𝑥))

𝑔11(1+𝑟11(Π(𝑡)+𝑇(0,𝑥))

𝐿𝑇)(Ψ(t) +𝐵(0, 𝑥))

𝑔21(1+𝑟21(Π(𝑡)+𝑇(0,𝑥))

𝐿𝑇)

− 𝛽1(Θ(t) + 𝐶(0, 𝑥)) (A.53)

𝐷0 𝑡𝛼Ψ(𝑡) = 𝜎2

𝑑2

𝑑𝑡2(Ψ(t) + 𝐵(0,𝑥)) + 𝛼2(Θ(t) + 𝐶(0, 𝑥))

𝑔12 (1+𝑟12(Π(𝑡)+𝑇(0,𝑥)

0)

𝐿𝑇)⁄(Ψ(t) +𝐵(0, 𝑥))

(𝑔22−𝑟22(Π(𝑡)+𝑇(0,𝑥))

𝐿𝑇)

− 𝛽2(Ψ(t) + 𝐵(0,𝑥))

(A.54)

𝐷0 𝑡𝛼Π(𝑡) = 𝜎4

𝑑2

𝑑𝑡2(Π(𝑡) + 𝑇(0, 𝑥)) + 𝑟𝑇(Π(𝑡) + 𝑇(0, 𝑥))

𝑎 (1 −(Π(𝑡) + 𝑇(0, 𝑥))

𝐿𝑇)

𝑏

(A.55)

𝐷0 𝑡𝛼Φ(𝑡) = 𝜎3

𝑑2

𝑑𝑡2(Φ(𝑡) + 𝑧(0, 𝑥)) − 𝑘1max{0, (Θ(t) + 𝐶(0, 𝑥)) − 𝐶𝑠𝑠}

+ 𝑘2max{0, (Ψ(t) + 𝐵(0, 𝑥)) − 𝐵𝑠𝑠}

(A.56)


𝐷0 𝑡𝛾(𝑑

𝑑𝑡Θ(t, x) + 𝜌(𝑥))

= 𝜎1𝑑2

𝑑𝑡2(Θ(t) + 𝐶(0, 𝑥))

+ 𝛼1(Θ(t) + 𝐶(0, 𝑥))𝑔11(1+𝑟11

(Π(𝑡)+𝑇(0,𝑥))𝐿𝑇

)(Ψ(t) + 𝐵(0, 𝑥))

𝑔21(1+𝑟21(Π(𝑡)+𝑇(0,𝑥))

𝐿𝑇)

− 𝛽1(Θ(t) + 𝐶(0, 𝑥))

(A.57)

85

𝐷0 𝑡𝛾(𝑑

𝑑𝑡Ψ(t, x) + 𝜌(𝑥))

= 𝜎2𝑑2

𝑑𝑡2(Ψ(t) + 𝐵(0, 𝑥))

+ 𝛼2(Θ(t) + 𝐶(0, 𝑥))𝑔12 (1+𝑟12

(Π(𝑡)+𝑇(0,𝑥)0)𝐿𝑇

)⁄(Ψ(t) + 𝐵(0, 𝑥))

(𝑔22−𝑟22(Π(𝑡)+𝑇(0,𝑥))

𝐿𝑇)

− 𝛽2(Ψ(t) + 𝐵(0, 𝑥))

(A.58)

𝐷0 𝑡𝛾(𝑑

𝑑𝑡Π(t, x) + 𝜌(𝑥)) = 𝜎4

𝑑2

𝑑𝑡2(Π(𝑡) + 𝑇(0, 𝑥)) + 𝑟𝑇(Π(𝑡) + 𝑇(0, 𝑥))

𝑎 (1 −(Π(𝑡) + 𝑇(0, 𝑥))

𝐿𝑇)

𝑏

(A.59)

𝐷0 𝑡𝛾(𝑑

𝑑𝑡Φ(t, x) + 𝜌(𝑥))

= 𝜎3𝑑2

𝑑𝑡2(Φ(𝑡) + 𝑧(0, 𝑥)) − 𝑘1max{0, (Θ(t) + 𝐶(0, 𝑥)) − 𝐶𝑠𝑠}

+ 𝑘2max{0, (Ψ(t) + 𝐵(0, 𝑥)) − 𝐵𝑠𝑠}

(A.60)

onde 𝛾 = 𝛼 − 1.

Documents

Modelação da Metástase Óssea utilizando Derivadas ... · fracionárias em equações diferenciais, a análise do corpotamento dinâmico do remodelamento do osso na presença ou