# Laboratório 5: Câncer e Tratamento Químico 

### Referente ao capítulo 10

Queremos minimizar a densidade de um tumor em um organismo e os efeitos colaterais das drogas para o tratamento de câncer por quimioterapia em um período de tempo fixo. É assumido que o tumor tenha um crescimento Gompertzian. A hipótese *log-kill* de Skipper é utilizada para modelar a ocisão de células cancerosas. Ela diz que a morte de células devido aos químicos é proporcional à população de tumor. 

Considere $N(t)$ a densidade normalizada do tumor no tempo $t$. Assim:

$$N'(t) = rN(t)\ln\left(\frac{1}{N(t)}\right) - u(t)\delta N(t)$$

onde $r$ é a taxa de crescimento do tumor, $\delta$ a magnitude da dose e $u(t)$ descreve a força do efeito da droga (famacocinética). 

Escolhemos o funcional da seguinte maneira para minimizar tanto a densidade de tumor quanto os efeitos da droga. Supomos aqui que quanto maior a força do efeito da droga, maior seu efeito negativo, o que é razoável. 

$$\min_u \int_0^T aN(t)^2 + u(t)^2 dt$$

Além disso, $u(t) \geq 0$ e $N(0) = N_0$. 

## Condições Necessárias 

### Hamiltoniano

$$
H = aN^2 + u^2 + \lambda\left[rN\ln\left(\frac{1}{N}\right) - u\delta N\right]
$$

### Equação adjunta 

$$
\lambda '(t) = - H_N = -2aN - \lambda\left(r\ln\left(\frac{1}{N}\right) - rN\frac{1}{1/N}\frac{1}{N^2} - u\delta\right) = -2aN + \lambda\left[r-r\ln\left(\frac{1}{N}\right) + u\delta\right]
$$

### Condição de transversalidade 

$$
\lambda(T) = 0
$$

### Condição de otimalidade

$$
H_u = 2u - \lambda \delta N 
$$

Nesse caso temos apenas uma desigualdade, a inferior e, portanto, $H_u < 0$ não é possível. Temos que $\delta > 0$ dada sua interpretação e $N(t) > 0$, afinal sua derivada não é nem definida caso contrário. Como o problema é de minimização, temos que lembrar que as desigualdades ficam invertidas em relação a $H_u$.

$$
H_u > 0 \implies u^*(t) = 0 \implies \lambda \delta N < 0 \implies \lambda(t) < 0 
$$
$$
H_u = 0 \implies 0 \le u^*(t) = \frac{\delta}{2}\lambda(t)N(t)   
$$

Assim $u^*(t) = \max\left\{0, \frac{\delta}{2}\lambda(t)N(t)\right\}$

### Importanto as bibliotecas


```python
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import sympy as sp

import sys  
sys.path.insert(0, '../pyscripts/')

from optimal_control_class import OptimalControl
```

### Usando a biblitoca sympy 


```python
x_sp,u_sp,lambda_sp, r_sp, a_sp, delta_sp = sp.symbols('N u lambda r a delta') 
H = a_sp*x_sp**2 + u_sp**2 + lambda_sp*(r_sp*x_sp*sp.log(1/x_sp) - u_sp*delta_sp*x_sp)
H
```




$\displaystyle N^{2} a + \lambda \left(- N \delta u + N r \log{\left(\frac{1}{N} \right)}\right) + u^{2}$




```python
print('H_x = {}'.format(sp.diff(H,x_sp)))
print('H_u = {}'.format(sp.diff(H,u_sp)))
print('H_lambda = {}'.format(sp.diff(H,lambda_sp)))
```

    H_x = 2*N*a + lambda*(-delta*u + r*log(1/N) - r)
    H_u = -N*delta*lambda + 2*u
    H_lambda = -N*delta*u + N*r*log(1/N)


Resolvendo para $H_u = 0$


```python
eq = sp.Eq(sp.diff(H,u_sp), 0)
sp.solve(eq,u_sp)
```




    [N*delta*lambda/2]



Aqui podemos descrever as funções necessárias para a classe. 


```python
parameters = {'a': None, 'delta': None, 'r': None}

diff_state = lambda t, N, u, par: -N*par['delta']*u + N*par['r']*np.log(1/N)
diff_lambda = lambda t, N, u, lambda_, par: -2*N*par['a'] + lambda_*(par['delta']*u + par['r']*(1 - np.log(1/N)))
update_u = lambda t, N, lambda_, par: np.maximum(0, 0.5*par['delta']*lambda_*N)
```

## Aplicando a classe ao exemplo 

Vamos fazer algumas exeperimentações. Sinta-se livre para variar os parâmetros. Nesse caso passaremos o limite inferior como padrão inicial. Note que `np.inf` deve ser escrito quando não há limite. 


```python
problem = OptimalControl(diff_state, diff_lambda, update_u, bounds = [(0,np.inf)])
```


```python
N0 = 0.975
T = 20
parameters['a'] = 3
parameters['delta'] = 0.45
parameters['r'] = 0.3
```


```python
t,x,u,lambda_ = problem.solve(N0, T, parameters)
ax = problem.plotting(t,x,u,lambda_)
```

