4.4
Monte Carlo Dinâmico Aplicado ao Modelo de
Crescimento
Para simular o crescimento de um tumor, discretizamos o espaço representando-o através de uma rede bidimensional com M = L × L sítios, e assumimos que cada sítio pode assumir apenas os status tumoral T ou vazio V . Definimos a probabilidade p0 de
que um sítio vazio q seja preenchido dado que ele tenha apenas um vizinho tumoral ao seu lado, sendo, portanto, (1 − p0) a chance de não haver o preenchimento. Dessa forma,
em uma vizinhança contendo ηq tumores, a chance de preenchimento do sítio q por um
de seus vizinhos será dada por
p0+ p0(1− p0) +· · · + p0(1− p0)ηq−1 = 1− (1 − p0)ηq. (4.7)
Com isso, a probabilidade por unidade de tempo de q ser ocupado será porporcional a Eq. (4.7), ou seja, gq ∝ 1 − (1 − p0)ηq. Essa probabilidade por unidade de tempo, gq,
carrega tanto as informações locais quanto as globais do sistema a cada instante; o termo ηq relaciona-se com as informações locais em um dado instante; as informações globais,
por outro lado, serão representadas por p0 e nesse contexto assumimos que p0 = p0(t) =
ω(t)/α, sendo ω(t) a taxa experimental obtida através da Eq. (4.1); consequentemente teremos a taxa de transição para cada sítio vazio com a forma
gq(t) = b 1− β α 1 1 + exp[γ(t− tc)] ηq , (4.8)
sendo b a frequência, no regime de tempos longos, com que um novo sítio T é adicionado ao sistema. Aqui, ηq pode variar entre um e ηmax, com ηmax = 8, ou seja, uma vizinhança
de Moore [80]. A taxa de transição, na literatura sobre propagação epidemica, comumente leva em conta as contribuições locais e globais separadamente [36, 42].
Consideramos que em cada intervalo de tempo ∆t apenas um evento pode ocorrer, e então a distância entre os estados será dada por a = |∆nT| = |∆nV| = 1. As grandezas
nT e nV representam o número de sítios4 com status de tumor, e nV o número de sítios
com status vazio. Podemos escrever uma equação estocástica de movimento, seguindo a
4
O número de sítios é diferente do número de células, portanto, denotamos nT como o número de sítios
na simulação e NT como número de células reais. Alguns trabalhos, como o de Jiang et al. (2005) [35]
4.4 - Monte Carlo Dinâmico Aplicado ao Modelo de Crescimento 33 Eq. (2.6), para esse problema como
d dtnT(t) = X j hg(t)in(j) 0 Pj(t)n (j) 0 , (4.9) em que X
j(·) representa a soma sobre todas as possíveis configurações j do sistema,
acessadas no tempo t, após várias simulações. hg(t)in(j) 0 =
Xn(j)0
q=1gq(t)/n (j)
0 representa
a taxa mesoscópica de aumento do número de sítios T ; com a mesma condição inicial, a probabilidade Pj(t) é gerada após rodarmos M simulações independentes; e n0 é o
número de sítios vazios na interface colônia-meio (alguns desses sítios podem estar no interior da colônia). Omitimos a partir daqui o sobrescrito j, para simplificarmos nossa notação, e quando nos referirmos a uma gradeza específica, estaremos tratando ela em uma configuração j genérica.
Aqui, M = nT + nV, e nV = n0+ ˜nV, ou seja, a quantidade de sítios vazios pode
ser dividida em n0, i.e., sítios vazios que contribuem para o crescimento de nT; e ˜nV,
que são os sítios vazios que não contribuem para o crescimento de nT. Note que assim,
nossos elementos interagentes serão os sítios de {n0}, isto é, N = n0 na Eq. (2.10). A
característica local dos elementos com status V é que para os sítios do conjunto {n0}, o
número de vizinhos tumores é ηq> 0, e para os sítios em {˜nV}, ηq = 0. nV é interpretada
com fonte de nT, ou seja, se um aumenta, o outro diminui e vice-versa [16], no entanto,
como gq(t) = 0 para os sítios do conjunto{˜nV}, apenas os sítios n0participam da dinâmica
do sistema, e, portanto, podem ser escritos como A†
j na Eq. (4.9). Negligenciamos a
transição T → V (morte celular), pois assumimos que as células dispõem de suprimento adequado de nutrientes, e o decréscimo na taxa de crescimento seria pequeno, uma vez que as células adjacentes logo tomariam o espaço deixado pela célula morta [53].
O tempo de espera entre dois eventos pode ser estimado através de uma equação com a forma da Eq. (2.14); para o nosso problema:
∆t(nT)= 1
Pn0
q=1gq(t)
. (4.10)
Como apenas um tipo de evento pode ocorrer nesse sistema, que é o evento de aparecimento de um novo sítio tumoral na colônia, temos fe = 1. O sobrescrito (nT)
denotará o tempo de espera entre o (nT)-ésimo e o (nT + 1)-ésimo sítio. E, finalmente,
4.4 - Monte Carlo Dinâmico Aplicado ao Modelo de Crescimento 34 forma da Eq. (2.11): Hq(t) = gq(t) max [gq(t)] = 1− (1 − p0) ηq 1− (1 − p0)ηmax , (4.11)
em que max [gq(t)] denota o maior valor possível de {gq(t)}. Operacionalmente escolhe-se
com probabilidade (1/n0) um sítio do conjunto {n0}, e tenta-se mudar seu status com a
probabilidade Hq(t). Essa tentativa é feita comparando-se Hq(t) com um número aleatório
ξ, uniformemente distribuído entre [0, 1). Se Hq(t) < ξ a nova configuração é aceita, e de
outra forma ela é rejeitada e o processo recomeça; em caso de aceite da nova configuração, o tempo avança, t → t + ∆t(nT); as variáveis n
T, nV, n0 e {ηq}Moore são atualizadas:
nT → nT + 1, nV → nV − 1 e n0 → n0+ ∆n0. O incremento ∆n0 é uma variável que
depende da localização do sítio q. Quanto maior ηq, menor ∆n0. Uma das possíveis
configurações é quando ηq = ηmax: nesse caso ∆n0 = −1 pois o sítio está em um buraco
e após a atualização de nT ele é retirado do conjunto {n0}. O termo {ηq}Moore representa
uma atualização do número de vizinhos de cada um dos elementos em volta do sítio q; Por fim um novo raio médio é calculado e a taxa X
qgq(t) é atualizada.
A probabilidade de ocupação de um sítio vazio é de Hq/n0 quando t é pequeno,
e para tempos grandes, essa probabilidade é de 1/n0, recobrando um dos procedimentos
clássicos dos modelos tipo Eden [81]. Na Fig. 4.8 mostramos a transição entre duas configurações possíveis. Acima temos a primeira configuração, com um sítio tumoral no centro da rede, e após ∆t = [8b(1−Γ)]−1 o sistema poderá adquirir uma das configurações
mostradas a direita. Nesse caso, como a chance de escolha de um sítio do conjunto {n0} é uniforme, então as probabilidades Pj(∆t(1)) de encontrar o sistema em algum
desses estados após o tempo ∆t será de 1/2. As configurações passam a se complicar e eventualmente o sistema pode adquirir a configuração mostrada abaixo, com um buraco rodeado de sítios cheios. Nesse caso, a chance de escolha do buraco será de 1/17 e Hq = 1,
uma vez que ηq = 8. Na Fig. 4.9 mostramos resumidamente um fluxograma da simulação
4.4 - Monte Carlo Dinâmico Aplicado ao Modelo de Crescimento 35
Figura 4.8: Acima mostramos a configuração inicial do sistema de tumores (quadrados pretos) e sua vizinhança (círculos) em rede. Após ∆t, ele pode assumir algum dos estados indicados. No primeiro caso, teremos ∆n0= 4 e no segunda caso, ∆n0= 2. Abaixo mostramos uma configuração
em que nT = 8 sítios e n0 = 17, essa configuração tem um buraco ao centro. O sistema pode
então assumir a nova configuração mostrada a direita, atualizando ∆n0=−1.
4.4.1
O Cálculo do Raio Médio e da Taxa Mesoscópica
O raio médio de uma colônia, em um instante, é dado por hri =√2d0rg, em que d0é
o diâmetro de uma única célula, e rgé o raio de giro, rg =
r (1/nT) XnT T =1[~rT − ~rcm(t)] 2, ~r T
a localização do sítio T na rede, e ~rcm(t) = (1/nT)
XnT
T =1~rT o centro de massa. Definimos
então S(t) ≡ XnT
T =1[~rT − ~rcm(t)]
4.4 - Monte Carlo Dinâmico Aplicado ao Modelo de Crescimento 36 Início da simulação Escolha de um q Hq(t) > ξ? Atualizar variáveis Fim Valores Médios não sim tempo suficiente após M simulações
Figura 4.9: Fluxograma descrevendo o procedimento computacional adotado.
intervalo de tempo ∆t, podemos calcular a grandeza ~δ, definida por: ~δ ≡ ∆~rcm
= ~rcm(t)− ~rcm(t− ∆t).
(4.12) Para otimizar o cálculo de rg, ao invés de passar por todos os sítios T , vamos considerar
a função S′ (t), dada por S′(t) = S(t − ∆t) + [~rnT − ~rcm(t)] 2. (4.13) Neste contexto, S(t − ∆t) =X(nT−1) T =1 [~rT − ~rcm(t− ∆t)]
2. Usando a Eq. (4.12), podemos
reescrever a Eq. (4.13) como S′(t) = nT−1 X T =1 {~rT − [~rcm(t)− ~δ ]}2+ [~rnT − ~rcm(t)] 2. (4.14)
Após alguma manipulação algébrica, temos S′(t) = nT−1 X T =1 [~rT − ~rcm(t)]2+ [~rnT − ~rcm(t)] 2 + 2~δ· nT−1 X T =1 [~rT − ~rcm(t)] + nT−1 X T =1 ~δ 2. (4.15)
4.5 - Os Ingredientes Básicos da Simulação Computacional 37