Modelo de Keller-Segel para relação população-economia

Grupo: Leonardo Barcelos, Luana Bianchi e Rubens Borrasca

O objetivo deste trabalho é implementar o modelo de Keller-Segel ^[1], que originalmente descreve o movimento de células em direção a um sinal químico, para um sistema englobando população e atividade econômica, problema proposto no livro Introduction to the Modeling and Analysis of Complex Sistems de Hiroki Sayama ^[2]. O método computacional utilizado para resolver o problema e implementar o modelo foi o FTCS (Forward Time Centered Space), método de resolução de equações diferenciais parciais, como a equação do calor e a da difusão, que é o caso deste trabalho.

Modelo de Keller-Segel

Proposto por Evelyn Fox Keller, física norte-americana, e Lee Aaron Segel, matemático também norte-americano, o modelo de Keller-Segel foi historicamente utilizado para descrever o movimento de bactérias. Introduzido primeiramente em 1970 para descrever a agregação de uma espécie de bolor limoso (ou slime mold) ameboide, Dictyostelium discoideum, o modelo se tornou um dos mais usados nos estudos biológicos-matemáticos. As células deste slime mold se comportam como amoebas individuais, e se alimentam de bactérias, mas quando a quantidade de comida fica pequena, elas se difundem pelo espaço e então se agregam em formato mais alongado, como o formato das lesmas, para uma migração de longa distância. Keller e Segel desenvolveram um modelo matemático para o processo de agregação, em que a chemotaxis, a taxa com que as células se movem em direção ao sinal químico, tem papel crítico na auto-organização das células. ^[3]

Baseados no que já era conhecido sobre esses organismos, Keller e Segel utilizaram as seguintes premissas ^[2]:

As células estão inicialmente distribuídas sobre o espaço de maneira mais ou menos homogênea, com algumas flutuações aleatótias;
As células apresentam chemotaxis em direção ao sinal químico denominado cAMP (cyclic adenosine monophosphate);
As células produzem moléculas cAMP;
As células e as moléculas cAMP difundem pelo espaço;
As células não morrem e não se dividem

De forma simplificada, ocultando alguns detalhes biológicos mais complicados a equação de Keller-Segel é a seguinte:

{\frac {\partial a}{\partial t}}=\mu \nabla ^{2}a-\chi \nabla \cdot (a\nabla c)

{\frac {\partial c}{\partial t}}=D\nabla ^{2}c+fa-kc

em que $a$ e $c$ são respectivamente as variáveis de estado para a concentração de células e a concentração de cMAP. $\mu$ é o parâmetro de mobilidade das células, $\chi$ é o parâmetro da chemotaxis celular, $D$ é a constante de difusão das moléculas cAMP, $f$ é a taxa de secreção de cMAP pelas células, e $k$ é a taxa de decaimento das moléculas cMAP.

Aplicação população-economia

De forma parecida com as premissas de Keller e Segel, os seguintes pontos são assumidos para modelar a relação entre a população e a atividade econômica^[2]:

A população não cresce e não decresce ao longo do tempo;
A economia é ativada por existir mais pessoas em uma região;
Sem pessoas a atividade econômica diminui;
População e atividade econômica difundem gradualmente;
As pessoas são atraídas por regiões com maior atividade econômica

O trabalho agora é traduzir essas premissas em equações. O primeiro ponto diz que que padrão não há mudança na população, dessa forma:

{\frac {\partial p}{\partial t}}=0

O segundo ponto nos diz que a economia tem um aumento que depende da quantidade de pessoas, e assim:

{\frac {\partial m}{\partial t}}=\alpha p

em que $\alpha$ é a taxa de aumento da atividade econômica per capita.

O terceiro ponto diz que a atividade econômica decai quando não há pessoas, e para isso adicionamos um termo de decaimento na equação:

{\frac {\partial m}{\partial t}}=\alpha p-\beta m

em que $\beta$ é a taxa de decaimento da atividade econômica.

O quarto ponto fala que tanto população quanto a atividade econômica difundem, o que uma premissa sobre o movimento espacial. Para isso são adicionados termos de difusão, os Laplacianos, as equações:

{\frac {\partial p}{\partial t}}=D_{p}\nabla ^{2}p

{\frac {\partial m}{\partial t}}=D_{m}\nabla ^{2}m+\alpha p-\beta m

onde $D_{p}$ e $D_{m}$ são as constantes de difusão da população e da atividade econômica respectivamente.

O quinto e último ponto fala que as pessoas são atraídas pela atividade econômica e assim se movem para áreas que possuem mais dinheiro. Para cumprir essa premissa, a equação do transporte é utilizada:

-\nabla \cdot (p\gamma \nabla m)

em que o gradiente de $m$ é utilizado para obter a velocidade média do movimento da população, sendo $\gamma$ uma constante que afeta essa velocidade.

O sistema de equação final fica da seguinte forma:

{\frac {\partial p}{\partial t}}=D_{p}\nabla ^{2}p-\gamma \nabla \cdot (p\nabla m)

{\frac {\partial m}{\partial t}}=D_{m}\nabla ^{2}m+\alpha p-\beta m

Comparando o sistema obtido com o problema original de Keller-Segel, percebe-se que se trocarmos células por pessoas e cMAP por atividade econômica os problemas ficam iguais, e até se poderia denominar como moneytaxis a migração das pessoas em direção a atividade econômica, como a chemotaxis descreve o movimento das células em direção ao cAMP.

Método FTCS