Note o formato do controle: Inicia com valor alto e depois reduz conforme se aproxima do dia 20. Isto é consistente com as práticas médicas atuais. Todavia, também observamos que o tumor tende a crescer após o dia 12. Por isso, vamos aumentar a importância de reduzir as células cancerosas, em troca de mais danos pelos remédios. Vamos notar que a força da droga é muito maior quando reduzimos a importância de seu efeito negativo. O formato das curvas é, entretanto, muito similar. 


```python
parameters['a'] = 10
t,x,u,lambda_ = problem.solve(N0, T, parameters)
ax = problem.plotting(t,x,u,lambda_)
```

Vamos comparar o controle escolhido conforme diferentes valor iniciais da densidade $N_0$. 


```python
parameters['a'] = 3
N0_values = [0.1, 0.5, 0.7, 0.9]
u_values = []
for N0 in N0_values: 
    _,_,u,_ = problem.solve(N0, T, parameters)
    u_values.append(u)
    
fig = plt.figure(figsize = (10,5))
plt.xlabel("Tempo")
plt.ylabel("Químico")
plt.title("Quimioterapia utilizada")
for i, N0 in enumerate(N0_values):
    plt.plot(t, u_values[i],label = r'$N_0$ = {}'.format(N0))
plt.legend()
plt.grid(alpha = 0.5)
```

Podemos observar acima que apenas o estágio inicial do controle é afetado pela densidade inicial, até que atinge uma espécie de equilíbrio. 

Outro ponto importante é o quanto a magnitude da dose influencia na dosagem do químico ao longo do tempo. Observe que quanto maior a magnitude, o valor inicial é mais alto, porém o decaimento é mais rápido, usando menos químico no total, como era esperado. 


```python
N0 = 0.8
delta_values = [0.25, 0.5, 0.75]
u_values = []
for delta in delta_values: 
    parameters['delta'] = delta
    _,_,u,_ = problem.solve(N0, T, parameters)
    u_values.append(u)
    
fig = plt.figure(figsize = (10,5))
plt.xlabel("Tempo")
plt.ylabel("Químico")
plt.title("Quimioterapia utilizada")
for i, delta in enumerate(delta_values):
    plt.plot(t, u_values[i],label = r'$\delta$ = {}'.format(delta))
plt.legend()
plt.grid(alpha = 0.5)
```

Podemos ver também diferentes valores temporais $T$ de duração do tratamento alteram o comportamento da quantidade de células cancerosas. Na verdade, o formato fica muito similar em todos os casos, variando o tempo de estabilidade no meio do caminho. 


```python
parameters['delta'] = 0.4
T_values = [10, 20, 40, 80]
x_values = []
u_values = []
t_values = []
for T in T_values: 
    t,x,u,_ = problem.solve(N0, T, parameters)
    t_values.append(t)
    x_values.append(x)
    u_values.append(u)
    
fig, ax = plt.subplots(2,2,figsize = (10,8))
fig.suptitle(r'Comparação temporal de $N(t)$')
for k in range(4):
    i = k//2
    j = k%2
    ax[i][j].set_ylabel(r"$N(t)$")
    ax[i][j].set_title("Densisade de células cancerosas")
    ax[i][j].plot(t_values[k], x_values[k])
    ax[i][j].grid(alpha = 0.5)

fig, ax = plt.subplots(2,2,figsize = (10,8))
fig.suptitle(r'Comparação temporal de $u(t)$')
for k in range(4):
    i = k//2
    j = k%2
    ax[i][j].set_ylabel(r"$u(t)$")
    ax[i][j].set_title("Força do efeito da quimioterapia")
    ax[i][j].plot(t_values[k], u_values[k], color = 'green')
    ax[i][j].grid(alpha = 0.5)
```

Você deve ter observado em todos os exemplos que $u(t)$ atinge 0 apenas pontualmente. É necessária a restrição $u(t) \ge 0$? 

## Experimentação 


```python
#N0 = 1
#T = 5
#parameters['r'] = 0.3
#parameters['a'] = 10
#parameters['delta'] = 0.4
#
#t,x,u,lambda_ = problem.solve(N0, T, parameters)
#roblem.plotting(t,x,u,lambda_)
```

### Este é o final do notebook

5a776214ccedf4b2ee94d7f932cfc25f8dbdb0bc
ipynb
Jupyter Notebook
notebooks/.ipynb_checkpoints/Laboratory5-checkpoint.ipynb
lucasmoschen/optimal-control-biological
642a12b6a3cb351429018120e564b31c320c44c5
MIT
2021-11-03T16:27:39.000Z
2021-11-03T16:27:39.000Z
notebooks/.ipynb_checkpoints/Laboratory5-checkpoint.ipynb
lucasmoschen/optimal-control-biological
642a12b6a3cb351429018120e564b31c320c44c5
MIT
notebooks/.ipynb_checkpoints/Laboratory5-checkpoint.ipynb
lucasmoschen/optimal-control-biological
642a12b6a3cb351429018120e564b31c320c44c5
MIT
Qwen/Qwen-72B
1. YES
2. YES