Análise da Estabilidade de um Problema em Imuno-oncologia: uma Abordagem Teórica Ampliada ^*

RESUMO

O câncer é uma questão de prioridade pública que aflige o mundo e muitos esforços através da pesquisa científica estão sendo realizados para o seu combate. Nesse sentido, a imunoterapia, como tratamento em imuno-oncologia, é considerada como modalidade terapeutica praticada nas duas últimas décadas. Neste trabalho, estuda-se o crescimento das células tumorais levando em consideração o microambiente determinado pela interação que háentre as células tumorais com as células efetoras, citocinas anti-inflamatórias e um fator imuno-supressivo. Apresentam-se duas variantes dos modelos matemáticos de⁽¹⁾ com a inserção de um termo switching^{(11), (12), (15)} nesses modelos. Faz-se o estudo qualitativo dos modelos e com a análise de estabilidade e as simulações numéricas dos mesmos ilustram-se de uma forma ampliada os resultados desta pesquisa.

Palavras-chave:
modelagem matemática; imuno-oncologia; imunoterapia; existência e unicidade; análise de estabilidade

ABSTRACT

This work studies the growing of tumoral cells taking into account the microenvironmental between effector cell, anti inflammatory cytokines and the immunosuppressive factor. It is presented two variants for the math models of⁽¹⁾ by inserting a switching term on these. Also, it is done a qualitative study of the models and by the stability analysis and the numerical development work, we showed the results of this research.

Keywords:
mathematical modelling; imunotherapy; existence and uniqueness; stabilty analysis

1 INTRODUÇÃO

O câncer é uma doença que, na última década, tem dado muita preocupação nos governos, sendo uma questão prioritária na saúde pública. A Agência Internacional para pesquisa em câncer da OMS - Organização Mundial da Saúde, estimou que houve 14,2 milhões de casos novos e 8,2 milhões de mortes anuais no mundo todo em 2012. Estimou-se que em 2015 esse número chegou a 9 milhões de casos novos e 13,2 milhões de mortes e para o ano 2025 o número de casos novos chegaráa 20 milhões. Nos Estados Unidos, as estatísticas apontam que 1600 pessoas morrem por dia por causa do câncer (representou 23% do total de mortes em 2010) e na Europa esses números não são muito diferentes. Dados fornecidos pelo Instituto Nacional do Câncer - INCA apontam que no Brasil, em 2016, se estima a ocorrência de 595 mil casos novos da doença (quase 300 mil casos do sexo feminino contra 295 mil do masculino). O câncer de próstata é o de maior incidência no sexo masculino (28,6% das neoplasias em homens) e no sexo feminino é o câncer de mama (28,1% das neoplasias em mulheres)⁴.

Denomina-se como câncer um conjunto de quase 100 doenças que afetam várias partes e órgãos do corpo, tendo como característica principal o crescimento descontrolado de células anormais e pouco diferenciadas com capacidade invasiva, incluindo à distancia (metástase). Nos dias atuais, o arsenal terapêutico contra o câncer tem-se ampliado. Além de se contar com a quimioterapia clássica, hormonoterapia, radioterapia e a cirurgia oncológica^{3), (14}; existem também as terapias biológicas dirigidas (ou terapia target)¹⁴, e, recentemente, a imuno-oncologia⁸. Esses avanços, na maior parte dos casos, devem-se à melhor compreensão do câncer no ponto de vista genético e molecular. Uma vez estabelecida a doença e sabendo que as interações da região afetada podem ser estudadas em nível molecular, celular ou tecidual^{3), (6} é possível propor uma abordagem matemática com uso de diversos modelos dependendo do tipo de interação e do nível que quer ser estudado⁹ sendo que os focos de estudo estipulam abordagens que, entre muitos outros, vão desde uma interação de competição⁷ até uma modelagem hospedeiro-parasita como variante de uma interação presa-predador com encontros aleatórios fazendo uso de duas ou mais equações diferenciais ordinárias^{10), (11), (12), (15}.

Nas últimas duas décadas surgiram outros tipos de tratamento levando em consideração a resposta imunológica do corpo frente a qualquer fator invasivo, pois háevidências que o organismo pode reconhecer e eliminar os tumores malignos⁵. Desde 1980, após o primeiro trabalho de Stepanova¹ têm sido publicados muitos outros trabalhos nesta direção. As pesquisas clínicas desse tipo fazem uso da imunoterapia (como tratamento em imuno-oncologia), como outra frente de combate ao câncer, que age estimulando o sistema imune via injeção direta de citocinas¹³, ou com uso de vacinas como a do Papilomavírus humano e a da hepatite como medida preventiva ao câncer de colo de útero e de fígado, respectivamente⁴. Nesta frente de pesquisa, ultimamente estão sendo aplicados ambos os tratamentos: imunoterapia e quimioterapia. Um exemplo de modelagem matemática nessa direção são os trabalhos de López et al.¹⁰ e as referências nele inclusas.

Seja qual for a modelagem matemática, e para um suficiente entendimento da imunoterapia a ser usada, é importante levar sempre em consideração as várias componentes do microambiente em que o tumor se desenvolve: as células imunes, fibroblastos, moléculas sinalizadoras e os fatores de crescimento, entre outros¹; fazer as respectivas interpretações e possíveis mecanismos de previsão do modelo a usar^{6), (13}. O intuito neste trabalho é estudar matematicamente o tratamento imuno-oncológico que leve em consideração a questão do microambiente tumoral, sujeito a um comportamento switching² como aplicado em; trabalhos anteriores^{11), (12), (15} assim como fazer o estudo qualitativo e análise de estabilidade dos sistemas respectivos, não discutidos nos trabalhos^{1), (6), (10}.

Após essa introdução, na seção 2 é proposto dois modelos do problema (com e sem tratamento) usando sistemas de equações diferenciais ordinárias e na seção 3 o estudo qualitativo dos sistemas: existência e unicidade, existência dos pontos de equilíbrio, assim como a análise de es tabilidade. Na seção 4, mostram-se os resultados e comentários das simulações numéricas dos sistemas e na última seção apresentam-se as conclusões deste trabalho.

2 O MODELO PROPOSTO

Tratamentos imuno-terapêuticos que consideram as interações das células tumorais com o microambiente tumoral tem o objetivo de eliminar ou diminuir o crescimento tumoral, como bem destaca o artigo de Arciero et al. (2004) e suas referências¹. Na modelagem, esse trabalho utiliza o rol que as citocinas anti-inflamatórias apresentam no crescimento tumoral aliado aos fatores imuno-supressivos (citocinas do tipo interleucinas, a prostaglandina e os fatores de crescimento β). Nessa linha de pesquisa, os autores fazem uma modelagem com duas situações: uma sem tratamento e outra com tratamento. Usando esse princípio e conforme a explicação mais adiante, neste trabalho se propõe inserir um comportamento switching^{11), (12), (15}. Então, o modelo aqui proposto, denominado modelo sem tratamento, é o seguinte:

$$\begin{gathered} \dot{T}=rT\left(1-\frac{T}{K}\right)-\frac{aET}{g_2+T}+\frac{p_{2}ST}{g_3+S} \\ \dot{E}=-{\mu }_{1}E+\frac{cT}{1+\gamma S}+\left(\frac{EI}{g_1+I}\right)\left(p_1-\frac{q_{1}S}{q_2+S}\right) \\ \dot{I}=\frac{p_{3}ET}{(g_4+T)(1+\alpha S)}-{\mu }_{2}I \\ \dot{S}=h_1(T,n)-{\mu }_{3}S \end{gathered}$$ (2.1)

com condições iniciais

$T\left(0\right)=T_0,E\left(0\right)=E_0,I\left(0\right)=I_0,S\left(0\right)=S_0.$ (2.2)

O modelo (2.1) descreve a taxa de variação no tempo das variáveis que representam, respectivamente, a densidade populacional das células tumorais T, células efetoras E, citocinas I (do tipo IL-2) e as protéınas S (fator de crescimento TGF-β). Todos os parâmetros colocados no modelo são considerados como constantes positivas descritas logo a seguir. O parâmetro r representa a taxa de crescimento intrínseco e o parâmetro K a capacidade de suporte do crescimento das células tumorais; a é a força da resposta imune; p ₂ é a máxima taxa de crescimento proliferativo de TGF-β; μ₁, μ₂, μ₃ são, respectivamente, as taxas de decaimento das células efetoras E, citocinas I e protéınas S. O fator antigênico c, representa a habilidade que o sistema imunológico tem de reconhecer as células tumorais; γ é o parâmetro inibidor que reduz a expressão antigênica; g ₁, g ₂, g ₃, são as taxas de saturação média das variáveis I, T, S, respectivamente. O parâmetro p ₁ representa a máxima taxa da atividade entre as populações E e I na ausência das protéınas; q ₁ é o parâmetro que representa a máxima taxa do efeito anti-proliferativo e q ₂ é a taxa de saturação média de S; p ₃ é a máxima taxa de citocinas produzidas pela interação entre as populações E e T. O parâmetro g ₄ mede a saturação média das células tumorais na ausência de S e α é a taxa de inhibição do crescimento de I na presença de S.

