NOVA MODELAGEM FRACIONÁRIA APLICADA À DINÂMICA TUMORAL (HPV 16) Lucas Kenjy Bazaglia Kuroda Tese apresentada à Universidade Estadual Paulista “Júlio de Mesquita Filho” para a ob- tenção do t́ıtulo de Doutor em Biometria. BOTUCATU São Paulo - Brasil Fevereiro – 2020 NOVA MODELAGEM FRACIONÁRIA APLICADA À DINÂMICA TUMORAL (HPV 16) Lucas Kenjy Bazaglia Kuroda Orientador: Prof. Dr. Rubens de Figueiredo Camargo Tese apresentada à Universidade Estadual Paulista “Júlio de Mesquita Filho” para a ob- tenção do t́ıtulo de Doutor em Biometria. BOTUCATU São Paulo - Brasil Fevereiro – 2020 ii iii Dedicatória Dedico este trabalho aos meus pais e irmãos que sempre apoiaram e incentivaram a minha vida acadêmica. Agradecimentos • Agradeço a Deus pelas oportunidades oferecidas e por me dar forças para en- carar todos os desafios; • Ao meu grande orientador Rubens de Figueiredo Camargo, obrigado por pas- sar um pouco dos seus conhecimentos, pelas sugestões propostas na elaboração deste trabalho, pelos ensinamentos, reflexões e mensagens não só no ramo acadêmico, mas também sobre como nós, pesquisadores e educadores, devemos encarar o Mundo. Refaço, mais uma vez, as grandes palavras desse grande professor que se tornou também um grande amigo: “Foi uma grande satisfação tê-lo como orientador”; • Aos professores Paulo Fernando de Arruda Mancera e Alexys Bruno Alfonso, que tiveram uma enorme contribuição nas sugestões, elaboração e resultados referententes à métodos computacionais, modelagem matemática e dinâmica tumoral. Muito obrigado! • À minha banca de doutorado, professor Rubens de Figueiredo Camargo, pro- fessor Fernando Luiz Pio dos Santos, professora Daniela dos Santos de Oliveira, professora Eliana Contharteze Grigoletto e o professor Junior Soares pelas su- gestões e discussões referentes ao presente trabalho; • À Coordenação de Aperfeiçoamento Pessoal de Ńıvel Superior (CAPES), pela concessaão da bolsa de doutorado; • Aos meus pais por me incentivarem nos estudos desde os primórdios de minha educação, tornando-me uma pessoa persistente na busca dos meus objetivos; vi • Ao meu irmão mais velho Tiago, pela paciência e exemplo de como devemos seguir na vida. Aos meus irmãos Paula, Matheus e Pedro, que estudaram comigo dia a dia até o ensino médio; • Ao meu amigo, inicialmente de doutorado e agora para a vida, Robinson Ta- voni que sempre esteve disposto a me ajudar na implementação de métodos computacionais e discussões referentes aos cálculo fracionário. • Aos meus grandes amigos de graduação: Eduardo Luppi, Evandro Tortora, Leandro Feitoza, Fernando Leandro e Juliana Finato, que puderam presenciar o ińıcio da minha jornada acadêmica e que incentivam até hoje o meu crescimento profissional. • Agradeço o meu grande amigo João Giglioti que me apoiou e acompanhou a finalização do doutorado. • Aos professores e funcionários do Departamento de Bioestat́ıstica; • Aos funcionários da Seção de Pós-Graduação do Instituto de Biociências. vii . “As fronteiras do cálculo fracionário são as da humanidade e temos uma infinidade de fenômenos fracionários no mundo em que vivemos”. Tenreiro Machado Siglas DMBA - Dimetilbenz(a)antraceno (Dimethylbenz(a)anthracene). DTM - Differential transform method; Método da transformada dife- rencial. GDTM - Generalized differential transform method; Método da trans- formada diferencial generalizada. HPV - Human Papillomavirus; Papilomav́ırus Humano. MF1 - Modelagem fracionária 1. MF2 - Modelagem fracionária 2. MF3 - Modelagem fracionária 3. MSDTM - Multi-step differential transform method; Método multi- passos com transformada diferencial. MSE - Mean square error; Erro quadrático médio. MSGDTM - Multi-step generalized differential transform method; Método multi-passos com transformada diferencial generalizada. RMSE - Root mean square error; Raiz do erro quadrático médio. TPA - Ant́ıgeno polipept́ıdeo tecidual. Sumário Página LISTA DE FIGURAS xii LISTA DE TABELAS xv RESUMO xvi SUMMARY xvii 1 INTRODUÇÃO 1 2 FUNDAMENTOS DO CÁLCULO FRACIONÁRIO 7 2.1 Integral Fracionária . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2 Derivada Fracionária de Riemann-Liouville . . . . . . . . . . . . . . . . . 10 2.3 Derivada Fracionária de Caputo . . . . . . . . . . . . . . . . . . . . . . . 13 2.4 Derivada Fracionária de Grünwald-Letnikov . . . . . . . . . . . . . . . . 16 3 CONTRIBUIÇÕES TEÓRICAS 18 3.1 Análise Dimensional na Modelagem Fracionária . . . . . . . . . . . . . . 18 3.1.1 Modelo Malthusiano . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.1.2 Oscilador Harmônico . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.1.3 Equação Loǵıstica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.1.4 Modelo de Gompertz . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.1.5 A Escolha do Parâmetro “τ” . . . . . . . . . . . . . . . . . . . . . . . 33 3.2 Análise do Método Multi-Passos com Transformada Diferencial Genera- lizada na Modelagem Fracionária . . . . . . . . . . . . . . . . . . . . . . 37 x 3.2.1 Método da Transformada Diferencial: Extensão e Generalização . . . . 39 3.2.2 Método da Transformada Diferencial Generalizada . . . . . . . . . . . 40 3.2.3 Método Multi-passos com Transformada Diferencial Generalizada . . . 45 3.2.4 Modelos Matemáticos . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.2.5 Modelagem Fracionária . . . . . . . . . . . . . . . . . . . . . . . . . . 46 4 MÉTODO DE GRÜNWALD-LETNIKOV 53 4.1 Discretização do Modelo de Malthus e Riccati . . . . . . . . . . . . . . . 56 5 MODELO DE CRESCIMENTO TUMORAL (HPV 16) 59 5.1 Conjunto de Dados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 5.2 Modelo de Gompertz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 5.3 Estimação dos parâmetros . . . . . . . . . . . . . . . . . . . . . . . . . . 64 5.4 Modelagem Fracionária . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 5.5 Programação / Softwares . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 5.6 Discretização . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 5.7 Resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 5.8 Observações Teóricas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 6 CONCLUSÕES 83 REFERÊNCIAS BIBLIOGRÁFICAS 85 7 APÊNDICES 103 7.1 Apêndice A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 7.1.1 Transformada de Laplace e Funções Especiais . . . . . . . . . . . . . . 103 7.1.2 Convolução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 7.1.3 Função Gama . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 7.1.4 Função Beta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 7.1.5 Função de Mittag-Leffler . . . . . . . . . . . . . . . . . . . . . . . . . . 108 7.1.6 Função Gel’fand-Shilov . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 xi 7.2 Apêndice B . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 7.3 Apêndice C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 7.4 Apêndice D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 7.4.1 Desenvolvimento Acadêmico . . . . . . . . . . . . . . . . . . . . . . . . 116 Lista de Figuras Página 1 Solução da equação fracionária de Malthus (25) para vários valores de β. 23 2 Solução da equação fracionária de Malthus (25) para 0 ≤ t ≤ 2 e 0 ≤ t ≤ 15, respectivamente. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3 Solução da equação fracionária do oscilador harmônico (32) até o instante t = 6 s e t = 15 s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 4 Solução da equação loǵıstica fracionária (43) para 0 ≤ t ≤ 10 e 0 ≤ t ≤ 50, respectivamente. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 5 Solução do modelo de Gompertz fracionária (50) para 0 ≤ t ≤ 10 e 0 ≤ t ≤ 50, respectivamente. . . . . . . . . . . . . . . . . . . . . . . . . . 32 6 Solução das diferentes modelagens no modelo Malthusiano (53). Sendo as linhas cont́ınuas referentes ao parâmetro rβ e as tracejadas quando assumimos que τ = 1 r . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 7 Método MSGDTM e solução anaĺıtica: As linhas tracejadas correspondem à solução anaĺıtica, enquanto as linhas cont́ınuas correspondem ao método MSGDTM para β = 1, 0, β = 0, 8, β = 0, 6, β = 0, 4 e β = 0, 2, respectiva- mente, de baixo para cima. . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 8 Método MSGDTM e solução anaĺıtica: As linhas tracejadas correspondem à solução anaĺıtica, enquanto as linhas cont́ınuas correspondem ao método MSGDTM para β = 1, 0, β = 0, 8, β = 0, 6, β = 0, 4 e β = 0, 2, respectiva- mente, de baixo para cima. . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 9 Solução da equação (66) utilizando o MSGDTM com um e dois passos. . 49 xiii 10 Solução da equação (66) utilizando o método MSGDTM com um, dois e quatro passos. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 11 Solução da equação (66) utilizando o método MSGDTM com 300 passos no intervalo de [0; 3]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 12 Método de Grünwald-Letnikov e solução anaĺıtica do modelo de Malthus (64) com r = 1, h = 0, 001: As linhas tracejadas correspondem à solução anaĺıtica, enquanto as linhas cont́ınuas correspondem ao método GL para β = 1, 0, β = 0, 8, β = 0, 6, β = 0, 4 e β = 0, 2, respectivamente, de baixo para cima. . . . . 57 13 Solução da equação de Riccati (66) utilizando o método de Grünwald- Letnikov para β = 0, 7, τ = 1 e h = 0, 001. . . . . . . . . . . . . . . . . . 58 14 Imagem A: As medições do tumor foram iniciadas quando o primeiro tumor tornou-se viśıvel, atingindo 1-2 mm de diâmetro. (i) Quando o primeiro tumor atingiu 3-4 mm de diâmetro, foi iniciado o tratamento com droga, (ii) até o tumor atingir cerca de 10 mm de diâmetro, (iii) altura em que o rato foi sacrificado. Imagem B: mostra uma patologia tumoral t́ıpica de papilomas de camundongo induzidos por DMBA / TPA (Imagem fornecida no documento auxiliar e também utilizada no artigo de Loizides et al. (2015)). . . . . . . . . . . . . . . . . . . . . . . . . . . 62 15 Modelagem fracionária 1 para diferentes valores de h (passo). . . . . . . . . . 68 16 Modelagem fracionária 2 para diferentes valores de h (passo). . . . . . . . . . 69 17 Modelagem fracionária 3 para diferentes valores de h (passo). . . . . . . . . . 70 18 Curvas de crescimento tumoral do Camundongo 01. . . . . . . . . . . . . 73 19 Curvas de crescimento tumoral do Camundongo 02. . . . . . . . . . . . . 73 20 Curvas de crescimento tumoral do Camundongo 03. . . . . . . . . . . . . 74 21 Curvas de crescimento tumoral do Camundongo 04. . . . . . . . . . . . . 74 22 Curvas de crescimento tumoral do Camundongo 05. . . . . . . . . . . . . 74 23 Curvas de crescimento tumoral do Camundongo 06. . . . . . . . . . . . . 75 24 Curvas de crescimento tumoral do Camundongo 07. . . . . . . . . . . . . 75 25 Curvas de crescimento tumoral do Camundongo 08. . . . . . . . . . . . . 75 xiv 26 Curvas com menores RMSE a partir do modelo clássico até β = 0, 50 e τ = 16, 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 27 Curvas com menores RMSE a partir do modelo clássico até β = 0, 989. . 78 28 Comportamento gráfico das três modelagens quando variado τ (modela- gem 1) e o valor da ordem da derivada fracionária β. . . . . . . . . . . . 79 29 Gráfico da função Γ(t), para valores reais. . . . . . . . . . . . . . . . . . 107 30 Gráficos da função de Mittag-Leffler, Eα(x) e Eα(x α), respectivamente, para 0 < α ≤ 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 31 Gráficos da função de Mittag-Leffler, Eα(x) e Eα(x α), respectivamente, para α = 1, 2, 3, 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 32 Gráficos da função φν , para 0 < ν < 1. . . . . . . . . . . . . . . . . . . . 112 33 Gráficos da função φν , para 1 < ν < 2. . . . . . . . . . . . . . . . . . . . 112 34 Gráfico da função φν , para 2 < ν < 3. . . . . . . . . . . . . . . . . . . . . 113 Lista de Tabelas Página 1 Modelagem para diferentes valores de τ . . . . . . . . . . . . . . . . . . . 37 2 Propriedades da transformada diferencial generalizada. Considera-se que r é uma constante e que as transformadas de x(t) e y(t) são X(k) e Y (k), respectivamente, com ordem β e centro em t0. . . . . . . . . . . . . . . . 44 3 Tempo de duplicação tumoral, µg, de cada camundongo em seus respec- tivos tumores. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 4 Menor erro quadátrico médio para a modelagem de Gompertz usual (RMSE - Usual), modelagem fracionária 1 (RMSE - MF1), modelagem fracionária 2 (RMSE - MF2) e modelagem fracionária 3 (RMSE - MF3) referente ao volume de dois tumores presentes em cada camundongo. . . 72 5 Desenvolvimento do programa informando o menor RMSE referente ao camundongo 07/T2 (Modelagem 1). . . . . . . . . . . . . . . . . . . . . . 81 6 Valor do parâmetro τ na Modelagem 1 para cada camundongo com seu respectivo tumor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 7 Desenvolvimento do programa informando o menor RMSE referente ao camundongo 06/T2 (Modelagem 3). . . . . . . . . . . . . . . . . . . . . . 82 8 Tabela de Transformadas integrais. . . . . . . . . . . . . . . . . . . . . . 121 NOVA MODELAGEM FRACIONÁRIA APLICADA À DINÂMICA TUMORAL (HPV 16) Autor: LUCAS KENJY BAZAGLIA KURODA Orientador: Prof. Dr. RUBENS DE FIGUEIREDO CAMARGO RESUMO O presente trabalho apresenta a nova modelagem fracionária, que con- sidera propriedades hereditárias e efeitos de memória, no modelo de Gompertz, para descrever a evolução do câncer causado pela infecção do HPV 16. Devido a varia- bilidade do desenvolvimento do câncer em humanos, utiliza-se o crescimento in vivo do tumor em camundongo transgênico que expressam os oncogenes E6 e E7 tratados com DMBA / TPA (inicializador e promotor do HPV 16) para capturar as carac- teŕısticas gerais dessa variabilidade. Resultados mostram que a inserção de um novo parâmetro na correção dimensional da modelagem fracionária, descreve, em com- paração ao modelo clássico, o progresso do volume tumoral em maior conformidade com os conjuntos de dados reais. Palavras-chave: Cálculo Fracionário, Modelagem Fracionária, Câncer, Correção Dimensional, Métodos Computacionais. A NEW FRACTIONAL MODELING APPLIED TO TUMOUR DYNAMICS (HPV 16) Author: LUCAS KENJY BAZAGLIA KURODA Advisor: Prof. Dr. RUBENS DE FIGUEIREDO CAMARGO SUMMARY The present work presents the fractional modeling, which considers hereditary properties and memory effects, to describe through the Gompertz model, the evolution of cancer caused by HPV 16 infection. Due to the variability of the development of cancer in humans, we used the in vivo growth of the transgenic mouse tumor expressing DMBA / TPA-treated E6 and E7 oncogenes (HPV 16 initiator and promoter) to capture the general characteristics of this variability. Results show that the insertion of a new parameter in the dimensional correction of fractional modeling describes, compared to the classical model, the progress of tumor volume in greater concorda with the actual data sets. Keywords: Fractional Calculus, Fractional Modeling, Cancer, Dimen- sional Correction, Computational Methods. 1 INTRODUÇÃO O cálculo fracionário está se tornando o cálculo da modernidade, onde a limitação das ordens de integração e derivação aos números inteiros é superada e se passa a trabalhar com as ordens desses operadores em valores arbitrários, reais ou complexos (de Carvalho & Ottoni, 2018). O significado da derivada dnf(x) dxn quando n não é inteiro, teve origem em 30 de setembro de 1695 quando L’Hôpital envia uma carta a seu amigo Leibniz perguntando o significado da derivada quando n fosse igual a 1/2. Leibniz afirmou que d 1 2f(x) dx 1 2 era igual a x √ dx : x. Não se sabe ao certo o sentido desta resposta, pois Leibniz foi sucinto no desenvolvimento desta questão, no entanto, o estudo por outros cientistas anos depois, mostram que o resultado obtido por Leibniz se assemelha a teoria em torno do cálculo fracionário mais recente (Camargo, 2009). Segundo Matlob & Jamali (2017), diante de diversas questões em re- ferência ao cálculo de ordem não inteira, cientistas, com o avanço da teoria, conse- guiram resultados intrigantes no ramo matemático e no campo de aplicações, em que exibem a importância desta teoria (Zhang et al., 2018). Pode-se citar: • Na f́ısica (Hifler, 2000; Caputo & Fabrizio, 2017); • Engenharia (Sabatier et al., 2007); • Sistemas financeiros (Chen, 2008); • Sistemas de controle (Baleanu et al., 2012); • Difusão (Metzler & Klafter, 2010); 2 • Comportamentos dinâmicos; – Sincronização de caos (Matouk, 2016); – Controle de consenso (Liu et al., 2017); – Problema de estabilização (Lu & Chen, 2009; Wang et al., 2016); – Bifurcação de Hopf (Deshpande et al., 2017; Huang et al., 2017). • Epidemiologia (Ahmed & Elgazzar, 2007); • Sistemas biológicos (Magin, 2010); • Dinâmica tumoral (Silva et al., 2019; Baba, 2019). Na natureza, para Sólis-Péres et al. (2019), muitos fenômenos f́ısicos têm descrições de ordens fracionárias “intŕınsecas”. Portanto, o cálculo de ordem não inteira em conjunto com a modelagem matemática são ferramentas necessárias para explicá-las (Zhang et al., 2018). A modelagem matemática é vista como o processo em que se trans- formam os problemas e eventos reais, em problemas matemáticos e resolvem os in- terpretando seus resultados no âmbito real (Bertolotto & Camargo, 2015). Neste contexto, a modelagem fracionária tem ganhado espaço em aplicações de situações reais, pois permite que os modelos de ordem fracionária contemplem comportamen- tos, situações que o cálculo usual não é capaz de captar. Segundo Area et al. (2016), nas últimas duas ou três décadas, um grande interesse tem sido dedicado ao estudo de equações diferenciais fracionárias. Esta modelagem mostra-se útil para descre- ver e explicar certos tipos de problemas e processos f́ısicos nos quais os efeitos de memória devem ser considerados. Esses modelos carregam informações de seus dife- rentes estados anteriores e podem ser uma maneira útil de incluir memória em um processo dinâmico, pois a ordem derivativa fracionária pode ser interpretada como um ı́ndice da memória (Du et al., 2013; Zhang et al., 2018). 3 Uma grande dificuldade em lidar com operadores fracionários, é a ex- trema dificuldade em resolver problemas analiticamente. Em alguns casos, podemos utilizar a metodologia das transformadas integrais (Camargo & Oliveira, 2015) para encontrar a solução de equações diferenciais fracionárias. No entanto, esta metodolo- gia possui uma limitação, não pode ser aplicada em modelos não lineares por se tratar de uma transformação linear. Assim, na maioria dos casos, não sabemos a solução exata para o problema e é preciso buscar uma aproximação numérica (Tavares et al., 2016). Segundo MacDonald et al. (2015), a computação no cálculo de soluções numéricas para equações diferenciais fracionárias pode ser intensiva devido ao efeito da não localidade (Akman Yildiz et al., 2018; Veeresha et al., 2019a; Owolab, 2019). Em geral, abordagens numéricas que dependem de truncamento, podem sofrer com altos graus de erros e imprecisões, como no Caṕıtulo 3.2 deste trabalho e em Kuroda et al. (2019), concluindo que o método multi-passos com transformada diferencial generalizada, dispońıvel na literatura, produz soluções equivocadas do modelo de estudo a partir do segundo passo. O problema é explicado em termos da não locali- dade da derivada de Caputo. Assim, MacDonald et al. (2015); Scherer et al. (2011); Ongun et al. (2013) apresentam o método fracionário de Grünwald-Letnikov, compu- tacionalmente eficiente, resulta em erros menores durante as simulações numéricas, além de levar em consideração as derivadas em tempos anteriores, o que preserva o efeito de memória (Scherer et al., 2011). Além do método de Grünwald-Letnikov, podemos citar os trabalhos (Diethelm & Luchko, 2004) que discutem um algoritmo para a solução numérica de problemas de valor inicial em equações diferenciais lineares de ordem fracionária, com coeficientes constantes, e derivadas fracionárias no sentido de Caputo. Em He et al. (2017), as soluções numéricas de equações lineares e não lineares de ordem fracionária são obtidas empregando o método “constructed conformable Adomian decomposition method”. El-Sayed et al. (2007), utilizam o método preditor-corretor do tipo Adams para a solução numérica de uma equação integral fracionária. Em 4 Bruno-Alfonso et al. (2019) a equação loǵıstica com derivada fracionária de Caputo é resolvida mediante um tratamento anaĺıtico-numérico, ou seja, num intervalo de tempo inicial, é usado o método da transformada diferencial generalizada, em que a série de potências fracionárias correspondente tem raio de convergência finito. Assim, este método é usado até metade do raio de convergência, a evolução posterior é obtida mediante generalização do método de Euler. Verificaram que a solução satisfaz numericamente a equação loǵıstica. A biologia é uma fonte rica para desenvolver e construir ideias ma- temáticas (Atici & Sengul, 2010; Bolton et al., 2014; Frunzo et al., 2019). Estudamos as aplicações de modelagem matemática e técnicas matemáticas para analisar os pro- blemas das biociências. Uma vez que um modelo é formulado, sua solução pode ser utilizada para confrontar a realidade usando técnicas matemáticas e o resultado pode ser comparado com observações. Se houver discrepâncias entre conclusão teórica e observações, o processo é repetido até que um modelo realmente satisfatório seja obtido (Salahshour et al., 2018). Neste trabalho, analisaremos como o modelo de ordem fracionária, com correção dimensional, descreve o crescimento tumoral conhe- cido como papilomav́ırus humano (HPV 16) (Loizides et al., 2015). A infecção do papilomav́ırus humano de alto risco (HPV) representa atualmente o fator de risco mais importante na gênese do carcinoma de colo uterino. O HPV 16, de alto risco oncológico, foi considerado pela IARC (Internacional Agency for Research on Cancer) como altamente carcinogênico para a raça humana (Rama et al., 2006). A relação entre o tumor e seu desenvolvimento temporal é de especial interesse, uma vez que a estimativa do crescimento é muito cŕıtica na prática cĺınica (Atici & Sengul, 2010). Em 1825, Benjamin Gompertz introduziu o que hoje co- nhecemos como função de Gompertz, função sigmóide aplicável a vários fenômenos de crescimento, em particular o crescimento de tumores, como visto nos trabalhos de Bolton et al. (2014), Frunzo et al. (2019) e Loizides et al. (2015). A equação diferencial de Gompertz descreve modelos de crescimento e esses modelos podem ser 5 estudados com base nos parâmetros, podendo ser a taxa de crescimento ou taxa de desaceleração do crescimento (Atici & Sengul, 2010). Muitos estudos têm focado na entrega de medicamentos no tratamento do câncer, a fim de obter uma melhor taxa de sobrevivência (Javadi et al., 2018). Entender a cinética do crescimento do tumor permite, por exemplo, que os médicos determinem o melhor tratamento dispońıvel. Há uma infinidade de dados experi- mentais dispońıveis para análise sistemática e, portanto, a modelagem matemática pode contribuir significativamente em muitas áreas de investigação experimental do câncer (Sengul, 2010). Neste trabalho, fizemos o estudo referente ao câncer causado pela in- fecção do HPV 16. Para descrever a evolução temporal do volume tumoral, utili- zamos a modelagem de ordem não inteira no modelo de Gompertz, considerando a correta dimensão do operador. As soluções foram exibidas através da aproximação do método computacional de Grünwald-Letnikov, que preserva o efeito de mémoria, do cálculo fracionário. O texto esta organizado da seguinte forma: no Caṕıtulo 2, apresen- tamos os fundamentos de cálculo fracionário. No Caṕıtulo 3, mencionamos algu- mas contribuições teóricas divididas em duas partes: a primeira parte é referente a dimensão das equações diferenciais quando realizada a modelagem fracionária. A segunda parte é relacionada ao Método multi-passos com transformada diferencial generalizada (MSGDTM), utilizado para resolver equações e sistemas de ordem fra- cionárias, que trata o cálculo fracionário como um operador local. No Caṕıtulo 4, estudamos o método de Grüwald-Letnikov, este método trata a derivada fracionária como um operador não local, que, aplicado em equações diferenciais fracionárias, é posśıvel descrever fenômenos com efeito memória. Por fim, no Caṕıtulo 5, realizamos três diferentes modelagens fracionárias (vistas no Caṕıtulo 3) referente ao volume tu- moral em camundongos transgênicos que expressam os oncogenes1 E6 e E7 do HPV 1Oncogenes são genes relacionados com o aparecimento e crescimento de tumores, malignos ou benignos. 6 16. 2 FUNDAMENTOS DO CÁLCULO FRA- CIONÁRIO O objetivo deste Caṕıtulo é apresentar as teorias advindas do cálculo fracionário como a integração e a derivação fracionária que serão definidas com base em (Camargo & Oliveira, 2015) e em (Kuroda, 2016). A teoria referente às funções especiais e transformada de Laplace, utilizadas neste caṕıtulo, podem ser vista no Apêndice A. Vale pontuar que diante a existência de diversas transformadas in- tegrais, este trabalho utiliza apenas transformada de Laplace por ser aplicada em problemas que possuem dependência temporal, como por exemplo, modelos tumo- rais. 2.1 Integral Fracionária A formalização da integração é tratada antes da derivação, por ser a melhor maneira de realizar a conexão entre integral e a derivada de ordem arbitrária, mesmo que esta ordenação seja pouco frequente no cálculo usual (Camargo & Oli- veira, 2015). Antes de apresentar o conceito de integração fracionária, mostra-se que uma integral de ordem n, com n ∈ N, de uma função f(t) com t ∈ R, pode ser vista como um produto de convolução de Laplace entre f(t) e a função Gel’fand-Shilov de ordem n (Camargo & Oliveira, 2015). 8 Define-se o operador integral I de ordem inteira2, sobre f(t) como, If(t) = ∫ t 0 f(t1)dt1, com isso I2f(t) = I[If(t)] ∫ t 0 ∫ t1 0 f(t2)dt2dt1. Admitindo a generalização, definimos o operador In, Inf(t) = ∫ t 0 ∫ t1 0 ∫ t2 0 . . . ∫ tn−2 0 ∫ tn−1 0 f(tn)dtndtn−1 . . . dt3dt2dt1. Teorema 1 Sejam n ∈ N, t ∈ R+ e f(t) uma função integravel, então Inf(t) = φn(t) ∗ f(t) := ∫ t 0 φn(t− τ)f(τ)dτ = ∫ t 0 (t− τ)n−1 (n− 1)! f(τ)dτ, (1) sendo φn(t) a função de Gel’fand-Shilov de parâmetro n ∈ N e ∗ denotando o produto de convolução de Laplace. A prova deste teorema pode ser vista em (Camargo & Oliveira, 2015). Definição 1 (Integral fracionária de Riemann-Liouville em intervalos finitos) Seja Ω = [a, b] com −∞ < a < b < ∞ um intervalo finito no eixo real R. As integrais de ordem arbitrárias de Riemann-Liouville, denotadas por, (Iνa+f)(t) ≡ aI ν t f(t) e (Iνb−f)(t) ≡ tI ν b f(t), ambas de ordem ν ∈ C, com Re(ν) > 0, são definidas, respectivamente, por aI ν t f(t) := 1 Γ(ν) ∫ t a f(τ)(t− τ)ν−1dτ, (2) sendo t > a. tI ν b f(t) := 1 Γ(ν) ∫ b t f(τ)(τ − t)ν−1dτ, (3) sendo t < b. 2Definimos I0f(t) = f(t). 9 Neste trabalho, com exceção do Caṕıtulo 3, utilizamos a definição (2) quando a = 0, uma vez que a modelagem fracionária será realizada em modelos tumorais no instante em que t = 0 e simplificaremos a notação para Iνf(t). Assim segue que, 0I ν t f(t) = Iνf(t) := φν(t) ∗ f(t) = 1 Γ(ν) ∫ t 0 f(τ)(t− τ)ν−1dτ, (4) sendo t > 0. Proposição 1 A partir da definição de integral fracionária de Riemann-Liouville dada pela expressão (4) e, sendo f(t) = tµ, temos Iνtµ = Γ(µ+ 1) Γ(µ+ ν + 1) tµ+ν . (5) Vejamos, Iνtµ = 1 Γ(ν) ∫ t 0 (t− τ)ν−1τµdτ = 1 Γ(ν) ∫ t 0 [ t ( 1− τ t )]ν−1 τµdτ. Com a mudança de variável u = τ t , segue que Iνtµ = 1 Γ(ν) ∫ 1 0 t1tν−1(1− u)ν−1(ut)µdu = tν+µ Γ(ν) ∫ 1 0 (1− u)ν−1uµdu. Pela definição da função Beta e pela relação existente entre as funções Beta e Gama, vide Apêndice A, temos Iνtµ = tν+µ Γ(ν) B(µ+ 1, ν) = Γ(µ+ 1) Γ(µ+ ν + 1) tµ+ν . Transformada de Laplace da Integral Fracionária de Riemann-Liouville A transfomada de Laplace da integral fracionária de Riemann-Liouville de ordem ν de uma função f é L [Iνf(t)] = L [f(t)] sν . (6) A obtenção desse resultado pode ser visto em (Camargo & Oliveira, 2015). 10 2.2 Derivada Fracionária de Riemann-Liouville Segundo Atangana (2016), o conceito de derivada foi introduzido para descrever a taxa de variação por uma função e posteriormente usado para construir equações matemáticas que descrevem o comportamento dos problemas do mundo real. No entanto, devido às complexidades dos problemas realistas, o conceito de derivada fracionária pode ser, em alguns casos, mais adequado para modelar o pro- blema do mundo real do que a derivada usual (Kilbas et al., 2006; Chechkin et al., 2005; Chen & Moore, 2002). Assim, muitos pesquisadores têm dedicado sua atenção no desenvolvimento de uma nova definição de derivada fracionária. O uso dessas derivadas fracionárias e sua aplicação, na última década, obtiveram uma melhoria notável na descrição de fenômenos reais, conforme demonstrado por muitos artigos cient́ıficos, conferências, monografias, dentre outros. Pode-se encontrar mais de uma definição de derivada fracionária na literatura, eles variam de Riemann-Liouville ao operador fracionário de Caputo-Fabrizio (Caputo & Frabrizio, 2015; Losada & Nieto, 2015; Atangana & Badr, 2015), recentemente proposta (Atangana, 2016). Consideramos a definição com o nome de Riemann e Liouville, que é a generalização natural da fórmula de Cauchy de uma função f(t) (Huang et al., 2017; Atangana & Gómez-Aguiar, 2017; Zhang et al., 2018; Kuroda, 2016). Definição 2 Sejam β ∈ C com Re(β) > 0 e n − 1 < Re(β) ≤ n, com n ∈ N. A derivada fracionária de Riemann-Liouville de ordem β de f(t), para t > 0, é definida como RLD βf(t) = Dn[In−βf(t)] = Dn[φn−β(t) ∗ f(t)], (7) sendo Dn a derivada usual de ordem inteira n e In−β a integral fracionária de Riemann-Liouville de ordem n− β definida em (4). Note que, se β ∈ N, a derivada fracionária de Riemann-Liouville equi- vale a derivada do cálculo usual, RLD βf(t) = Dn[In−βf(t)] = Dn[I0f(t)] = Dnf(t). 11 Sendo n ∈ N, temos que RLD nInf(t) = f(t). No entanto, RLD n não é o operador inverso à direita de In. De fato (Gorenflo & Mainardi, 2008), In [RLD nf(t)] = f(t)− n−1 ∑ k=0 f (k)(0+) tk k! , t > 0. Desta maneira, a derivada fracionária de ordem β > 0 segundo Riemann-Liouville também foi definida como sendo a operação inversa à esquerda da integral fracionária, ou seja, sendo Id o operador identidade, RLD βIβ = Id. Proposição 2 A derivada fracionária de Riemann-Liouville, de ordem β ∈ C com Re(β) > 0, da função f(t) = tµ, com µ > −1, é da forma RLD βtµ = Γ(µ+ 1) Γ(µ− β + 1) tµ−β . (8) Vejamos, pela definição da derivada fracionária na equação (7), sendo ν = n− β temos RLD βtµ = Dn[Iνtµ]. Vimos na Proposição 1 de integral fracionária, equação (5), Iνtµ = Γ(µ+ 1) Γ(µ+ ν + 1) tµ+ν . Assim, segue que 12 RLD βtµ = Dn[Iνtµ] = Dn [ Γ(µ+ 1) Γ(µ+ ν + 1) tµ+ν ] = Γ(µ+ 1) Γ(µ+ ν + 1) Γ(µ+ ν + 1) Γ(µ+ ν − n + 1) tµ+ν−n = Γ(µ+ 1) Γ(µ+ ν − n + 1) tµ+ν−n· Realizando a substituição ν = n − β na equação acima, temos que a derivada fra- cionária de Riemann-Liouville da função tµ é RLD βtµ = Γ(µ+ 1) Γ(µ− β + 1) tµ−β . Proposição 3 A derivada fracionária de Riemann-Liouville, de ordem β ∈ C com Re(β) > 0, da função constante f(t) = c, é da forma RLD βc = c tβΓ(1− β) · Pela definição da derivada fracionária (7) e com ν = n− β temos RLD βc = Dn[Iνc] = Dn[In−βc]. Temos que, Iνc = tνc Γ(ν + 1) . Logo, RLD βc = Dn[In−βc] = Dn [ tn−βc Γ(n− β + 1) ] = c Γ(n− β + 1) Γ(1− β) t−β Γ(n− β + 1) = c tβΓ(1− β) · 13 Transformada de Laplace da Derivada Fracionária de Riemann-Liouville A transformada de Laplace da derivada fracionária de Riemann-Liouville é L [RLD βf(t)] = sβF (s)− n−1 ∑ k=0 DkI(n−β)f(0+)sn−1−k. (9) Podemos escrever a equação (9) também na forma L [RLD βf(t)] = sβF (s)− n−1 ∑ k=0 sn−1−kg(k)(0+), (10) com F (s) = L [f(t)] e g(0+) = In−βf(0+). A obtenção desse resultado pode ser visto em (Camargo & Oliveira, 2015). 2.3 Derivada Fracionária de Caputo Caputo (1967) em seu artigo reformulou a definição da derivada fra- cionária de Riemann-Liouville, trocando a ordem da derivada ordinária com o ope- rador integral fracionário. Ao fazer isso, a transformada de Laplace dessa nova de- rivada depende das condições iniciais de ordem inteira, diferentemente da derivada fracionária Riemann-Liouville, que envolve condições iniciais de ordem não inteira (10) (Almeida, 2017). Motivado por este conceito, apresentamos a seguinte definição (Fernandez-Anaya et al., 2017; Xiao et al., 2015). Definição 3 Sejam β ∈ C com Re(β) > 0, e n o menor inteiro maior que ou igual a Re(β), isto é, n − 1 < Re(β) ≤ n. A derivada fracionária de Caputo de ordem β de f(t), para t > 0, é definida como CD βf(t) = In−β[Dnf(t)] = φn−β(t) ∗ [Dnf(t)]. (11) Assim como na derivada fracionária de Riemann-Liouville, quando β = n a derivada fracionária de Caputo corresponde a derivada usual, CD βf(t) = In−β[Dnf(t)] = In−n[Dnf(t)] = Dnf(t). 14 Proposição 4 A derivada fracionária de Caputo de ordem β ∈ C, com Re(β) > 0, da função f(t) = tµ com µ > −1, é da forma CD βtµ = Γ(µ+ 1) Γ(µ− β + 1) tµ−β . (12) Pela definição da derivada fracionária de Caputo na equação (11), te- mos CD βf(t) = Iν [Dnf(t)], como f(t) = tµ, CD βtµ = Iν [Dntµ]. Calculando a derivada de ordem n da função tµ, temos Dntµ = Γ(µ+ 1) Γ(µ− n+ 1) tµ−n. Assim, segue que CD βtµ = Iν [Dntµ] = Iν [ Γ(µ+ 1) Γ(µ− n+ 1) tµ−n ] = Γ(µ+ 1) Γ(µ− n + 1) Iν [tµ−n]. Pela equação (5), Iνtµ = Γ(µ+ 1) Γ(µ+ ν + 1) tµ+ν . Voltando à equação acima, CD βtµ = Γ(µ+ 1) Γ(µ− n + 1) Γ(µ− n+ 1) Γ(µ− n+ ν + 1) tµ−n+ν = Γ(µ+ 1) Γ(µ− n + ν + 1) tµ−n+ν . Realizando a substituição ν = n − β na equação acima, temos que a derivada fra- cionária de Caputo da função tµ é CD βtµ = Γ(µ+ 1) Γ(µ− β + 1) tµ−β . 15 Proposição 5 A derivada fracionária de Caputo de ordem β ∈ C, com Re(β) > 0, da função f(t) = c, é da forma CD βc = 0. Pela definição da derivada fracionária de Caputo na equação (11), te- mos CD βc = In−β[Dnc] = In−β[0]. Podemos observar que a derivada fracionária de Caputo aplicada em uma função constante é sempre zero, o que não acontece que a derivada fracionária de Riemann-Liouville, como vimos na Proposição 3. Proposição 6 A função de Mittag-Lefller pode ser considerada uma generalização “fracionária” da função exponencial. Pela definição da função de Mittag-Leffler, temos que Eβ(t β) = ∞ ∑ n=0 (tβ)n Γ(βn+ 1) , = 1 + tβ Γ(β + 1) + t2β Γ(2β + 1) + t3β Γ(3β + 1) + . . . . Aplicando a derivada fracionária de Caputo na função de Mittag-Lefller e utilizando o resultado da Proposição 4, temos Dβ [ Eβ(t β) ] = Dβ [ 1 + tβ Γ(β + 1) + t2β Γ(2β + 1) + t3β Γ(3β + 1) + . . . ] ; = 0 + 1 Γ(β + 1) Γ(β + 1) Γ(1) t0 + 1 Γ(2β + 1) Γ(2β + 1) Γ(β + 1) tβ + . . . , = 1 + tβ Γ(β + 1) + t2β Γ(2β + 1) + t3β Γ(3β + 1) + . . . , = Eβ(t β). Portanto, a derivada fracionária de Caputo da função de Mittag-Leffler resulta na própria função. Por este motivo, tal função é considerada uma generalização para a função exponencial. Assim, as soluções de equações diferenciais fracionárias com coeficientes constantes serão dadas em termos das funções de Mittag-Leffler. 16 Transformada de Laplace da Derivada Fracionária de Caputo A transformada de Laplace da derivada fracionária de Caputo é L [CD βf(t)] = sβF (s)− n−1 ∑ k=0 sβ−1−kf (k)(0+). (13) com F (s) = L [f(t)]. A obtenção desse resultado pode ser visto em (Camargo & Oliveira, 2015). Notemos que a transformada de Laplace da derivada fracionária de Riemann-Liouville necessita das condições iniciais dadas em termos da integral In−β e de suas derivadas de ordem k = 1, 2, . . . , n−1. Por outro lado, a derivada fracionária segundo Caputo parece ser mais conveniente quando utilizamos a transformada inte- gral, uma vez que requer o conhecimento de condições iniciais dadas na função e em suas derivadas de ordem inteira, que são fisicamente interpretáveis (Camargo, 2009). 2.4 Derivada Fracionária de Grünwald-Letnikov A derivada fracionária de Grünwald-Letnikov é uma extensão da de- rivada usual, que foi introduzido por Anton Karl Grünwald em 1867, e depois por Aleksey Vasilievich Letnikov em 1868 (Diethelm, 2010; Matlob & Jamali, 2017). Temos que a derivada de primeira ordem de f(t) é definida como: D1f(t) = df(t) dt = lim h→0 f(t)− f(t− h) h . (14) Aplicando novamente a derivada, D2f(t) = d2f(t) dt2 = lim h→0 f(t)− 2f(t− h) + f(t− 2h) h . (15) Pelo prinćıpio de indução, Dnf(t) = dnf(t) dtn = lim h→0 1 hn n ∑ k=0 (−1)k ( n k ) f(t− kh), (16) sendo ( n k ) = n(n− 1)(n− 2) . . . (n− k + 1) k! . 17 Expandindo a expressão (16) válida para todos os n ∈ N, obtemos assim a importante equação fracionária conhecida como derivada fracionária de Grünwald-Letnikov, com β ∈ R, dada por GLD βf(t) = dβf(t) dtβ = lim h→0 1 hβ ∞ ∑ k=0 (−1)k ( β k ) f(t− kh). (17) Com n = t h , obtemos GLD βf(t) = lim h→0 h−β Γ(−β) n ∑ k=0 Γ(k − β) k + 1 f(t− kh). (18) 3 CONTRIBUIÇÕES TEÓRICAS Esta seção foi elaborada com o objetivo de mostrar algumas contri- buições teóricas adquiridas neste trabalho quando analisávamos a viabilidade de alguns métodos numéricos que seriam utilizados para resolver o problema espećıfico do Caṕıtulo 5, bem como a mais adequada forma de tratar a modelagem fracionária sem erro dimensional. Para a exibição dos gráficos deste Caṕıtulo, foi utilizado o software Wolfram Mathematica 11.2, uma das grandes vantagens deste programa são as funções de Mittag-Leffler inseridas no seu código. O processador da máquina utili- zada foi da Intel(R) CoreTM i7-7500U CPU @ 2.70GHz 2.90GHz. 3.1 Análise Dimensional na Modelagem Fracionária Na literatura, há várias definições para as derivadas fracionárias, como por exemplo, as representações de Riemann-Liouville, Riesz, Weyl, Grünwald- Letnikov, Caputo, dentre outras (Gómez-Aguilar et al., 2015). Uma análise minu- ciosa dos sistemas dinâmicos fracionários é necessária para alcançar uma definição apropriada da derivada fracionária. Por exemplo, a definição de Riemann-Liouville envolve condições iniciais fisicamente inaceitáveis (condições iniciais de ordem fra- cionária) (Diethelm, 2010); ao contrário da representação de Caputo, as condições iniciais são expressas em termos de derivadas de ordem inteira com significado f́ısico (Diethelm et al., 2005). Esta definição é usada, principalmente, para incluir efeitos de memória (Gómez-Aguilar et al., 2015). Diante de tantos cuidados referentes a uti- lização do cálculo fracionário, nesta seção será realizada a modelagem fracionária de 19 modelos cujas soluções foram exibidas sem considerar a correta dimensão da equação quando há a troca da ordem da derivada inteira pela não inteira. A modelagem fracionária realizada em diversos trabalhos consiste em substituir a ordem da derivada inteira pela derivada fracionária, sem considerar a questão dimensional dos modelos, como por exemplo, os trabalhos de Mainardi (1996); El-Sayed et al. (2007); Camargo (2009); Mainardi (2010); Grigoletto & Oli- veira (2013); Camargo & Oliveira (2015); Kuroda (2016); Kuroda et al. (2017b); D’Ovidio et al. (2018); Veeresha et al. (2019b); Oliveira et al. (2019); Sólis-Péres et al. (2019). Pode-se ver, no decorrer deste caṕıtulo, que existem diferentes maneiras de se realizar a modelagem fracionária tratando a questão dimensional, coincidindo com os resultados obtidos em modelos cuja dimensão não foi considerada, havendo apenas alterações nos valores dos parâmetros do modelo. Note que o operador diferencial, d dt tem dimensão3 t−1 e quando rea- lizada a modelagem fracionária d dt → dβ dtβ sua dimensão passa a ser [ dβ dtβ ] = t−β, o que pode tornar imprecisa a descrição do modelo, caso as constantes envolvidas já tenham sua unidade de medida fixa. Assim, Gómez-Aguilar et al. (2012) propõem um procedimento alter- nativo simples para a construção da equação diferencial fracionária, introduzindo um novo parâmetro τ com 0 < β ≤ 1 da seguinte maneira, d dt → 1 τ 1−β dβ dtβ . (19) Desta forma, a equação (19) é dimensionalmente consistente, se e so- mente se, o novo parâmetro τ , tem dimensão de tempo [τ ] = t. Assim, [ 1 τ 1−β dβ dtβ ] passa a ser uma derivada no tempo no sentido usual, cuja dimensão é t−1. Assim, para n − 1 < β ≤ n com n ∈ N, temos a seguinte modelagem fracionária com correção dimensional, 3t simboliza o tempo referente a unidade de medida utilizada em modelos, podendo ser segundos, dias, horas, dentre outros. 20 dn dtn → 1 τn−β dβ dtβ . (20) Para analisar o comportamento dinâmico de um sistema fracionário é necessário usar uma definição apropriada de derivada fracionária. De fato, a definição da derivada de ordem fracionária não é única como já foi visto. No caso de Caputo, a derivada de uma constante é zero e a transformada de Laplace da derivada de Caputo depende de condições iniciais dadas em termos das derivadas de ordem inteira da função, que são fisicamente interpretáveis, o mesmo não ocorre com a derivada de Riemann-Liouville. Por estes motivos, a derivada de Caputo será utilizada nos modelos a seguir, levando-se em conta a correção dimensional dada em (20) (Gómez- Aguilar et al., 2012). 3.1.1 Modelo Malthusiano A forma mais simples de representar o processo de dinâmica populaci- onal é pelo modelo proposto por Thomas Robert Malthus, no qual é assumido que a probabilidade da população se reproduzir ou morrer permanece constante (Seti et al., 1999). dN(t) dt = (γ − δ)N(t), (21) sendo N(t) a população de indiv́ıduos no instante t, γ taxa de natalidade e δ a taxa de mortalidade. Em nossas aplicações utilizamos modelos matemáticos a fim de des- crever o crescimento de células tumorais, neste caso, as constantes utilizadas nos modelos de Malthus, loǵıstico e Gompertz serão positivas, ou seja, a taxa de natali- dade será maior que a taxa de mortalidade. Assim, assumindo r = (γ − δ) > 0, dN(t) dt = rN(t). (22) A solução do modelo (22), obtida via separação de variáveis, é dada 21 por N(t) = ertN(0), (23) sendo N(0) o tamanho inicial do tumor. A partir do modelo de crescimento tumoral de Malthus na equação (22), iremos utilizar a derivada fracionária segundo Caputo nesta equação, sendo 0 < β ≤ 1 a ordem da derivada (consequentemente n = 1) e τ o parâmetro de correção referente a dimensão da equação. Assim, 1 τ 1−β dβN(t) dtβ = 1 τ 1−β [ CD βN(t) ] = rN(t). (24) Sendo F (s) a transformada de Laplace de N(t), temos L [ 1 τ 1−β [ CD βN(t) ] ] = L [rN(t)] 1 τ 1−β [ sβF (s)− sβ−1N(0) ] = rF (s) sβF (s)− sβ−1N(0) = rτ 1−βF (s) F (s) = N(0) sβ−1 sβ − rτ 1−β . Utilizando a transformada de Laplace inversa e os resultados já obtidos das transformadas de Laplace das funções especiais, temos N(t) = L −1[F (s)] = L −1 [ N(0) sβ−1 sβ − rτ 1−β ] = N(0)Eβ(rτ 1−βtβ). Note que a solução anaĺıtica do modelo de Malthus (22) depende de dois parâmetros, β e τ . Assim Gómez-Aguilar et al. (2012, 2013, 2015), Tavoni (2019) e Cardoso (2019) propõem a seguinte substituição τ = β r , com o intuito da 22 solução anaĺıtica do modelo apresentar apenas um parâmetro, a ordem da derivada fracionária (β). Como 0 < β ≤ 1, temos que 0 < τ ≤ 1 r . Dáı, N(t) = N(0)Eβ(rτ 1−βtβ) = N(0)Eβ( r1−β r−β τ 1−βtβ) = N(0)Eβ(r 1−βτ 1−β tβ r−β ) = N(0)Eβ((rτ) 1−βrβtβ) = N(0)Eβ(β 1−β(rt)β). Assim, uma solução fracionária alternativa, com dimensão corrigida do modelo de Malthus (24) dependendo apenas do valor da derivada β é N(t) = N(0)Eβ(β 1−β(rt)β). (25) Quando β = 1, recuperamos a solução de ordem inteira do modelo, ou seja, N(t) = N(0)Eβ(β 1−β(rt)β) = N(0)ert. Note que o comportamento descrito na Figura4 1 é controverso com trabalhos vistos referente ao modelo fracionário de Malthus, o que seria, até então, um comportamento inesperado (menor ordem da derivada, implica no maior o cres- cimento da curva) (Camargo & Oliveira, 2015; Kuroda, 2016; Kuroda et al., 2017b), se torna condizente com as pesquisas feitas com a maior parte dos modelos do cálculo fracionário, ou seja, quanto menor a ordem de derivada, menor será o crescimento. Veja que para β = 0, 9 ainda está com crescimento mais elevado do que a ordem inteira β = 1, 0, mas com o passar do tempo a curva de ordem inteira tende a crescer 4Vale destacar que todos gráficos foram criados pelo autor, com exceção da Figura 14 retirado de (Loizides et al., 2015). Maiores informações sobre os programas e métodos utilizados estão descritos na Seção 5.5. 23 0 1 2 3 t 10 20 30 N(t) �=1.0 �=0.9 �=0.8 �=0.7 �=0.6 Figura 1: Solução da equação fracionária de Malthus (25) para vários valores de β. mais rápido (Figura 2). Podemos notar também, que nas proximidades do instante inicial, todas as curvas de β < 1 tem crescimento mais elevado (Figura 2). 0.0 0.5 1.0 1.5 t 0.5 1.0 1.5 2.0 2.5 3.0 3.5 N(t) β=1.0 β=0.9 β=0.8 β=0.7 β=0.6 0 2 4 6 8 10 12 14 t 5.0×106 1.0×107 1.5×107 N(t) �=1.0 �=0.9 �=0.8 �=0.7 �=0.6 Figura 2: Solução da equação fracionária de Malthus (25) para 0 ≤ t ≤ 2 e 0 ≤ t ≤ 15, respectivamente. Podemos escrever a equação de Malthus (24) de tal maneira que o parâmetro de correção dimensional esteja inclúıdo implicitamente nas constantes do modelo, ou seja, 24 1 τ 1−β [ CD βN(t) ] = rN(t) ⇔ CD βN(t) = β1−βrβN(t). (26) Vejamos, L [ CD βN(t) ] = L [ β1−βrβN(t) ] sβF (s)− sβ−1N(0) = β1−βrβF (s) (sβ − β1−βrβ)F (s) = sβ−1N(0) F (s) = sβ−1N(0) sβ − β1−βrβ L −1[F (s)] = N(0)L −1 [ sβ−1 sβ − β1−βrβ ] N(t) = N(0)Eβ(β 1−βrβtβ) N(t) = N(0)Eβ(β 1−β(rt)β). Teorema 2 Seja 1 τn−β [ CD βN(t) ] = ±anN(t), com n ∈ N, tal que n− 1 < β ≤ n e a ∈ N ∗, sujeito às condições iniciais, N(0) = N0 e N (m−1)(0) = 0, com m = 2, ..., n. A solução da equação fracionária é da forma N(t) = N0Enβ ( ±βn(1−β)(at)nβ ) , com 0 < τ ≤ 1 a e 0 < β ≤ 1. De fato, aplicando a transformada de Laplace, 25 1 τn−β [ CD βN(t) ] = ±anN(t) 1 τn(1−β) [ CD nβN(t) ] = ±anN(t), 0 < β ≤ 1, 1 τn(1−β) L [ CD nβN(t) ] = ±L [anN(t)] 1 τn(1−β) [ snβF (s)− snβ−1N(0) ] = ±anF (s) snβF (s)− snβ−1N0 = ±τn(1−β)anF (s) [ snβ ∓ anτn(1−β) ] F (s) = snβ−1N0 F (s) = snβ−1N0 snβ ∓ anτn(1−β) A partir da transformada de Laplace inversa, L −1[F (s)] = N0L −1 [ snβ−1 snβ ∓ anτn(1−β) ] N(t) = N0Enβ(±anτn(1−β)tnβ), β = aτ, 0 < τ ≤ 1 a N(t) = N0Enβ ( ± an(1−β) a−nβ τn(1−β)tnβ ) N(t) = N0Enβ ( ±(aτ)n(1−β)anβtnβ ) N(t) = N0Enβ ( ±βn(1−β)(at)nβ ) . Consequentemente, podemos reescrever a equação da seguinte maneira, 1 τn−β [ CD βN(t) ] = ±anN(t) ⇔ CD nβN(t) = ±βn(1−β)anβN(t). (27) Veremos que o Teorema 2 pode ser utilizado no modelo do oscilador harmônico a seguir. 26 3.1.2 Oscilador Harmônico Segundo Aguiar & Guedes (2013), dentre todos os sistemas f́ısicos ca- racterizados por algum tipo de oscilação, o oscilador harmônico simples (OHS) é, sem dúvida, o caso mais estudado, pois em primeira aproximação, considerando pequenas amplitudes de oscilação, todo sistema oscilatório pode ser aproximado por um OHS, o qual, segundo evidências históricas, é estudado desde o século XVIII (Bassalo & Cattani, 2009). Provavelmente, o f́ısico e matemático súıço Leonhard Euler foi o primeiro cientista a estudar o OHS, quando resolveu a equação diferencial de movi- mento, utilizando um método numérico conhecido como Método das Quadraturas, em um trabalho estritamente matemático. Após os estudos puramente matemáticos de Euler, o OHS se tornou um dos maiores objetos de estudo dos pontos de vista f́ısico e matemático, embora saibamos que todo sistema f́ısico real tem dissipação (sistemas não conservativos), como por exemplo, o oscilador harmônico amortecido (OHA) (Aguiar & Guedes, 2013). Seja ω2 constante positiva igual a k m (constante elástica sobre a massa) e considerando que não há atritos nem forças externas atuando sobre o sistema, temos (Rodrigues & Oliveira, 2015; Gómez-Aguilar et al., 2012; Kuroda et al., 2015) d2x(t) dt2 + ω2x(t) = 0. (28) Temos que a solução para este modelo, sendo x(0) a posição inicial, é dada por (Gómez-Aguilar et al., 2012) x(t) = x(0) cos(ωt). (29) A partir da equação (28), iremos utilizar a derivada fracionária de Caputo para 1 < β ≤ 2 e τ o parâmetro de correção referente a dimensão da equação. Desta forma, τ tem dimensão de tempo [τ ] = t e [ dβ dtβ ] = t−β . Assim, [ 1 τ 2−β dβ dtβ ] passa a ser uma derivada no tempo no sentido usual, cuja dimensão é t−2. Podemos então realizar a seguinte modelagem fracionária, 27 d2 dt2 → 1 τ 2−β dβ dtβ . (30) Logo, pela equação (20), a equação do oscilador harmônico fracionário será, 1 τ 2−β dβx(t) dtβ + ω2x(t) = 0 ⇔ 1 τ 2−β [ CD βx(t) ] + ω2x(t) = 0. (31) Considerando 0 < β ≤ 1, podemos reescrever a equação (31) da se- guinte maneira, 1 τ 2(1−β) [ CD 2βx(t) ] + ω2x(t) = 0. (32) Sujeito as condições iniciais x(0) = x0 e x′(0) = 0, pelo Teorema 2, temos que a solução da equação (32) é: x(t) = x0E2β ( −β2(1−β)(ωt)2β ) . (33) Consequentemente, podemos reescrever a equação do oscilador harmônico fracionário da seguinte maneira, sem erro dimensional, 1 τ 2(1−β) [ CD 2βx(t) ] = −ω2x(t) ⇔ CD 2βx(t) = −β2(1−β)ω2βx(t). (34) Quando β = 1, recuperamos a solução de ordem inteira do modelo (Gómez-Aguilar et al., 2012), ou seja, x(t) = x0E2β ( −β2(1−β)(ωt)2β ) = x0 cos(ωt). A solução completa da equação (33) é apresentada no APÊNDICE B. Portanto, percebe-se que ao utilizarmos a derivada fracionária na equação do oscilador harmônico simples (sem atritos ou forças externas) obtemos uma equação que, de acordo com a ordem da derivada, retratará alterações no amor- tecimento e frequência do oscilador, permitindo descrever este fenômeno de maneira mais realista (Kuroda et al., 2015). 28 1 2 3 4 5 t -1.0 -0.5 0.5 x(t) β=1.0 β=0.95 β=0.9 β=0.85 β=0.8 2 4 6 8 10 12 14 t -1.0 -0.5 0.5 x(t) β=1.0 β=0.95 β=0.9 β=0.85 β=0.8 Figura 3: Solução da equação fracionária do oscilador harmônico (32) até o instante t = 6 s e t = 15 s. 3.1.3 Equação Loǵıstica Considere a equação loǵıstica obtida por P. F. Verhulst (Seti et al., 1999): dN(t) dt = rN(t) [ 1− N(t) K ] , (35) sendo N(t) o número de indiv́ıduos de uma população, K a capacidade de suporte e r > 0 a taxa de crescimento. A solução para esta equação quando K = 1 (sem perda de generali- dade) é da forma N(t) = 1 1 + [ 1 N(0) − 1 ] e−rt · (36) Antes de realizarmos a modelagem fracionária na equação loǵıstica, va- mos considerar a seguinte mudança de variável v(t) = 1 N(t) ⇒ v′(t) = −N−2(t)N ′(t), uma vez que não é posśıvel encontrar a solução anaĺıtica, via transformada de La- place, em modelos não lineares. Multiplicando a equação (35) por −N−2(t) temos, −N−2(t)N ′(t) = r [ 1− 1 N(t) ] , 29 ou seja, dv(t) dt = r[1− v(t)]. (37) Realizando a modelagem fracionária (20) segundo Caputo, com 0 < β ≤ 1, na equação (37), 1 τ 1−β [ CD βv(t) ] = r[1− v(t)]. (38) Sendo V (s) a transformada de Laplace de v(t), ou seja V (s) = L [v(t)], temos L [ 1 τ 1−β [ CD βv(t) ] ] = L [r[1− v(t)]] L [CD βv(t)] = rτ 1−β L [1]− rτ 1−β L [v(t)] sβV (s)− sβ−1v(0) = rτ 1−β 1 s − rτ 1−βV (s) (sβ + rτ 1−β)V (s) = rτ 1−βs−1 + v(0)sβ−1 V (s) = rτ 1−β [ s−1 sβ + rτ 1−β ] + v(0) [ sβ−1 sβ + rτ 1−β ] L −1[V (s)] = rτ 1−β L −1 [ s−1 sβ + rτ 1−β ] + v(0)L −1 [ sβ−1 sβ + rτ 1−β ] . Sabe-se que, vide Apêndice A, L [Eα(λt α)] = sα−1 sα − λ , (39) e L [tβ−1Eα,β(λt α)] = sα−β sα − λ . (40) Assim, v(t) = rτ 1−βtβEβ,β+1(−rτ 1−βtβ) + v(0)Eβ(−rτ 1−βtβ). Utilizando a relação, Eα,α+1(−zα) = 1−Eα(−zα) zα . (41) Temos que, 30 v(t) = rτ 1−βtβ 1− Eβ(−rτ 1−βtβ) rτ 1−βtβ + v(0)Eβ(−rτ 1−βtβ). v(t) = 1−Eβ(−rτ 1−βtβ) + v(0)Eβ(−rτ 1−βtβ) v(t) = 1 + Eβ(−rτ 1−βtβ)[v(0)− 1]. Como v(t) = 1/N(t), a solução fracionária da equação loǵıstica é dada por N(t) = 1 1 + [ 1 N(0) − 1 ] Eβ(−rτ 1−βtβ) , β = rτ, 0 < τ ≤ 1 r . (42) N(t) = 1 1 + [ 1 N(0) − 1 ] Eβ(−β1−β(rt)β) · (43) Consequentemente, podemos reescrever a equação loǵıstica fracionária da seguinte maneira sem erro dimensional. 1 τ 1−β [ CD βN(t) ] = rN(t)[1−N(t)] ⇔ CD βN(t) = β1−βrβN(t)[1−N(t)]. (44) É importante mencionar que a equação (43) não representa a solução da equação loǵıstica fracionária, 1 τ 1−β [ CD βN(t) ] = rN(t) [1−N(t)] , e sim a equação loǵıstica linearizada, equação (38). Embora sejam similares, elas não coincidem sempre, uma vez que em geral não é válida a regra da cadeia para a derivada fracionária (Cardoso & Camargo, 2015). Por outro lado, observa-se que lim β→1 N(t) = 1 1 + [ 1 N(0) − 1 ] Eβ(−β1−β(rt)β) = 1 1 + [ 1 N(0) − 1 ] e−rt , recuperando a solução da equação loǵıstica de ordem inteira. Note pela Figura 4 quando β = 1, isto é, derivada de ordem inteira, temos a população convergindo rapidamente para seu tamanho máximo, no entanto, quando consideramos a derivada fracionária no modelo, ao diminuir a ordem da derivada o crescimento populacional converge para a capacidade de suporte mais lentamente. Com o passar do tempo, observa-se que a as populações tendem a atingir a capacidade de suporte. 31 2 4 6 8 t 0.4 0.6 0.8 N(t) β=1 β=0.9 β=0.8 β=0.7 β=0.6 10 20 30 40 t 0.4 0.6 0.8 N(t) β=1 β=0.9 β=0.8 β=0.7 β=0.6 Figura 4: Solução da equação loǵıstica fracionária (43) para 0 ≤ t ≤ 10 e 0 ≤ t ≤ 50, respectivamente. 3.1.4 Modelo de Gompertz A equação diferencial de Gompertz é dada por (Tjorve & Tjorve, 2017; Sólis-Péres et al., 2019), dN(t) dt = rN(t) ln ( k N(t) ) , (45) sendo N(t) o tamanho populacional, r > 0 a taxa intŕınseca de crescimento e k a capacidade de suporte. Sujeito a condição inicial N(0) = N0, a solução para o modelo de Gompertz é dada por N(t) = ke − ln ( k N0 ) e−rt . (46) Assim como a equação loǵıstica, o modelo de Gompertz é não linear. Logo, não é posśıvel encontrar a solução via transformada de Laplace quando reali- zada a modelagem fracionária. Seja G(t) = ln(k)− ln(N), temos dG(t) dt = −rG(t), (47) com G(0) = ln ( k N0 ) . 32 Realizando a modelagem fracionária (20), com 0 < β ≤ 1, na equação (47), 1 τ 1−β [ CD βG(t) ] = −rG(t), (48) com G(0) = ln ( k N0 ) . Usando o Teorema 2, a solução do modelo (48) é da forma G(t) = G(0)Eβ(−β1−β(rt)β). (49) Fazendo a mudança de variável G(t) = ln(k)− ln(N), segue que N(t) = ke − ln ( k N0 ) Eβ(−β1−β(rt)β) . (50) Note que quando β = 1, recuperamos a solução anaĺıtica clássica do modelo (46). Consequentemente, podemos reescrever o modelo de Gompertz fra- cionário da seguinte maneira sem erro de dimensão. 1 τ 1−β [ CD βN(t) ] = rN(t) ln ( k N(t) ) ↔ CD βN(t) = β1−βrβN(t) ln ( k N(t) ) . (51) 2 4 6 8 t 0.4 0.6 0.8 N(t) β=1 β=0.9 β=0.8 β=0.7 β=0.6 10 20 30 40 t 0.70 0.75 0.80 0.85 0.90 0.95 N(t) β=1 β=0.9 β=0.8 β=0.7 β=0.6 Figura 5: Solução do modelo de Gompertz fracionária (50) para 0 ≤ t ≤ 10 e 0 ≤ t ≤ 50, respectivamente. 33 O comportamento do modelo de Gompertz, Figura 5, se assemelha com a solução obtida por Malthus, Figura 4, quando β = 1 a população converge rapidamente para seu tamanho máximo, no entanto, quando consideramos a derivada fracionária no modelo, ao diminuir a ordem da derivada o crescimento populacional converge para a capacidade de suporte mais lentamente. Assim como no modelo loǵıstico, é importante mencionar que a equação (50) não representa a solução do modelo de Gompertz fracionária, 1 τ 1−β [ CD βN(t) ] = rN(t) ln ( k N(t) ) , e sim o modelo de Gompertz linearizado, equação (48). Embora sejam similares, elas não coincidem sempre, uma vez que em geral não é válida a regra da cadeia para a derivada fracionária (Cardoso & Camargo, 2015). A vantagem de linearizar é poder encontrar um posśıvel valor de τ para utilizar a modelagem fracionária de alguns modelos que não conseguimos encontrar a solução via transformada de Laplace. 3.1.5 A Escolha do Parâmetro “τ” Na seção anterior realizamos a modelagem fracionária de Caputo, com o parâmetro τ em quatro modelos: modelo de Malthus (22), oscilador harmônico (28), modelo loǵıstico (35) e no modelo de Gompertz (45). Em cada modelo foi sugerido um valor de τ , em função de β, a partir do momento em que encontramos a solução da equação via transformada de Laplace. No modelo loǵıstico e de Gompertz, foi sugerida a linearização das equações, para assim encontrar uma solução da equação linearizada e um posśıvel valor para τ . Note que em alguns casos, o desafio passa ser em encontrar o valor do parâmetro τ . No entanto, esse novo parâmetro na modelagem fracionária nos abre novos horizontes sobre diferentes modelagens sem erros de dimensão se tratando do mesmo modelo. Por exemplo, referente ao modelo malthusiano (24), 34 1 τ 1−β [ CD βN(t) ] = rN(t), (52) quando tomado 0 < τ ≤ 1 r podemos reescrever o modelo da seguinte maneira, como já foi visto: CD βN(t) = β1−βrβN(t). Note que a nova modelagem vai depender apenas do parâmetro β uma vez que tomamos β = rτ . Agora, seja τ = 1 r (como a dimensão de [τ ] = t, assume-se que τ = 1 r , uma vez que [r] tem dimensão t−1) no modelo (52), podemos reescrever esta equação da seguinte maneira, 1 τ 1−β [ CD βN(t) ] = rN(t) CD βN(t) = τ 1−βrN(t) CD βN(t) = ( 1 r )1−β rN(t) CD βN(t) = ( r r1−β ) N(t) CD βN(t) = ( r r.r−β ) N(t) CD βN(t) = rβN(t). Esta modelagem foi utilizada para analisar o método da transformada diferencial generalizada com multi-passos no modelo malthusiano (64), como pode ser visto na Seção 3.2. Para τ = 1 r , temos que 1 τ 1−β [ CD βN(t) ] = rN(t) ⇔C DβN(t) = rβN(t). (53) Note pela Figura 6, que ambas soluções de (53) coincidem. Além disso, o comportamento inesperado visto em Kuroda et al. (2017b) é descrito, ou seja, 35 0 1 2 3 t 10 20 30 40 50 N(t) β=1.0 β=0��� (r0.9) β=0��� (r0.8) β=0���(r�� ) β=0� � (r���) β=0��� (τ=r-1) β=0��� (τ=r-1) β=0��� (τ=r-1) β=0� � (τ=r-1) Figura 6: Solução das diferentes modelagens no modelo Malthusiano (53). Sendo as linhas cont́ınuas referentes ao parâmetro rβ e as tracejadas quando assumimos que τ = 1 r . quanto menor a ordem da derivada maior será o crescimento populacional. Situação contrária quando tomamos β = rτ , Figura 2. A modelagem feita anteriormente, equação (53), é a mesma utilizada por Kuroda et al. (2019) para analisar o método MSGDTM e na proposta por Ro- drigues & Oliveira (2015) na equação do oscilador harmônico (31), quando substitui ω2 por ωβ, com 1 < β ≤ 2 e τ = 1 ω . Vejamos, 1 τ 2−β [ CD βx(t) ] = ω2x(t) CD βx(t) = τ 2−βω2x(t) CD βx(t) = ( 1 ω )2−β ω2x(t) CD βx(t) = ωβx(t). Nota-se a extrema importância desta nova modelagem fracionária de- pendendo do valor de τ (imposta para a correção dimensional, não analisada até o momento), que, apesar de ser mais um parâmetro introduzido no modelo, podemos descrever diversas soluções e aproximar essas curvas ao dados reais (como será visto 36 no Caṕıtulo 5, onde apresentamos as aplicações). A grande dificuldade está em es- colher qual valor de τ será utilizado no modelo. Em alguns casos se faz necessário escolher aquele em que a curva se ajusta melhor aos dados reais. Para modelos mais complexos se faz necessário a utilização de métodos computacionais para encontrar a solução fracionária, como por exemplo, o método de Grünwald-Letnikov descrito no próximo Caṕıtulo. Podemos fazer algumas sugestões referentes ao valor de τ , Tabela 1, de tal maneira que esta equação dependa apenas do valor da ordem da derivada fracionária. Assim, uma maneira de se corrigir o erro dimensional é elevar a β nos parâmetros do modelo, como visto em (53). A Tabela 1 contém os modelos de Malthus (M), oscilador harmônico (OH), loǵıstico (L), Gompertz (G) e o modelo de von Bertalanffy (V) juntamente com algumas sugestões para o valor de τ em cada equação e sua respectiva modelagem para o valor do parâmetro escolhido (τ). Encontra-se também alguns trabalhos em que o valor do parâmetro foi utilizado, mesmo que implicitamente5. Com o objetivo de encontrar soluções de equações diferenciais fra- cionárias via métodos computacionais, vale ressaltar a importância de se escolher qual o método utilizar quando se refere a modelagem fracionária. A seguir será des- crito o método MSGDTM, que se faz inaplicável por tratar a modelagem fracionária como um operador local. No Caṕıtulo referente as infecções por papilomav́ırus humano (HPV 16) administrados em camundongos, utilizaremos os três tipos de modelagem fra- cionária de Caputo vistos nessa seção, e exibiremos a melhor aproximação aos dados reais. 5Modelagem utilizada sem citar o valor do parâmetro τ . 37 Tabela 1: Modelagem para diferentes valores de τ . Modelo Equação τ Modelagem Referência M 1 τ1−β [ CDβN(t) ] = rN(t) τ = β r CDβN(t) = β1−βrβN(t) - M 1 τ1−β [ CDβN(t) ] = rN(t) τ = 1 r CDβN(t) = rβN(t) (West, 2015; Kuroda et al., 2019) OH 1 τ2−β [ CDβx(t) ] = −ω2x(t) τ = β ω CD2βx(t) = −β2(1−β)ω2βx(t) (Gómez-Aguilar et al., 2012) OH 1 τ2−β [ CDβx(t) ] = −ω2x(t) τ = 1 ω CDβx(t) = −ωβx(t) (Rodrigues & Oliveira, 2015) L 1 τ1−β [ CDβN(t) ] = rN(t)[1 − N(t)] τ = β r CDβN(t) = β1−βrβN(t)[1 − N(t)] - L 1 τ1−β [ CDβN(t) ] = rN(t)[1 − N(t)] τ = 1 r CDβN(t) = rβN(t)[1 − N(t)] (Ortigueira & Bengochea, 2017) G 1 τ1−β [ CDβN(t) ] = rN(t) ln ( k N(t) ) τ = β r CDβN(t) = β1−βrβN(t) ln ( k N(t) ) - G 1 τ1−β [ CDβN(t) ] = rN(t) ln ( k N(t) ) τ = 1 r CDβN(t) = rβN(t) ln ( k N(t) ) - V 1 τ1−β [ CDβV (t) ] = aV (t) 2 3 − bV (t) τ = 3β b - (Tavoni, 2019) 3.2 Análise do Método Multi-Passos com Transformada Di- ferencial Generalizada na Modelagem Fracionária Esta seção é baseada no artigo de Kuroda et al. (2019), em que é apre- sentada uma análise cŕıtica de uma técnica numérica que tem sido usada na resolução de equações diferenciais de ordem fracionária com derivadas de Caputo. Trata-se do método multi-passos com transformada diferencial generalizada. Verifica-se que a versão do método dispońıvel na literatura produz soluções equivocadas a partir do segundo passo, isto é mostrado em aplicações aos modelos de Malthus e de Riccati. O problema é explicado em termos da não localidade da derivada de Caputo e das propriedades da transformada diferencial generalizada. 38 Considerando uma equação diferencial que descreve um fenômeno, uma maneira comum de usar a modelagem fracionária é substituir as derivadas de ordem inteira por derivadas de ordem não inteira, geralmente com ordem menor que ou igual à ordem das derivadas originais, de modo que a solução do problema original possa ser recuperada como um caso particular (Camargo & Oliveira, 2015; Kuroda et al., 2017b). Assim, o principal desafio passa ser buscar técnicas de resolução das equações, por exemplo, se tratando de equações lineares de ordem fracionárias, podemos utilizar a metodologia da transformada de Laplace. Deve-se notar que há limitações na utilização dessa metodologia, para modelos são não lineares, que são resolvidas via métodos numéricos, em geral. Entre estes está o método multi-passos com transformada diferencial generalizada. O método MSGDTM foi utilizado para investigar um modelo de aban- dono do tabagismo, incluindo comparações com método de Runge-Kutta clássico no caso de derivadas de ordem inteira (Erturk et al., 2012). Aplicado também num modelo epidemiológico para v́ırus de computador de (Handam & Freihat, 2015), no estudo da co-infecção por HIV e malária (Ebenezer et al., 2015) e no tratamento de dinâmica tumoral (Arshad et al., 2015). Em relação com as soluções do sistema Chua, os autores destacaram que o método tem a vantagem de fornecer uma forma anaĺıtica aproximada da solução dentro de cada sub-intervalo de tempo (Freihat & Momani, 2012). É importante salientar que uma série de aplicações do método têm sido reportadas, antes da sua validade ter sido rigorosamente verificada. Neste sentido, Momenzadeh & Sarv Ahrabi (2017) apontaram problemas no MSGDTM. A questão é que existe grande risco ao estender a validade, à modelagem de ordem não inteira, dos métodos numéricos demonstrados apenas para equações diferenciais de ordem inteira. De um lado, o método da transformada diferencial pode ser generalizado para derivadas de ordem fracionária, do outro lado, o método multi-passos, válido para ordem inteira, mas não se estende de forma automática para o caso fracionário. A não localidade da derivada fracionária e as propriedades da transformada diferencial 39 generalizada devem ser consideradas, uma vez que a derivada fracionária de Caputo é definida em termos de uma integral (Camargo & Oliveira, 2015). Com o intuito de esclarecer a situação, o referido trabalho traz uma análise detalhada do MSGDTM. A análise utiliza uma série de resultados numéricos para os modelos de Malthus e de Riccati (Momenzadeh & Sarv Ahrabi, 2017). 3.2.1 Método da Transformada Diferencial: Extensão e Generalização O DTM (método da transformada diferencial) permite resolver uma grande variedade de equações diferenciais de ordem inteira Kanth & Aruna (2009); Odibat (2008); Soltanalizadeh (2011). A transformada diferencial converte cada função (derivável infinitas vezes) na sequência dos seus coeficientes de Taylor. O DTM converte a equação diferencial num sistema de equações algébricas para os termos da transformada. Então, a solução é dada por uma série de Taylor, cujo raio de convergência pode ser finito ou infinito. A avaliação numérica da série também pode apresentar instabilidade numérica, limitando mais o intervalo de resolução da equação. Uma maneira de ampliar o intervalo de resolução de equações diferenci- ais de ordem inteira, é aplicar o DTM e utilizar multi-passos, ou seja, uma sequência de intervalos. Isto define o método MSDTM (método multi-passos da transformada diferencial). Em cada passo, a condição inicial é dada pelos valores finais da solução do passo anterior. Além disso, o tamanho de cada intervalo deve ser menor que o raio de convergência da série de Taylor. Para resolver equações diferenciais fracionárias com derivada de Ca- puto de ordem β, o DTM é generalizado, dando lugar ao GDTM (método da trans- formada diferencial generalizada). A transformada diferencial generalizada de cada função é dada pela sequência dos coeficientes de uma série de potências da mesma, em que os expoentes são múltiplos inteiros não negativos da ordem da derivada β. Novamente, o intervalo de resolução pode ser limitado pelo raio de convergência da série ou por instabilidades numéricas. 40 Para estender o intervalo de resolução do GDTM, tem sido utilizado o MSGDTM Kuroda et al. (2017a); Kuroda (2016); Odibat et al. (2010). A ideia é a mesma do MSDTM, mas cada passo utiliza o GDTM. Como comentado antes, esta extensão deve ser verificada e analisada rigorosamente. A seguir, com o objetivo de contribuir para a modelagem fracionária, são descritos detalhadamente apenas o GDTM e o MSGDTM. 3.2.2 Método da Transformada Diferencial Generalizada Segundo Kuroda (2016); Odibat et al. (2010, 2008), a técnica da trans- formada diferencial é um dos métodos numérico-anaĺıticos para equações diferenci- ais ordinárias e parciais que utiliza a forma de polinômios como aproximações das soluções exatas que são suficientemente diferenciáveis. A transformada diferencial, F (k), da k-ésima derivada da função f(t) é definida como, F (k) = 1 k! [ dkf(t) dtk ] t=t0 . (54) Já a transformada inversa de F (t) é dada por Odibat (2008); Kuroda (2016); Odibat et al. (2010, 2008); Mirzaee (2011): f(t) = ∞ ∑ k=0 F (k)(t− t0) k. (55) A correspondência entre as equações (54) e (55) pode ser evidenciada ao aplicar o resultado da (54) na equação (55). De fato, seguindo este procedimento temos: f(t) = ∞ ∑ k=0 (t− t0) k k! [ dkf(t) dtk ] t=t0 . Assim, a transformada diferencial fornece os coeficientes da expansão em série de Taylor. No entanto, as correspondentes derivadas são calculadas ite- rativamente pelas equações transformadas da função original. Para fins de imple- mentação, f é expressa por uma série finita. Logo, a equação (55) pode ser escrita como (Mirzaee, 2011), 41 f(t) ≈ S ∑ k=0 F (k)(t− t0) k, em que S é suficientemente grande. Para resolver problemas não-lineares utili- zando o DTM, precisamos, primeiramente, aplicar a transformação diferencial (54) na equação de estudo, resultando em uma relação de recorrência, depois resolvendo esta equação utilizando a transformada diferencial inversa (55), obtemos a solução do problema. O método de transformada diferencial generalizada (GDTM) é utili- zado para exibir a solução numérica de equações diferenciais de ordem fracionária. O método proposto é baseado na fórmula generalizada de Taylor. Definimos a trans- formada diferencial generalizada da k-ésima derivada da função f(t) como (Odibat et al., 2010), F (k) = 1 Γ(βk + 1) [ (CD β)kf(t) ] t=t0 , (56) sendo 0 < β ≤ 1, (CD β)k a derivada CD β k-vezes e CD β a derivada de Caputo de ordem β como visto na sessão anterior, com ińıcio em t = t0 (Camargo & Oliveira, 2015; Caputo & Frabrizio, 2015; Caputo, 1967), CD βf(t) = 1 Γ(n− β) ∫ t t0 f (n)(τ)(t− τ)n−β−1 dτ, (57) com n sendo o menor dos inteiros maiores ou iguais que β. Seja CD β a derivada de Caputo (57), f(t) = (t − t0) α e u = τ − t0 t− t0 , temos 42 CDβ(t− t0) α = Γ(α+ 1) Γ(α− n+ 1)Γ(n − β) ∫ t t0 (τ − t0) α−n(t− τ)n−β−1dτ, CDβ(t− t0) α = Γ(α+ 1) Γ(α− n+ 1)Γ(n − β) ∫ 1 0 [u(t− t0)] α−n [−u(t− t0) + (t− t0)] n−β−1 (t − t0)du, CDβ(t− t0) α = Γ(α+ 1) Γ(α− n+ 1)Γ(n − β) ∫ 1 0 uα−n(t− t0) α−n(t − t0) n−β−1(t − t0)(1 − u)n−β−1du, CDβ(t− t0) α = Γ(α+ 1)(t − t0)α−β Γ(α− n+ 1)Γ(n − β) ∫ 1 0 uα−n(1− u)n−β−1du, CDβ(t− t0) α = Γ(α+ 1)(t − t0)α−β Γ(α− n+ 1)Γ(n − β) B(α − n+ 1, n− β), CDβ(t− t0) α = Γ(α+ 1)(t − t0)α−β Γ(α− n+ 1)Γ(n − β) Γ(α− n+ 1)Γ(n − β) Γ(α − n+ 1 + n− β) , CDβ(t− t0) α = Γ(α + 1) Γ(α− β + 1) (t− t0) α−β , desde que α > n− 1. Temos que a transformada inversa de (56) é dada por (Odibat et al., 2010), f(t) = ∞ ∑ k=0 F (k)(t− t0) βk. (58) Vejamos, (CDβ)0f(t) = f(t) = F (0) + F (1)(t − t0) β + F (2)(t − t0) 2β + . . .+ F (k)(t − t0) kβ + . . . . Quando t = t0, temos que (CD β)0f(t0) = F (0). Como visto em (3.2.2), (CDβ)1f(t) = F (1)Γ(β + 1) Γ(1) + F (2)Γ(2β + 1) Γ(β + 1) (t − t0) β + . . .+ F (k)Γ(kβ + 1) Γ((k − 1)β + 1) (t− t0) (k−1)β + . . . . Quando t = t0, segue que 43 (CD β)1f(t0) = F (1)Γ(β + 1). (CD β)2f(t) = F (2)Γ(2β + 1) Γ(β + 1) Γ(β + 1) 1 + F (3)Γ(3β + 1) Γ(2β + 1) Γ(2Γ + 1) Γ + 1 (t− t0) β + . . .+ = + F (k)Γ(kβ + 1) Γ((k − 1)β + 1) Γ((k − 1)β + 1) Γ((k − 2)β + 1) (t− t0) (k−1)β + . . . . Quando t = t0, (CD β)2f(t0) = F (2)Γ(2β + 1). Assim, (CD β)kf(t0) = F (k)Γ(kβ + 1). Como visto em (56), F (k) = [ (CD β)kf(t) Γ(kβ + 1) ] t=t0 . Substituindo (56) em (58) temos, f(t) = ∞ ∑ k=0 (t− t0) βk Γ(βk + 1) [ (CD β)kf(t) ] t=t0 . Na Tabela 2, apresentamos algumas propriedades da transformada di- ferencial generalizada. As duas primeiras exprimem a linearidade da transformada e a terceira expressa a transformada do produto como convolução. A quarta pro- priedade é fundamental para o GDTM, uma vez que fornece uma fórmula simples da transformada da derivada de Caputo. A demonstração dos resultados tabelados estão no APÊNDICE C ou podem ser vistos em (Odibat et al., 2010; Freihat & Momani, 2012; Momenzadeh & Sarv Ahrabi, 2017). Consideremos um problema de valor inicial de equação diferencial de ordem fracionária, β, em relação à x(t), da forma Dβ t0x(t) = f(t, x(t)), 44 Tabela 2: Propriedades da transformada diferencial generalizada. Considera-se que r é uma constante e que as transformadas de x(t) e y(t) são X(k) e Y (k), respecti- vamente, com ordem β e centro em t0. Função Original Transformada Diferencial Generalizada f(t) = r x(t) F (k) = r X(k) f(t) = x(t)± y(t) F (k) = X(k)± Y (k) f(t) = x(t)y(t) F (k) = k ∑ l=0 X(l)Y (k − l) f(t) = Dβ t0x(t) F (k) = Γ[(k + 1)β + 1] Γ(kβ + 1) X(k + 1) em que f(t, x) depende de t e x suavemente, não necessariamente de forma linear. A equação deve ser resolvida para t > t0, com x(t0) = x0. Ao aplicar a transformada à condição inicial temos que X(0) = x0, enquanto a aplicação à equação diferencial produz X(k + 1) = Γ(kβ + 1) Γ[(k + 1)β + 1] F (k,X), (59) em que F (k,X) é a transformada da composta f(t, x(t)). Aqui o śımbolo X repre- senta todos os valores da transformada de x(t). Quando f(t, x) = r x, a Tabela 2 fornece F (k,X) = r X(k), levando à relação de recorrência X(k + 1) = r Γ(kβ + 1) Γ[(k + 1)β + 1] X(k). (60) Um outro exemplo é f(t, x) = 2x− x2 + 1, donde F (k,X) = 2X(k)− k ∑ l=0 X(l)X(k − l) e X(k + 1) = Γ(kβ + 1) Γ[(k + 1)β + 1] [2X(k)− k ∑ l=0 X(l)X(k − l)]. (61) 45 3.2.3 Método Multi-passos com Transformada Diferencial Generalizada Nesta subseção é apresentado o MSGDTM, de acordo com Arshad et al. (2015). Para resolver o problema de valor inicial num intervalo [t0, T ], a ideia do método multi-passos é dividir o intervalo em M subintervalos da forma [tm−1, tm], m = 1, 2, ...,M , todos de tamanho h = T−t0 M . A solução tem a forma x(t) =                x1(t) se t ∈ [t0, t1], x2(t) se t ∈ [t1, t2], ... xM(t) se t ∈ [tM−1, tM ], (62) sendo x1(t) a solução obtida pelo GDTM no intervalo [t0, t1]. Para m ≥ 2, aplica-se o GDTM no intervalo [tm−1, tm], usando a condição inicial um(tm−1) = um−1(tm−1). Deve-se notar que Arshad et al. (2015) expressam um(t) como uma série generalizada de potências de t − tm−1, e não de t − t0. Ao mesmo tempo, a derivada de Caputo não é Dβ tm−1 , senão Dβ t0 . Acontece que, a derivada de Caputo leva em conta toda a memória do sistema, o que está relacionado com a sua caracteŕıstica não local. Isto faz com que a aplicação da equação (59) seja questionável para m ≥ 2. A seguir, veremos o que acontece ao aplicar o método MSGDTM em modelos mais simples. 3.2.4 Modelos Matemáticos Como visto anteriormente, o modelo de Malthus é utilizado na mo- delagem de crescimento populacional. Alguns autores utilizam este modelo para descrever o crescimento de células tumorais (Kuroda, 2016); Supõe-se que a taxa segundo qual a população cresce em um determinado instante é proporcional ao tamanho da população naquele instante, ou seja dN(t) dt = r N(t), (63) sendo N(t) a população no instante t e r > 0, a taxa de crescimento intŕınseca. 46 3.2.5 Modelagem Fracionária O modelo fracionário de Malthus é uma generalização da equação (63). No lugar da primeira derivada, usa-se a derivada de Caputo de ordem β com ińıcio em t = 0, sendo 0 < β ≤ 1. Portanto, vale (53) CD βN(t) = rβN(t). (64) Neste caso, a aplicação do GDTM leva à relação de recorrência na equação (60). Através da metodologia da transformada de Laplace, temos que a solução anaĺıtica deste modelo (Camargo & Oliveira, 2015; Kuroda, 2016), é dada por N(t) = N(0)Eβ((r t) β). (65) sendo Eβ(r t β) a função de Mittag-Leffler (Camargo & Oliveira, 2015). A seguir, são comparadas as soluções anaĺıticas com os resultados com- putacionais, utilizando o método MSGDTM, sendo r = 1 e N(0) = 1. Nas Figuras 7 e 8, pode-se notar que o GDTM está aparentemente em bom acordo com a solução anaĺıtica da equação (65). Ao utilizar o MSGDTM, a solução numérica apresenta singularidades nos pontos da tm, com 1 ≤ m ≤ M−1. A inexistência de tais singularidades na solução anaĺıtica prova que o MSGDTM, como apresentado na literatura, não está correto. De fato, Momenzadeh & Sarv Ahrabi (2017) apontaram a inconsistência da versão atual do método. Com o objetivo de aprofundar mais a análise dos casos considerados por eles, aplicaremos o MSGDTM a um problema de valor inicial da equação de Riccati com derivada de Caputo de ordem β, sendo 0 < β ≤ 1, com ińıcio em t0 = 0. Com a condição inicial x(0) = 0, a equação diferencial tem a forma 1 τ 1−β CD β 0x(t) = 2x(t)− x2(t) + 1. (66) Nestas condições, o GDTM leva em X(0) = 0 e na relação de recorrência na equação (61). 47 0 1 2 3 4 0 20 40 60 80 100 120 Tempo P o p u la ç ã o MSGDTMcomumpasso 0 1 2 3 4 0 20 40 60 80 100 120 Tempo P o p u la ç ã o MSGDTMcomdois passos Figura 7: Método MSGDTM e solução anaĺıtica: As linhas tracejadas correspondem à solução anaĺıtica, enquanto as linhas cont́ınuas correspondem ao método MSGDTM para β = 1, 0, β = 0, 8, β = 0, 6, β = 0, 4 e β = 0, 2, respectivamente, de baixo para cima. Os resultados numéricos para a equação de Riccati (66) de ordem β = 0, 7 foram obtidos no intervalo de [0; 0, 4]. Neste intervalo existe a convergência do método no modelo, no entanto, se estender esse intervalo de solução nota-se que o método deixa de convergir em t = 0, 6 aproximadamente. Por isso, a utilização no método MSDTM é utilizado em modelos com ordem de derivação inteira para resolver o problema de convergência. Aplicado o GDTM, temos x(t) ≈ 1, 10 t0,7 + 1, 61 t1,4 + 1, 14 t2,1 − 0, 60 t2,8 − 2, 54 t3,5, para 0 ≤ t ≤ 0, 4. Ao contrário do modelo de Malthus, a equação de Riccati é não linear, o que dificulta o uso da transformada de Laplace. Assim, neste segundo caso não dispomos de solução alternativa para comparação. 48 0 1 2 3 4 0 20 40 60 80 100 120 Tempo P o p u la ç ã o MSGDTMcomquatro passos 0 1 2 3 4 0 20 40 60 80 100 120 Tempo P o p u la ç ã o MSGDTMcom o� �� � passos Figura 8: Método MSGDTM e solução anaĺıtica: As linhas tracejadas correspondem à solução anaĺıtica, enquanto as linhas cont́ınuas correspondem ao método MSGDTM para β = 1, 0, β = 0, 8, β = 0, 6, β = 0, 4 e β = 0, 2, respectivamente, de baixo para cima. Dividindo o intervalo em dois, ou seja, o MSGDTM com dois passos de [0; 0, 2] e [0, 2; 0, 4] teremos como condições iniciais, y(0) = 0 e y(0, 2) ≈ 0, 55. Desta vez, a solução tem a forma x(t) ≈ 1, 10 t0,7 + 1, 61 t1,4 + 1, 14 t2,1 − 0, 60 t2,8 − 2, 54 t3,5, para 0 ≤ t ≤ 0, 2 e x(t) ≈ 0, 55 + 1, 98(t− 0, 2)0,7 + 1, 30(t− 0, 2)1,4 − 1, 54(t− 0, 2)2,1− −3, 07(t− 0, 2)2,8 + 0, 66(t− 0, 2)3,5, 49 para 0, 2 ≤ t ≤ 0, 4. Na Figura 9 são apresentadas as soluções com um passo e dois passos do método MSGDTM para o problema da equação (66). Observa-se que o MSGDTM com dois passos produz uma singularidade em t = 0, 2. Novamente, isto exibe a inconsistência do método. 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 t 0 0.2 0.4 0.6 0.8 1 1.2 x( t) MSGDTM: 1 Passo MSGDTM: 2 Passos Figura 9: Solução da equação (66) utilizando o MSGDTM com um e dois passos. Resta explicar por que um método inconsistente tem sido amplamente aplicado como se fosse válido. De um lado está a dificuldade de verificar numeri- camente a qualidade das soluções. Do outro lado está, ao que parece, a evolução do gráfico produzido pelo MSGDTM quando o número de passos, M , aumenta suficientemente. Para deixar isto mais claro, divideremos o intervalo em quatro partes, ou seja, aplicaremos o método MSGDTM em 4 passos, [0; 0, 1], [0, 1; 0, 2], [0, 2; 0, 3] e [0, 3; 0, 4]. As respectivas condições iniciais são x(0) = 0, x(0, 1) ≈ 0, 29, x(0, 2) ≈ 0, 68 e x(0, 3) ≈ 1, 12. As soluções obtidas são, x(t) ≈ 1, 10 t0,7 + 1, 61 t1,4 + 1, 14 t2,1 − 0, 60 t2,8 − 2, 54 t3,5, para 0 ≤ t ≤ 0, 1; 50 x(t) ≈ 0, 29 + 1, 65(t− 0, 1)0,7 + 1, 71(t− 0, 1)1,4 − 0, 16(t− 0, 1)2,1− −2, 75(t− 0, 1)2,8 − 2, 53(t− 0, 1)3,5, para 0, 1 ≤ t ≤ 0, 2; x(t) ≈ 0, 68 + 2, 09(t− 0, 2)0,7 + 0, 97(t− 0, 2)1,4 − 2, 12(t− 0, 2)2,1− −2, 54(t− 0, 2)2,8 + 2, 54(t− 0, 2)3,5, para 0, 2 ≤ t ≤ 0, 3; x(t) ≈ 1, 12 + 2, 19(t− 0, 3)0,7 − 0, 37(t− 0, 3)1,4 − 2, 65(t− 0, 3)2,1+ +1, 06(t− 0, 3)2,8 + 4, 52(t− 0, 3)3,5, para 0, 3 ≤ t ≤ 0, 4. Na Figura 10 representamos graficamente as soluções com um, dois e quatro passos do MSGDTM, aplicado ao problema da equação (66). Pode-se observar a presença de novas singularidades espúrias, ao comparar a curva de quatro passos com a do método GDTM. O resultado para M = 300, no intervalo 0 ≤ t ≤ 3, é mostrado na Figura 11. Desta vez há 299 singularidades espúrias, mas elas não são percept́ıveis no gráfico. Aparentemente, o limite da sequência de soluções produzidas pelo MSGDTM converge para uma função suave. Isto pode ter sido interpretado como ind́ıcio dela ser uma solução do problema. Assim, nesta Seção, foi apresentado e analisado o método multi-passos com transformada diferencial generalizada Arshad et al. (2015) para a resolução de equações diferenciais de ordem fracionária com derivada de Caputo. Mediante aplicações aos modelos de Malthus e de Riccati, verificou-se que o método é incon- sistente na sua versão multi-passos, comparando com a solução obtida pelo método da transformada diferencial generalizada com apenas um passo. A inconsistência é devido à extensão descuidada de técnicas válidas apenas para derivadas de ordem 51 0 0.1 0.2 0.3 0.4 t 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 x( t) MSGDTM: 1 Passo MSGDTM: 2 Passos MSGDTM: 4 Passos Figura 10: Solução da equação (66) utilizando o método MSGDTM com um, dois e quatro passos. 0 0.5 1 1.5 2 2.5 3 t 0 0.5 1 1.5 2 2.5 x( t) MSGDTM: 300 Passos Figura 11: Solução da equação (66) utilizando o método MSGDTM com 300 passos no intervalo de [0; 3]. inteira. Tal extensão desconsidera a não localidade da derivada de Caputo e detalhes da propriedade da transformada diferencial generalizada dessa derivada. A aparente convergência das soluções para uma função suave e a carência de soluções exatas para os modelos investigados podem ter dificultado a detecção da inconsistência. 52 A análise cŕıtica aqui apresentada contribuirá para reforçar o necessário rigor ma- temático nas diversas aplicações do cálculo fracionário. Visando um aporte mais construtivo, correções e alternativas para o método investigado estão em andamento (Kuroda et al., 2019). Utizando o método GDTM, o raio de convergência pode ser finito ou infinito. A avaliação numérica da série também pode apresentar instabilidade numérica, limitando mais o intervalo de resolução da equação. Assim, Bruno-Alfonso et al. (2019) resolveram a equação loǵıstica com derivada fracionária de Caputo me- diante um tratamento anaĺıtico-numérico. Inicialmente utilizaram o GDTM. A série de potências fracionárias correspondente tem raio de convergência finito. Foi usada até metade do raio de convergência. A evolução posterior foi obtida mediante gene- ralização do método de Euler. Verificaram que a solução satisfaz numericamente a equação loǵıstica. No próximo Caṕıtulo iremos apresentar o método de Grünwald- Letnikov que leva em consideração as derivadas em tempos anteriores, o que reflete que preserva o efeito da memória (Scherer et al., 2011). 4 MÉTODO DE GRÜNWALD-LETNIKOV Segundo Ongun et al. (2013), embora a discussão sobre derivadas de ordem não inteira remonte quase ao desenvolvimento da teoria clássica do cálculo diferencial de ordem inteira, somente no final do século XIX se percebeu as gran- des melhorias que poderiam ser alcançadas pela exploração do poder de cálculo fracionário; por meio de equações diferenciais fracionárias é posśıvel descrever, de maneira natural, fenômenos da vida real com efeito de memória (Kempfle et al., 2002). Na maioria dos casos, a solução de uma equação diferencial fracionária não existe em termos de um número finito de funções elementares. Portanto, é fundamen- tal utilizar métodos numéricos para avaliar soluções aproximadas (Diethelm et al., 2005; Garrappa & Popolizio, 2011; He et al., 2017). Neste caṕıtulo abordamos o método numérico de Grünwald-Letnikov (GL) baseado nos trabalhos de Scherer et al. (2011); Ongun et al. (2013); MacDonald et al. (2015); Tavoni (2019). Aplicado em equações diferenciais da ordem fracionária β, este método é uma extensão do método de Euler para equações diferenciais or- dinárias (Scherer et al., 2011). A derivada de Caputo é aproximada pela abordagem de Grünwald- Letnikov usando diferenças finitas de ordem fracionária. Para isso, vamos estabelecer a relação entre as derivadas de Riemann-Liouville, Caputo e de Grünwald-Letnikov, vistas no primeiro caṕıtulo, RLD βf(t) = Dn[In−βf(t)] = 1 Γ(n− β) ( dn dtn ) ∫ t 0 f(τ)dτ (t− τ)β−n+1 , CD βf(t) = In−β[Dnf(t)] = 1 Γ(n− β) ∫ t 0 f (n)(τ)dτ (t− τ)β−n+1 , 54 GLD βf(t) = lim h→0 1 hβ ∞ ∑ k=0 (−1)k ( β k ) f(t− kh). Primeiramente, mencionaremos a relação existente entre a derivada de Caputo com de Riemann-Liouville para m − 1 < β ≤ m, sendo m o menor inteiro positivo tal que β < m (Scherer et al., 2011; Tavoni, 2019), CD βf(t) = RLD βf(t)− m−1 ∑ ν=0 rβν (t)f (ν)(0), sendo rβν (t) = tν−β Γ(ν + 1− β) . Utilizaremos a aplicação do método de Grünwald-Letnikov em modelos fracionários com 0 < β ≤ 1, assim, com m = 1 temos CD βf(t) = RLD βf(t)− rβ0 (t)f(0), sendo rβ0 (t) = t−β Γ(1− β) . Como visto na Proposição 3, referente a derivada fracionária de Riemann-Liouville de uma constante c, RLD βc = ct−β Γ(1− β) · Neste caso c = f(0), assim CD βf(t) = RLD β (f(t)− f(0)) . Como as derivadas de Riemann-Liouville e de Grünwald-Letnikov coin- cidem (Ongun et al., 2013), CD βf(t) = GLD β (f(t)− f(0)) , com 0 < β ≤ 1. (67) Além disso, CD βf(t) = GLD βf(t)− m−1 ∑ ν=0 rβν (t)f (ν)(0), com m− 1 < β ≤ m, (68) sendo rβν (t) = tν−β Γ(ν + 1− β) . 55 Para implementar o método de Grünwald-Letnikov, vamos considerar (Scherer et al., 2011; Ongun et al., 2013; Tavoni, 2019), (1− z)β = ∞ ∑ j=0 ω (β) j zj , (69) em que ω (β) j = (−1)j ( β j ) . Os pesos do operador podem ser avaliados por meio da seguinte fórmula de recorrência (Momani & Obibat, 2007; Jiang et al., 2017), ω (β) 0 = 1, ω (β) j = ( 1− 1 + β j ) ω (β) j−1, j = 1, 2, 3 . . . , e ω (β−1) 0 = 1, ω (β−1) j = ( 1− β j ) ω (β−1) j−1 , j = 1, 2, 3 . . . . Proposição 7 Considere o peso ω (β) j , com j = 0, . . . , n, temos que n ∑ j=0 ω (β) j = ω(β−1) n . A prova da Proposição 7 pode ser vista em Tavoni (2019). Considere uma equação diferencial fracionária com derivada de Caputo de ordem 0 < β ≤ 1, sujeita à condição inicial y(t0) = y0, CD βy(t) = f(t, y(t)). (70) Sabe-se, pela equação (67), que CD βf(t) = GLD β (f(t)− f(0)) , utilizando a derivada fracionária de Grünwald-Letnikov e a relação (69), GLD βf(t) = lim h→0 1 hβ ∞ ∑ k=0 ω (β) k f(t− kh). Assim, a aproximação de Grünwald-Letnikov da equação (70) em [t0, t], 0 < t < T , com passo h → 0 e n = t− t0 h , 56 CD βy(t) = f(t, y(t)), GLD β (y(t)− y(0)) = f(t, y(t)), 1 hβ n ∑ j=0 ω (β) j ( y(t− jh)− y ( t0 − j t0 − t0 n )) = f(t, y(t− (j + 1))h), 1 hβ n ∑ j=0 ω (β) j (y(t− jh)− y0) = f(t, yt−1), n ∑ j=0 ω (β) j (yn−j − y0) = hβf(t, yn−1), ω (β) 0 (yn − y0) + ω (β) 1 (yn−1 − y0) + . . .+ ω(β) n (y0 − y0) = hβf(t, yn−1), ω (β) 0 yn + ω (β) 1 yn−1 + . . .+ ω(β) n y0 − ω (β) 0 y0 − . . .− ω(β) n y0 = hβf(t, yn−1), yn + ( n ∑ j=1 ω (β) j yn−j ) − y0 ( ω (β) 0 + ω (β) 1 + . . .+ ω(β) n ) = hβf(t, yn−1), yn + ( n ∑ j=1 ω (β) j yn−j ) − y0ω (β−1) n = hβf(t, yn−1). Portanto, o método expĺıcito de Grünwald-Letnikov para uma equação diferencial fracionária é dado pela equação, yn = y0ω (β−1) n − ( n ∑ j=1 ω (β) j yn−j ) + hβf(t, yn−1). (71) Outros resultados como a ordem de aproximação, ordem de con- sistência, estabilidade e erro de truncamento pode ser vistos em Podlubny (1999); Scherer et al. (2011); Ongun et al. (2013); Tavoni (2019). 4.1 Discretização do Modelo de Malthus e Riccati Para a implementação do método de Grünwald-Letnikov, foi utilizado o Code::Blocks em linguagem C. As soluções numéricas obtidas foram descritas em 57 arquivo .dat e para a exibição gráfica foi utilizado o programa Xmgrace. Como foi feito no caṕıtulo referente ao método MSGDTM, exibiremos as soluções do mesmo modelo de Malthus (64) pelo método de Grünwald-Letnikov juntamente com sua solução anaĺıtica. Quando ao modelo de Riccati (66), exibiremos apenas sua solução via método GL por se tratar de uma equação não linear. Assim, a discretização expĺıcita do modelo de Malthus pelo método de Grünwald-Letnikov é da forma, xn = ω(β−1) n x0 − ( n ∑ j=1 ω (β) j xn−j ) + hβrβxn−1. (72) A seguir na Figura 12, está a solução anaĺıtica (65) e a solução do modelo de Malthus via método de GL (72), 0 1 2 3 4 Tempo 0 20 40 60 80 100 120 P o p u la çã o Figura 12: Método de Grünwald-Letnikov e solução anaĺıtica do modelo de Malthus (64) com r = 1, h = 0, 001: As linhas tracejadas correspondem à solução anaĺıtica, enquanto as linhas cont́ınuas correspondem ao método GL para β = 1, 0, β = 0, 8, β = 0, 6, β = 0, 4 e β = 0, 2, respectivamente, de baixo para cima. A discretização expĺıcita do modelo de Riccati pelo método de Grünwald-Letnikov é da forma, xn = ω(β−1) n x0 − ( n ∑ j=1 ωβ j xn−j ) + hβτ 1−β [ 2xn−1 − x2 n−1 + 1 ] . (73) 58 A seguir na Figura 13, está a solução do modelo de Riccati via método de GL (72), 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 Tempo 0 0.2 0.4 0.6 0.8 1 1.2 x( t) Figura 13: Solução da equação de Riccati (66) utilizando o método de Grünwald- Letnikov para β = 0, 7, τ = 1 e h = 0, 001. Nota-se que existe uma boa aproximação da solução anaĺıtica com a solução via GL, ao contrário do método MSGDTM que não leva em consideração a não localidade das derivadas fracionárias. O método de Grünwald-Letnikov leva em consideração as derivadas em tempos anteriores, ou seja, preserva o efeito da memória (Scherer et al., 2011). Este método é computacionalmente eficiente e resulta em erros menores durante as simulações numéricas (MacDonald et al., 2015; Tavoni, 2019). No próximo caṕıtulo, utilizaremos as diferentes modelagens citadas no Caṕıpulo 2 e exibiremos suas respectivas curvas de soluções de ordem fracionárias em modelo de câncer conhecido como HPV 16 via método de Grünwald-Letnikov. 5 MODELO DE CRESCIMENTO TUMORAL (HPV 16) Segundo Mirabello et al. (2016), cerca de 200 tipos de infecções por papilomav́ırus humano (HPV) não causam patologia evidente ou apenas papilomas benignos (verrugas). No entanto, os tipos de HPV 16, 18, 31, 33, 35, 39, 45, 51, 52, 56, 58, 59 e 68 (Bouvard et al., 2009) podem causar todos os casos de câncer cervical e outros cânceres anogenitais6 e orofaŕıngeos7. Mais de meio milhão de pessoas a cada ano são acometidas por um câncer relacionado ao HPV em todo o mundo (Forman et al., 2012) e mais de 200.000 mortes ocorrem a cada ano (Mirabello et al., 2017). Em Santos et al. (2011), o Brasil pode ser considerado um páıs com alta incidência de infecções pelo HPV, e o Ministério da Saúde registra a cada ano aproximadamente 137 mil novos casos. Diferentes estudos epidemiológicos indicam que os tipos de HPV prevalentes são diferentes de acordo com a região de origem, sendo que, no Brasil, o HPV 16 foi o tipo oncogênico mais encontrado (Giulliano et al., 2008). Essas infecções carcinogênicas8 por HPV são extremamente comuns e facilmente transmitidas aos epitélios9 suscet́ıveis pelo contato f́ısico direto (Hans- Ulrich et al., 2010). Apenas uma pequena proporção de infecções por alto risco 6Relativo ao ânus e aos órgãos genitais. 7Plural de “orofaŕınge”, relativo à orofaringe. Pequeno espaço da cavidade bucal que compreende a raiz da ĺıngua, o palato mole e a epiglote. 8Carcinógeno, também chamado de canceŕıgeno ou carcinogênico é a qualidade daquilo capaz de provocar ou estimular o aparecimento de carcinomas (câncer que afeta o tecido epitelial, a camada superior da pele (epitélio) ou o revestimento dos órgãos internos de um organismo, podendo se espalhar através de metástase) ou câncer em um organismo. 9Epitélio ou tecido epitelial é um dos quatro tipos de tecidos básicos do organismo humano. 60 persiste e progride para o pré-câncer, que por sua vez pode levar ao câncer (McCredie et al., 2008; Schiffman et al., 2007). Não sabe-se por que essas infecções comuns e tipicamente benignas às vezes causam câncer. Estudos anteriores sugerem que a variação genética viral está associada ao risco de câncer. O câncer do colo do útero humano é causado por papilomav́ırus hu- manos de alto risco que codificam as oncoprotéınas E6 e E7, cada uma das quais altera a função de alvos distintos que regulam o ciclo celular, a apoptose10 e a dife- renciação11 (Riley et al., 2003). Dentro os diferentes tipos de HPV, o HPV 16 causa metade dos cânceres cervicais em todo o mundo e é um dos carcinógenos humanos mais importantes (Mirabello et al., 2016). A história natural do carcinoma do colo uterino pode ser dividida em três fases: a primeira quando está presente a infecção por HPV, sem outras mani- festações detectáveis; a segunda quando já estão presentes alterações morfológicas das células do epitélio do colo uterino, que caracterizam as lesões intra-epiteliais; e a terceira com a presença de lesão atravessando a membrana basal do epitélio, carac- terizando o carcinoma invasor, fase esta irreverśıvel e que se não tratada levará ao óbito (Rama et al., 2006). Segundo Giulliano et al. (2008), dentre os métodos de tratamento não invasivos dispońıveis atualmente, existe uma toxina com propriedades citotóxicas12 chamada podofilina (Scheinfeld & Lehman, 2006) e o ácido tricloroacético13, um 10Apoptose ou Morte Celular Programada é um tipo de “autodestruição celular” que requer energia e śıntese protéica para a sua execução. Está relacionado com a homeostase na regulação fisiológica do tamanho dos tecidos, exercendo um papel oposto ao da mitose. 11A diferenciação celular consiste em um conjunto de processos que transformam e especializam as células embrionárias. Após estas transformações, sua morfologia e fisiologia são definidas, o que as tornam capazes de realizar determinada função. 12A citotoxicidade é a propriedade nociva de uma substância em relação às células. Células citotóxicas são células que têm a capacidade de destruir outras células através da libertação de certas substâncias nocivas. 13O ácido tricloroacético faz parte do tratamento de algumas doenças, por ser um cauterizante. No geral, o ácido tricloroacético é usado para tratamentos de pele. 61 agente cáustico14 que combate verrugas anogenitais por meio do contato (Tchernev, 2009). É descrito também o uso de interferons, protéınas endógenas intracelulares que apresentam atividade antitumoral e antiviral (Bharti et al., 2009) e imiquimod, um medicamento usado no tratamento de verrugas anogenitais, cujo mecanismo de ação consiste na potencialização da produção de interferons, com efeito antiviral, antiproliferativo e antiangiogênico (Yang et al., 2009; Garland et al., 2006). Além destes, cidofovir é um importante antiviral usado em tumores induzidos pelo HPV, com a capacidade de provocar morte celular programada entre as células tumorais (Tchernev, 2009). Um dos métodos invasivos mais confiáveis para a remoção de lesões únicas ou múltiplas na região anogenital é o método cirúrgico (Scheinfeld & Lehman, 2006). O tratamento com laser apresenta eficácia entre 60% e 77%, no entanto é um método de custo elevado e que pode liberar DNA viral na vaporização das lesões durante o procedimento (Bharti et al., 2009). A criocirurgia é outra medida utilizada em verrugas genitais, levando à necrose da epiderme e, ocasionalmente, da parte superior da derme. 5.1 Conjunto de Dados Esta seção é referente ao estudo do artigo “Model-Based Tumor Growth Dynamics and Therapy Response in a Mouse Model of De Novo Carci- nogenesis” de Loizides et al. (Loizides et al., 2015). Entramos em contato com um dos autores, Georgios D. Mitsis, o qual forneceu o conjunto de dados e o material auxiliar utilizados em sua pesquisa na descrição do modelo referente ao volume tumoral em camundongos transgênicos que expressam os oncogenes15 E6 e E7 do HPV 16 administrados em seu dorso. Os oncogenes E6 e E7 contribuem de maneira diferente para a ca