O FTCS (Forward Time Centered Space, em tradução livre significa "avançado no tempo, centrado no espaço), é um método de discretização de Equações Diferenciais Parciais(EDP). Para a derivada temporal teremos, ^[4]

{\frac {\partial f}{\partial t}}\rightarrow {\frac {f^{n+1}-f^{n}}{\Delta t}}

e para a parte espacial,

{\frac {\partial ^{2}f}{\partial s^{2}}}\rightarrow {\frac {f_{i-1}-2f_{i}+f_{i+1}}{(\Delta s)^{2}}}

onde $s$ é uma variável espacial qualquer $(x,y,z,...)$ e $t$ é o tempo.

Discretização do Modelo de Keller-Segel em 1D

Em 1D o sistema de equações diferenciais parciais será:

{\frac {\partial p}{\partial t}}=D_{p}{\frac {\partial ^{2}p}{\partial x^{2}}}-\gamma \left[{\frac {\partial p}{\partial x}}{\frac {\partial m}{\partial x}}+p{\frac {\partial ^{2}m}{\partial x^{2}}}\right]

{\frac {\partial m}{\partial t}}=D_{m}{\frac {\partial ^{2}m}{\partial x^{2}}}+\alpha p-\beta m

Agora utilizando a discretização FTCS teremos:

{\frac {p_{i}^{n+1}-p_{i}^{n}}{\Delta t}}={\frac {D_{p}}{(\Delta x)^{2}}}\left[p_{i-1}^{n}-2p_{i}^{n}+p_{i+1}^{n}\right]-{\frac {\gamma }{(\Delta x)^{2}}}\left[(p_{i+1}^{n}-p_{i}^{n})(m_{i+1}^{n}-m_{i}^{n})+p_{i}^{n}(m_{i-1}^{n}-2m_{i}^{n}+m_{i+1}^{n})\right]

{\frac {m_{i}^{n+1}-m_{i}^{n}}{\Delta t}}={\frac {D_{m}}{(\Delta x)^{2}}}\left[m_{i-1}^{n}-2m_{i}^{n}+m_{i+1}^{n}\right]+\alpha p_{i}^{n}-\beta m_{i}^{n}

onde o sub-índice $i$ se refere à coordenada $x$ ; e o superíndice $n$ se refere ao tempo. Reorganizando as equações e agrupando alguns termos teremos:

p_{i}^{n+1}=p_{i}^{n}\left[1-2K_{1}-K_{2}\left(m_{i-1}^{n}-m_{i}^{n}\right)\right]+K_{1}\left[p_{i-1}^{n}+p_{i+1}^{n}\right]-K_{2}\left[p_{i+1}^{n}(m_{i+1}^{n}-m_{i}^{n})\right]

m_{i}^{n+1}=m_{i}^{n}\left[1-K_{3}-\lambda \right]+K_{3}\left[m_{i-1}^{n}+m_{i+1}^{n}\right]+Vp_{i}^{n}

onde os termos agrupados são: $K_{1}={\frac {D_{p}\Delta t}{(\Delta x)^{2}}}$ , $K_{2}={\frac {\gamma \Delta t}{(\Delta x)^{2}}}$ , $K_{3}={\frac {D_{m}\Delta t}{(\Delta x)^{2}}}$ , $V=\alpha \Delta t$ , $\lambda =\beta \Delta t$

Discretização do Modelo de Keller-Segel em 2D

Em 2D o sistema de equações diferenciais parciais será:

{\frac {\partial p}{\partial t}}=D_{p}\left[{\frac {\partial ^{2}p}{\partial x^{2}}}+{\frac {\partial ^{2}p}{\partial y^{2}}}\right]-\gamma \left[{\frac {\partial p}{\partial x}}{\frac {\partial m}{\partial x}}+{\frac {\partial p}{\partial y}}{\frac {\partial m}{\partial y}}+p\left({\frac {\partial ^{2}m}{\partial x^{2}}}+{\frac {\partial ^{2}m}{\partial y^{2}}}\right)\right]

{\frac {\partial m}{\partial t}}=D_{m}\left[{\frac {\partial ^{2}m}{\partial x^{2}}}+{\frac {\partial ^{2}m}{\partial y^{2}}}\right]+\alpha p-\beta m

Agora utilizando a discretização FTCS e assumindo que $\Delta x=\Delta y=\Delta s$ teremos:

{\frac {p_{i,j}^{n+1}-p_{i,j}^{n}}{\Delta t}}={\frac {D_{p}}{(\Delta s)^{2}}}\left[(p_{i-1,j}^{n}-2p_{i,j}^{n}+p_{i+1,j}^{n})+(p_{i,j-1}^{n}-2p_{i,j}^{n}+p_{i,j+1}^{n})\right]-{\frac {\gamma }{(\Delta s)^{2}}}\left[(p_{i+1,j}^{n}-p_{i,j}^{n})(m_{i+1,j}^{n}-m_{i,j}^{n})+(p_{i,j+1}^{n}-p_{i,j}^{n})(m_{i,j+1}^{n}-m_{i,j}^{n})+p_{i,j}^{n}\left((m_{i-1,j}^{n}-2m_{i,j}^{n}+m_{i+1,j}^{n})+(m_{i,j-1}^{n}-2m_{i,j}^{n}+m_{i,j+1}^{n})\right)\right]

{\frac {m_{i,j}^{n+1}-m_{i,j}^{n}}{\Delta t}}={\frac {D_{m}}{(\Delta s)^{2}}}\left[(m_{i-1,j}^{n}-2m_{i,j}^{n}+m_{i+1,j}^{n})+(m_{i,j-1}^{n}-2m_{i,j}^{n}+m_{i,j+1}^{n})\right]+\alpha p_{i,j}^{n}-\beta m_{i,j}^{n}

onde os sub-índices $i$ e $j$ se referem às coordenadas $x$ e $y$ respectivamente; e o superíndice $n$ se refere ao tempo. Reorganizando as equações e agrupando alguns termos teremos:

p_{i,j}^{n+1}=p_{i,j}^{n}\left[1-4K_{1}-K_{2}\left(m_{i-1,j}^{n}-2m_{1,j}^{n}+m_{i,j-1}^{n}\right)\right]+K_{1}\left[p_{i-1,j}^{n}+p_{i,j-1}^{n}+p_{i+1,j}^{n}+p_{i,j+1}^{n}\right]-K_{2}\left[p_{i+1,j}^{n}(m_{i+1,j}^{n}-m_{i,j}^{n})+p_{i,j+1}^{n}(m_{i,j+1}^{n}-m_{i,j}^{n})\right]

m_{i,j}^{n+1}=m_{i,j}^{n}\left[1-4K_{3}-\lambda \right]+K_{3}\left[m_{i-1,j}^{n}+m_{i,j-1}^{n}+m_{i+1,j}^{n}+m_{i,j+1}^{n}\right]+Vp_{i,j}^{n}

Resultados

1D

Com o intuito de explorar a equação e suas consequências, os resultados foram divididos em várias simulações diferentes.

Para todas as simulações realizadas, exceto onde indicado, os parâmetros utilizados foram os seguintes:

$N=100$

N indica o limite do eixo x. Quanto maior o seu valor, mais demorada fica a simulação. Porém, valores pequenos de N podem gerar resultados inconclusivos.

$\Delta x=1$

Indica a unidade infinitesimal do eixo x (neste caso, o único eixo do problema).

$\Delta t=0.3$

Indica quanto tempo cada passo da simulação percorre.

$D_{m}=1$

Constante de difusão da função de renda.

$D_{p}=1$

Constante de difusão da função de densidade populacional.

$\gamma =1$

$\gamma$ é um parâmetro intimamente ligado com o movimento da população para regiões com alta concentração de dinheiro. Nas palavras de Sayama,^[2] indica a "habilidade que as pessoas da rede tem de 'sentir o cheiro' do dinheiro".

$\alpha =1$

O parâmetro $\alpha$ é utilizado no cálculo de $v$ , como descrito acima, seguindo a fórmula $v=\alpha \Delta t$ . Está relacionado com a velocidade do progresso da economia, ou seja, quão rápido cresce a renda em uma posição da malha através da influência da população.

$\beta =1$

$\beta$ é utilizado na fórmula $\lambda =\beta \Delta t$ , que indica a taxa de decaimento da renda em regiões da malha com baixa população.

Estes parâmetros foram utilizados pois representam bem a característica do sistema de Keller-Segel que queremos representar: a formação de clusters de população na malha ^[1]. Uma nova combinação de parâmetros poderia resultar em resultados diferentes, mas que fogem do escopo deste projeto.

Além disso, foram utilizadas condições periódicas de contorno (PBC) para a solução das equações diferenciais parciais. Deste modo, pode-se pensar no eixo x como um anel, onde, considerando um sistema de tamanho $N$ , os pontos $x=0$ e $x=N$ estão conectados.

População e Dinheiro em pontos separados

Para esta simulação, considera-se que no tempo 0, toda a população está concentrada em 1 ponto $x={\mathcal {C}}_{1}$ , enquanto todo o dinheiro está em um outro ponto, distante deste, $x={\mathcal {C}}_{2}$ . Deste modo, temos as seguintes equações para as condições iniciais:

$p(x,t=0)=\left\{{\begin{array}{lc}1,\quad {\text{p/}}\quad x={\mathcal {C}}_{1},{\mathcal {C}}_{1}\in [0,N]\\0,\quad {\text{caso contrario}}\end{array}}\right.$

$m(x,t=0)=\left\{{\begin{array}{lc}1,\quad {\text{p/}}\quad x={\mathcal {C}}_{2},{\mathcal {C}}_{2}\in [0,N],{\mathcal {C}}_{2}\neq {\mathcal {C}}_{1}\\0,\quad {\text{caso contrario}}\end{array}}\right.$

Na figura abaixo, consegue-se observar o resultado da construção do sistema desta maneira:

Resultados da simulação para o caso de população e dinheiro em pontos separados e distantes na malha

Com toda a população concentrada em 1 ponto ( ${\mathcal {C}}_{1}=20$ ), a atividade econômica cresce consideravelmente neste intervalo ao longo do tempo. Em contrapartida, o local que continha todo o dinheiro no começo da simulação ( ${\mathcal {C}}_{2}=80$ ), em pouco tempo tem a sua renda líquida reduzida, pois o dinheiro decai a uma taxa $\beta$ , como já falado. Essa tendência indica, portanto, que o sistema é construído de tal forma que a atração da população por regiões de alta renda líquida é menor que a taxa de evolução do sistema monetário em pontos de alta densidade populacional.

Além disso, outra observação interessante é que nota-se para $t\rightarrow \infty$ uma tendência inerente da densidade populacional em seguir uma distribuição de formato gaussiano sob a malha. Por isso, é interessante analisar a solução estacionária (se existir).

Como as equações diferenciais do sistema são dependentes umas das outras (isto é, a equação de m(x,t) depende de p(x,t) e a de p depende de m), a solução não é trivial. Utilizando o Método de Coeficientes a Determinar para solução de equações diferenciais em estado estacionário, conseguimos a seguinte relação genérica, que se mostra útil:

$m(x,t)=A_{1}e^{-B_{1}x^{2}}+A_{2}p(x,t)+A_{3}$

Onde $A_{1},A_{2},A_{3}$ e $B_{1}$ são constantes que dependem das condições iniciais do problema. Como o objetivo do trabalho não é resolver analiticamente esta EDP, estas constantes não foram calculadas.

Mas a equação acima mostra que, para o estado estacionário, o formato gaussiano obtido era previsto na solução analítica também. Portanto, o resultado obtido apresenta coerência.

Mais um ponto a ser aferido é que, depois de se desfazer de seu formato inicial, o total de dinheiro sob a malha tende a seguir a distribuição populacional, porém com um desvio padrão maior (maior largura na Gaussiana). Essa observação indica que, para centros econômicos (regiões com alto $m$ ) a tendência é que suas periferias também possuam valores altos de renda, apesar da população consideravelmente menor. Além disso, para regiões fora do contorno de centros econômicos (distância maior do que 3 vezes o desvio padrão da gaussiana) a atividade econômica é basicamente nula, assim como a densidade populacional. Este último fato descreve de forma genérica e simplista o comportamento atual observado em metrópoles nos dias de hoje: uma cidade grande possui alto número de habitantes, alta renda, seus contornos também apresentam atividade econômica forte (porém menor que o centro), mas para um raio suficientemente grande, tanto dinheiro quanto população caem exponencialmente.

População uniforme e sem dinheiro no sistema para t = 0

Nesta simulação, considera-se que, para t=0, não há dinheiro sob a malha. Deste modo, a equação que descreve o dinheiro no sistema ao longo do tempo pode ser escrita como:

$m(x,t=0)=0,\forall x\in [0,N]$

Além disso, a população é iniciada de forma aleatória sob a malha. Deste modo, não há tendência inicial à formação de centros com alta densidade de população.

Um cuidado importante na iniciação deste sistema é que, como estamos modelando a densidade populacional, e consideramos que não há perda de população com o tempo, a integral da função p(x,t) deve ser sempre igual a 1. Em t=0, na iniciação do sistema, é necessário implementar esta condição ao problema, para que ela se mantenha durante a simulação. Para isso, basta sortear valores e normalizá-los depois. Um exemplo de código simples em Python que realiza esta operação está indicado abaixo.

import numpy as np
N = 100   #Limite no eixo x
p_auxiliar = np.random.rand(N)   #Vetor auxiliar
p = p_auxiliar / np.sum(p_auxiliar)   #Normalização, para que a soma de p em todo x seja 1

Esta condição inicial da população se assemelha muito à proposição inicial que Keller e Segel fizeram para um sistema celular, como descrito acima. Temos uma concentração homogênea de $p$ com pequenas flutuações ao longo do eixo.

Resultados da simulação para o caso de m(x,t=0)=0, e população iniciada aleatoriamente.

Na imagem acima, para t = 0 (início da simulação) compreende-se melhor as condições iniciais do sistema. Enquanto que a população, aleatoriamente distribuída sob o eixo x, se assemelha a um ruído branco, o dinheiro não existe na malha.

Na segunda coluna de imagens, nota-se um ponto interessante: a formação de clusters de população (e consequentemente, de dinheiro). Estes clusters são, na verdade, picos que aparecem no gráfico de p(x,t), e indicam alta concentração da população em pontos específicos. Além dos picos claramente visíveis (um deles próximo a x=50, e outro próximo a x=90) pode-se notar, também, sub-picos nas bases destes picos de população. Para t=24.9 e dt=0.3, deduz-se que o sistema, nesta representação, havia passado por 83 iterações até então, o que indica que, durante estas 83 iterações, haviam mais clusters de população em tempos passados, menores porém definidos. E estes "mini-clusters" se agruparam até formar os 2 picos que vemos.

Com o passar do tempo na simulação, nota-se que o comportamento continua, de modo que para a última coluna de figuras é visível que apenas 1 dos picos iniciais se manteve, enquanto o outro foi praticamente inteiro "engolido" pela cauda do maior. Este comportamento de formação de 1 único cluster de população, com formato gaussiano, já havia sido observado na simulação anterior, para $t$ suficientemente grande.

Simulação do mesmo sistema anterior, para t suficientemente grande (t = 270 neste caso) a ponto de chegar em um estado próximo ao estacionário, onde as funções que descrevem população e renda do sistema praticamente não se alteram mais com o tempo.

A figura acima mostra o que acontece caso deixemos o mesmo sistema apresentado antes evoluir até um estado de equilíbrio, onde não há alterações para a população ou renda do sistema. Neste caso, observa-se com mais clareza uma curva de formato gaussiano, em localização bem próxima àquela que vimos para t=124.8 na figura anterior, tanto para a distribuição da população quanto da renda. E mais uma vez, mesmo que não muito perceptível pois as 2 curvas apresentadas são bem largas, o desvio padrão da curva que descreve a renda aparenta ser maior que o desvio padrão da curva que descreve a população.

População e renda iniciados aleatoriamente na rede

Sob a perspectiva de testar a formação de clusters de população, foi criada uma rede completamente aleatória, sem nenhum viés, seja populacional ou econômico. Para isto, tanto população quanto renda possuem valores aleatórios para $t=0$

Mais uma vez, como iniciamos aleatoriamente a população, é importante normalizá-la para que a integral seja igual a 1, como dito acima.

Resultados da simulação para rede iniciada com população e renda aleatórias em t = 0.

Enquanto que a primeira coluna de gráficos mostra o estado inicial aleatório do sistema, a segunda coluna (t = 9.9) indica exatamente o esperado: 5 clusters bem visíveis são formados para a população. Nota-se, entretanto, que os mesmos 5 picos não são tão visíveis no gráfico de rendas. Uma explicação para isto seria que, mais uma vez, com os valores de desvio padrão da distribuição mais elevados, um pico acaba se sobrepondo a outro, tornando as distribuições difícies de se distinguir.

Conforme a simulação avança no tempo, é perceptível mais uma vez a tendência entre os picos de se mergirem. E, como visto anteriormente, para $t\rightarrow \infty$ , a tendência é que o sistema entre em um estado com somente um pico de formato gaussiano na população, e um pico na renda, com largura maior.

Com isso, conclui-se que a formação destes clusters é, de fato, inerente ao sistema, e consequência do modelo utilizado.

2D

Para o caso em duas dimensões, foi utilizada uma distribuição populacional uniforme em todo o espaço. Já a distribuição econômica, no instante t=0 começou da seguinte forma: Em cada canto do espaço foi atribuído um valor de 0.125, no centro 2 e ao redor do centro em 4 pontos 1. A seguir, é confirmado um comportamento que foi observado no caso unidimensional, em que os picos concentrados, após a evolução do sistema, tomam a forma de gaussianas. É possível notar também que a população tende a "clusterizar" em torno dos locais em que a atividade econômica tinha valores altos no início da simulação. Isso se deve principalmente ao termo que é influenciado pela constante $\gamma$ do sistema de EDPs que modela o sistema. Além disso, foi utilizado PBC.

Evolução da atividade econômica e da população para população inicial uniformemente distribuída

A seguir, foi gerada uma animação com a evolução do sistema até a estabilização da atividade econômica. A estabilidade da atividade econômica foi entendida como $max(m_{i,j}^{n+10}-m_{i,j}^{n})\leqslant \epsilon$ , onde $\epsilon$ é o valor que regula o erro. Para este caso $\epsilon =1\times 10^{-9}$ .

Animação da evolução do sistema

Os parâmetros utilizados para gerar as imagens foram os seguintes: $L=100$ , $ds=1$ , $dt=0.3$ , $D_{m}=0.5$ , $D_{p}=0.5$ , $\alpha =1.2$ , $\beta =0.03$ , $\gamma =1$

Para a geração da evolução do sistema e do gif, o código utilizado foi o seguinte:

L = 100 # Tamanho do Grid
D_p = 0.5 # Coeficiente de difusão de pessoas
D_m = 0.5 # Coeficiente de difusão de dinheiro
ds = 1 # Diferencial Espacial
dt = 0.3 # Diferencial Temporal
alfa = 1.2 # Taxa de produção de economia per capita
beta = 0.03 # Taxa de decaimente da economia
gamma = 1 # Taxa de velocidade com que as pessoas migram em direção ao dinheiro

# CRIA OBJETO DE PARÂMETROS PARA O MODELO #
parametros = ParametrosKellerSegelModel(L, L, D_p, D_m, ds, dt, alfa, beta, gamma)

# INICIALIZA O MODELO COM OS PARÂMETROS #
modelo = KellerSegelModel(parametros)

# GERA CONDIÇÕES INICIAIS A SEREM ESTUDADAS #
condicao_inicial_populacao = np.matrix(np.full((parametros.N_x, parametros.N_x), 1 / (L ** 2)))
condicao_inicial_dinheiro = np.matrix(np.zeros((parametros.N_y, parametros.N_y)))

condicao_inicial_dinheiro[(0, 0)] = 0.125
condicao_inicial_dinheiro[(0, 99)] = 0.125
condicao_inicial_dinheiro[(99, 99)] = 0.125
condicao_inicial_dinheiro[(99, 0)] = 0.125
condicao_inicial_dinheiro[(49, 49)] = 2
condicao_inicial_dinheiro[(49, 69)] = 1
condicao_inicial_dinheiro[(49, 29)] = 1
condicao_inicial_dinheiro[(69, 49)] = 1
condicao_inicial_dinheiro[(29, 49)] = 1

# SETA AS CONDIÇÕES INICIAIS DO MODELO #
modelo.setEstadoInicial(condicao_inicial_populacao, condicao_inicial_dinheiro)

# GERA O OBJETO DE GRÁFICO ESTÁTICO
jpeg = JpegTool(nome_imagem = 'comparacao_tempos_pop_eco_2d')

# GERA GRÁFICOS
jpeg.geraJpeg(modelo)

# REINICIALIZA MODELO
modelo = KellerSegelModel(parametros)

# SETA AS CONDIÇÕES INICIAIS DO MODELO #
modelo.setEstadoInicial(condicao_inicial_populacao, condicao_inicial_dinheiro)

# GERA O OBJETO DE ANIMACAO #
ani = AnimacaoTool(nome_gif = 'dinamica_pop_eco_2d')

# GERA GIF COM A EVOLUCAO DO MODELO #
ani.geraGif(modelo)

Código completo se encontra no Github.

Discussão

Nota-se que, independente do caso 1D ou 2D, as conclusões tiradas dos modelos construídos se assemelham bastante. Assim como previsto pela teoria de Keller e Segel, base deste estudo, observa-se a formação de clusters de densidade populacional e de renda a partir do momento em que se inicia a simulação, para um intervalo de tempo suficientemente grande a ponto de que estes sejam formados. E essa tendência independe das condições iniciais apresentadas.

Nas figuras abaixo, sendo a primeira uma figura retirada do livro do Sayama ^[2], e a segunda gerada utilizando o modelo construído, a fim de comparar com a primeira, observa-se, partindo de uma população uniforme na malha e de uma renda aleatoriamente distribuída, a formação de grupos (pontos mais escuros), deste modo provando que o modelo construído chega nos mesmos resultados que a referência.

Fonte: Sayama

Clusters obtidos para os casos estudados, sob condições iniciais semelhantes às aplicadas por Sayama na imagem acima. Na imagem, observa-se, assim como na referência, o principal objetivo da aplicação do modelo: a presença de clusters de população sob a malha.

Também foi observado, para ambas as análises, a tendência da função que representa a renda da rede em seguir o comportamento da função que descreve a densidade populacional, algo que também é previsto pelo modelo natural de Keller e Segel, visto que cada foco de população é considerado também um foco de criação de renda por definição.

Por último, apesar deste modelo não ser ideal na simulação de sistemas reais população-renda, pois deixa de fora várias variáveis importantes para manter um modelo realista, algumas conclusões obtidas aqui podem ser observadas no mundo real. É o caso da distribuição populacional gaussiana por exemplo, que modela convincentemente a população, a nível macro-geográfico, em centros urbanos, assim como a concentração monetária ser intimamente ligada com a concentração da população.

Um estudo mais aplicado sob a estabilidade do método, algo que não é ponto de foco deste trabalho, pode possibilitar a exploração de novas combinações de parâmetros do modelo, vindo a surgir análises sob pontos de vista diferentes.

Programas

Código disponível da plataforma GitHub[2]

Modelagem unidimensional

import numpy as np
import matplotlib.pyplot as plt

def FTCS(p,m,k1,k2,l,k3,v):
    N = p.shape[0]
    
    p_ = np.zeros(N)
    m_ = np.zeros(N)
    
    for j in np.arange(0,N):
        back = (j-1)%N
        forward = (j+1)%N
        p_[j] = (1 - 2*k1 - k2*(m[back] - m[j]))*p[j] + k1*(p[back]+p[forward]) - k2*(m[forward] - m[j])*p[forward]
        m_[j] = (1 - l - k3)*m[j] + k3*(m[back] + m[forward]) + v*p[j]
    return p_,m_

def plot(p,m,filename):
    #Plota o estado final do sistema em um arquivo
    fig,ax = plt.subplots(2,1)
    
    #População
    ax[0].set_title("População")
    ax[0].plot(p)

    #Renda
    ax[1].set_title("Renda")
    ax[1].plot(m)

    plt.savefig(filename + ".png")

##########################################################################################
"""DECLARAÇÃO DE CONSTANTES"""
N = 100  #Tamanho do eixo x
dx = 1
dt = 0.3
Dm = 1.0  #Constante de difusão p/ renda.
Dp = 1.0  #Constante de difusão p/ população.
gamma = 1.0
alpha = 1.0
beta = 1.0

k1 = Dp*dt/(dx**2)
k2 = gamma*k1/Dp
k3 = Dm*dt/(dx**2)
v = alpha*dt
l = beta*dt
##########################################################################################

print("=================================================================================")
print("Caso 1: População e dinheiro em pontos separados")
print("População começa totalmente concentrada em 1 ponto da rede.")
print("Dinheiro começa totalmente focado em 1 ponto da rede longe do ponto da população.")
print("=================================================================================")

N_simulations = 500  #Número de simulações desejado.
T = N_simulations * dt  #Tempo máximo de duração da simulação.

#Setando população para t = 0
p = np.zeros(N)
p[20] = 1
    
#Setando renda para t = 0
m = np.zeros(N)
m[80] = 1

for t in np.arange(0,T,dt):
    p,m = FTCS(p,m,k1,k2,l,k3,v)

#Criação da figura com o estado final
plot(p,m,"caso_2")

##########################################################################################
print("=================================================================================")
print("Caso 2: Sem dinheiro p/ t = 0")
print("População começa totalmente concentrada em 1 ponto da rede.")
print("Dinheiro começa zerado")
print("=================================================================================")

N_simulations = 500  #Número de simulações desejado.
T = N_simulations * dt  #Tempo máximo de duração da simulação.

#Setando população para t = 0
p1 = np.random.rand(N)
p = p / np.sum(p)

#Setando renda para t = 0
m = np.zeros(N)

for t in np.arange(0,T,dt):
    p,m = FTCS(p,m,k1,k2,l,k3,v)

#Criação da figura com o estado final
plot(p,m,"caso_2")

##########################################################################################
print("=================================================================================")
print("Caso 3: Rede aleatória")
print("População começa distribuída aleatoriamente sobre a rede.")
print("Dinheiro começa totalmente focado em um ponto no meio da rede.")
print("=================================================================================")

N_simulations = 500  #Número de simulações desejado.
T = N_simulations * dt  #Tempo máximo de duração da simulação.

#Setando população para t = 0
p1 = np.random.rand(N)
p = p / np.sum(p)

#Setando renda para t = 0
m = np.random.rand(N)

for t in np.arange(0,T,dt):
    p,m = FTCS(p,m,k1,k2,l,k3,v)

#Criação da figura com o estado final
plot(p,m,"caso_3")

Script para gerar os gifs

import numpy as np
import matplotlib.pyplot as plt
import imageio  """Biblioteca p/ gerar os gifs. Para instalar: pip install imageio"""
import os

def FTCS(p,m,k1,k2,l,k3,v):
    N = p.shape[0]
    
    p_ = np.zeros(N)
    m_ = np.zeros(N)
    
    for j in np.arange(0,N):
        back = (j-1)%N
        forward = (j+1)%N
        p_[j] = (1 - 2*k1 - k2*(m[back] - m[j]))*p[j] + k1*(p[back]+p[forward]) - k2*(m[forward] - m[j])*p[forward]
        m_[j] = (1 - l - k3)*m[j] + k3*(m[back] + m[forward]) + v*p[j]
    return p_,m_

def plot_gif(p,m,t,images):
    fig,ax = plt.subplots(1,2,figsize=(10,5))

    #População
    ax[0].set_title("População, t = {}".format(round(t,1)))
    ax[0].plot(100*p,color='purple')
    ax[0].set_ylabel("%")
    ax[0].grid(False)

    #Renda
    ax[1].set_title("Renda, t = {}".format(round(t,1)))
    ax[1].plot(m,color='darkgreen')
    ax[1].fill_between(np.arange(m.shape[0]),m,color='green',alpha=0.5,label="Total money: $ {}".format(round(np.sum(m),2)))
    ax[1].set_ylabel("$")
    ax[1].grid(False)
    ax[1].legend()
    
    plt.tight_layout()
    plt.savefig(imgname) #Salva o frame do sistema
    plt.close()
    
    images.append(imageio.imread(imgname)) #Adiciona a imagem na lista indicada
    os.remove(imgname) #Remove a imagem recém criada
    time.sleep(0.5)


##########################################################################################
"""DECLARAÇÃO DE CONSTANTES"""
N = 100  #Tamanho do eixo x
dx = 1
dt = 0.3
Dm = 1.0  #Constante de difusão p/ renda.
Dp = 1.0  #Constante de difusão p/ população.
gamma = 1.0
alpha = 1.0
beta = 1.0

k1 = Dp*dt/(dx**2)
k2 = gamma*k1/Dp
k3 = Dm*dt/(dx**2)
v = alpha*dt
l = beta*dt
##########################################################################################

N_simulations = 500  #Número de simulações desejado.
T = N_simulations * dt  #Tempo máximo de duração da simulação.

#Setando população para t = 0. Varia conforme a condição inicial desejada.
p = np.zeros(N)
p[20] = 1
    
#Setando renda para t = 0. Varia conforme a condição inicial desejada.
m = np.zeros(N)
m[80] = 1

images = [] #Vetor para alocar as imagens que formarão o gif.
for t in np.arange(0,T,dt):
    if i % 10 == 0:  #Gera 1 imagem a cada 10 iterações, para não gerar um .gif muito pesado.
        plot_gif(p,m,t,images)
    p,m = FTCS(p,m,k1,k2,l,k3,v)

imageio.mimsave("foo.gif", images,fps=3) #Ajustar fps conforme desejar. Quanto maior, + imagens/segundo

Script para gerar os snapshots do sistema

import numpy as np
import matplotlib.pyplot as plt

def FTCS(p,m,k1,k2,l,k3,v):
    N = p.shape[0]
    
    p_ = np.zeros(N)
    m_ = np.zeros(N)
    
    for j in np.arange(0,N):
        back = (j-1)%N
        forward = (j+1)%N
        p_[j] = (1 - 2*k1 - k2*(m[back] - m[j]))*p[j] + k1*(p[back]+p[forward]) - k2*(m[forward] - m[j])*p[forward]
        m_[j] = (1 - l - k3)*m[j] + k3*(m[back] + m[forward]) + v*p[j]
    return p_,m_

def plot_grid(fig,p,m,i,t,ncols,nrows=2):
    #Population
    ax_p = fig.add_subplot(nrows,ncols,i+1)
    ax_p.set_title("t = {}".format(t),fontsize=18)
    ax_p.plot(100*p,color='purple')
    ax_p.fill_between(np.arange(N),100*p,color='purple',alpha=0.5)
    
    #Money
    ax_m = fig.add_subplot(nrows,ncols,i+1+ncols)
    ax_m.set_title("t = {}".format(t),fontsize=18)
    ax_m.plot(m,color='darkgreen')
    ax_m.fill_between(np.arange(N),m,color='green',alpha=0.5,label="Renda líquida: $ {}".format(round(np.sum(m),2)))
    ax_m.legend(fontsize=10)
    
    if i == 0:
        ax_p.set_ylabel("%",fontsize=15)
        ax_m.set_ylabel("$", fontsize=15)

##########################################################################################
"""DECLARAÇÃO DE CONSTANTES"""
N = 100  #Tamanho do eixo x
dx = 1
dt = 0.3
Dm = 1.0  #Constante de difusão p/ renda.
Dp = 1.0  #Constante de difusão p/ população.
gamma = 1.0
alpha = 1.0
beta = 1.0

k1 = Dp*dt/(dx**2)
k2 = gamma*k1/Dp
k3 = Dm*dt/(dx**2)
v = alpha*dt
l = beta*dt
##########################################################################################

N_simulations = 500  #Número de simulações desejado.
T = N_simulations * dt  #Tempo máximo de duração da simulação.

#Setando população para t = 0. Varia conforme a condição inicial desejada.
p = np.zeros(N)
p[20] = 1
    
#Setando renda para t = 0. Varia conforme a condição inicial desejada.
m = np.zeros(N)
m[80] = 1

ncols = 5  #Número de colunas da imagem // Número de snapshots
snapshots = np.array([0*dt,5*dt,10*dt,15*dt,20*dt])  #Vetor que aloca os tempos onde serão tirados os snapshots do sistema

fig, big_axes = plt.subplots(nrows=2,ncols=1,figsize=(20,6),sharey=True)
titles = ["População (%)","Renda ($)"]
for row,big_ax in enumerate(big_axes):
    big_ax.set_title(titles[row],fontsize=22,y=1.15)
    big_ax.axis('off')

count = 0
for t in np.arange(0,T,dt):
    if t in snapshots:  #Gera 1 imagem a cada 10 iterações, para não gerar um .gif muito pesado.
        plot_grid(fig,p,m,count,t,ncols)
        count += 1
    p,m = FTCS(p,m,k1,k2,l,k3,v)

fig.subplots_adjust(hspace=0.55)
plt.savefig("foo.png")

Modelagem bidimensional

class KellerSegelModel():
    def __init__(self, parametros):
        self.parametros = parametros

        self.zeros = np.matrix(np.zeros((parametros.N_x, parametros.N_y)))

        self.tempo = 0

    def contagemPopulacao(self):
        return self.estado_populacao.sum()

    def contagemDinheiro(self):
        return self.estado_dinheiro.sum()

    def setEstadoInicial(self, matriz_populacao, matriz_dinheiro):
        self.estado_populacao = np.matrix(matriz_populacao)
        self.estado_dinheiro = np.matrix(matriz_dinheiro)

    def getEstado(self):
        return (self.estado_populacao, self.estado_dinheiro, self.tempo)

    def atualizaEstado(self):
        # Pega o estado atual da população e dinheiro
        pn = self.estado_populacao
        mn = self.estado_dinheiro

        # Inicializa as variáveis que receberão o estado seguinte
        pn1 = self.zeros.copy()
        mn1 = self.zeros.copy()

        k1 = self.parametros.k1
        k2 = self.parametros.k2
        k3 = self.parametros.k3
        lamb = self.parametros.lamb
        v = self.parametros.v

        # Realiza o FTCS
        for i in np.arange(0, self.parametros.N_x):
            i_previous = (i - 1) % self.parametros.N_x # Garante que em i == 0 o item anterior seja o último da lista
            i_next = (i + 1) % self.parametros.N_x # Garante que em i == N_x o item posterior seja o primeiro da lista

            for j in np.arange(0, self.parametros.N_y):
                j_previous = (j - 1) % self.parametros.N_y # Garante que em j == 0 o item anterior seja o último da lista
                j_next = (j + 1) % self.parametros.N_y # Garante que em j == N_y o item posterior seja o primeiro da lista

                pn1[(i, j)] = pn[(i, j)] * (1 - 4 * k1 - k2 * (mn[(i_previous, j)] - 2 * mn[(i, j)] + mn[(i, j_previous)])) \
                            + k1 * (pn[(i_previous, j)] + pn[(i, j_previous)] + pn[(i_next, j)] + pn[(i, j_next)]) \
                            - k2 * (pn[(i_next, j)] * (mn[(i_next, j)] - mn[(i, j)]) + pn[(i, j_next)] * (mn[(i, j_next)] - mn[(i, j)]))

                mn1[(i, j)] = mn[(i, j)] * (1 - 4 * k3 - lamb) + k3 * (mn[(i_previous, j)] + mn[(i, j_previous)] + mn[(i_next, j)] + mn[(i, j_next)]) + v * pn[(i, j)]

        self.estado_populacao = pn1 # Atualiza o estado da população
        self.estado_dinheiro = mn1 # Atualiza o estado do dinheiro
        self.tempo += self.parametros.dt # Atualiza o tempo decorrido no modelo

    def atualizaEstadoMultiplasVezes(self, n = 1):
        for _ in range(0, n):
            self.atualizaEstado()

Referências

↑ ^1,0 ^1,1 Evelyn F. Keller, Lee A. Segel, Initiation of slime mold aggregation viewed as an instability, Journal of Theoretical Biology, Volume 26, Issue 3, 1970, Pages 399-415
↑ ^2,0 ^2,1 ^2,2 ^2,3 ^2,4 Sayama, H. Introduction to the Modeling and Analysis of Complex Systems. 2015
↑ Hoffmann, F. C. O. Keller-Segel-Type Models and Kinetic Equations for Interacting Particles: Long-Time Asymptotic Analysis , Tese de doutorado. 2017 [1]
↑ Scherer, C. Métodos Computacionais da Física. 2010

[keller-segel-original-1] 1,0 ^1,1 Evelyn F. Keller, Lee A. Segel, Initiation of slime mold aggregation viewed as an instability, Journal of Theoretical Biology, Volume 26, Issue 3, 1970, Pages 399-415

[sayama-2] 2,0 ^2,1 ^2,2 ^2,3 ^2,4 Sayama, H. Introduction to the Modeling and Analysis of Complex Systems. 2015

[thesis-3] Hoffmann, F. C. O. Keller-Segel-Type Models and Kinetic Equations for Interacting Particles: Long-Time Asymptotic Analysis , Tese de doutorado. 2017 [1]

[scherer-4] Scherer, C. Métodos Computacionais da Física. 2010

[1]

[2]

[3]

[4]

Modelo de Keller-Segel para relação população-economia

Índice