Conforme usado em trabalhos anteriores^{11), (12), (15}, o termo h ₁(T, n) dado por

$h_1(T,n)=\frac{p_4{T}^n}{{\tau }_c^n+T^n},n\ge 2$ (2.3)

descreve o comportamento switching, no modelo aqui proposto, e mede o efeito estabilizador da população de células tumorais³. Quando falamos em switching neste trabalho, dá-se uma conotação similar (e diferente em forma) à fornecida em^{11), (12), (15}. Esses trabalhos, com aplicação biológica-ecológica, referem-se a populações num nível macroscópico e aqui, no microambiente tumoral, é considerado num nível microscópico. Arciero et al.¹ e Lopez et al.¹⁰ usam essas funções com conotação bioqúımica, cujos gráficos resultam ser sigmoides como mostrado em^{11 4}.

Na última equação do modelo (2.1), o termo h ₁ representa a influência que as células tumorais têm no crescimento das protéınas S, dessa forma τ_c é o parâmetro que indica o valor crítico de células tumorais e p ₄ é o parâmetro que representa a máxima taxa de produção de TGF-β pela influência das células tumorais. Por último, n é a intensidade do switching.

2.1 Modelagem com tratamento imuno-terapêutico

Como segunda variante, aqui se considera a inserção de um tratamento no modelo. O tratamento fornece um modelo pré-clínico do tipo siRNA⁵ (small interfering RNA) em que se limita ou suprime a produção da expressão TGF-β nas células tumorais. Essa metodologia terapêutica usa as citocinas do tipo IL-2 como estratégia de neutralização para inibir a produção e os efeitos imuno-supressivos do TGF-β. Nos últimos anos, estão sendo publicados trabalhos que relacionan o siRNA com os ligantes 1 e 2 do receptor PD-1 (ou Programmed Cell Death-1), outro pathway muito relevante em terapias imuno-oncolôgicas clínicas⁸.

Ao incorporar o tratamento imunoterapêutico siRNA no modelo, faz-se necessária a inserção de uma nova variável que a denominamos como A que representa a densidade populacional de ligantes vinculados ao siRNA. Dessa forma, insere-se uma equação adicional ao modelo (2.1), dando como origem um modelo de cinco equações em que se troca o termo h ₁(T, n) por

$h_2(T,n)=\frac{p_4{T}^n}{{\tau }_c^n+(1+\frac{fA}{k_1})T^n},n\ge 2.$

Assim, o modelo é representado por

$$\begin{gathered} \dot{T}=rT\left(1-\frac{T}{K}\right)-\frac{aET}{g_2+T}+\frac{p_{2}ST}{g_3+S} \\ \dot{E}=-{\mu }_{1}E+\frac{cT}{1+\gamma S}+\left(\frac{EI}{g_1+I}\right)\left(p_1-\frac{q_{1}S}{q_2+S}\right) \\ \dot{I}=\frac{p_{3}ET}{(g_4+T)(1+\alpha S)}-{\mu }_{2}I \\ \dot{S}=\frac{p_4{T}^n}{{\tau }_c^n+(1+\frac{fA}{k_1})T^n}-{\mu }_{3}S,n\ge 2 \\ \dot{A}=D_1\left(t\right)-{\mu }_{A}A \end{gathered}$$ (2.4)

com condições iniciais

$T\left(0\right)=T_0,E\left(0\right)=E_0,I\left(0\right)=I_0,S\left(0\right)=S_0,A\left(0\right)=A_0.$ (2.5)

O que diferencia o nosso modelo dos modelos apresentados em¹ e¹⁰ é a modificação do primeiro termo na quarta equação por eles propostos. Neste caso, h₂(T, n) deixa de ser um comportamento switching (como explicado no modelo (2.1)) e os valores de n, neste caso, melhoran os resultados quando comparados com o modelo de Arciero et al.

O modelo (2.4) apresenta três novos parâmetros positivos: k ₁, f e μ_A . O siRNA age como um inibidor não-competitivo do supressor com taxa constante de inibição k ₁ e f é a proporção de ligantes A, dessa forma f A representa a quantidade de fios siRNA ligados. Na última equação, D ₁ descreve a dose de infusão contínua que neste trabalho a consideramos como sendo constante. Por último, o parâmetro μ_A representa a taxa de decaimento da população A.

3 ESTUDO QUALITATIVO E ANÁLISE DE ESTABILIDADE

Com o intuito de estudar a existência e unicidade dos sistemas de equações diferenciais ordinárias (2.1) e (2.4) com condições iniciais dadas respectivamente por (2.2) e (2.5), representamos os sistemas nas formas dadas a seguir.

