Karen Pieri Um Modelo Espacial para a Resposta Imunológica Adaptativa da Hepatite B Botucatu 2010 Karen Pieri Um Modelo Espacial para a Resposta Imunológica Adaptativa da Hepatite B Monografia apresentada ao Instituto de Biociências da Universidade Estadual Paulista "Júlio de Mesquita Filho", Campus de Botucatu, para obtenção do título de Bacharel em Física Médica. Orientadora: Profa. Dra. Cláudia Pio Ferreira Botucatu 2010 Agradecimentos Agradeço primeiramente à Deus, por sempre mostrar Sua presença na minha vida, mesmo nas coisas mais pequeninas (como é bom caminhar debaixo de Sua graça!). Ao meu pai, o homem que sem dúvida eu mais admiro, pois sempre investiu nos meus sonhos sem nunca duvidar, sempre confiando e acreditando no que eu sou (mesmo nas horas em que eu mesma não confiei). Por ser exemplo de uma das poucas pessoas que conheço que com muito esforço e integridade foi capaz de mudar sua história de vida. (Eu deveria receber dois diplomas pai: um para mim e outro para você, pois, com certeza você trilhou esse caminhou todo, semestre a semestre e se formou comigo). À minha mãe e à minha irmã, Karina, pelo aconchego de todo final de semana e o carinho com que sempre me recebiam de volta. Pelas inúmeras conversas, pelas risadas. Por se fazerem tão presentes mesmo quando estavamos tão longe. Ao Adriel, por ocupar um lugar tão importante na minha vida, sempre me apoiando em cada decisão, em cada palavra, e por nas horas mais difíceis ter me ajudado a prosseguir. E claro, por sempre me mostrar o lado engraçado das coisas. Aos meu avôs, Wanda e Pedro, que agora tem alegrado nossa casa, e mesmo antes, sempre me ajudaram infinitas vezes. Pelas esperas no ponto de ônibus, pelos lanchinhos para levar no estágio, pela disposição de quem já vive uma fase sossegada e aceita voltar à correria por um tempo para fazer parte do sonho de outra pessoa. (Como prometido farei uma cópia do meu diploma e darei para vocês num quadro, pois um pedacinho dele pertence a vocês também). Aos amigos de Botucatu que me deixaram tanta saudade, foram minha família nestes anos, sempre com muito companheirismo e alegria, dividiram comigo cada momento dessa nossa pequena luta. Em especial à Carol, Tost, Lana, Horofote, Nhonho, Katatau e Berinjela (insisti em manter os apelidos, pois quero me lembrar de vocês exatamente assim). Dificilmente encontrarei amigos tão bons quanto vocês souberam ser para mim. A todos da IV turma de Física Médica da UNESP. À minha orientadora, professora Cláudia Pio, pelo seu tão grande comprometimento e esforço com o meu aprendizado e trabalho. Mais que isso, pelo carinho com que ensina seus alunos, sendo sempre acessível para auxiliar e para dar inúmeros conselhos profissionais. Não poderia esquecer da minha colega de curso e veterana, Priscila Muniz, quem com tanto empenho iniciou este trabalho estudando a resposta imunológica adaptativa, me dando a oportunidade de continuá-lo e desenvolvê-lo. A todos do Departamento de Bioestatística e desta universidade que contribuiram para este projeto ou simplesmente fizeram parte destes anos, tornando meus dias aqui muito mais felizes. Este trabalho foi financiado pela Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP). “...Seja bendito o nome de Deus de eternidade a eternidade, porque Dele são a sabedoria e a força; e ele muda os tempos e as estações; Ele remove os reis e estabelece os reis; Ele dá sabedoria aos sábios e conhecimento aos entendidos. Ele revela o profundo e o escondido; conhece o que está em trevas, e com Ele mora a luz. Ò Deus de meus pais, eu Te dou graças e Te louvo, porque me deste sabedoria e força...” Dn(2.20-23) Resumo Este trabalho propôs o desenvolvimento de um modelo espacial bidimensional que descrevesse a resposta imunológica adaptativa para a infecção viral da hepatite B. Este modelo foi construído baseando-se na dinâmica de interação entre seis populações: he- patócitos saudáveis T , hepatócitos infectados Y , células virais V (representando o vírus da hepatite B, i.e., HBV ), sistema imune inativo I, sistema imune ativo X e células de memória X. Inicialmente construiu-se um modelo compartimental adimensional e estudou-se as soluções de equilíbrio e os valores limiares que determinam a estabilidade de cada uma delas. Com este modelo foi possível reproduzir as diversas evoluções observadas para a doença, são elas: indivíduos que eliminam a infecção sem formar resposta imunológica, portadores agudos e portadores crônicos. Ao incluirmos difusão das células de defesa do sistema imunológico e do vírus (modelo espacial), analisamos duas situações: modelo homogêneo, no qual os parâmetros do mo- delo são os mesmos em todos os pontos da rede, e modelo heterogêneo, o qual caracteriza células mais permeáveis ou menos permeáveis à invasão do vírus. Para os dois modelos espaciais observou-se que, os tempos de eliminação viral e/ou invasão do vírus tornam-se menores em relação ao modelo compartimental. Os resultados mostraram ainda que, para o conjunto de valores utilizados nas simulações, o modelo é sensível às variações na taxa de difusão viral e não depende da difusão das células de memória, considerando o caso em que as duas taxas de difusão são diferentes de zero. Finalmente, o modelo heterogêneo quando comparado ao modelo homogêneo, mostra que podemos ter a infecção limitada espacialmente dependendo do tipo de célula envolvida no processo de infecção. Palavras-chave: hepatite B; resposta imunológica adaptativa; modelagem matemática; modelos espaciais. Abstract This paper proposed a two-dimensional spatial model to describe the adaptive immune response for viral hepatitis B. This model considered six populations: healthy hepatocytes T , infected hepatocytes Y , hepatitis B virus V , innate immune system I, active immune system X and memory cells, X. First, a compartmental model was constructed and its equilibrium solutions and also the threshold values related to the stability of each solution were obtained. Using this model, we was able to reproduce the different trends observed for the disease, which are: individuals that eliminate the infection without forming immune response, patients with acute and chronic carriers. By including dispersion of defense cells of the immune system and virus (spatial mo- del), we analyze two situations: homogeneous model, in which the model parameters are the same at all points of the network, and heterogeneous model, which characterizes cells more permeable and less permeable to virus invasion. For the two spatial models (ho- mogeneous and heterogeneous) the times relatead to the viral erradication and/or virus invasion and persistence becoming smaller in relation to the compartmental model. The results also showed that for the set of values used in the simulations and if the two diffu- sion rates are different from zero, the model is sensitive to variations in the rate of viral spread and not dependent on the dispersion of memory cells. Finally, the heterogeneous model when compared to the homogeneous model shows that the infection can be spatially limited depending on the type of the cell involved in the infection process. Key-words: hepatitis B; adaptive immune response; mathematical modeling, spatial models. 7 Sumário Resumo 5 Abstract 6 1 Introdução 9 1.1 Hepatite B (HBV ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.2 Curso sorológico da hepatite B . . . . . . . . . . . . . . . . . . . . . . . . 10 2 Objetivos 14 3 Modelo Dinâmico Viral Sem Acoplamento Espacial 14 3.1 Modelo Matemático da Dinâmica Viral do HBV . . . . . . . . . . . . . . . 14 3.2 Modelo Adimensonal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.3 Soluções de equilíbrio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.4 Estudo da Estabilidade dos Pontos de Equilíbrio . . . . . . . . . . . . . . . 18 3.4.1 Estabilidade de P0 . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.4.2 Estabilidade de P1 . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.4.3 Estabilidade de P2 . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.5 Resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 4 Modelo Dinâmico Viral Espacial 26 4.1 Modelo Matemático Espacial da Dinâmica Viral do HBV . . . . . . . . . . 26 4.2 Metodologia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4.3 Resultados e Discussão . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 4.3.1 Acoplamento Espacial Homogêneo . . . . . . . . . . . . . . . . . . 29 4.3.2 Acoplamento Espacial Heterogêneo . . . . . . . . . . . . . . . . . . 34 5 Conclusão 36 8 Lista de Figuras 1 Curso sorológico da hepatite B para um portador agudo. . . . . . . . . . . 12 2 Curso sorológico da infecção crônica pelo vírus da hepatite B para um portador crônico. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3 Evolução temporal das células para o modelo compartimental convergindo para o ponto de equilíbrio livre da doença. . . . . . . . . . . . . . . . . . . 23 4 Evolução temporal das células para o modelo compartimental convergindo para o ponto de equilíbrio livre da doença com formação de células de memória. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 5 Evolução temporal das células para o modelo compartimental convergindo para o ponto de equilíbrio endêmico, onde as populações de células coexis- tem com a população viral. . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 6 Dispersão espacial da população viral com acoplamento espacial e parâme- tros relativos à solução de eliminação viral. . . . . . . . . . . . . . . . . . . 30 7 Evolução temporal das células segundo o modelo espacial e conjunto depa- râmetros relativo ao ponto de equilíbrio livre da doença. . . . . . . . . . . 31 8 Dispersão da população viral com acoplamento espacial e conjunto de pa- râmetros relativos a coexistência das populações. . . . . . . . . . . . . . . . 32 9 Evolução temporal das células no modelo espacial e conjunto de parâmetros relativos a coexistência das populações. . . . . . . . . . . . . . . . . . . . . 32 10 Tempo de eliminação do vírus em função da taxa de difusão do vírus. . . . 33 11 Tempo de invasão do vírus em função da taxa de difusão do vírus. . . . . . 33 12 Dispersão espacial da população viral em malha heterogênea. . . . . . . . . 34 13 Evolução temporal das células segundo o modelo espacial, em malha hete- rogênea. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 9 1 Introdução Nesta seção discutiremos brevemente os principais aspectos envolvidos na infecção pelo HBV (vírus da hepatite B) com o objetivo de destacar pontos importantes para a for- mulação matemática deste problema biológico. Entre eles: estudo da cinética viral e, identificação dos cursos sorológicos de portadores crônicos e agudos. 1.1 Hepatite B (HBV ) A doença hepática devida ao HBV é um problema global com uma taxa mundial estimada de 300 milhões de individuos infectados [8]. No Brasil, sua distribuição geográfica não é uniforme e apresenta áreas de baixa prevalecência no sul e sudeste e, alta prevalecência na Região Amazônica [9]. O vírus da hepatite B (HBV ) pode evoluir de diversas formas após o contágio, po- dendo produzir hepatite aguda; hepatite crônica não progressiva; hepatite crônica progres- siva culminando em cirrose; hepatite fulminante com necrose hepática maciça; um estado de portador assintomático, com ou sem doença subclínica progressiva e, ainda servir de pano de fundo para a hepatite D (HDV ). Ele é 100 vezes mais infeccioso que o vírus da imunodeficiência humana (HIV ) e 10 vezes mais que o da hepatite C (HCV ). O HBV é um membro dos Hepadnaviridae, uma família de vírus contendo DNA que causa hepatite em múltiplas espécies de animais [6]. Sua partícula viral é constituída por uma molécula de ácido nucléico, envolvida por uma estrutura protéica chamada capsídeo. Além disso, ele é revestido por duas camadas: uma externa (envelope) e, outra interna (core). A camada externa é constituída por HBsAg, o antígeno de superfície do vírus da hepatite B. Este antígeno é a substância que o organismo identificará como estranha e que induzirá a produção de proteínas específicas pelo sistema imune. A interna é constituída por HBcAg, o antígeno c do vírus da hepatite B que assim como o antígeno s também desencadeará resposta imunológica. O HBV possui ainda uma outra proteína, HBeAg [3]. Estudos de cinética viral realizados in vivo indicam alta taxa de replicação do HBV em torno de 1011 a 1012 partículas virais por dia. Estima-se que cerca de 109 células infectadas sejam destruídas por dia, indicando que aproximadamente de 1% a 7% da 10 massa hepatocitária do fígado seja reciclada diariamente nos pacientes infectados. No entanto tem sido demontrado que a meia-vida dos hepatócitos que contém o HBV é al- tamente variável, dependendo principamente do grau de ativação do sistema imunológico do hospedeiro. Modelos matemáticos indicam que pacientes com hepatite crônica B que possuem meia vida das células infectadas ao redor de 10 dias necessitariam aproximada- mente de 1 ano de supressão da viremia para obter a erradicação completa do HBV do organismo, enquanto que pacientes com meia vida das células infectadas da ordem de 100 dias necessitariam de mais de 10 anos de supressão para obter a cura [5]. É interessante ressaltar que o HBV não causa diretamente a patologia observada, já que a lesão hepatocelular ocorre devido a ação do próprio sistema imune do hospedeiro contra os hepatócitos infectados, ou seja, o vírus faz com que as células de defesa ataquem as células do fígado, sendo esse fato essencial para a resolução da doença hepática. 1.2 Curso sorológico da hepatite B A seguir, estudaremos separadamente os cursos sorológicos das hepatites aguda e crônica resultantes da infecção pelo HBV . a. Hepatite B Aguda: A Figura 1 apresenta a evolução temporal dos marcadores sorológicos da hepatite B para um portador agudo. Observamos que – o HBsAg está presente no soro antes mesmo do início das alterações da alanina transferase (ALT ) e do aparecimento dos sintomas. Este é portanto, o primeiro a aparecer e sua concentração máxima coincide com o início dos sintomas. Por isso, é um dos marcadores pesquisados tanto na triagem de doadores de sangue em unidades hemoterápicas quanto no diagnóstico em laboratórios de saúde pública. Declina a níveis indetectáveis no prazo de 3 a 6 meses [1]; – OHBeAg surge no início da infecção, logo após oHBsAg, e pode ser detectado durante várias semanas. Este marcador indica intensa multiplicação viral e, portanto, um maior potencial de transmissão. Nota-se que quando a replicação viral é máxima a curva de ALT também apresenta seu máximo, uma vez que 11 a última é uma enzima que é liberada quando ocorre algum tipo de agressão à membrana dos hepatócitos; – Já o anti-HBcIgM surge na fase aguda, atingindo quantidades detectáveis no soro pouco antes do aparecimento dos sintomas, podendo permanecer até oito meses após o início da infecção. É o único capaz de ser detectado na fase de janela imunológica, sendo de extrema importância na detecção de casos agudos. Geralmente, não é detectado em portadores crônicos; – O anti-HBcIgG substitui gradativamente o anti-HBcIgM e, em geral, persiste por toda vida, sendo útil para indicar que houve em algum momento infecção pelo vírus da hepatite B; – O anti-HBcTotal(IgM + IgG) é o único marcador que define a etiologia da doença, detectável entre o desaparecimento dos antígenos e o aparecimento dos anticorpos. Como contém o anti-HBcIgG também persiste por toda vida sem conferir imunidade, servindo de parâmetro para saber se em algum momento da vida do indivíduo houve a infecção; – Por fim, o anti-HBe pode ser detectado após o desaparecimento do HbeAg in- dicando diminuição da multiplicação viral e a provável evolução para a cura da infecção, enquanto que o anti-HBs pode ser detectado após o desaparecimento do HBsAg. A presença de anti-HBs indica imunidade em relação à infecção, sendo o princípio das estratégias atuais de vacinação, que empregam HBsAg não infeccioso [1, 3]. Em pessoas vacinadas contra a hepatite B, espera-se portanto encontrar anti-HBs. Em cerca de 5% a 10% dos casos de infecção pelo vírus da hepatite B não há desenvolvimento de imunidade, configurando a evolução para a forma crônica [1]. b. Hepatite B crônica: Na Figura 2, podemos conferir como se caracteriza sorologicamente a evolução para a forma crônica da infecção pelo vírus da hepatite B. A permanência do HBsAg por mais de seis meses caracteriza este tipo de portador, sendo que podemos observar 12 Figura 1: Curso sorológico da hepatite B para um portador agudo [3]. que este marcador está sempre presente ao longo do seu curso sorológico. No curso sorológico da infecção crônica podem ainda ser detectados o HBeAg ou o anti-HBe (não mostrados na Figura 2): o primeiro sinaliza uma forma mais grave de infecção enquanto que o segundo indica um estágio com diminuição da multiplicação viral. Em geral nos pacientes que desenvolvem a forma crônica o Anti-HBcIgM deixa de ser detectado (como mostra a Figura 2, a partir de oito meses de infecção ), a não ser em casos de alta replicação viral. A presença do anti-HBcIgG, da mesma forma que na infecção aguda, persiste por toda a vida do paciente sem trazer qualquer ganho em termos de imunidade. Figura 2: Curso sorológico da infecção crônica pelo vírus da hepatite B para um portador crônico [3]. 13 Baseado na argumentação biológica anterior desenvolveu-se um modelo matemático para o estudo da evolução temporal da infecção pelo vírus da hepatite B enfatizando sua resposta imunológica adaptativa. As principais referências estudadas foram o artigo de Funk [2] e a monografia de Muniz [7]. A seção 3, Modelo Dinâmico Viral Sem Aco- plamento Espacial, traz a formulação do modelo matemático, a adimensionalização do mesmo, a obtenção dos pontos de equilíbrio e estudo da estabilidade para cada um deles e, por fim, os resultados. Na seção 4, Modelo Dinâmico Viral Espacial, estão descritas as alterações necessárias para adaptarmos o modelo compartimental a um modelo espacial, a metodologia seguida e, resultados obtidos para acoplamentos espaciais homogêneos e heterogêneo. Por fim, a seção 5 encerra esta monografia com as Conclusões. Esta monografia é baseada no trabalho de iniciação científica da aluna Karen Pieri, a qual realizou estágio no Departamento de Bioestatística do Instituto de Biociências da UNESP/Botucatu de Janeiro de 2007 à Dezembro de 2009, com a Profa. Dra Claudia P. Ferreira, no tema de modelagem matemática. Durante o estágio a aluna se familiarizou com o ambiente Linux, aprimorando seus conhecimentos em Linguagem de Programação C, processador de texto latex e outros pacotes do Linux, como por exemplo, os pacotes gráficos xmgrace e gnuplot. Foi bolsista FAPESP de Fevereiro de 2008 à Janeiro de 2010 sob responsabilidade da Profa Dra Claudia Pio Ferreira. Resultados parciais do trabalho foram apresentados em congressos, como o CONFIAM- Congresso de física aplicada a medicina - em 2009 (Botucatu - SP), Workshop da Pós - Graduação "ciência, mercado e sociedade - em 2009 (Botucatu - SP) e Encontro regional de matemática aplicada e computacional (Petrópolis - RJ). No ano de 2010 a aluna realizou estágio de janeiro à novembro no Centro Infantil de Investigações Hematológica Dr. Domingos A. Boldrini, no setor de física da radioterapia, e também desenvolveu o trabalho de pesquisa "Implementação do Sistema de Dosimetria in Vivo no Hospital Infantil Boldrini”, trabalhando com um programa de controle de qualidade adicional baseado em dosímetros semicondutores, com o objetivo de utilizá-lo em procedimentos de T.B.I (total body irradiation), tratamento de radioterapia aplicado comumente após transplantes de medula. Este trabalho está inscrito para ser apresentado no XII Congresso Brasileiro de Radioterapia em outubro de 2010. 14 2 Objetivos Desenvolver um modelo espacial que descreva a dinâmica da resposta imunológica adap- tativa para a infecção do vírus da hepatite B. Para isso, utilizou-se a idéia apresentada no artigo de Funk [2], isto é, considerou-se acoplamento espacial entre as variáveis do modelo compartimental original, o qual considera seis populações de células: hepatócitos saudáveis, hepatócitos infectados, células virais, sistema imune inativo, sistema imune ativo e células de memória. Como a maioria das infecções tem uma localização espacial e as interações entre as populações de vírus e células são locais, modelos com estrutura espacial são um importante passo para entendermos o que acontece quando a suposição de populações bem misturadas é relaxada. A comparação entre os dois modelos (com e sem estrutura espacial) foi feita através da medida da intensidade e da velocidade de espalhamento da infecção. 3 Modelo Dinâmico Viral Sem Acoplamento Espacial Esta seção apresenta o modelo matemático da dinâmica viral da hepatite B, sem aco- plamentos espaciais, seguido pelo estudo de seus pontos de equilíbrio e estabilidade do mesmos, analiticamente e/ou numericamente. 3.1 Modelo Matemático da Dinâmica Viral do HBV Este modelo foi construído baseando-se na dinâmica de interação entre seis populações: hepatócitos saudáveis T , hepatócitos infectados Y , células virais V (representando o vírus da hepatite B, i.e., HBV ), sistema imune inativo I, sistema imune ativo X e células de memória X. Sabemos que as células do fígado se regeneram e que o vírus ao sair da célula infectada não a destrói, de maneira que células saudáveis e infectadas crescem de acordo com um termo logiśtico, o que nos leva a incluir um termo de suporte K, relativo ao número de células presentes no fígado. No modelo, a população de hepatócitos saudáveis T varia no tempo de acordo com as seguintes condições: as células são produzidas a uma taxa r1 (crescimento logístico) e morrem a uma taxa u1 e, ainda têm uma fração transformada em hepatócitos infectados Y 15 devido à taxa de encontro e entre o vírus e os hepatócitos sadios. Por sua vez, hepatócitos infectados, Y , são produzidos a uma taxa e, devido ao encontro entre vírus HBV e hepatócitos normais, e a uma taxa r2 (crescimento logístico). Sua morte pode estar associada a dois fatores: a taxa u2, de degradação celular natural, ou devido ao encontro com probabilidade α com uma célula do sistema imune ativo. O vírus V , varia apenas com a taxa de crescimento regida pelo parâmetro p e, pela taxa de mortalidade dada por c. As células imunes inativas I são produzidas no timo de acordo com a taxa f e caso encontrem hepatócitos infectados são ativadas em a uma taxa a (sendo também propor- cionais à quantidade de células infectadas). Células imunológicas ativadas, geram duas células imunes ativas (reprodução binária) X, evento expresso por 2aIY , e têm taxa de mortalidade natural g. O sistema imune ativo, X, segundo a taxa k é capaz de eliminar os hepatócitos infecta- dos. Células ativas também têm reprodução binária, originando duas células de memória M , sendo este evento controlado pelo termo q Xl l+X , em que l é a taxa de saturação e q é taxa de conversão de células imunes ativas em imunes de memória (capacidade de memorização limitada do sistema imunológico). As células de memória podem retornar ao estado de células imunes ativas, em duas situações: quando o organismo necessitar da apresentação rápida, pois a taxa de formação de células imunes ativas a partir das células de memória depende de um fator z1a, sendo z1 > 1, enquanto que, sua formação a partir do sistema imune ativo está relacionado apenas ao fator a, o que torna o primero método mais veloz; ou quando uma taxa de ativação de background n estiver presente. Finalmente, as células de memória morrem a uma taxa s. Essas condições biológicas nos levam ao seguinte conjunto de equações diferenciais ordinárias, o qual descreve a evolução temporal de cada uma das populações: Ṫ = r1T (1 − H K ) − eTV − u1T Ẏ = r2Y (1 − H K ) + eTV − αXY − u2Y V̇ = pY − cV (1) İ = f − gI − aIY 16 Ẋ = 2aIY − qX l +X − gX + (z1aY + n)M Ṁ = 2qX l +X − sM − (z1aY + n)M em que T + Y = H, sendo T, Y, V, I,X,M , respectivamente, hepatócitos não infectados, hepatócitos infecta- dos, vírus HBV livres, células imunes inativas, células imunes ativas e células de memória. A Tabela 1 resume a discussão anterior, associando cada um dos parâmetros ao seu respectivo sentido biológico. Tabela 1: Descrição dos parâmetros do modelo viral-imunológico da HBV Parâmetro Sentido Biológico r1 Taxa de crescimento dos hepatócitos saudáveis (T ) r2 Taxa de crescimento dos hepatócitos infectados (T ) K Capacidade de suporte do Meio e Taxa de Encontro entre hepatócitos (T ) e vírus HBV (V ) u1 Taxa de Mortalidade dos hepatócitos saudáveis u2 Taxa de Mortalidade dos hepatócitos infectados α Taxa de Encontro entre células imunes ativas (X) e hepatócitos infectados p Taxa de Crescimento do vírus HBV c Taxa de Mortalidade dos vírus HBV f Taxa de Crescimento das células imunes inativas (I) g Taxa de Mortalidade das células imunológicas (I e X) a Taxa de Encontro entre hepatócitos infectados e células imunes inativas q Taxa de Transição de I para células de memória M l Termo de Saturação da produção de M z1 Taxa de Ativação das células de memória n Taxa de Ativação de background das células M s Taxa de Mortalidade das células de memória 17 3.2 Modelo Adimensonal Utilizaremos o processo de adimensionalização para reduzir o número de parâmetros atra- vés da seguinte mudança de variáveis: ϕ = r1θ, t = T K , y = Y K , v = θ−1 1 V, i = I l , x = X l , e m = M l (2) onde θ representa a variável tempo e ϕ a variável tempo após a adimensionalização. Substituindo (2) em (1) e após algumas manipulações algébricas, temos: ṫ = t(1 − (t+ y)) − tv − η1t ẏ = ry(1− (t+ y)) + tv − γxv − η2y v̇ = βy − ηcv (3) i̇ = Ψ − ηgi− σiy ẋ = 2σiy − θ x 1 + x − ηgx+ (zy + ηn)m ṁ = 2θ x x+ 1 − ηsm− (zy + ηn)m, sendo η1 = u1 r1 , r = r2 r1 , η2 = u2 r1 , ηc = c r1 , Ψ = f θ2r1 , ηg = g r1 , σ = ak r1 , θ = q lr1 , z = z1ak r1 , ηn = n r1 , e ηs = s r1 . Escolhendo eθ1 r1 = 1 obtemos que θ1 = r1 e , γ = αl ke e β = pe r2 1 . (4) A admensionalização proposta reduz o sistema de equações (1) o qual possui dezenove parâmetros ao sistema de equações (3) o qual possui doze parâmetros. 3.3 Soluções de equilíbrio Para a obtenção das soluções de equilíbrio basta igualarmos o conjunto de equações (3) à zero. As soluções de equilíbrio são escritas da forma Pn=(t̂,ŷ,v̂,̂i,x̂,m̂), assim: 1. Se v̂ = 0 obtemos duas soluções: 18 P0 = (1 − η1, 0, 0, Ψ ηg , 0, 0) e, P1 = (1 − η1, 0, 0, Ψ ηg , θ ηg (−1 + 2ηn ηs+ηn ) − 1, 2θ (ηs+ηn) θ ηg (−1+ 2ηn ηs+ηn )−1 θ ηg (−1+ 2ηn ηs+ηn ) ) 2. Se v̂ �= 0 P2 = (1 − ηc β v̂ − v̂ − η1, ηc β v̂, v̂, Ψ (ηg+σ ηc β v̂) , x̂, 2θ ηs+z ηcv̂ β +ηn x̂ x̂+1 ). 3.4 Estudo da Estabilidade dos Pontos de Equilíbrio Tendo encontrado os pontos de equilíbrio estudaremos separadamente a estabilidade de cada um deles através da análise de seus autovalores (soluções de seus polinômios ca- racterísticos). Com este propósito construiremos a matriz de derivadas parciais de (3), também conhecida como matriz jacobiana J . 3.4.1 Estabilidade de P0 Substituindo o ponto P0 em J obtemos a matriz jacobiana calculada em P0: J0 = ⎡ ⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ −1 + η1 −1 + η1 −1 + η1 0 0 0 0 r − r(1 − η1) − η2 1 − η1 0 0 0 0 β −ηc 0 0 0 0 −σi 0 −ηg 0 0 0 2σi 0 0 −θ − ηg ηn 0 0 0 0 2θ −ηs − ηn ⎤ ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ . (5) O polinômio característico é dado por det|J0 − λI|, em que I é a matriz identidade de mesmas dimensões que J0. Devido a grande dificuldade de manipular os termos, (trata-se de uma matriz 6× 6), reescrevemos J0 − λI como uma matriz de dimensões 4× 4 através do método dos cofatores e, então aplicamos a propriedade: M = ⎡ ⎣ A B C D ⎤ ⎦ = ⎡ ⎣ A 0 C I ⎤ ⎦ ⎡ ⎣ I A−1B 0 D − CA−1B ⎤ ⎦ (6) portanto, detM = detAdet(D − CA−1B), em que A, B, C e D são matrizes de ordem 2 × 2. A resolução do determinante deste modo, nos permitiu encontrar um polinômio característico fatorado, com termos no qual λ tem grau de no máximo 2, dado por (−1)1+1(−1 + η1 − λ)(−1)3+3(−ηg − λ)Δ1Δ2, (7) 19 sendo Δ1 = (r − r(1 − η1) − η2 − λ)(−ηc − λ) − β(1 − η1) e Δ2 = (−θ − ηg − λ)(−ηs − ηn − λ) − (ηn2θ). (8) Resolvendo cada termo do polinômio encontraremos seus autovalores, dados por λ1 = −1 + η1, λ2 = −ηg, λ3 = −1 2 (rη1 − η2 − ηc ± √ r2η2 1 − 2rη1η2 + 2rη1ηc + η2 2 − 2η2ηc + η2 c + 4β − 4βη1) e, λ4 = 1 2 (−θ − ηg − ηs − ηn ± √ θ2 + 2θηg − 2θηs + 6θηn + η2 g − 2ηgηs − 2ηgηn + η2 s + 2ηsηn + η2 n), onde as soluções λ3 e λ4 foram econtradas através do software de manipulação algébrica Maple. A estabilidade dos pontos de equilíbrio P0 será avaliada partindo do princípio que um ponto estável é aquele que possui todos os seus autovalores com parte real negativa ou parte real nula (neste caso, desde que a multiplicidade algébrica e geometrica destes au- tovalores sejam iguais). Nossos autovalores são distintos, o que implica em multiplicidade algébrica m = 1 e, sendo 1 ≤ μ ≤ m, a multiplicidade geométrica de cada autovalor é μ = 1. Uma vez que todos os parâmetros são positivos e, consequentemente λ2 < 0, as restri- ções incluem apenas λ1, λ3 e λ4. • Para que λ1 tenha parte real negativa ou nula: −1 + η1 ≤ 0 =⇒ η1 ≤ 1. Os autovalores λ3 e λ4 são compostos por duas raízes reais sendo uma positiva e outra negativa. Analisaremos somente a raiz positva de cada um deles, pois, caso ela satisfaça o critério λ < 0, sua raiz negativa consequentemente também atenderá essa condição. Assim, • Para que λ3 tenha parte real negativa ou nula:(√ r2η2 1 − 2rη1η2 + 2rη1ηc + η2 c − 2η2ηc + η2 2 − 4βη1 + 4β )2 ≤ (−rη1 + η2 + ηc) 2, logo β(1 − η1) ηc(η2 − rη1) ≤ 1. 20 • Por fim, para que λ4 tenha sua parte real negativa ou nula: √ θ2 + 2θηg − 2θηs + 6θηn + η2 g − 2ηgηs − 2ηgηn + η2 s + 2ηsηn + η2 n ≤ (θ + ηg + ηs + ηn)2, teremos 1 ≤ θηs + ηgηs + ηgηn θηn . Definindo R0 = β(1−η1) ηc(η2−rη1 e, R1 = θηs+ηgηs+ηgηn θηn , obtemos que, se R0 ≤ 1, R1 ≤ 1 e η1 < 1 a solução convergirá ao ponto de equilíbrio P0. 3.4.2 Estabilidade de P1 Do ponto de vista biológico, para que as células imunes ativas existam devemos considerar x̂ > 0 =⇒ θ ηg ( ηn − ηs ηs + ηn ) > 1. Definindo R2 = θ ηg ηn−ηs ηs+ηn , temos que x̂ = R2 − 1 e, portanto P1 pode ser escrito como P1 = (1 − η1, 0, 0, Ψ ηg , R2 − 1, 2θ (ηs + ηn) R2 − 1 R2 ). Utilizando os mesmos conceitos empregados para estudar a estabilidade de P0, a estabi- lidade do ponto de equilíbrio P1 será analisada. A matriz Jacobiana calculada em P1 é dada por: J1 = ⎡ ⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ −1 + η1 −1 + η1 −1 + η1 0 0 0 0 r − r(1 − η1) − γ(R2 − 1) − η2 1 − η1 0 0 0 0 β −ηc 0 0 0 0 −σ Ψ ηg 0 −ηg 0 0 0 2σ Ψ ηg + 2zθ(R2−1) (ηs+ηn)R2 0 0 −θ R2 + θ(R2−1) R2 2 − ηg ηn 0 − 2zθ(R2−1) (ηs+ηn)R2 0 0 2θ R2 − 2θ(R2−1) R2 2 −ηs − ηn ⎤ ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ , (9 sendo o polinômio característico dado por det|J1 − λI|: (−1)1+1(−1 + η1 − λ)(−1)3+3(−ηg − λ)Δ3Δ4, onde Δ3 = (r − r(1 − η1) − γ(R2 − 1) − η2 − λ)(−ηc − λ) − β(1 − η1) e, Δ4 = ( −θ R2 + θ(R2 − 1) R2 2 − ηg − λ)(−ηs − ηn − λ) − ηn( 2θ R2 − 2θ(R2 − 1) R2 2 ). 21 Resolvendo o polinômio através do software de manipulação algébrica Maple encontrare- mos seus autovalores, dados por: λ5 = η1 − 1, λ6 = −ηg, λ7 = 1 2 (rη1 − γR2 + γ − η2 − ηc ± √ Υ1), λ8 = 1 2 1 R2 2 (−θ − ηgR 2 2 −R2 2ηs −R2 2ηn ± √ Υ2), em que Υ1 = −2rη1γR2 + r2η2 1 + 2rη1γ − 2rη1η2 + 2rη1ηc + γ2R2 2 − 2γ2R2 + 2γR2η2 −2γR2ηc + γ2 − 2γη2 + +2γηc + +η2 2 − 2η2ηc + η2 c + 4β(1 − η1) e, Υ2 = θ2 + 2θηgR 2 2 − 2θR2 2ηs + 6θR2 2ηn + η2 gR 4 2 − 2ηgR 4 2ηn + +R4 2η 2 s + 2R4 2ηsηn +R4 2η 2 n. Todos os parâmetros são positivos, portanto λ6 < 0. Com relação aos outros autovalores temos que: • Para que λ5 tenha sua parte real negativa ou nula: η1 − 1 ≤ 0 =⇒ η1 ≤ 1; (10) • Para que λ7 tenha sua parte real negativa ou nula: rη1 − γR2 + γ − η2 − ηc + √ Υ1 ≤ 0 então temos 1 ≤ ηc(η2 − γ − rη1 + γR2) β(1 − η1) ; • Para que λ8 tenha sua parte real negativa ou nula: −θ − ηgR 2 2 − R2 2ηs −R2 2ηn + √ Υ2 ≤ 0 tal que, θ(ηn − ηs) R2 2ηg(ηs + ηn) ≤ 1. Definindo R3 = ηc(η2−γ−rη1+γR2) β(1−η1) e R4 = θ(ηn−ηs) R2 2 ηg(ηs+ηn) , se, η1 ≤ 1, R3 ≥ 1 e R4 ≤ 1, o sistema convergirá ao ponto de equilíbrio P1. 3.4.3 Estabilidade de P2 A análise da estabilidade de P2 só é viável se realizada numericamente, portanto o resul- tado obtido para este ponto está descrito na proxima secção. 22 3.5 Resultados O resumo da análise da estabilidade dos pontos de equilíbrio do modelo matemático adimensional da dinâmica viral do HBV está colocado na Tabela 2. Tabela 2: Condições de estabilidade dos três pontos de equilíbrio referentes ao modelo dinâmico viral-imunológico do HBV Pontos de Equilíbrio η1 R0 R1 R2 R3 R4 P0 < 1 ≤ 1 ≥ 1 ≤ 1 − − P1 < 1 − < 1 > 1 ≥ 1 ≤ 1 P2 < 1 − < 1 > 1 < 1 ≤ 1 Pela Tabela 2 observamos que se R0 = β(1−η1) ηc(η2−rη1) ≤ 1 e R1 = θηs+ηgηs+ηgηn θηn ≥ 1 a solu- ção converge para o equilíbrio trivial onde não há população viral V , nem hepatócitos infectados Y , sistema imune ativo X ou células de memória M . O parâmetro R0 pode ser interpretado da seguinte forma: β e r relacionam-se respec- tivamente com o crescimento das populações de vírus HBV e dos hepatócitos infectados; enquanto que, ηc, η1 e η2 são termos associados respectivamente à mortalidade do vírus HBV , dos hepatócitos saudáveis e dos hepatócitos infectados. O termo 1 − η1 é o res- ponsável pela existência ou não das células saudáveis no fígado. Esse termo aparece como autovalor para os três pontos de equilíbrio estudados η1 < 1, e reflete o fato de que não temos paciente vivo com o fígado totalmente composto por células infectadas. Já o termo η2 − rη1 traz uma perspectiva em torno da eliminação da infecção. Se η2 > rη1, siginifica que a infecção decresce com o tempo, sendo que a taxa de mortalidade dos hepatócitos infectados é maior do que o produto da mortalidade dos hepatócitos saudáveis pelo cres- cimento dos hepatócitos infectados. Por fim o termo β ηc reflete o descrescimento com o tempo da população de HBV quando ηc > β, a taxa de mortalidade do vírus é maior que a de seu crescimento. Por sua vez, em R1 temos ηn representando a taxa de ativação das células ativas em células de memória; enquanto que, ηs e ηg são respectivamente, as taxas de mortalidade das células de memória e imune ativa. Uma vez que R2 é a condição biológica que determina a existência do sistema imunológico a partir de células imunes ativasX, podemos estabelecer uma relação entre R2 e R1. Rearranjando R2 como R2 (ηs−ηn) = θ ηg 1 (ηs+ηn) e, invertendo os dois lados da equação teremos (ηs−ηn) R2 = ηg θ (ηs + ηg). Podemos reescrever R1 como 23 R1 = ηs ηn + ηg θηn (ηs +ηn) e substituir nele a nova expressão de R2, R1 = ηs ηn + ηs−ηn R2 . Fazendo algumas manipulações algébricas, chegamos à R1 = 2ηs ηn − 1. Logo, para que R1 satisfaça o ponto P0 devemos ter 2ηs ηn − 1 ≥ 1, que resulta em ηs ≥ ηn. Este resultado siginifica que a taxa de ativação não é suficiente para manter a população das células memória. Consequentemente, no ponto P0 não há formação de resposta imune. Figura 3 mostra a evolução temporal das populações, dado o conjunto de valores que obedece às quatro situações discutidas até aqui, i.e., η1 < 1, η2 > rη1, ηc > β e, ηs > ηn. A evolução temporal das populações de células e vírus foram obtidas através da implementação em linguagem C do modelo adimensional via um runge-kutta de quarta ordem. Biologicamente o ponto P0 está associado à indivíduos que nunca entraram em contato com o vírus ou indivíduos que entraram em contato com o vírus e foram capaz de eliminá-lo sem ativar a resposta imunológica adaptativa. Utilizou-se r = 0.2, γ = 10.0, ηs = 0.3, ηc = 0.6, ηg = 0.14, ηn = 0.2, η1 = 0.4, η2 = 2.0, β = 0.2, ψ = 0.06 e θ = 0.3. 0 0,05 0,1 0,15 t 0,3 0,4 0,5 0,6 0,7 P( t) Figura 3: Evolução temporal das células segundo o modelo adimensional convergindo para o ponto de equilíbrio P0 com R0 ≈ 0.10 e R1 ≈ 2.67. Em vermelho temos T̂ = 0.6, em azul escuro temos Î ≈ 0.43 e as outras curvas, não representadas no gráfico, vão a zero. Novamente da Tabela 2 observamos que se R1 = θηs+ηgηs+ηgηn θηn < 1, R2 = θ ηg ηs−ηn ηs+ηn > 1, R3 = ηc(η2−γ−rη1+γR2) β(1−η1) ≥ 1 e R4 = θ(ηn−ηs) R2 2 ηg(ηs+ηn) ≤ 1 a solução converge para o equilíbrio no qual não há população viral HBV (V ) nem hepatócitos infectados Y mas há sistema imune ativo X e células de memória M . Partindo da mesma análise que fizemos para P0, dado que R1 = 2ηs ηn − 1, para que a condição dada por R1 < 1 seja satisfeita, devemos ter 24 ηs < ηn, caracterizando o aparecimento e a sustentação no equilíbrio de populações de células imunes ativas X e criação de memória imunológica M , mesmo sem a existência de vírus HBV e de hepatócitos infectados Y . Se a ativação de background das células de memória for maior que uma combinação das taxas de mortalidade dessas células M e das células imunes efetoras X, essas células X são ativadas mesmo na ausência do vírus específico. Esta situação surge quando as células de memória M são convertidas em linfócitos efetores, células imunes ativas X (parâmetro ηn). Em particular, se a taxa de ativação de background for nula, encontramos: P1 = (1 − η1, 0, 0, Ψ ηg , −θ ηg − 1, 2θ 1 ηs −θ ηg − 1 −θ ηg ), (11) (população negativa) e, portanto, podemos concluir que sob essa condição (ηn = 0) só é possível obter os pontos de equilíbrio P0 e P2 [7]. Quanto à R4, podemos reescrevê-lo em função de R2, obtendo R4 = 1 R2 , ou seja, este critério também está intimamente relacionado ao surgimento de resposta imunológica: basta que R2 > 1, para que a condição R4 ≤ 1 seja atendida. O Ponto P1 com ηn > 0 representa duas típicas situações relacionadas à hepatite B : • Um indivíduo que entrou em contato com o vírus HBV , desenvolvendo a doença mas conseguiu eliminá-la, adquirindo imunidade pelo resto de sua vida; caracteriza o portador agudo da doença ou, • Imunidade induzida pela vacinação; Este quadro pode ser observado na Figura 4, que mostra a evolução temporal das popula- ções com eliminação do vírus da hepatite e com ganho de imunidade. Utilizou-se r = 0.2, γ = 0.3, ηs = 0.04, ηc = 0.4, ηg = 0.16, ηn = 0.8, η1 = 0.46, η2 = 0.23, β = 0.12, ψ = 0.27 e θ = 0.27. O critério dado por R3 determina se após o desenvolvimento da infecção ha- verá eliminação das populações de vírus HBV e hepatócitos infectados Y (características do portador agudo) ou a solução tenderá à coexistência: quando R3 ≥ 1 a solução se es- tabiliza em P1 enquanto que se R3 < 1 tenderá ao ponto P2 (portador crônico). Portanto, se R1 < 1, R2 > 1, R3 < 1 e R4 ≤ 1 a solução converge ao ponto de equilíbrio no qual existem todas as populações. Neste caso, o que mantém a população de células imunes efetoras e das células de memória não é a taxa de ativação de background ηn, mas sim a 25 0 0,1 0,2 0,3 0,4 t 0 0,5 1 1,5 2 2,5 3 P( t) 0 0,1 0,2 t 0 0,1 0,2 P( t) Figura 4: Evolução temporal das células segundo o modelo adimensional convergindo para o ponto de equilíbrio P1 com R1 ≈ 0.43, R2 ≈ 1.51, R3 ≈ 1.79 e R4 ≈ 0.663. Em vermelho temos T̂ = 0.54, em azul escuro temos Î ≈ 1.67, em verde temos X̂ ≈ 0.51, em azul claro temos M̂ ≈ 0.21 e as outras duas curvas vão para zero, sendo mostradas na figura menor. própria presença dos vírus HBV e das células infectadas que permanecem no organismo [7]. O ponto P2 representa o portador crônico, o qual nunca consegue curar-se da infecção. Esta situação pode ser observada na Figura 5. Utilizou-se r = 0.71, γ = 0.57, ηs = 0.03, ηc = 0.29, ηg = 0.1144, ηn = 0.57, η1 = 0.328, η2 = 0.23, β = 0.67, ψ = 0.07 e θ = 0.29. 0 0,2 0,4 0,6 0,8 1 t 0 0,5 1 1,5 2 P( t) Figura 5: Evolução temporal das células segundo o modelo adimensional convergindo para o ponto de equilíbrio P2 com R1 ≈ 0.47, R2 ≈ 2.26, R3 ≈ 0.45 e R4 ≈ 0.44. Em vermelho temos T̂ ≈ 0.50, em preto temo Ŷ ≈ 0.05, em laranja temos V̂ ≈ 0.12, em azul escuro temos Î ≈ 0.31, em verde temos X̂ ≈ 2.21 e, em azul claro temos M̂ ≈ 0.48. 26 4 Modelo Dinâmico Viral Espacial Nesta seção desenvolveremos um modelo espacial bidimensional, a partir do modelo apre- sentada na seção anterior, seguindo o desenvolvimento apresentado no artigo de Funk [2], isto é, considerando acoplamentos espaciais entre as variáveis do modelo não espacial. Como a maioria das infecções tem uma localização espacial e as interações entre as po- pulações de vírus e células são locais, modelos com estrutura espacial são um importante passo para entendermos o que acontece quando a suposição de populações bem misturadas é relaxada. 4.1 Modelo Matemático Espacial da Dinâmica Viral do HBV Seguindo a idéia proposta por Funk [2], construímos a versão espacial do modelo descrito na seção anterior acrescentando os índices i e j a cada uma das populações (T, Y, V, I,X,M), de maneira que as mesmas tenham uma posição de coordenada (i,j ) definidas numa rede quadrada. Sendo assim, Ṫi,j = Ti,j(1 − (Ti,j + Yi,j)) − Ti,jVi,j − η1Ti,j Ẏi,j = rYi,j(1 − (Ti,j + Yi,j)) + Ti,jVi,j − γXi,jVi,j − η2Yi,j V̇i,j = βYi,j − ηcVi,j − mv 8 i+1∑ i0=i−1 j+1∑ j0=j−1 (Vi,j − Vi0,j0) İi,j = Ψ − ηgIi,j − σIi,jYi,j (12) Ẋi,j = 2σIi,jYi,j − θ Xi,j 1 +Xi,j − ηgXi,j + (zYi,j + ηn)Mi,j − mx 8 i+1∑ i0=i−1 j+1∑ j0=j−1 (Xi,j −Xi0,j0) Ṁi,j = 2θ Xi,j 1 +Xi,j − ηsMi,j − (zYi,j + ηn)Mi,j. Deve ser ressaltado que com relação ao modelo não espacial alteramos a notação das variáveis t, y, v, i, x,m de minúsculas para maiúsculas, evitando confundir com os novos parâmetros de locação espacial da matriz i, j e com o termo de difusão dado por mv e mx. A vizinhança de um sítio é uma coleção N de elementos usada para determinar o estado σi,j ou, neste caso, a densidade de população do sítio alvo (i, j). Ela pode ser local (elementos adjacentes) ou não-local, pode ser previamente definida ou arbitrária (randômica), pode conter o sítio alvo ou não. Uma vizinhança pode assumir qualquer 27 configuração espacial e sua escolha depende do processo a ser modelado. A vizinhança de Moore de raio 1 é composta pelos vizinhos adjacentes ao sítio, i. e., N ≡ {σi−1,j−1, σi−1,j , σi−1,j+1, σi,j−1, σi,j+1, σi+1,j−1, σi+1,jσi+1,j+1}. Os novos termos acrescentados a equação do vírus e das células imunes ativas de- terminam o tipo de acoplamento espacial considerado no modelo. Os parâmetros mν e mx denotam as taxas de difusão, respectivamente, do vírus (HBV ) e das células imunes ativas para sítios adjacentes. Estamos supondo que as demais células não se movem, o que é uma simplificação do modelo. A dinâmica da infecção será estudada numa rede quadrada L × L, representando o local onde o processo de infecção acontece, ou seja, no caso da hepatite B representa o fígado. Segundo Funk [2], instabilidades devido a acoplamentos espaciais não acontecem, de maneira que as diferenças observadas entre os dois modelos (com e sem acoplamento espacial) devem ser perceptíveis nos primeiros instantes da evolução da infecção (compor- tamento transiente). A dinâmica do processo de infecção com e sem estrutura espacial, isto é, a intensidade e velocidade de espalhamento da infecção serão comparadas. Estu- daremos também como a intensidade dos acoplamentos (mν e mx) afeta a dinâmica do processo infeccioso. 4.2 Metodologia Para o modelo espacial utilizaremos o mesmo método de resolução e análise do modelo não adimensional (Runge-kutta), no entanto com algumas diferenças relativas as imple- mentações das rotinas computacionais: no modelo anterior resolviamos um sistema de Equações Diferencias Ordinárias (EDO), agora precisamos resolver este mesmo sistema em cada sítio da rede lembrando que eles interagem entre si. Escolhemos uma malha (matriz) bidimensional 301 × 301, onde cada sítio (i, j) da matriz é visto como um vetor de seis elementos: T, Y, V, I,X e M , os quais correspondem as populações estudadas. Utilizamos condições de contorno abertas e a atualização é feita de modo assíncrona. A cada iteração, sorteia-se aleatoriamente um sítio qualquer da rede (i, j) para atualiza- 28 lo. No caso das populações de vírus e sistema imune ativo precisamos analisar primeiro o processo de difusão, o qual ocorre entre os sítios vizinhos: a quantidade de população que se difunde depende do gradiente de concentração entre o sítio analisado e seus vizinho, por exemplo, se o sítio (i, j) tiver concentração de vírus igual a 0.5 e o sítio (i, j+1) tiver 0.1, então a quantidade 0.4/8×mv deverá ser somada a (i, j+1) e subtraída do sítio que descendeu (i, j), para respeitarmos a condição de conservação de massas dada pela Lei de Lavoisier. Depois a evolução temporal da população é feita pelo Runge-kutta. Esse processo ocorrerá até que tenhamos sorteado uma quantidade suficiente de posi- ções, dada pelo número de sítos (L×L). Passamos para uma segunda iteração onde todo processo descrito anteriormente é retomado, só que agora para um passo de tempo t+ 1. Com relação aos resultados apresentados, mostramos a dispersão viral na matriz e a evolução temporal em função do tempo. A discussão terá como alicerce a análise da dinâmica viral em torno da estabilidade dos três ponto de equilíbrio obtidos para o modelo sem acoplamento espacial: P0, P1 e P2. 29 4.3 Resultados e Discussão 4.3.1 Acoplamento Espacial Homogêneo Nesta subseção apresentaremos os resultados obtidos para os pontos P1 e P2 com difusão. A análise referente ao ponto de equiíbrio P0 é trivial e portanto optamos por não apresentá- la. Os resultados serão apresentados na forma de dois tipos de gráficos. O primeiro avaliará a dispersão viral na matriz, ou seja a capacidade do vírus ser difundido. Para esta análise da dispersão consideramos a intesidade viral de acordo com as cores apresentadas na Tabela 3 : Tabela 3: Código de Cores para a intensidade viral durante a dispersão. Cor Intensidade Viral Correspondente preto V ≥ 0.20 vermelho 0.13 ≤ V < 0.20 amarelo 0.08 ≤ V < 0.10 verde 0.05 ≤ V < 0.08 branco V ≤ 0.05 (13) O segundo gráfico avaliará a evolução temporal das populações em função do tempo no sítio central i, j, em um segundo sítio i+ 60; j + 60 e um terceiro sítio i+ 140; j + 140 (a matriz de trabalho é uma 301 × 301). Assumiremos que o ambiente é homogêneo, de maneira que os parâmetros do modelo são os mesmos em todos os pontos da rede e que somente o vírus e as células imunes ativas se dispersam (dispersão local). No caso do ponto de estabilidade P1, a infeção viral é introduzida no centro da matriz, nos vetores V dos sítios de (i, j) até (i− 50, j+ 50), (i, j+ 50), (i+ 50, j+ 50); (i− 50, j), (i+50, j) e (i−50, j−50), (i, j−50), (i+50, j−50); de forma que conseguimos obter uma grande parcela da malha com vírus inicialmente, simulando um extenso processo de infec- ção. A intenção simples é notarmos que independente do tamanho da infecção, quando trabalhamos com P1 ela será eliminada, ou seja, a eliminação viral ocorrerá certamente e independente das condições iniciais. Todos os outros sítios contém inicialmente em seus 30 Figura 6: Dispersão da população viral com acoplamento espacial (difusão) correspon- dente ao estado estacionário P1, para mv = mx = 0.1. Da esquerda para a direita temos as iterações t = 2, 3, 4, 5 e 6. vetores T = I = 1 e Y = V = X = M = 0. A Figura 6 e a Figura 7 mostram a evolução temporal referente ao ponto de equilíbrio P1 para o modelo com acoplamento (mv = mx = 0.1). A solução foi construída conver- gindo para o ponto de equilíbrio P1 com R1 ≈ 0.43, R2 ≈ 1.51, R3 ≈ 1.79 e R4 ≈ 0.663 e com os parâmetros r = 0.2, γ = 0.3, ηs = 0.04, ηc = 0.4, ηg = 0.16, ηn = 0.8, η1 = 0.46, η2 = 0.23, β = 0.12, φ = 3.0, ψ = 0.27, θ = 0.27 e z = 6.0. Na Figura 6 a intensidade viral ≤ 0.05 ocorre após seis passos de tempo, havendo eliminação total ao termínio das iterações. Como já era esperado, uma vez que utilizamos um modelo homogêneo em que todos os sítios recebem idênticos valores de parâmetros que os levarão ao mesmo ponto de equilíbrio, a eliminação da infecção torna-se também homogênea, o que pode ser observado pelos padrões de cores bem distribuídas que foram formados. Enquanto isso, através da análise da Figura 7 observamos que a eliminação da infeção ocorre antes que o vírus chegue à borda da malha e, quanto mais afastado o sítio analisado estiver do ponto central da matriz menor a intensidade viral o atinge. Em vermelho temos T̂ = 0.54, em azul escuro temos Î ≈ 1.67, em verde temos X̂ ≈ 0.51, em azul claro temos M̂ ≈ 0.21 e as outras duas curvas vão para zero. Como já era previsto por Funk [2] as diferenças no comportamento das curvas são perceptíveis somente no início da evolução temporal (comportamento transiente) e quando comparamos os resultados com o gráfico de P1 obtido para o modelo não espacial concluimos que o valor de equilíbrio de cada população não foi alterado. Para a análise do ponto de equilíbrio P2 a infecção viral é introduzida apenas no vetor V do sítio cetral (i = L 2 , j = L 2 ) da matriz, em que L × L é o tamanho da malha, sendo 31 0 100 200 300 0 0,5 1 1,5 2 2,5 3 0 100 200 300 0 0,2 0,4 0,6 0,8 1 0 100 200 300 0 0,5 1 1,5 2 P( t) 0 100 200 300 0 0,001 0,002 0 100 200 300 t 0 0,5 1 1,5 2 0 100 200 300 t 0 0,5 1 Figura 7: Evolução temporal das células segundo o modelo adimensional espacial, convergindo para o ponto de equilíbrio P1. De cima para baixo temos os sítios 150×150, 210×210 e 290×290 L = 301. A idéia neste caso é iniciar com uma quantidade aparentemente insignificante de sítios infectados, para que possamos notar que independente do tamanho da infecção, quando trabalhamos com P2, ela será difundida por toda a malha, ou seja, os parâmetros adotados irão caracterizá-la como persistente. Todos os outros sítios contém inicialmente em seus vetores T = I = 1 e Y = V = X = M = 0. A Figura 8 e a Figura 9 apresentam a evolução temporal das células segundo o modelo adimensional convergindo para o ponto de equilíbrio P2 com R1 ≈ 0.47, R2 ≈ 2.26, R3 ≈ 0.45 e R4 ≈ 0.44. A solução foi construída para r = 0.71, γ = 0.57, ηs = 0.03, ηc = 0.29, ηg = 0.1144, ηn = 0.57, η1 = 0.328, η2 = 0.23, β = 0.67, φ = 2.14, ψ = 0.07, θ = 0.29 e z = 4.29. Em vermelho temos T̂ ≈ 0.50, em preto temo Ŷ ≈ 0.05, em laranja temos V̂ ≈ 0.12, em azul escuro temos Î ≈ 0.31, em verde temos X̂ ≈ 2.21 e, em azul claro temos M̂ ≈ 0.48. Em P2 também observamos que após as primeiras iterações (Figura 9) as soluções convergem para os valores de equilíbrio obtido no modelo compartimental e, a dispersão da infeção pela malha ocorre homogeneamente, atingindo todos os sítios da matriz incluindo os mais próximos das bordas (Figura 8). O efeito da difusão sobre o modelo, pode ser melhor estudado nas Figura 10 e a Figura 32 Figura 8: Disperão da população viral com acoplamento espacial (difusão) para o estado estacionário P2. Da esquerda para a direita temos as iterações de número t = 10, 30, 60, 110 e 180. 0 100 200 300 0 0,5 1 1,5 2 2,5 3 0 100 200 300 0 0,05 0,1 0,15 0,2 0,25 0 100 200 300 0 0,5 1 1,5 2 2,5 3 P( t) 0 100 200 300 0 0,1 0,2 0,3 0,4 0 100 200 300 t 0 0,5 1 1,5 2 2,5 3 0 100 200 300 0 0,1 0,2 0,3 0,4 Figura 9: Evolução temporal das células segundo o modelo adimensional espacial, convergindo para o ponto de equilíbrio P2. De cima para baixo temos os sítios 150×150, 210×210 e 290×290. 11. A Figura 10 apresenta o tempo de eliminação do vírus no sítio central (L 2 , L 2 ) e a Figura 11 o tempo de chegada do vírus aos sítios da borda (L 2 + 140, L 2 + 140), ambos em função da taxa de difusão viral mv igual a taxa de difusão de células de imunes ativas mx. Para estes dois gráficos notamos que com o aumento da difusão a velocidade de eliminação viral no sítio central e a velocidade de chegada do vírus na borda da malha aumentam, mas de maneira não proporcional, uma vez que quanto maior a difusão mais suave a inclinação da curva se torna. Em seguida, geramos mais duas curvas, nas quais mantivemos o valor demx fixo em 0.5 e variamos apenas mv, mantendo todas as outras condições iguais as das Figuras 10 e 11. Deste modo poderíamos estabelecer alguma relação de comportamento das velocidades 33 0 0,2 0,4 0,6 0,8 1 Taxa de difus. 0 10 20 30 40 50 60 70 80 90 100 110 T em po d e el m in . v ir al mv Figura 10: Análise da difusão para o ponto de equiíbrio P1. Tempo de eliminação viral no sítio central em função da taxa de difusão mv (mv = mx). 0 0,2 0,4 0,6 0,8 1 Taxa de difus. 0 20 40 60 80 100 120 140 160 180 200 T em po d e ch eg ad a vi ra l mv Figura 11: Análise da difusão para o ponto de equiíbrio P2. Tempo de chegada do vírus no sítio da borda (L 2 + 140, L 2 + 140) em função da taxa de difusão mv (mv = mx). de eliminação e de chegada viral com os valores assumidos pelas taxas de difusão quando mv ≤ mx e mv > mx. Os resultados obtidos para esta situação foram exatamente os mesmos, de modo que as curvas obtidas eram idênticas as Figuras 10 e 11, o que nos leva a concluir que para a ordem de grandeza das taxas de difuão com que trabalhamos, o tempo de eliminação viral e o tempo de chegada do vírus não dependem da taxa de difusão de células imunes. Em resumo, quando a escolha dos parâmetros é tal que a situação estudada estiver associada á eliminação viral (P1), isto ocorrerá de maneira mais rápida em relação ao modelo compartimental e, da mesma forma se a escolha dos parâmetros estiver associada ao equilíbrio endêmico, i.e., onde as polulações de células e virus coexistem (P2), a ve- 34 locidade com que este evento ocorrerá também aumentará. Quanto as taxas de difusão, podemos afirmar que o modelo matemático é sensível apenas á variação de mv (difusão viral). 4.3.2 Acoplamento Espacial Heterogêneo Até a seção anterior estavamos trabalhando com uma malha de parâmetros homogêneos, na qual cada sítio, ao fim das iterações terminaria no mesmo ponto de equilíbrio. Nesta seção trabalharemos com uma malha de parâmetros heterogêneos, repartindo nossa matriz em quatro quadrantes: o superior esquerdo e o inferior direito com seus sítios contendo os parâmetros relativos ao ponto de equilíbrio P1 e, o superior direito e o inferior esquerdo, com os sítios com os parâmetros realtivos a P2, de modo que cada referido quadrante converge para seu designado ponto de equilíbrio, atribuindo-se no runge-kutta valores que atendessem as respectivas condições para R0, R1 e R2. Inicialmente inserimos na malha uma infecção extensa exatamente como havíamos realizado para o ponto de equilíbrio P1 para a malha homogênea. O padrão de cores apresentado no gráfico de dispersão viral respeita os limites estabelecidos pela Tabela 3 de Código de Cores para a Intensidade Viral Durante a Dispersão. Na Figura 12 observamos que apesar da malha ser heterôgenea não houveram dife- renças significativas na dispersão viral, cada quadrante resolveu individualmente a infec- ção para seus próprios parâmetros sem qualquer interferência dos quadrantes adjacentes. Portanto, como já esperavamos somente os parâmetros atribuídos a cada sítio podem determinar a evolução ou a eliminação da infecção. Figura 12: Dispersão da população viral com acoplamento espacial (difusão) em malha heterogênea. Da esquerda para a direita temos as iterações de número t = 0, 5, 18, 53 e 119. 35 0 100 0 0,5 1 1,5 2 2,5 0 100 0 0,2 0,4 0 100 0 0,5 1 1,5 2 2,5 P( t) 0 100 0 0,2 0,4 0 100 0 0,5 1 1,5 2 2,5 0 100 0 0,4 0,8 0 100 0 0,5 1 1,5 2 2,5 0 100 0 0,001 0,002 0 100 t 0 0,5 1 1,5 2 P( t) Figura 13: Evolução temporal das células segundo o modelo adimensional espacial. De cima para baixo temos os sítios 10 × 290, 90 × 210, com parâmetros de E2, sítio central, 210 × 210 e 290 × 290, com parâmetros de E1 . Através da Figura 13 notamos que no modelo heterogêneo , cada quadrante também converge para seus respectivos pontos de equilíbrio em relação ao modelo não espacial e as populações estabilizam-se nos mesmos valores. Comparando-se este gráfico com as Figura 8 e 9 notamos que não houve variação do tempo de eliminação (P1) ou tempo de chegada viral (P2). Na Figura 13 deve ser ressaltado que como a infecção não atingiu o sítio de 290 × 290, sendo eliminada antes disso no seu quadrante, o gráfico que deveria conter as poluações de vírus e hepatócitos infectados correpondentes a este sítio não foi construído. 36 5 Conclusão Os resultados obtidos para modelo compartimental o qual descreve a dinâmica da infecção pelo vírus da hepatite B e sua resposta imunológica adaptativa podem ser resumidos na Tabela 4 onde indicamos cada uma das soluções de equilíbrio e os valores limiares que determinam a estabilidade de cada uma delas. Tabela 4: Condições de estabilidade relevantes Pontos de Equilíbrio R0 R1 R3 P0 ≤ 1 ≥ 1 − P1 − < 1 ≥ 1 P2 − < 1 <1 As soluções P0, P1 e P2 indicam, respectivamente, indivíduos que eliminam a infecção sem formar resposta imunológica, portadores agudos e portadores crônicos, de forma que o modelo reproduz as possíveis evoluções da doença. Ao incluirmos difusão (modelo espacial), tanto no modelo heterogêneo como no homo- gêneo, os tempos de eliminação viral para P1 e chegada viral para P2 tornam-se menores em relação ao modelo compartimental, ou seja, a velocidade de eliminação ou invasão do vírus é aumentada. Para valores da ordem de grandeza das taxas de difusãomv emx com que trabalhamos, somente a taxa de difusão viral mv influi na dinâmica do modelo, mostrando que este modelo matemático não é sensível à mx (difusão associada ás células de memória). O modelo heterogêneo quando comparado ao modelo homogêneo, mostra que podemos ter a infecção limitada espacialmente dependendo do tipo de célula envolvida no processo de infecção, i.e., células mais permeáveis ou menos permeáveis a invasão do vírus. 37 Referências [1] ABBAS,K.,FAUSTO, N. Robbins Patologia-Bases Patológicas. 8ed. Rio de Ja- neiro: Guanabara Koogan, 2008. [2] FUNK, G.A.; et al. Spatial models in virus-immune dynamics. Journal of Theore- tical Biology, v.233, p.221-236, nov.2004. [3] GASPAR, A.M.C. et al. Hepatites Virais-Triagem e Diagnóstico sorológico em unidades hemoterápicas.2ed.Brasília: Ministério da Saúde, Coordenação Nacional de Doenças Sexualmente Transmissíveis e Aids. 1998. Série TELELAB. Acessado em 08 de fevereiro de 2009. Disponível em bvsms.saude.gov.br/bvs/publicacoes/cd0710.pdf [4] JUNG, M.C., PAPE,G.R. Immunology of Hepatitis B infection. The Lancet: In- fectious Diseases, v.2, p.43-49,jan.2002. [5] KALIL, A.N., COELHO, J., STRAUSS, E. Fígado e Vias Biliares-Clínica e Cirurgia. 1ed. Rio de Janeiro: Livraria e Editora RevinteR Ltda, 2001. [6] MIGUEL, J.C.Avaliação da Imunogenicidade de Vacina Recombinante Bra- sileira (Butag) cotra Hepatite B entre crianças do Rio de Janeiro. 2005. 100f. Tese de Mestrado (Mestre em Biologia Parasitária). Instituto Oswald Cruz, Fundação Oswaldo Cruz - Ministério da Saúde.Rio de Janeiro. [7] MUNIZ, P.G. Estudo da Eficiêncida da Resposta Imunológica Adaptativa. 2007. 47f. Trabalho de Conclusão de Curso (Bacharelado em Física Médica). Intituto de Biociências, Universidade Estadual Paulista Júlio Mesquita Filho, Botucatu. [8] NOWAK,M. A., MAY,R.M.Vyrus Dynamics: Mathematical Principles of Im- munology and Virology.1.ed.Oxford: Oxford University Press, 2000. v.1. 44-51p. [9] PAIXÃO,J.B.A.Gastroenterologia: Hepatites- Aspectos epidemiológicos das hepatites virais.1.ed. Rio de Janeiro: Sociedade de Gastroenterologia (livraria Ru- bio), 2001. [10] RUGGIERO,Márcia.A.Gomes. Cálculo Numérico: Aspectos Teóricos e Com- putacionais.1ed. São Paulo:McGraw-Hill,1988. Capa Folha de Rosto Agradecimentos Epígrafe Resumo Abstract Sumário Lista de Figuras 1 Introdução 2 Objetivos 3 Modelo Dinâmico Viral Sem Acoplamento Espacial 4 Modelo Dinâmico Viral Espacial 5 Conclusão Referências