$$\begin{gathered} \left\{\begin{array}{l}\dot{X}=F\left(X\right)\\ X\left(0\right)=X_0\end{array}\right., \\ F\left(X\right)=\left[\begin{array}{c}{f}_1\left(X\right)\\ f_2\left(X\right)\\ f_3\left(X\right)\\ f_4\left(X\right)\end{array}\right],X=X\left(t\right)=\left[\begin{array}{c}T\left(t\right)\\ E\left(t\right)\\ I\left(t\right)\\ S\left(t\right)\end{array}\right],t\ge 0 \end{gathered}$$ (3.1)

e para k = 1, ..., 4, as funções f _k: R ⁴ → R são dadas por:

$$\begin{gathered} f_1\left(X\right)=rT\left(1-\frac{T}{K}\right)-\frac{aET}{g_2+T}+\frac{p_{2}ST}{g_3+S} \\ {f}_2\left(X\right)=-{\mu }_{1}E+\frac{cT}{1+\gamma S}+\left(\frac{EI}{g_1+I}\right)\left(p_1-\frac{q_{1}S}{q_2+S}\right) \\ {f}_3\left(X\right)=\frac{p_{3}ET}{(g_4+T)(1+\alpha S)}-{\mu }_{2}I \\ {f}_4\left(X\right)=\frac{p_4{T}^n}{{\tau }_c^n+T^n}-{\mu }_{3}S,n\ge 2. \end{gathered}$$

No caso do sistema (2.4)-(2.5) que modela o problema com tratamento imuno-terapêutico:

$$\begin{gathered} \left\{\begin{array}{l}\dot{X}=G\left(X\right)\\ X\left(0\right)=X_0\end{array}\right. \\ G\left(X\right)=\left[\begin{array}{c}{f}_1\left(X\right)\\ f_2\left(X\right)\\ f_3\left(X\right)\\ g_4\left(X\right)\\ g_5\left(X\right)\end{array}\right],X=X\left(t\right)=\left[\begin{array}{c}T\left(t\right)\\ E\left(t\right)\\ I\left(t\right)\\ S\left(t\right)\\ A\left(t\right)\end{array}\right],t\ge 0, \end{gathered}$$ (3.2)

em que f _k: R ⁵ → R, k = 1, 2, 3 tem a mesma lei de formação das funções do problema sem tratamento dado em (3.1), e as funções g ₄, g ₅: R ⁵ → R são:

$$\begin{gathered} \begin{array}{l}{g}_4\left(X\right)=\frac{p_4{T}^n}{{\tau }_c^n+(1+{\displaystyle \frac{fA}{k_1})T^n}}-{\mu }_{3}S,n\ge 2\end{array} \\ {g}_5\left(X\right)=D_1-{\mu }_{A}A. \end{gathered}$$

Consideramos também os espaços normados R ⁿ e C¹ ([0, ∞), R ⁿ ). Neles usamos, respectivamente, as seguintes normas: ||X|| = max _i {x _i } para X ∈ R ⁿ e ||H|| = sup_{t ≥ 0} ||h _i || para H ∈ C ¹ ([0, ∞), R ⁿ ) com h _i: R ⁿ → R.

Seja o cone positivo de R ⁿ dado por: ${\mathrm{ℝ}}_{+}^n$ = {X = [x ₁, x ₂, ..., x _n ]^T ∈ R ⁿ /x ⁱ > 0} usaremos V ⊂ ${\mathrm{ℝ}}_{+}^n$ da forma V = Π_{i = 1, ..., n} [a _i , b _i ] com [a _i , b _i ] intervalos fechados de R, e, pelo caráter biológico dos problemas modelados, as condições iniciais dos sistemas serão consideradas todas como sendo positivas, isto é, X ₀ ∈ ${\mathrm{ℝ}}_{+}^n$, n = 4, 5.

3.1 Existência e Unicidade

Antes de propor o Teorema de Existência e Unicidade de cada sistema precisaremos demonstrar dois lemas que serão dados nesta seção.

Lema 3.1. Para cada k, ||∇f _k || do sistema (3.1) é limitado.

Demonstração. Por ser cada função f _k diferenciável aplicamos o teorema do valor médio^6,7:

$f_k\left(X_2\right)-f_k\left(X_1\right)=\nabla f_k\left(\overline{X}\right)·(X_2-X_1)\text{para algum}\overline{X}\in \mathcal{L}(X_1,X_2)\subset U$

$\left|f_k\right(X_2)-f_k(X_1\left)\right|\le \left|\right|\nabla f_k\left(\overline{X}\right)\left|\right|\left|\right|X_2-X_1\left|\right|.$

sendo que

$\nabla f_k\left(X\right)=\left(\frac{\partial f_k}{\partial T}\left(X\right),\frac{\partial f_k}{\partial E}\left(X\right),\frac{\partial f_k}{\partial I}\left(X\right),\frac{\partial f_k}{\partial S}\left(X\right)\right)\text{para}k=1,\dots ,4.$

Então,

$\begin{array}{l}\nabla f_1\left(X\right)=\left(r-\frac{2r}{k}T-\frac{a{g}_{2}E}{(g_2+T{)}^2}+\frac{p_{2}S}{g_3+S},-\frac{aT}{g_2+T},0,\frac{p_2{g}_{3}T}{(g_3+S{)}^2}\right)\\ \nabla f_2\left(X\right)=\begin{array}{l}\left(\frac{c}{1+\gamma S},-{\mu }_1+\frac{p_{1}I}{g_1+I}-\frac{q_{1}S}{(g_2+S)(g_1+I)},\frac{p_1{g}_{1}E}{(g_1+I{)}^2}-\frac{q_1{g}_{1}ES}{(q_2+S)(g_1+I{)}^2},\right.\left.-\frac{c\gamma T}{(1+\gamma S{)}^2}-\frac{q_1{q}_{2}EI}{(g_1+I)(q_2+S{)}^2}\right)\end{array}\\ \nabla f_3\left(X\right)=\left(\frac{g_4{p}_{3}E}{(g_4+T{)}^2(1+\alpha S)},\frac{p_{3}T}{(g_4+T)(1+\alpha S)},-{\mu }_2,-\frac{\alpha p_{3}ET}{(g_4+T)(1+\alpha S{)}^2}\right)\\ \nabla f_4\left(X\right)=\left(\frac{n{p}_4{\tau }^n}{{\tau }^{2n}{T}^{1-n}+2{\tau }^{n}T+T^{n+1}},0,0,-{\mu }_3\right).\end{array}$

Considerando T ≥ T ₁ > 0, E ≥ E ₁ > 0, I ≥ I ₁ > 0, S ≥ S ₁ > 0, e T ≤ T _max, E ≤ E _max, I ≤ I _max, S ≤ S _max, e lembrando que todos os parámetros dos modelos são positivos, temos, se k = 1:

$$\begin{gathered} \left|\frac{\partial f_1}{\partial T}\right|<r+\frac{2r}{k}{T}_{\max }+\frac{a{g}_2{E}_{\max }}{(g_2+T_1{)}^2}+\frac{p_{2}S}{{\displaystyle \frac{g_3}{S_{\max }}+1}}=valo{r}_1, \\ \left|\frac{\partial f_1}{\partial E}\right|<\frac{a}{{\displaystyle \frac{g_2}{T_{\max }}+1}}=valo{r}_2, \\ \left|\frac{\partial f_1}{\partial I}\right|=0, \\ \left|\frac{\partial f_1}{\partial S}\right|<\frac{p_2{g}_3{T}_{\max }}{(g_3+S_1{)}^2}=valo{r}_3. \end{gathered}$$

Para cada j o valor _j > 0, então encontramos L ₁ > 0 tal que

$\left|\right|\nabla f_k\left(X\right)\left|\right|\le L_1\text{com}{L}_1=\max \{valo{r}_j,j=1,2,3,4\}$

Seguindo um procedimento similar para k = 2, 3, 4, determinamos valores positivos que limitam cada derivada parcial. Dai que:

$\left|\right|\nabla f_k\left(X\right)\left|\right|\le L_k\text{com}{L}_k=\max \{valo{r}_{kj},j=1,2,3,4\}\text{e}valo{r}_{kj}>0.$

Dessa forma terminamos a demonstração do lema. □

Teorema 3.1 (Existência e unicidade do sistema sem tratamento). Se n é par, existe uma única solução do problema de Cauchy dado por (2.1) - (2.3) no intervalo [0, ∞) com E(t) > 0, T(t) > 0, I (t) > 0, S(t) > 0 para todo t ≥ 0.

Demonstração. Usamos a representação (3.1) do sistema (2.1)-(2.3). Seja t₀ ≥ 0 e dados (t₀, X ₁), (t ₀, X ₂) ∈ [0, ∞) × U ⊂ [0, ∞) × R ⁿ ,⁸

$\left|\right|F(t_0,X_2)-F(t_0,X_1)\left|\right|=\left|\right|F\left(X_2\right)-F\left(X_1\right)\left|\right|.$

F é diferenciável pois cada função f _k de (3.1) é diferenciável. Consideremos V = [T ₁, T _max]× [E ₁, E _max]× [I ₁, I _max]× [S ₁, S _max] ⊂ R ⁴ e U uma vizinhança de V. Sem perda de generalidade consideramos a unitário⁹, então pelo teorema do valor médio existe $\overline{X}$ ∈ tal que:

$\Vert F\left(X_2\right)-F\left(X_1\right)\Vert \le \Vert DF\left(\overline{X}\right)\Vert \Vert (X_2-X_1)\Vert$

Com

$DF\left(\overline{X}\right)=\left[\begin{array}{c}\nabla f_1\left(\overline{X}\right)\\ \nabla f_2\left(\overline{X}\right)\\ \nabla f_3\left(\overline{X}\right)\\ \nabla f_4\left(\overline{X}\right)\end{array}\right].$ (3.3)

Aplicando o resultado do Lema (3.1) a cada componente de (3.3), obtemos

$\left|\right|DF\left(\overline{X}\right)\left|\right|\le {\Lambda }_1\text{com}{\Lambda }_1=\max \{L_k,k=1,\dots ,4\},\overline{X}\in U.$

Assim,

$\Vert F\left(X_2\right)-F\left(X_1\right)\Vert \le {\Lambda }_1\Vert X_2-X_1\Vert$

concluimos que F é uma função Lipchitz contínua local e uniformemente em X, a constante de Lipchitz é Λ₁ > 0, i.e., o sistema (3.1) tem uma única solução.

Vejamos agora que X(t) ∈ ${\mathrm{ℝ}}_{+}^4$, ou seja, E(t) > 0, T (t) > 0, I(t) > 0, S(t) > 0 para todo t ≥ 0.

Na última equação do sistema dado em (2.1): $\dot{S}=\frac{p_4{T}^n}{{\tau }_c^n+T^n}$ - μ₃ S, como n ≥ 2. Como n é par todos os parâmetros são positivos, temos: $\dot{S}>-{\mu }_{3}S\Rightarrow S\left(t\right)>S_0{e}^{-{\mu }_{3}t}$ ⇒ S(t) >, t ≥ 0. Da primeira equação do mesmo sistema:

$$\begin{gathered} \dot{T}=T\left(r\left(1-\frac{T}{K}\right)-\frac{aE}{g_2+T}+\frac{p_{2}S}{g_3+S}\right),\text{ou seja}, \\ T\left(t\right)=T_0{e}^{{\int }_0^t\left(r\left(1-\frac{T}{K}\right)-\frac{aE}{g_2+T}+\frac{p_{2}S}{g_3+S}\right)ds}\Rightarrow T\left(t\right)>0,t\ge 0. \end{gathered}$$

Na segunda equação:

$$\begin{gathered} \dot{E}=E\left(-{\mu }_1+\left(\frac{I}{g_1+I}\right)\left(p_1-\frac{q_{1}S}{q_2+S}\right)\right)+\frac{cT}{1+\gamma S},\text{então} \\ E\left(t\right)>E_0{e}^{{\int }_0^t\left(-{\mu }_1+\left(\frac{I}{g_1+I}\right)\left(p_1-\frac{q_{1}S}{q_2+S}\right)\right)ds}\Rightarrow E\left(t\right)>0,t\ge 0. \end{gathered}$$

Por último, na equação: $\dot{I}=\frac{p_{3}ET}{(g_4+T)(1+\alpha S)}$ -μ₂ I, e como T, E, S são positivas para todo t ≥ 0 temos que $\dot{I}>-{\mu }_{2}I\Rightarrow I\left(t\right)>S_0{e}^{-{\mu }_{2}t}$ ⇒ I(t) > 0, t ≥ 0.

Concluindo assim a demonstração do teorema.

O seguinte Lema será útil na demonstração da existência e unicidade do sistema que modela o problema com tratamento.

Lema 3.2. Dado o sistema (3.2) , cada gradiente de ||g ₄ ||, ||g ₅ || e ||∇f _k || é limitado para k = 1, 2, 3.

Demonstração. Considerando a diferenciabilidade das funções f _k e por elas serem as mesmas do Lema (3.1), utilizamos as mesmas limitantes L _k , então

$\left|\right|\nabla f_k\left(X\right)\left|\right|\le L_k\text{com}{L}_k=\max \{valo{r}_{kj},j=1,2,3,\}\text{e}k=1,2,3.$

As funções g ₄, g ₅ também são diferenciáveis, então aplicando o teorema do valor médio: |g _j (X ₂) - g _j (X ₁)| ≤ ||∇g _j ($\overline{X}$)||||X ₂ - X ₁||j = 4, 5, $\overline{X}$ ∈ L(X ₁, X ₂) ⊂ U.

Neste caso, $\nabla g_j\left(X\right)=\left(\frac{\partial g_j}{\partial T}\left(X\right),\frac{\partial g_j}{\partial E}\left(X\right),\frac{\partial g_j}{\partial I}\left(X\right),\frac{\partial g_j}{\partial S}\left(X\right),\frac{\partial g_j}{\partial A}\left(X\right)\right)$ para j = 4, 5. Então,

$$\begin{gathered} \nabla g_4\left(X\right)=\left(\frac{n{p}_4{\tau }^n}{{\tau }^{2n}{T}^{1-n}+2(1+\frac{f}{k_1}A){\tau }^{n}T+(1+\frac{f}{k_1}A{)}^2{T}^{n+1}},,,0,,,0,,,-,{\mu }_3,,,,\frac{\frac{f{p}_4}{k_1}}{{\left({\left(\frac{\tau }{T}\right)}^n+1+\frac{f}{k_1}A\right)}^2}\right), \\ \nabla g_5\left(X\right)=\left(0,0,0,0,-{\mu }_A\right). \end{gathered}$$

Considerando 0 < T ₁ ≤ T ≤ T _max, 0 < E ₁ ≤ E ≤ E _max, 0 < I ₁ ≤ I ≤ I _max, 0 < S ₁ ≤ S ≤ S _max, 0 < A ₁ ≤ A ≤ Amax e com um procedimento similar ao demonstrado no Lema (3.1) para encontrar os valores limitantes de cada componente do gradiente, encontramos que:

$$\begin{gathered} \begin{array}{l}\left|\frac{\partial g_4}{\partial T}\right|\le gvalo{r}_1,gvalo{r}_1=\frac{n{p}_4{\tau }^n}{{\tau }^{2n}{T_1}^{1-n}+2\tau (1+\frac{f}{k_1}{A}_1)T_1+(1+\frac{f}{k_1}{A}_1{)}^2(T_1{)}^{n+1}}>0\end{array} \\ \left|\frac{\partial g_4}{\partial A}\right|\le gvalo{r}_2,gvalo{r}_2=\frac{\frac{f{p}_4}{k_1}}{{\left({\left(\frac{\tau }{T_{max}}\right)}^n+1+\frac{f}{k_1}{A}_1\right)}^2}>0 \end{gathered}$$

Escolhendo M ₁ = max{gvalor ₁, gvalor ₂, μ₃} concluímos que ||∇g ₄(X)|| ≤ M ₁. Também, é imediato mostrar que ||∇g ₅(X)|| = μA. Concluimos a prova.

Na demonstração do seguinte teorema, alguns ítens de prova do Teorema (3.1) serão re-aproveitados.

Teorema 3.2 (Existência e unicidade do sistema com tratamento). Se n é par, existe uma única solução do problema de Cauchy dado pelo sistema (2.4) - (2.5) . No intervalo [0, ∞) com E(t) > 0, T(t) > 0, I(t) > 0, S(t) > 0, A(t) > 0 para todo t ≥ 0.

Demonstração. Usaremos a representação (3.2) do sistema (2.4)-(2.5). G é diferenciável pois cada função f _k , g ₄, g ₅ de (3.2) é diferenciável. Então, considerando U vizinhança de V, existe $\overline{X}$ ∈ L(X ₁, X ₂) ⊂ U ⊂ V tal que:

$\Vert G\left(X_2\right)-G\left(X_1\right)\Vert \le \Vert DG\left(\overline{X}\right)\Vert \Vert (X_2-X_1)\Vert$

com

$DG\left(\overline{X}\right)=\left[\begin{array}{c}\nabla f_1\left(\overline{X}\right)\\ \nabla f_2\left(\overline{X}\right)\\ \nabla f_3\left(\overline{X}\right)\\ \nabla g_4\left(\overline{X}\right),\\ \nabla g_5\left(\overline{X}\right)\end{array}\right]$ (3.4)

Usando V = [T ₁, T _max] × [E ₁, E _max] × [I ₁, I _max] × [S ₁, S _max] × [A ₁, A _max]. Aplicando o resultado do Lema (3.2) a cada componente de (3.4), obtemos ||DG(X)|| ≤ Λ₂ com Λ₂ = max{L _k , μ₃, M ₁} k = 1, ..., 4, $\overline{X}$ ∈ U. Assim,

$\Vert G\left(X_2\right)-G\left(X_1\right)\Vert \le {\Lambda }_2\Vert X_2-X_1\Vert$

concluimos que G é uma função Lipchitz contínua local e uniformemente em X com constante de Lipchitz Λ₁ > 0, i.e., o sistema (3.2) e em consequência o sistema (2.4)-(2.5) tem uma única solução.

Vejamos agora E(t) > 0, T(t) > 0, I(t) > 0, S(t) > 0, A(t) > 0 para todo t ≥ 0.

Da última equação do sistema dado em (3.2), $\dot{A}$ = D ₁ - μ_A A, como D ₁ é um parâmetro positivo, segue que:

$\dot{A}=-{\mu }_{A}A\Rightarrow A\left(t\right)>A_0{e}^{-{\mu }_{A}t}\Rightarrow A\left(t\right)>0,t\ge 0.$

Na penúltima equação do mesmo sistema: $\dot{S}=\frac{p_4{T}^n}{{\tau }_c^n+(1+\frac{fA}{k_1})T^n}-{\mu }_{3}S$, e como todos os parâmetros são positivos, n par e A > 0 temos:

$\dot{S}>-{\mu }_{3}S\Rightarrow S\left(t\right)>S_0{e}^{-{\mu }_{3}t}\Rightarrow S\left(t\right)>0,t\ge 0.$

Sendo que as três primeiras equações dos sistemas (3.1) e (3.2) são iguais, então para t ≥ 0 a positividade das soluções componentes T(t), E(t), I(t) estão garantidas, como visto no Teorema (3.1). Concluimos a demonstração do teorema.

Teorema 3.3. Se p ₁ < (μ ₁ + q ₁ ) e n par, as soluções dos sistemas (2.1) e (2.4) com respectivas condições iniciais (2.2) e (2.5) , sã o limitadas para todo t > 0.

Demonstração. Como as três primeiras equações dos sistemas (2.1) e (2.4) são iguais, a prova das soluções T, S, I serem limitadas será dada apenas uma vez. Também, em todo o processo usaremos o fato, já provado, de E(t), T(t), I(t), S(t), A(t) serem positivas para todo t ≥ 0.

Para a solução T, usamos a equação

$\dot{T}=rT\left(1-\frac{T}{K}\right)-\frac{aET}{g_2+T}+\frac{p_{2}ST}{g_3+S}.$

Então, da positividade das soluções e parâmetros temos $\dot{T}\le rT\left(1-\frac{T}{K}\right)+p_{2}T$ e como

$\dot{u}=ru\left(1-\frac{u}{K}\right)+p_{2}u,u\left(0\right)=u_0>0$

Tem solução

$u\left(t\right)=\frac{(r+p_2)u_0}{\frac{u_0}{K}+(1-\frac{u_0}{K})e^{-t}}$

limitada por $K(1+\frac{p_2}{r})$ para todo t ≥ 0, isto é, u(t) ≤ $K(1+\frac{p_2}{r})$, então por², segue que T(t) é limitada para todo t ≥ 0. No próximo passo, consideramos M _T = $K(1+\frac{p_2}{r})$ > 0 como sendo o valor que limita a solução T.

Para a solução E, usamos a equação

$\dot{E}=-{\mu }_{1}E+\frac{cT}{1+\gamma S}+\left(\frac{EI}{g_1+I}\right)\left(p_1-\frac{q_{1}S}{q_2+S}\right).$

De forma similar ao processo anterior temos que

$\dot{E}\le cT+\left[p_1-({\mu }_1+q_1)\right]E\le c{M}_T+\left[(p_1-({\mu }_1+q_1)\right]E,$

isto é, $\dot{E}$ ≤ c _M T + [(p ₁ - (μ₁ + q ₁)]E.

Como a solução de $\dot{u}$c _M T + [(p ₁ - (μ₁ + q ₁)]u, u(0) = u ₀ > 0 é dada por

$u\left(t\right)=\frac{c{M}_T}{-p_1+({\mu }_1+q_1)}+(u_0-\frac{c{M}_T}{-p_1+({\mu }_1+q_1)})e^{-\left(\right({\mu }_1+q_1)-p_1)t}$

está limitada por N = $\frac{c{M}_T}{-p_1+({\mu }_1+q_1)}$ > 0, para todo t ≥ 0 pois por hipótese (μ₁ + q ₁) − p ₁ > 0, então segue que E é uma solução limitada por N = $\frac{c{M}_T}{-p_1+({\mu }_1+q_1)}$ > 0, para todo t ≥ 0.

Antes de mostrar que a solução I(t) é limitada provemos que a solução S(t) do sistema (2.1) é limitada. Esta solução vem da equação

$\dot{S}=\frac{p_4{T}^n}{{\tau }_c^n+T^n}-{\mu }_{3}S,n\ge 2.$

Como T ≥ T ₀ > 0 e T é limitada por M _T = $K(1+\frac{p_2}{r})$ > 0, é imediato verificar que $\frac{p_4{T}^n}{{\tau }_c^n+T^n}$ está limitada por $\frac{p_4{M_T}^n}{{\tau }_c^n+{M_T}^n}$. Então, $\dot{S}$ ≤ L - μ₃ S com L = $\frac{p_4{M_T}^n}{{\tau }_c^n+{M_T}^n}$ ≥ 2 (n par).

Como $\dot{u}$ = L - μ₃ u, u(0) = u ₀ > 0 tem como solução

$u\left(t\right)=\frac{L}{{\mu }_3}+(u_0-\frac{L}{{\mu }_3})e^{-{\mu }_{3}t}$

então u(t) é limitada e em consequência a solução S(t) é limitada por $\frac{L}{{\mu }_3}$ para todo t ≥ 0.

Seguindo um procedimento similar, prova-se facilmente que tanto as soluções I(t), S(t), A(t) (essas últimas soluções do sistema (2.4)), são limitadas para todo t ≥ 0. Prova concluída.

Os Teoremas provados acima fornecem às soluções dos sistemas respectivos a propriedade de invariança nos cones positivos de R ⁴ e R ⁵, mostrando dessa forma que o problema de conotação biológica e médica é matematicamente viável de ser tratado numericamente.

3.2 Análise de estabilidade

Nesta subseção estudaremos a existência dos pontos de equilíbrio assim como a a estabilidade dos sistemas envolvidos no problema. Para tal, e sem perda de generalidade, estudaremos os sistemas adimensionalizados respectivos.

Para adimensionalizar os modelos consideramos as seguintes mudanças de variáveis e parâmetros.

No caso do modelo (2.1):

$\begin{array}{l}w=E/g_2,x=T/g_2,y=I/g_1,z=S/g_3,\overline{t}={\mu }_{2}t,\overline{c}=c/{\mu }_2,\overline{{\mu }_i}={\mu }_i/{\mu }_2,(i=1,3),\\ \overline{q_1}=q_1/{\mu }_2,\overline{q_2}=q_2/g_3,\overline{\gamma }=\gamma g_3,\overline{p_1}=p_1/{\mu }_2,\overline{p_2}=p_2{\mu }_2,\overline{r}=r/{\mu }_2,\overline{k}=K/{\mu }_2,\\ \overline{a}=a/{\mu }_2,\overline{p_3}=\frac{p_3{g}_2}{{\mu }_2{g}_1},\overline{\alpha }=\alpha g_3,\overline{g_4}=g_4/g_2,\overline{p_4}=\frac{p_4}{{\mu }_2{g}_3},{\overline{\tau }}_c={\tau }_c/g_2.\end{array}$

Para o modelo (2.4), além das variáveis e parâmetros acima, consideramos:

$v=\frac{f}{k_1}A,{\overline{D}}_1=D_{1}f/{\mu }_2{k}_1,{\overline{\mu }}_A={\mu }_A/{\mu }_2.$

Ao desconsiderar as barras das variáveis e parâmetros, obtemos os seguintes modelos adimensionalizados.^10,11

Para o modelo sem tratamento:

$\dot{x}=rx\left(1-\frac{x}{k}\right)-\frac{axw}{1+x}+\frac{p_{2}xz}{1+z}$ (3.5)

$\dot{w}=-{\mu }_{1}w+\frac{cx}{1+\gamma z}+\left(\frac{yw}{1+w}\right)\left(p_1-\frac{q_{1}z}{q_2+z}\right)$ (3.6)

$\dot{y}=\frac{p_{3}wx}{(g_4+x)(1+\alpha z)}-y$ (3.7)

$$\begin{gathered} \dot{z}=\frac{p_4{x}^n}{{\tau }_c^n+x^n}-{\mu }_{3}z,n\ge 2 \\ \text{c.i.:}w\left(0\right)=w_0,x\left(0\right)=x_0,y\left(0\right)=y_0,z\left(0\right)=z_0. \end{gathered}$$ (3.8)

Jáo modelo com tratamento é dado por:

$\dot{x}=rx\left(1-\frac{x}{k}\right)-\frac{axw}{1+x}+\frac{p_{2}xz}{1+z}$ (3.9)

$\dot{w}=-{\mu }_{1}w+\frac{cx}{1+\gamma z}+\left(\frac{yw}{1+w}\right)\left(p_1-\frac{q_{1}z}{q_2+z}\right)$ (3.10)

$\dot{y}=\frac{p_{3}wx}{(g_4+x)(1+\alpha z)}-y$ (3.11)

$\dot{z}=\frac{p_4{x}^n}{{\tau }_c^n+(1+v)x^n}-{\mu }_{3}z,n\ge 2$ (3.12)

$$\begin{gathered} \dot{v}=D_1-{\mu }_{A}v \\ \text{c.i.:}w\left(0\right)=w_0,x\left(0\right)=x_0,y\left(0\right)=y_0,z\left(0\right)=z_0,v\left(0\right)=v_0. \end{gathered}$$ (3.13)

3.3 Existência dos pontos de equilíbrio

Ao igualar a zero cada equação do lado direito do sistema (3.5)-(3.8) obtemos o que segue.

Da equação (3.5), x = 0 ou r$\left(1-\frac{x}{k}\right)-\frac{aw}{1+x}+\frac{p_{2}z}{1+z}=0$. Se x = 0, das equações (3.8) e (3.7) obtemos z = 0, y = 0 e em (3.6): w = 0. Assim, obtemos o primer ponto de equilíbrio (o trivial): E ₀ = (0, 0, 0, 0)¹².

Se x ≠ 0 $\tilde{\alpha }=k\left(1+\frac{\frac{p_2}{r}z}{1+z}\right)$, da equação $\left(1-\frac{x}{k}\right)-\frac{aw}{1+x}+\frac{p_{2}z}{1+z}=0$ obtemos dois valores de x:

$x_{1,2}=\frac{1}{2}\left[(1-\tilde{\alpha })\pm \sqrt{(1+\tilde{\alpha }{)}^2-4\frac{k}{r}aw}\right],$ (3.14)

em que $2\sqrt{\frac{k}{r}a}\sqrt{w}\le 1+\tilde{\alpha }$ o que implica

$w\ge \frac{(1+k)(1+\frac{2p_2}{r}k)+1}{\frac{4a}{r}(1+\frac{p_2}{r}k{)}^2}>0$

para x ser um valor real. Para x ₁, x ₂ ser positivos (considerando x ₁, x ₂ serem, respectivamente, os valores de x para a parte positiva e negativa da raiz) teremos as seguintes condições entre z e w:

$x_1>0\text{se}z>-\frac{1-\frac{1}{k}}{1-\frac{1}{k}+\frac{p_2}{r}}\text{e}z>-\frac{1-\frac{aw}{r}}{1-\frac{aw}{r}+\frac{p_2}{r}}$ (3.15)

$x_2>0\text{se}-\frac{1-\frac{1}{k}}{1-\frac{1}{k}+\frac{p_2}{r}}<z<-\frac{1-\frac{aw}{r}}{1+\frac{p_2}{r}}$ (3.16)

as condições acima descritas determinam uma regiões de pontos no plano WZ. No caso de (3.15) determina o conjunto de pontos (w, z) que se encontam acima do valor

$z=-\frac{1-\frac{1}{k}}{1-\frac{1}{k}+\frac{p_2}{r}}$

e acima da hipérbole trasladada (dada em foma implícita): $-\frac{a}{r}zw+(1+\frac{p_2}{r})z-\frac{a}{r}w+1=0$ e a região de pontos dada por (3.16) são aqueles acima de $z=-\frac{1-\frac{1}{k}}{1-\frac{1}{k}+\frac{p_2}{r}}$ e abaixo da reta

$z=\frac{1-\frac{aw}{r}}{1+\frac{p_2}{r}}.$

Considerando as condições acima relatadas e os valores x _1,2, calculamos duas possibilidades para z: z _1,2 ao zerar a equação (3.8), isto é;

$z_{1,2}=\frac{p_4{x}_{1,2}^n}{{\mu }_3({\tau }_c^n+x_{1,2}^n)},n\ge 2.$ (3.17)

Ao zerar as equações do lado direito de (3.6) e (3.7) encontramos o seguinte sistema para w e y por resolver:

$$\begin{gathered} w\left[-{\mu }_1+y(-{\mu }_1+p_1-\frac{q_{1}z}{q_2+z})\right]=-\frac{cx(1+y)}{1+\gamma z} \\ y=\frac{p_{3}wx}{(g_4+x)(1+\alpha z)} \end{gathered}$$ (3.18)

inserindo este último valor de y na penúltima equação e considerando os valores de x _1,2 e z _1,2 dados por (3.14) e (3.17) obtemos a equação quadrática s ₁ w ₂ + s ₂ w + s ₃ = 0 que tem por solução:

$w_{1,2}=\frac{-s_2\pm \sqrt{s_2^2-4s_1\frac{cx}{1+\gamma z}}}{2s_1}$ (3.19)

com

$s_1=\left(-{\mu }_1+p_1-\frac{q_{1}z}{q_2+z}\right)\frac{p_{3}x}{(g_4+x)(1+\alpha z)},s_2=-{\mu }_1+\frac{c{p}_3{x}^2}{(1+\gamma z)(g_4+x)(1+\alpha z)}.$

Ao igual que as outras variáveis, w deve ser real e positivo. Para ser real,

$\Delta =s_2^2-4s_1\frac{cx}{1+\gamma z}>0,$

colocando os valores de s ₁, s ₂ nesta expressão e considerando a positividade de z e dos parâmetros obtemos que 0 < $z\le \frac{(-{\mu }_1+p_1)q_2}{{\mu }_1+q_1-p_1}$ e a seguinte inequação descrita pelos pontos na parábola ou acima da parábola h(t) = [t − (μ₁ + 2u)]² - 4u(μ₁ + u) ≥ 0¹³ onde

$t=\frac{c{p}_3{x}^2}{(1+\gamma z)(g_4+x)(1+\alpha z)}$

e u = -μ₁ + p ₁ - $\frac{q_{1}z}{q_2+z}$ > 0.¹⁴ A positividade de w: como -μ₁ + p ₁ - $\frac{q_{1}z}{q_2+z}$ > 0 então s ₁ > 0 e na equação (3.19) analisaremos a positividade do numerador que nos darão as duas ráızes w ₁, w ₂, descartando a primeira por ser negativa. Então, ficamos apenas com w ₂ > 0 se s ₂ + s ₄ < 0. Por um tratamento algébrico simples, isto equivale a ter de ficar com os pontos na hipérbole ou dentro dos ramos da hipérbole de pontos (z, t) da forma

$\left(\sqrt{\frac{3}{2}}t+{\mu }_1-p_1\right)(\frac{q_2}{z}+1\mathrm{\backslash Big})+q_1\le 0$

onde t foi dado acima.

Obtidos o valor de w voltamos na equação (3.18) para termos o valor de y¹⁵ e da forma como estão definidos jásão positivos.

O procedimento analítico com ferramentas algébricas e geométricas descrito acima faz com que encontremos dois pontos de equilíbrio (de coexistência) E ₁, E ₂ para o sistema (3.5).

Na procura dos pontos de equilíbrio do modelo com tratamento (3.9)-(3.13), no que segue aproveitaremos parte da metodologia que garantiu a existência dos equilíbrios de coexistência do modelo sem tratamento.

Zerando o lado direito da equação (3.9) obtemos que x = 0 ou x ≠ 0. Considerando o primeiro valor nulo para x, e, ao zerar os lados direitos de (3.11), (3.12) e (3.13) e substituindo-os em (3.10) obtemos de forma imediata o equilíbrio quase trivial, $E_0^t$ = (0, 0, 0, 0, D ₁/μ_A ).

Considerando x ≠ 0 segue, como feito antes, o mesmo resultado que em (3.14) (e suas respectivas condições). Sabendo que neste caso também vale v = D ₁/μ_A , encontramos os valores de z:

$z_{1,2}^t=\frac{p_4{x}_{1,2}^n}{{\mu }_3({\tau }_c^n+(1+\frac{D_1}{{\mu }_A}\left)x_{1,2}^n\right)},n\ge 2.$ (3.20)

que ao ser substituido em (3.18) (3.19) (com as respectivas e similares condições) obtemos os valores para $w_{1,2}^t$ e $y_{1,2}^t$ e descrevendo dessa forma os pontos de coexistência $E_1^t$, $E_2^t$ para o sistema tratado. Isto finaliza a prova do seguinte Lema.

Lema 3.3 (Existência dos pontos de equilíbrio dos sistemas). Os sistemas dados em (3.5) - (3.8) e (3.9) - (3.13) têm:

Um ponto de equilíbrio livre da doença: o trivial E ₀ = 0 e o quase trivial $E_0^t$ = (0, 0, 0, 0, D ₁ /μ _A ), respectivamente.

Dois pontos de coexistência E ₁ , E ₂ e $E_1^t$ , $E_2^t$ para cada sistema, respectivamente ¹⁶ .

No seguinte Teorema apresentaremos a estabilidade dos pontos de equilíbrio garantidos no Lema 3.3 quando tratados pelos sistemas linearizados correspondentes aos modelos com e sem tratamento adimensionalizados.

Teorema 3.4 (Estabilidade dos Sistemas). Sejam o sistema (3.5) - (3.8) correspondente ao sistema sem tratamento e o sistema com tratamento (3.9) - (3.13) . Nos pontos de equilíbrio livre da doença os sistemas sã o instáveis e nos pontos de coexistência a estabilidade é assintótica.

Demonstração. Estudaremos a estabilidade dos respectivos sistemas linearizados da forma

$\dot{X}=JX.$

onde J é a matriz jacobiana cuja ordem estáassociada ao número de variáveis do sistema respectivo. No caso, usaremos as matrizes J _G e J _F , dos sistemas com e sem tratamento (adimensionalizado) (3.9)-(3.13) e (3.5)-(3.8), respectivamente. Sendo que os sistemas são muito parecidos, serásuficiente linearizar o sistema com tratamento cujo jacobiano colocaremos a seguir.

$J_G=\left[\begin{array}{cc}{J}_F& 0\\ 0& -{\mu }_A\end{array}\right],$

J_F é a matriz de ordem 4 (dada abaixo) e os zeros são matrizes linha e coluna¹⁷,

$J_F=\left[\begin{array}{cccc}{a}_{11}& -\frac{ax}{1+x}& 0& \frac{p_{2}x}{(1+z{)}^2}\\ \frac{c}{1+\gamma z}& a_{22}& a_{23}& a_{24}\\ a_{31}& a_{32}& -1& a_{34}\\ a_{41}& 0& 0& -{\mu }_3\end{array}\right]$

$$\begin{gathered} \begin{array}{l}{a}_{11}=r-\frac{2r}{k}x-\frac{aw}{1+x^2}+\frac{p_{2}z}{1+z},a_{22}=-{\mu }_1+\frac{p_{1}y}{1+y}-\frac{q_{1}yz}{(1+y)(q_2+z)}, \\ {a}_{23}=\frac{w}{(1+y{)}^2}\left[p_1-\frac{q_{1}z}{q_2+z}\right],a_{24}=-\frac{c\gamma x}{(1+\gamma z{)}^2}-\frac{q_1{q}_{2}wy}{(1+y)(q_2+z{)}^2}, \\ {a}_{31}=\frac{p_3{g}_{4}w}{(g_4+x{)}^2(1+\alpha z)},a_{32}=\frac{p_{3}x}{(g_4+x)(1+\alpha z)}, \\ {a}_{34}=-\frac{\alpha p_{3}wx}{(g_4+x)(1+\alpha z{)}^2},a_{41}=\frac{n{p}_4{\tau }_c^n{x}^{n-1}}{({\tau }_c^n+x^n{)}^2}\left(\text{ou}{a}_{41}^G=\frac{n{p}_4{\tau }_c^n{x}^{n-1}}{({\tau }_c^n+(1+\frac{D_1}{{\mu }_A})x^n{)}^2}\right).\end{array} \end{gathered}$$

Avaliamos as matrizes J _F e J _G nos pontos de equilíbrio caracterizados pelo Lema (3.3) e analisamos a estabilidade como segue.

Ao ser avaliada no ponto E ₀, a matriz J _F é diagonal, então em E ₀ o sistema sem tratamento é instável¹⁸ por ter um autovalor positivo(os outros são negativos). De forma similar, ao avaliar J _G em $E_0^t$ a matriz diagonal resultante tem um autovalor positivo e os demais negativos, então em $E_0^t$ o sistema tratado é instável.

Pela forma como construimos J _G podemos, de uma forma simultânea, encontrar as equações caracteristicas nos outros pontos de equilíbrio¹⁹:

$P_G\left(\lambda \right)=|J_G-\lambda I|=-({\mu }_A+\lambda )P_F\left(\lambda \right)=0$

com P _F (λ) = λ₄ + A ₃λ³ + A ₂λ² + A ₁λ + A ₀ em que

$$\begin{gathered} \begin{array}{l}{A}_0=({\mu }_3{a}_{11}+\frac{a_{41}{p}_{2}x}{(1+z{)}^2})(a_{22}+a_{23}{a}_{32})+\left[{\mu }_3{a}_{23}{a}_{31}+(a_{24}+a_{23}{a}_{34})a_{41}+\frac{{\mu }_{3}c}{1+\gamma z}\right]\frac{ax}{1+x}, \\ {A}_1={\mu }_3(1+a_{23}+a_{32})-(a_{22}+a_{23}{a}_{32})a_{11}-\left[a_{23}{a}_{31}+a_{24}{a}_{41}+\frac{c(1+{\mu }_3)}{1+\gamma z}\right]\frac{ax}{1+x}, \\ {A}_2=-{\mu }_3(a_{11}+a_{22})+(a_{11}-1)a_{22}+{\mu }_3-a_{11}-a_{23}{a}_{32}+\frac{acx}{(1+x)(1+\gamma z)}-\frac{a_{41}{p}_{2}x}{(1+z{)}^2}, \\ {A}_3=1+{\mu }_3-(a_{11}+a_{22}).\end{array} \end{gathered}$$

Dessa forma, o sistema tratado é assintóticamente estável se o sistema sem tratamento o for. Utilizando o método de Hurwitz para P _F (λ), no(s) ponto(s) de equilíbrio(s), surgem as seguintes condições para estabilidade assintótica dos sistemas:

$$\begin{gathered} \begin{array}{l}{A}_3>0,A_0>0, \\ {A}_3{A}_2>A_1\text{e} \\ {A}_3{A}_2{A}_1>A_1^2+A_3^2{A}_0\end{array} \end{gathered}$$

O próximo exemplo ilustra a aplicação do Lema (3.3) e o Teorema 3.4 no caso de um sistema sem tratamento cujos parâmetros são dados por r = 0, 18, K = 10⁹, a = 1, p ₂ = 0, 27, μ₁ = 0, 03, μ₂ = 10, μ₃ = 10 , c ∈ [0, 0, 035], γ = 10, α = 10^-3, n ∈ [2, 200], g ₁ = 2 × 107, g ₂ = 105, g ₃ = 2 × 107, p ₁ = 0, 1245, q ₁ = 0, 1121, q ₂ = 2 × 106, p ₃ = 5, g ₄ = 103, p ₄ ∈ [0, 3 × 108], τ_c = 106.

Exemplo. Seja o modelo sem tratamento com os parâmetros dados acima e n > 2, c = 0, 035 e p ₄ = 121. Dos dados fornecidos verifica-se que p ₁ < μ₁ + q ₁, então as condições do Teorema 3.3 ficam satisfeitas.

Após considerar as condições estabelecidas pelo Lema 3.3 obtemos o ponto de coexistência²⁰E ₁ = (0.2764, 0.018, 4.33(10)^-5, 9.75(10)^-15),

$J_F=\left[\begin{array}{rrrr}0,0169& -0,0217& 0& 0,7463\\ 0,0035& -0,003& 0,000224& 193480\\ 0,00000548& 0.0024& -1& 0.8686\\ 1,7655(10{)}^{-13}& 0& 0& -1\end{array}\right]$

$P_F\left(\lambda \right)={\lambda }^4+1.9861{\lambda }^3+0,9722{\lambda }^2+1,0025\lambda +2,5106(10{)}^{-5}.$

Como $a_{41}^G$ = 6,8(10)^-18, o polinômio P _G (λ) correspondente acaba sendo muito similar numericamente à P _F (λ), verificando-se as condições de Hurwitz em ambos casos. Então, pelo Teorema 3.4, os sistemas são estáveis assintóticamente.

Observação. Embora, fazendo-se uma aproximação numérica de um sistema não linear resolveria o problema da procura do ponto de equilíbrio lidando em forma aproximada e com um grande número de parâmetros fixados, neste trabalho, decidiu-se fazer a procura de forma analítica, permitindo dessa forma a modificação de alguns dos parâmetros como serávisto na próxima seção.

4 RESULTADOS DAS SIMULAÇÕES NUMÉRICAS

Nesta seção exibimos os resultados numéricos da evolução das células tumorais dos modelos adimensionalizados (3.5)-(3.8) e (3.9)-(3.13).

Nas simulações aqui mostradas os parâmetros usados, em quase todos casos, foram aqueles dos trabalhos¹ e¹⁰. Além dos parâmetros jáfornecidos no exemplo da seção anterior, utilizam-se D ₁ = 5 × 10¹⁰, μ_A = 0, 66, f = 0, 9, k ₁ = 1. Todo o processamento dos dados foi realizado usando um software matemático que usa o método de Runge-Kutta de ordem 4 para resolução das equações diferenciais ordinárias. A discretização foi realizada usando um centésimo como passo no tempo.

A Figura 1a exibe três resultados simultâneos obtidos para um valor de c pequeno (c = 0, 5 × 10^-6), mostrando em todos os casos atingir a capacidade de suporte para valores p ₄ ∈ {0, 5 × 107, 3 × 108}. Esse comportamento acontece em valores próximos a 1000 unidades e não muda ao modificar n entre 2 e 65. Quando c deixa de ser pequeno (c = 0, 035) e p ₄ = 121, observa-se que a população de células tumorais vai para estabilidade assintótica. Esse resultado é mostrado na Figura 1b.

Figura 1:
Resultados do modelo (2.1) (sem tratamento).

Figura 2:
Comparação de resultados: à esquerda o modelo de Arciero e à direita o modelo (2.4) aqui proposto (ambos com tratamento).

Na Figura 2 se comparam os resultados do modelo com tratamento aqui proposto (modelo (2.4)), com aquele de¹ que utiliza o valor fixo n = 2. Quando o valor da resposta imune aumenta de a = 1 para a = 1, 3, percebe-se uma resposta oscilatória na densidade das células tumorais (adimensionais) por conta da dose de infusão contínua. Observa-se que o valor dos picos decresce (valores próximos de 1,5 para 1) em quanto o número dos mesmos aumenta (de 7 para 8) como se exibe na Figura 2. Ao mesmo tempo acontece um atraso nos picos na escala adimensional do tempo ao compararmos as Figuras 2a e 2b, o que fornece uma forma de controlar o crescimento do tumor ao longo do tempo (adimensional) e com um valor máximo de menor intensidade a cada vez.

4 CONCLUSÕES

Neste artigo apresentamos duas variantes dos modelos, com e sem tratamento, dados em¹ que consideram a evolução no tempo das células tumorais no seu microambiente. Microambiente determinado pela interação entre as células tumorais com as células efetoras, citocinas anti-inflamatórias e um fator imuno-supressivo. Essas variantes consideram o comportamento switching como trabalhado em Palomino et al.^{11), (12), (15}.

Foi realizado o estudo qualitativo dos sistemas mostrando condições para a existência e unicidade assim como para a existência dos pontos de equilíbrio, assuntos não explorados nas referências¹ e¹⁰. De forma direta encontramos em cada sistema um ponto de equilíbrio livre da doença, que por não ter interpretação biológica do problema carece de importância. Os resultados do modelo sem tratamento são similares aos obtidos em Arciero et al.¹ em que ao considerarse uma fraca habilidade do sistema imunológico em reconhecer as células tumorais T e uma baixa produção de TGF-β, se observa que a quantidade de células tumorais cresce rapidamente até atingir a capacidade de suporte das mesmas. Mas, ao aumentar esses valores, consegue-se modificar o nível de células T e se obter estabilidade assintótica com apenas aumentar o valor de n em algumas unidades (n > 2, enquanto o modelo em¹ faz atingir a capacidade de suporte (crescimento acelerado do nível das células tumorais). Quando inserindo um tratamento, como mostra o modelo (2.4), e ao comparar os resultados, o modelo aqui proposto fornece um comportamento oscilatório com diminuição do nível dessas células e retardando o crescimento das mesmas. Dessa forma o modelo aqui proposto mostra que simulando com valores de n até 100 a densidade de células tumorais decresce mais²¹ que ao utilizar o modelo de Arciero et al.¹.

Para outro conjunto de parâmetros, e mudando os parâmetros c, p ₄ e a, observa-se também a presença de ciclos limites que mostra um estudo futuro de bifurcações com relação a esses parâmetros. Também, como trabalho futuro, e dado que sóhátrabalhos de cunho clínico⁸, propõe-se estudar matematicamente a terapia siRNA que inibe os ligantes 1 e 2 do receptor PD-1 para morte programada das células tumorais T.

REFERÊNCIAS

• ¹
  J.C. Arciero, T. Jackson & D. Kirschner. A mathematical model of tumor-immune evasion and siRNA treatment. Discrete and Cont. Dynam. Syst., 4(1) (2004), 39-58.
• ²
  B.G. Bachpatte. Comparised theorems related to a certain inequality used in the theory of differential equations. S. Journal of Math., 22(3) (1996), 383-394.
• ³
  N. Bellomo & L. Preziosi. Modelling and mathematical problems related to tumor evolution and its interaction with the immune system. Math. Comp. Model, 32 (2000), 413-452.
• ⁴
  Brasil. “Estimativa 2016: Incidência do Câncer no Brasil”, INCA, Rio de Janeiro, (2015).
• ⁵
  H. Byrne, S. Cox & C. Kelly. Macrophage-tumor interactions: in vivo dynamics. Discr. Cont. Dyn. Syst. B, 4(1) (2004), 81-98.
• ⁶
  L. De Pillis, A. Radunskaya & C. Wiseman. A validated mathematical model of cell mediated immune response to tumor growth. Cancer Res., 65(17) (2005), 7950-7958.
• ⁷
  B.I. Freedman. “Modeling cancer treatment using competition: a survey”. Mathematics for Life Science and Medicine, (Y. Takeuchi, Y. Iwasa e K. Sato, eds), Springer-Verlag, Berlin Heidelberg, (2007).
• ⁸
  K. Iwamura, T. Kato, Y. Miyahara, H. Naota, J. Mineno, H. Ikeda & H. Shiku. siRNA mediated silencing of PD-1 ligands enhances tumor-specific humam T-cell effector functions. Gene Ther., 19(10) (2012), 959-966.
• ⁹
  T. Kammertoes, T. Shuler & T. Blankestein. Immunotherapy targets the stroma to hit the tumor. Trends Mol. Med., 11(5) (2005), 225-231.
• ¹⁰
  A. López, J. Seoane & M. Sanjuán. A validated mathematical model of tumor growth including tumorhost interaction, Cell mediated immune response and chemotherapy. Bull. Math. Biol., 76 (2014), 2884-2906.
• ¹¹
  S. Palomino, A. Vilcarromero, O. Bonato & J. Fernandes. Co-existência de Espécies em Sistemas Presa-Predador com Switching. TEMA - Tendências em Matemática Aplicada e Computacional, 7(2) (2006), 317-326.
• ¹²
  S. Palomino, D. Coutinho & K. Barbosa. A convex cpproach for controlled Lotka-Volterra multiSpecies model, em Abstract Book CMPD2, “II Conf. Comp. Math. Population Dynamics”, pp. 121, Campinas, (2007).
• ¹³
  S. Rosenberg, Y. Yang & N. Restifo. Cancer immunotherapy moving beyond current vacines. Nat. Med., 10 (2004), 909-915.
• ¹⁴
  A. Sudhakar. History of cancer, ancient and modern treatment methods. J. Cancer Sci. Ther., 1(2) (2009), 1-7.
• ¹⁵
  M. Uzeda & S. Palomino. Análise de estabilidade em interações de tipo hospedeiro-parasita: o caso generalizado, em XXXV CNMAC, “Proceeding Series of the Brazilian Society of Applied and Computational Mathematics”, Natal, RN, 3 (2015), 1-2.

• *
  Trabalho apresentado no XXXVI CNMAC.

• 1
  Stepanova estudou a reação imune do desenvolvimento de um tumor maligno (Biophysics, 24, 917-923).

• 2
  Na modelagem do problema, o comportamento switching se refere ao uso de funções da forma h(w, z) =

• 3
  Na bioqúımica este tipo de função é conhecida como função ativadora (ou repressora) de Hill.

• 4
  Em(12) se generaliza o conceito da interação entre populações de n espécies com comportamento switching. Em(11) se exibe, para o caso de três espécies (duas presas e um predador), como são os gráficos das superfícies sigmoidais nos casos n = 1 e n = 10.

• 5
  Modelo pré-clínico, na terminologia médica indica um tratamento realizado in vitro.

• 6
  O Teorema do valormédio (TVM) para F função diferenciável de Rn em Rm é um resultado bem conhecido: se procura:

• 7
  U é uma a vizinhança de Rn.

• 8
  Por ser o sistema autônomo F(t, X) = F(X).

• 9
  Quando a é unitário o enunciado do TVM se restringe ao uso das normas em questão.

• 10
  Serão estes modelos os que serão sujeitos a a simulação numérica na próxima seção.

• 11
  Pela forma como foram definidos, w, x, y, z são variáveis positivas.

• 12
  Na verdade, E0 = (0, 0, 0, 0) não tem nenhuma interpretação biológica.

• 13
  Essa solução é dada na região dos pontos do plano da forma (t, h(t)).

• 14
  Note que x e z são conhecidos.

• 15
  Na verdade são dois valores para w e y: um por cada valor de x.

• 16
  Os pontos de coexistência são obtidos baixo condições geométricas mostradas no processo de prova.

• 17
  Com quatro elementos.

• 18
  Também, E0 carece de importância por não ter interpretação biológica do problema dado.

• 19
  Apenas um detalhe, ao calcular o polinômio caracteristico de JG , o elemento 41 da matriz JF será trocado pelo valor $a_{41}^G=\frac{n{p}_4{\tau }_c^n{x}^{n-1}}{{\left({\tau }_c^n+\left(1+\frac{D_1}{{\mu }_A}\right)x^n\right)}^2}.$

• 20
  Para encontrar o ponto, nas condições geométricas e algébricas do Lema, foi feita a escolha w = 0.13299, z = 0.0001399. O ponto E2 é desconsiderado pois, para aquela escolha, teve algumas das suas coordenadas negativas.

• 21
  Decresce no tempo (adimensional).

Datas de Publicação

• Publicação nesta coleção
  Sep-Dec 2017

Histórico

• Recebido
  12 Dez 2016
• Aceito
  15 Maio 2017