Uma abordagem geométrica da estabilidade local em dinâmicas discretas com aplicação em um modelo de cooperação

A geometric approach to local stability in discrete dynamics with a cooperation model application

Resumo

Na busca pela compreensão teórica de fenômenos naturais, pesquisadores muitas vezes recorrem a modelos matemáticos. Os modelos baseados em indivíduos (IBMs) geralmente são descritos por regras probabilísticas, dificultando a análise analítica de sua estabilidade. Entretanto, é possível escrever uma Equação Mestra que descreve a dinâmica da média das variáveis de interesse, viabilizando tal análise. Neste trabalho revisitamos a teoria de sistemas dinâmicos e apresentamos por uma argumentação geométrica como calcular a estabilidade local de pontos fixos de um sistema determinístico de equações de diferença. Em seguida, desenvolvemos um IBM estocástico inspirado no fenômeno de cooperação entre machos de aves que formam grupos cooperativos para atrair fêmeas. Mostramos como construir as Equações Mestras que descrevem de forma exata a média do tamanho populacional de cooperadores e não cooperadores. Analisamos então a estabilidade local de tais equações e discutimos os possíveis mecanismos de manutenção da cooperação. Assim, esse trabalho auxilia estudantes de diferentes formações nos seus primeiros contatos com a teoria de sistemas dinâmicos.

Palavras-chave:
Modelagem; sistemas dinâmicos; cooperação

Abstract

In the search for understanding natural phenomena with complex dynamics, researchers often resort to mathematical tools. Individual-based models (IBMs) are often described by probabilistic rules, which hampers the analytical description of their stability. However, it is possible to write a Master Equation that describes the average behavior of the variables of interest, enabling such analysis. In this paper, we revisit Dynamical Systems Theory through geometric arguments in order to evaluate the local stability of fixed points in discrete-time equations. After that, we present an IBM, inspired by the cooperation phenomenon observed in some male neotropical birds that form groups to attract females. From this example, we show how to construct the Master Equations that accurately describe the mean population size of cooperators and non-cooperators in the system. Finally, we analyze the local stability of these equations and discuss possible maintenance mechanisms of cooperation. Thus, this work assists students from different backgrounds in their first encounters with Dynamical Systems Theory.

Keywords:
Modeling; dynamical systems; cooperation

1. Introdução

Com frequência sistemas biológicos apresentam dinâmicas complexas. Interações ecológicas do tipo predador-presa [¹], a ascensão e declínio de epidemias [²], o piscar sincronizado de vaga-lumes [³] e o surgimento de novas espécies [⁴] são alguns exemplos desses fenômenos que fascinam os cientistas há muito tempo. Nesse contexto, a matemática aplicada à formulação de modelos fornece um solo fértil para a exploração de hipóteses teóricas.

Dentre estes modelos, estão aqueles descritos por equações determinísticas, sejam diferenciais ou de diferença. Nos modelos de equações diferenciais, como as equações de Lotka-Volterra (que descrevem interações do tipo predador-presa) [¹], o tempo é uma grandeza contínua, enquanto nos modelos de equações de diferença, como o mapa logístico (que, em ecologia, descreve a variação do tamanho de uma população) [¹], o tempo é uma grandeza discreta [⁵]. A escolha entre um ou outro depende das características do fenômeno e dos objetivos da pesquisa. Um aspecto importante de modelos determinísticos é a existência de técnicas analíticas para se estudar a estabilidade dos pontos fixos do sistema, a partir das quais se pode investigar suas propriedades dinâmicas [⁵].

Outra abordagem utilizada no estudo de processos biológicos é a Modelagem Baseada em Indivíduo, ou IBM (do inglês “Individual Based Modeling”), também chamada, em um contexto mais geral, Modelagem baseada em Agentes. Nesse caso, o sistema é descrito por meio de um algoritmo, isto é, uma sequência de regras, que dita o comportamento dos indivíduos que o compõem [⁶, ⁷]. Muitas vezes essas regras são probabilísticas, trazendo estocasticidade para os modelos. Esse algoritmo pode então ser replicado em simulações computacionais. Entretanto, o estudo das propriedades dinâmicas de um IBM a partir de simulações pode exigir um alto custo computacional, uma vez que geralmente são necessárias muitas réplicas para uma análise estatística apropriada. Uma alternativa é a descrição da dinâmica a partir de um conjunto de equações, denominadas Equações Mestras, que descrevem a evolução temporal da probabilidade do sistema se encontrar em um dado estado, ou ainda, a evolução da média das variáveis de interesse [⁸]. Uma das vantagens de se escrever as Equações Mestras é a possibilidade de retomar a análise de estabilidade de forma determinística. Há ainda a possibilidade de se interpretar um IBM utilizando teoria de jogos, escrevendo as matrizes de payoff, as quais descrevem quanta “vantagem” um indivíduo possui em relação a outro [⁹].

A idealização de um IBM está muito presente na biologia evolutiva, na qual os mecanismos de surgimento e fixação de caracteres observáveis, i.e. do fenótipo, em uma população são de bastante relevância [¹⁰]. Um caso interessante é o de comportamentos cooperativos, no qual indivíduos de uma comunidade promovem o aumento do sucesso reprodutivo – fitness, em inglês – de outros indivíduos, em detrimento do próprio fitness. Observada em diversas espécies, como de aves e peixes, a cooperação vem sendo estudada pela biologia e pela matemática de diferentes formas, mas ainda não é dada como uma questão resolvida [¹¹]. Um exemplo ainda não explorado pela modelagem matemática é a curiosa dinâmica reprodutiva das espécies de aves do gênero Chiroxiphia (e.g. C. caudata e C. linearis). Nela, os machos se reúnem em grupos para realizar exibições, compostas por uma série de movimentos e vocalizações coordenados, na expectativa de convencer uma fêmea da espécie a acasalar, não com qualquer indivíduo do grupo, mas sim apenas com um macho dominante previamente determinado¹[¹², ¹³]. Superficialmente, os machos subalternos parecem beneficiar o dominante ao mesmo tempo em que prejudicam a si mesmos, uma vez que a cooperação exige deles recursos importantes para sua sobrevivência [⁹, ¹¹, ¹², ¹³].

Neste trabalho, iniciamos expondo uma parte da teoria de sistemas dinâmicos para equações de diferença. Mais especificamente, apresentamos argumentos geométricos para a análise da estabilidade local em um sistema de uma e duas dimensões (subseções 2.1 e 2.2 respectivamente). Em seguida, como um exemplo de aplicação dos conceitos abordados, apresentamos um IBM que descreve, de forma simplificada, uma dinâmica de cooperação inspirada no comportamento das aves do gênero Chiroxiphia (subseção 3.1). Por fim, obtemos as Equações Mestras para a dinâmica média do sistema e avaliamos a estabilidade local dos diferentes pontos fixos (subseção 3.3), permitindo discutir, com base nos resultados analíticos, a manutenção dos comportamentos cooperativos nesse sistema (seção 4).

Uma vez que abordamos conceitos da teoria de sistemas dinâmicos, do formalismo teórico à aplicação, esse trabalho constitui-se em uma introdução bastante completa à modelagem matemática de fenômenos biológicos, sendo relevante para estudantes de física e matemática interessados em biologia quantitativa, ou ainda para estudantes de biologia com interesses em descrições analíticas. Ainda, professores de disciplinas como “Introdução a sistemas dinâmicos” e “Matemática para a biologia” podem se beneficiar do seu conteúdo.

2. Avaliando a Estabilidade Local em Modelos de Tempo Discreto

2.1. Caso unidimensional

Partimos de um sistema unidimensional representado por uma equação de tempo discreto com a forma geral

(1) $$x_{t+1}=F\left(x_t\right),$$

onde $t\in {\mathbb{Z}}_{+}$ (subconjunto dos inteiros positivos e o $0$) é o índice temporal, e $x_t\in ℝ$ representa o estado do sistema no respectivo instante. Iterando a (equação 1) ao longo do tempo, os valores de $x_t$ são atualizados, podendo atingir, eventualmente, um ponto fixo $x^{\ast }$. Um ponto fixo $x^{\ast }$ é um possível valor de $x_t$ tal que, se $x_t$ de fato assumi-lo, $x_{t+n}$ será igual a $x_t$ para qualquer $n\ge 0$:

(2) $$x_t=x^{\ast }\to x_{t+n}=x_t.$$

Dado um ponto fixo $x^{\ast }$, podemos determinar sua estabilidade local, isto é, calcular como o sistema se comporta caso esteja nas proximidades de $x^{\ast }$. No que segue, vamos apresentar com argumentos geométricos como isto pode ser feito.

Consideremos uma condição inicial $x_0$ próxima a um ponto fixo $x^{\ast }$, isto é,

(3) $$x_0=x^{\ast }+ϵ,$$

na qual $ϵ\in ℝ$. Pela (equação 1), calculamos o valor de $x_t$ no próximo instante de tempo

(4) $$x_1=F\left(x_0\right)=F\left(x^{\ast }+ϵ\right).$$

Qualquer que seja o valor de $x_1$, podemos expressá-lo como

(5) $$x_1=x^{\ast }+\mu .$$

onde $\mu \in ℝ$.

Se $ϵ=0$, temos $x_1=x_0=x^{\ast }$, já que $x^{\ast }$ é um ponto fixo. Para $ϵ\ne 0$, $x_1$ é, tal como representado na (equação 5), igual a $x^{\ast }$ mais uma quantidade $\mu$. O valor absoluto de $ϵ$, $\left|ϵ\right|$, pode ser menor, igual ou maior que $\left|\mu \right|$. Se $\left|ϵ\right|<\left|\mu \right|$, $x_0$ é mais próximo de $x^{\ast }$ do que $x_1$. Se $\left|ϵ\right|>\left|\mu \right|$, $x_1$ é mais próximo de $x^{\ast }$ do que $x_0$. Por fim, se $\left|ϵ\right|=\left|\mu \right|$ e $ϵ=\mu$, então $x_0=x_1$, de modo que $x_0$ também é ponto fixo do sistema; porém, se $ϵ=-\mu$, temos que $x_0$ e $x_1$ estão à mesma distância de $x^{\ast }$, mas $x_0\ne x_1$.

Consideremos a seguinte relação

(6) $$\mu =Jϵ,$$

onde $J\in ℝ$. Podemos então avaliar a evolução do sistema através de $J$, ao invés de $ϵ$ e $\mu$ diretamente, ou de $x_0$ e $x_1$. Isso se provará útil assim que encontrarmos uma forma alternativa de determinar $J$.

Na Figura 1 são representados dois pontos, $\left(x^{\ast },F\left(x^{\ast }\right)=x^{\ast }\right)$ e $\left(x_0,F\left(x_0\right)=x_1\right)$. A distância no eixo horizontal entre esses pontos é dada por $ϵ$ e, no eixo vertical, dada por $\mu$, formando o triângulo traçado na figura, com

(7) $$\tan \theta =\frac{F\left(x_0\right)-F\left(x^{\ast }\right)}{x_0-x^{\ast }}=\frac{F\left(x^{\ast }+ϵ\right)-F\left(x^{\ast }\right)}{ϵ}=\frac{\mu }{ϵ}=J.$$

Figura 1
Representação gráfica genérica da função $F\left(x_t\right)$. Aqui, representamos os pontos $(x^{\ast },x^{\ast })$ e $(x_0,x_1)$, a distância $ϵ$ ((equação 3)) e $\mu$ ((equação 5)). Observe que, apesar dos dois eixos, há apenas uma dimensão sendo representada no gráfico, pois todo ponto possui a forma $(x_t,x_{t+1})$.

Se aplicarmos o limite $ϵ\to 0$ na (equação 7), obtemos

(8) $$J=\underset{ϵ\to 0}{lim}\frac{F\left(x^{\ast }+ϵ\right)-F\left(x^{\ast }\right)}{ϵ}={\frac{dF\left(x_t\right)}{d{x}_t}|}_{x^{\ast }},$$

então

(9) $$\left|J\right|=\left|\left(\frac{dF\left(x_t\right)}{d{x}_t}{|}_{x^{\ast }}\right)\right|=\frac{\left|\mu \right|}{\left|ϵ\right|}.$$

Assim, para um sistema dinâmico de tempo discreto dado pela (equação 1), e que possui um ponto fixo $x^{\ast }$ (definido pela (equação 2)), podemos estudar como o sistema evolui nas proximidades de $x^{\ast }$ através de $J$, que se iguala à derivada de $F\left(x_t\right)$ em relação a $x_t$ no ponto $x^{\ast }$, dado um $ϵ$ suficientemente pequeno. Se $\left|J\right|>1$, classificamos $x^{\ast }$ como instável, ou repulsor, porque $x_t$ afasta-se desse ponto. Se $\left|J\right|<1$, classificamos $x^{\ast }$ como assintoticamente estável, ou atrator, porque $x_t$ aproxima-se desse ponto. No apêndice A.1, discutimos a dinâmica temporal em torno dos atratores e repulsores.

2.2. Caso bidimensional

Podemos generalizar esse método para $n\ge 1$ dimensões. Começaremos analisando o caso bidimensional ($n=2$), para então generalizar as equações obtidas para um $n$ qualquer.

Com $n=2$ temos duas variáveis independentes (duas dimensões), que chamamos de $x_t$ e $y_t$. A (equação 1) é substituída por uma equação vetorial

(10) $${\overrightarrow{x}}_{t+1}=\overrightarrow{F}\left({\overrightarrow{x}}_t\right),$$

na qual

(11) $$\overrightarrow{F}\left({\overrightarrow{x}}_t\right)=\left[\begin{array}{c}f(x_t,y_t)\\ g(x_t,y_t)\end{array}\right]=\left[\begin{array}{c}{x}_{t+1}\\ y_{t+1}\end{array}\right]={\overrightarrow{x}}_{t+1}.$$

Note que a função anteriormente dada por $F\left(x_t\right)$ se torna também um vetor dado por duas componentes, $f(x_t,y_t)$ e $g(x_t,y_t)$. Ainda, um ponto fixo ${\overrightarrow{x}}^{\ast }$ dessa dinâmica é definido como

(12) $${\overrightarrow{x}}_t={\overrightarrow{x}}^{\ast }\to {\overrightarrow{x}}_{t+n}={\overrightarrow{x}}_t.$$

De forma análoga ao caso unidimensional, seja uma condição inicial ${\overrightarrow{x}}_0$ próxima de um ponto fixo ${\overrightarrow{x}}^{\ast }$, temos

(13) $${\overrightarrow{x}}_0={\overrightarrow{x}}^{\ast }+\overrightarrow{ϵ},$$

(14) $${\overrightarrow{x}}_1=\overrightarrow{F}\left({\overrightarrow{x}}_0\right)=\overrightarrow{F}\left({\overrightarrow{x}}^{\ast }+\overrightarrow{ϵ}\right)={\overrightarrow{x}}^{\ast }+\overrightarrow{\mu },$$

nas quais introduzimos os vetores $\overrightarrow{ϵ}$ e $\overrightarrow{\mu }$, e definimos a matriz J tal que

(15) $$\overrightarrow{\mu }=\mathbf{\text{J}}\overrightarrow{ϵ}.$$

É importante destacar a diferença entre a (equação 6) e a (equação 15). Enquanto na primeira, $J$ é um escalar que relaciona os escalares $ϵ$ e $\mu$, na segunda J é uma matriz que relaciona os vetores $\overrightarrow{ϵ}={({ϵ}_x,{ϵ}_y)}^T$ e $\overrightarrow{\mu }={({\mu }_x,{\mu }_y)}^T$ (que podem ter módulo, direção e sentido diferentes). Vamos agora calcular essa matriz J, e mostrar como utilizá-la para avaliar a estabilidade local de um ponto fixo do sistema.

Em notação matricial, da definição de multiplicação de matrizes e da (equação 15), concluímos que J deve ser uma matriz $2\times 2$, pois a perturbação $\overrightarrow{ϵ}$ é escrita como uma matriz $2\times 1$ que, ao ser multiplicada por J, é igual a $\overrightarrow{\mu }$, também uma matriz $2\times 1$. Escrevendo

(16) $$\mathbf{\text{J}}=\left[\begin{array}{ccc}{J}_{11}& & J_{12}\\ J_{21}& & J_{22}\end{array}\right],$$

obtemos

(17) $$\mathbf{\text{J}}\overrightarrow{ϵ}=\left[\begin{array}{c}{J}_{11}{ϵ}_x+J_{12}{ϵ}_y\\ J_{21}{ϵ}_x+J_{22}{ϵ}_y\end{array}\right]=\left[\begin{array}{c}{\mu }_x\\ {\mu }_y\end{array}\right]=\overrightarrow{\mu }.$$

Portanto, os elementos de J devem satisfazer o sistema de equações

(18) $$\begin{array}{c}{J}_{11}{ϵ}_x+J_{12}{ϵ}_y={\mu }_x\\ J_{21}{ϵ}_x+J_{22}{ϵ}_y={\mu }_y\end{array}.$$

A Figura 2 é análoga à Figura 1, mas para o caso bidimensional, no qual diferentes representações da dinâmica podem ser feitas. Na Figura 2a temos um plano formado pelos eixos $x_t$ e $y_t$, também conhecido como espaço de fase. Cada ponto nesse plano corresponde a um vetor ${\overrightarrow{x}}_t$. Na Figura 2b temos um gráfico de três dimensões: duas delas formam o plano descrito anteriormente, enquanto a terceira representa a função $f(x_t,y_t)$. Note que não representamos a função $g(x_t,y_t)$ presente na (equação 11), já que essa seria descrita como uma quarta dimensão na figura.

Figura 2
Diferentes representações gráficas associadas à função $\overrightarrow{F}\left({\overrightarrow{x}}_t\right)$. O item a) é um plano formado pelos eixos $x_t$ e $y_t$; em b) há um gráfico com três dimensões, duas delas formam o plano descrito anteriormente, enquanto a terceira representa a função $f(x_t,y_t)$; em c) um gráfico de $f(x_t,y_t)$ com um eixo $S$ paralelo a $\overrightarrow{ϵ}$; em d) um gráfico de $g(x_t,y_t)$ com o mesmo eixo $S$.

No plano $x_t\times y_t$ (Figura 2a) temos três pontos: ${\overrightarrow{x}}^{\ast }=(x^{\ast },y^{\ast })$, ${\overrightarrow{x}}_0=(x_0,y_0)$ e ${\overrightarrow{x}}_1=(x_1,y_1)$. Estão representados também suas componentes e os vetores $\overrightarrow{ϵ}$ e $\overrightarrow{\mu }$. Na Figura 2b temos outros dois pontos: $(x^{\ast },y^{\ast },f\left({\overrightarrow{x}}^{\ast }\right))$ e $(x_0,y_0,f\left({\overrightarrow{x}}_0\right))$. Com esses últimos, podemos traçar um triângulo retângulo, sendo o segmento de reta que os liga a hipotenusa, $ϵ=\left|\overrightarrow{ϵ}\right|$ um cateto e, como $f\left({\overrightarrow{x}}_0\right)-f\left({\overrightarrow{x}}^{\ast }\right)=x_1-x^{\ast }={\mu }_x$, ${\mu }_x$ corresponde ao outro cateto. Além disso,

(19) $$\tan \theta =\frac{{\mu }_x}{ϵ}=\frac{f\left({\overrightarrow{x}}_0\right)-f\left({\overrightarrow{x}}^{\ast }\right)}{ϵ}=\frac{f\left({\overrightarrow{x}}^{\ast }+\overrightarrow{ϵ}\right)-f\left({\overrightarrow{x}}^{\ast }\right)}{ϵ},$$

sendo $\theta$ o ângulo representado na Figura 2b. Note que $\overrightarrow{ϵ}$ foi escolhido, sem perda de generalidade, de modo a cair sobre a reta formada entre a origem do gráfico e o ponto fixo $(x^{\ast },y^{\ast })$, que definimos como um novo eixo $S$. Na Figura 2c, observamos a dinâmica sobre o plano formado pelo eixo $f\left({\overrightarrow{x}}_t\right)$ e pelo eixo $S$. Dessa forma, podemos lidar com a função $f(x_t,y_t)$ como sendo igual a $f\left(S\right)$, com as variáveis $x_t$ e $y_t$ correspondendo a pontos sobre esse novo eixo. Por exemplo, o vetor ${\overrightarrow{x}}^{\ast }=(x^{\ast },y^{\ast })$ torna-se o escalar $S^{\ast }$. Também podemos deduzir, da Figura 2b,

(20) $$\begin{array}{c}x\left(S\right)=S\mathrm{}\cos \mathrm{}\alpha \\ y\left(S\right)=S\mathrm{}\sin \mathrm{}\alpha \end{array}.$$

Logo, a (equação 19) torna-se

(21) $$\tan \theta =\frac{{\mu }_x}{ϵ}=\frac{f\left(S_0\right)-f\left(S^{\ast }\right)}{ϵ}=\frac{f\left(S^{\ast }+ϵ\right)-f\left(S^{\ast }\right)}{ϵ},$$

e no limite onde $ϵ\to 0$ temos

(22) $$\tan \theta =\underset{ϵ\to 0}{lim}\frac{f\left(S^{\ast }+ϵ\right)-f\left(S^{\ast }\right)}{ϵ}={\frac{df\left(S\right)}{dS}|}_{S^{\ast }}.$$

Uma vez que $f$ é função de $x_t$ e $y_t$ que, por sua vez, são agora funções de $S$, a derivada é dada por

(23) $${\frac{df\left(S\right)}{dS}|}_{S^{\ast }}={\frac{\partial f}{\partial x_t}|}_{{\overrightarrow{x}}^{\ast }}\frac{d{x}_t}{dS}+{\frac{\partial f}{\partial y_t}|}_{{\overrightarrow{x}}^{\ast }}\frac{d{y}_t}{dS}.$$

Derivando as (equações 20) em relação a $S$, e substituindo o resultado na (equação 23), obtemos

(24) $${\frac{df\left(S\right)}{dS}|}_{S^{\ast }}=\frac{{\mu }_x}{ϵ}={\frac{\partial f}{\partial x_t}|}_{{\overrightarrow{x}}^{\ast }}\cos \alpha +{\frac{\partial f}{\partial y_t}|}_{{\overrightarrow{x}}^{\ast }}\sin \alpha ,$$

e, ainda,

(25) $$\begin{array}{l}{\mu }_x={\frac{\partial f}{\partial x_t}|}_{{\overrightarrow{x}}^{\ast }}ϵ\cos \alpha +{\frac{\partial f}{\partial y_t}|}_{{\overrightarrow{x}}^{\ast }}ϵ\sin \alpha \\ ={\frac{\partial f}{\partial x_t}|}_{{\overrightarrow{x}}^{\ast }}{ϵ}_x+{\frac{\partial f}{\partial y_t}|}_{{\overrightarrow{x}}^{\ast }}{ϵ}_y.\end{array}$$

Comparando as (equações 18) e (25), percebemos que $J_{11}={\frac{\partial f}{\partial x_t}|}_{{\overrightarrow{x}}^{\ast }}$, e $J_{12}={\frac{\partial f}{\partial y_t}|}_{{\overrightarrow{x}}^{\ast }}$. Aplicando o mesmo raciocínio para $g(x_t,y_t)$, o que pode ser feito com o auxílio da Figura 2d, encontramos

(26) $${\mu }_y={\frac{\partial g}{\partial x_t}|}_{{\overrightarrow{x}}^{\ast }}{ϵ}_x+{\frac{\partial g}{\partial y_t}|}_{{\overrightarrow{x}}^{\ast }}{ϵ}_y.$$

Enfim, para um sistema dinâmico bidimensional de tempo discreto, mostramos que

(27) $$\mathbf{\text{J}}=\left[\begin{array}{ccc}{\frac{\partial f}{\partial x_t}|}_{{\overrightarrow{x}}^{\ast }}& & {\frac{\partial f}{\partial y_t}|}_{{\overrightarrow{x}}^{\ast }}\\ {\frac{\partial g}{\partial x_t}|}_{{\overrightarrow{x}}^{\ast }}& & {\frac{\partial g}{\partial y_t}|}_{{\overrightarrow{x}}^{\ast }}\end{array}\right].$$

A matriz cujos elementos são derivadas parciais de primeira ordem das componentes de uma função $\overrightarrow{F}\left(\overrightarrow{x}\right)$ é chamada de Matriz Jacobiana de $\overrightarrow{F}$. Aqui temos a Jacobiana de $\overrightarrow{F}$ aplicada no ponto fixo ${\overrightarrow{x}}^{\ast }$, que descreve a transformação do vetor $\overrightarrow{ϵ}$ em $\overrightarrow{\mu }$.

Agora, devemos calcular os autovalores e autovetores de J, isto é, encontrar escalares ${\lambda }_i$ (autovalores) e vetores não nulos ${\overrightarrow{v}}_i$ (autovetores) que satisfaçam a equação

(28) $$\mathbf{\text{J}}{\overrightarrow{v}}_i={\lambda }_i{\overrightarrow{v}}_i.$$

Geometricamente, se for possível satisfazer a (equação 28), então existem vetores não nulos ${\overrightarrow{v}}_i$ que ao serem multiplicados pela matriz J, resultam em vetores paralelos a ${\overrightarrow{v}}_i$, e diferentes apenas por uma constante de proporcionalidade ${\lambda }_i$. Uma propriedade importante dos autovetores é que, se ${\overrightarrow{v}}_i$ é um autovetor de uma matriz J associado a um autovalor ${\lambda }_i$, e $r\in {ℝ}^{\ast }$ (subconjunto dos reais sem o $0$), então $r{\overrightarrow{v}}_i$ é também um autovetor de J associado ao mesmo autovalor ${\lambda }_i$.

Se for possível calcular um autovetor ${\overrightarrow{v}}_i$ de J, o autovalor ${\lambda }_i$ associado indicará a relação entre $\overrightarrow{ϵ}$ e $\overrightarrow{\mu }$ na direção de ${\overrightarrow{v}}_i$, e portanto, como o sistema se comportará (nessa direção) nas proximidades do ponto fixo ${\overrightarrow{x}}^{\ast }$. Assim, ao encontrar todos os autovetores de J, podemos compreender a estabilidade do sistema de equações ao redor de seu ponto fixo.

Se ${\overrightarrow{v}}_i$ e $\overrightarrow{ϵ}$ forem paralelos, sabemos que existe um autovetor ${\overrightarrow{v}}_i^{\prime }$ que não possui apenas direção, mas também módulo e sentido iguais a $\overrightarrow{ϵ}$. Os módulos de ${\overrightarrow{v}}_i$ e ${\overrightarrow{v}}_i^{\prime }$ podem ser iguais ou não, mas pela mesma propriedade ambos possuem o mesmo autovalor ${\lambda }_i$. Sendo assim, supondo que ${\overrightarrow{v}}_i$, $\overrightarrow{ϵ}$ e $\overrightarrow{\mu }$ são paralelos, pelas (equações 15) e (28)

(29) $$\overrightarrow{\mu }=\mathbf{\text{J}}\overrightarrow{ϵ}=\mathbf{\text{J}}{\overrightarrow{v}}_i^{\prime }={\lambda }_1{\overrightarrow{v}}_i^{\prime },$$

(30) $$\overrightarrow{\mu }={\lambda }_i\overrightarrow{ϵ}.$$

Se a suposição feita anteriormente for inválida, será necessário substituir $\overrightarrow{\mu }$ e/ou $\overrightarrow{ϵ}$ por suas componentes paralelas a ${\overrightarrow{v}}_i$ para obter um resultado semelhante

(31) $${\overrightarrow{\mu }}_i={\lambda }_i{\overrightarrow{ϵ}}_i.$$

Em geral, porém, pouco importam as direções de $\overrightarrow{ϵ}$ e $\overrightarrow{\mu }$, pois eles são necessários apenas para compreender a função dos autovalores e autovetores de J. De forma análoga ao caso unidimensional, encontramos que, se $\left|{\lambda }_i\right|>1$, então $\left|{\overrightarrow{\mu }}_i\right|>\left|{\overrightarrow{ϵ}}_i\right|$, enquanto se $\left|{\lambda }_i\right|<1$, $\left|{\overrightarrow{\mu }}_i\right|<\left|{\overrightarrow{ϵ}}_i\right|$.

A título de exemplo, podemos resgatar a Figura 2a e supor que do ponto fixo representado ($x^{\ast }$) obtemos infinitos autovetores em duas direções diferentes, como representado na Figura 3. Cada conjunto de autovetores em uma determinada direção possui um mesmo autovalor associado, e nessa figura (por exemplo) podemos deduzir que o módulo desses autovalores são maiores que $1$, porque em cada direção a trajetória está se afastando do ponto fixo. De outra forma podemos afirmar que ${\mu }_1>{ϵ}_1$, donde concluímos (pela equação 31 e discussão anterior) que o autovalor ${\lambda }_1$ é maior que $1$. O mesmo vale para ${\lambda }_2$. Aqui nos baseamos em suposições e no comportamento dinâmico do sistema representado na imagem para ilustrar o papel dos autovetores e autovalores, porém ao aplicar esse conhecimento no estudo de um sistema particular fazemos o inverso, ou seja, calculamos a Matriz Jacobiana em um ponto fixo, e em seguida os autovalores desta, a partir dos quais deduzimos o comportamento dinâmico nas proximidades do ponto.

Figura 3
Ponto fixo e autovetores. Baseando-se na Figura 2a, as setas pretas indicam os autovetores no espaço de fase, e a interseção entre elas, o ponto fixo estudado. Desse cenário podemos inversamente deduzir que os dois autovalores associados aos autovetores representados são maiores que $1$, dado que ${\mu }_1>{ϵ}_1$ e ${\mu }_2>{ϵ}_2$.

Assim, para cada ponto fixo em um sistema bidimensional é necessário encontrar dois autovalores (que, nesse trabalho, assumimos que são números reais) associados a direções diferentes, a partir dos quais classificamos o ponto fixo com relação à sua estabilidade. Temos então que, para $\left|{\lambda }_{\mathrm{1,2}}\right|>1$, o ponto fixo é instável, ou repulsor. Para $\left|{\lambda }_{\mathrm{1,2}}\right|<1$, o ponto fixo é assintoticamente estável, ou atrator. Por fim, para $\left|{\lambda }_1\right|<1$ e $\left|{\lambda }_2\right|>1$, ou $\left|{\lambda }_1\right|>1$ e $\left|{\lambda }_2\right|<1$, o ponto fixo é um ponto de sela hiperbólico, o que também faz dele um ponto instável. No apêndice A.2 ilustramos essas possíveis dinâmicas, e no apêndice A.3 comentamos sobre a estabilidade em sistemas com $n$ dimensões.

3. Aplicação a Sistemas Cooperativos

3.1. Apresentação

Para ilustrar a aplicação dos conceitos da seção anterior, escolhemos apresentar um modelo de um sistema biológico real envolvendo comportamentos cooperativos. O estudo da cooperação [⁹] é interessante não só para a biologia evolutiva [¹⁴], como também para a economia [¹⁵] e antropologia [¹⁶]. Em geral, modelos envolvendo cooperação, como o dilema do prisioneiro [¹⁷] e o modelo gavião-pombo [¹⁵], demonstram que ela não pode ser estável na ausência de mecanismos que a favoreçam [⁹].

Um exemplo interessante de comportamento cooperativo é o caso das aves do gênero Chiroxiphia (e.g. C. caudata e C. linearis), devido aos seus sistemas de acasalamento [¹², ¹³]. Neles, os machos se agrupam em coalizões cooperativas para realizar o cortejo das fêmeas através de exibições motoras e sonoras [¹², ¹³], em outras palavras, eles fazem sua “dança do acasalamento” em conjunto, coordenando seus movimentos uns com os outros. Apesar de comportamentos de exibição serem recorrentes em diferentes grupos animais [¹⁸], o que os diferencia é o fato desse cortejo ser feito em grupo – fator decisivo na escolha das fêmeas. Machos que realizam a exibição sozinhos têm uma taxa de sucesso de cópula drasticamente menor [¹², ¹³, ¹⁹]. O indivíduo dominante do grupo é geralmente o único que copula com a fêmea, após uma “finalização solo” da performance, executada apenas por ele [¹², ¹³]. Como esse comportamento aparentemente beneficia apenas um indivíduo do grupo, ainda que todos paguem o custo energético e de exposição a potenciais predadores, não é claro por que ele se mantém, do ponto de vista evolutivo. A forma como pressões de seleção sexual interagem com dinâmicas cooperativas é pouco explorada do ponto de vista teórico [¹¹], reforçando a particularidade desse sistema específico.

Para modelar essa dinâmica cooperativa, o primeiro passo é construir uma situação simplificada, mas que mantenha a capacidade de replicar os fenômenos de interesse. Vamos imaginar uma espécie hipotética na qual alguns indivíduos já realizam um comportamento cooperativo, associando-se a grupos e auxiliando o dominante no cortejo em detrimento das suas chances individuais de se reproduzir. Aqui, definimos grupo como um conjunto de um ou mais indivíduos associados por interações de cooperação. Ainda, definimos o macho dominante como aquele que realiza a cópula.

Em nossa situação hipotética, vamos considerar que: (i) a capacidade de cooperar é um fenótipo discreto e estritamente hereditário: um indivíduo ou coopera por toda a sua vida, ou nunca coopera ativamente, e isso é determinado pelo comportamento de seu progenitor; (ii) a chance da fêmea escolher qualquer grupo é diretamente proporcional ao tamanho do grupo, conferindo assim uma potencial vantagem ao comportamento de cooperar/formar grupos, com essa vantagem impactando primariamente o dominante do grupo; (iii) um tamanho populacional constante, que pode ser interpretado como a capacidade de suporte do ambiente; e (iv) no caso da morte do dominante de um grupo com mais de um indivíduo, outro integrante qualquer do grupo tomará a posição.

3.2. Modelo baseado em indivíduos

Com o contexto apresentado, a modelagem baseada em indivíduos consiste em modelar explicitamente cada indivíduo e associar regras de iteração para a evolução do sistema. Iniciamos considerando uma população de $N$ ($\ge 2$) machos, na qual certos indivíduos não cooperam, assumindo o fenótipo (0), e outros cooperam, assumindo o fenótipo (1). A sequência completa de todas as regras de iteração define um passo de tempo $t$ (Figura 4a), sendo estas regras detalhadas a seguir.

Figura 4
Dinâmica temporal do modelo de Cooperação. a) Diagrama dos eventos que acontecem a cada passo de tempo em uma dinâmica com regras baseadas em indivíduos (IBM); b) Exemplificação de uma possibilidade de mudança de estado em um sistema hipotético com 10 indivíduos, com as probabilidades de cada evento particular acontecer. Os círculos preenchidos representam indivíduos sem o fenótipo de cooperar (0), enquanto os círculos sem preenchimento representam indivíduos com o fenótipo de cooperar (1). No caso apresentado, a fêmea escolheu com probabilidade $P^s$ um grupo de Cooperadores, gerando assim um novo indivíduo também com o fenótipo de cooperar. No tempo $t+1$ esse novo indivíduo escolheu, com probabilidade $P^c$, associar-se a um Desertor, caracterizando-o como um Traidor. Uma situação especial resultante das regras de mortalidade também está representada: na eventual morte de um dominante Desertor, o novo dominante do grupo será um indivíduo com fenótipo (1), de modo que a morte de um Desertor implica no recrutamento de todos os Traidores associados para a população de Cooperadores. A probabilidade de um indivíduo morrer é representada por $P^m$.

As fêmeas não são modeladas explicitamente, mas assumimos que a cada iteração do modelo um grupo de $n$ machos tem probabilidade

(32) $$P^s\left(n\right)=\frac{n}{N}$$

de ser escolhido por uma fêmea, de forma que um novo indivíduo com o mesmo fenótipo do macho dominante desse grupo será gerado. Se o fenótipo do descendente for o de cooperar (1), este novo macho irá se associar a um grupo já existente, com probabilidade proporcional ao tamanho $n$ desse grupo (independente dos fenótipos do grupo),

(33) $$P^c\left(n\right)=\frac{n}{N}.$$

Caso tenha o fenótipo de não cooperar (0), ele irá fundar um novo grupo, inicialmente composto apenas por este indivíduo. Dessa forma, enfatizamos que: não é possível haver mais de um macho com fenótipo (0) em um mesmo grupo; esses machos sempre serão os dominantes em seus respectivos grupos; e esta é a única maneira pela qual novos grupos podem surgir.

Para que essa população se mantenha constante, também assume-se que cada nascimento implica em uma morte. Portanto, um indivíduo é sorteado aleatoriamente para morrer, com probabilidade uniforme

(34) $$P^m=\frac{1}{N}.$$

A classe de indivíduos que não cooperam (0) pode receber os benefícios da cooperação (tendo chances maiores de ser escolhida pela fêmea devido à formação de grupos), nunca abrindo mão da própria reprodução em benefício de outros. Chamaremos essa classe de Desertores (D). Já os machos com fenótipo (1) se beneficiam dependendo do dominante do grupo a que se associamml: quando o dominante ajudado também tem o fenótipo (1), estes fazem parte da classe que chamaremos de Cooperadores (C), já que contribuem para o crescimento dessa subpopulação (indivíduos dominantes com fenótipo de cooperar também são Cooperadores). Alternativamente, quando o dominante ajudado é um Desertor, estes fazem parte da classe que denominaremos Traidores (T), já que não contribuem para o crescimento da subpopulação de fenótipo (1) (Figura 4b). Grupos de mais de um indivíduo podem ser compostos por apenas Cooperadores, ou por um Desertor dominante e Traidores subordinados (com o segundo tipo se tornando o primeiro, caso o Desertor morra). No instante inicial, todos os $N$ machos são dominantes, portanto há apenas Cooperadores e Desertores.

A partir dessa estruturação, já temos uma proposta de um IBM que modela de forma simplificada a dinâmica dos Chiroxiphia. No caso, modelamos a prevalência do fenótipo cooperativo ao longo do tempo, o que nos permite explorar como esse comportamento pode ser mantido (ou não) na natureza. Para avaliar a frequência fenotípica ao longo do tempo, o IBM exige uma implementação computacional dos processos descritos (a Figura 5a mostra algumas curvas de evolução desse sistema, o script em $R$ está disponível em https://github.com/Natale57/Individual-Based-Cooperation-Model) ou pode ser estudado analiticamente a partir de Equações Mestras. A seguir, mostramos como obter tais equações.

Figura 5
Espaço de fase e evolução temporal (N = $500$). Em a) o eixo horizontal representa os passos de tempo $t$, e o eixo vertical os tamanhos das populações $D$, $T$ e $C$. Em destaque, estão as trajetórias previstas pelas equações ((equação 52)) enquanto, em cinza, estão as trajetórias de $10$ simulações computacionais da dinâmica baseada em indivíduos, utilizando os mesmos parâmetros e condições iniciais. Note que, apesar de começar mais próximo de $P_2^{\ast }$, o sistema se afasta deste ponto fixo, e evolui em direção a $P_1^{\ast }$. Em b) a área sombreada representa a região de interesse para o modelo, na qual todas as dimensões têm valores positivos e menores ou iguais a $N$. O eixo horizontal representa o número de Traidores, e o eixo vertical o número de Cooperadores. As setas em torno de cada ponto representam as direções dos respectivos autovetores. Em preto, estão várias trajetórias que o sistema pode percorrer a depender das condições iniciais. A trajetória em destaque é a mesma representada em a).

3.3. Tratamento analítico

O conjunto de regras apresentadas na seção anterior é simples o suficiente para podermos descrever de forma exata o tamanho médio das populações de cada classe (Desertores, Cooperadores e Traidores) ao longo do tempo. Uma vantagem de escrever essas equações é a possibilidade de identificar os pontos fixos e estudar suas propriedades de estabilidade.

Seguindo os eventos apresentados (Figura 4a), o número de Desertores em um tempo $t+1$ pode ser descrito como o número de indivíduos dessa classe no tempo $t$, mais o número de novos Desertores ($D^{+}$) e menos o número de Desertores que desapareceram ($D^{-}$). Consequentemente, o número médio de Desertores pode ser descrito pelo valor médio desses elementos:

(35) $$\begin{array}{l}{D}_{t+1}=D_t+D_t^{+}-D_t^{-}\to {⟨D⟩}_{t+1}\\ ={⟨D⟩}_t+{⟨D^{+}⟩}_t-{⟨D^{-}⟩}_t.\end{array}$$

sendo $⟨\cdot ⟩$ a média de ensemble no tempo $t$, isto é, a média sobre infinitas realizações da dinâmica no tempo $t$.

Novos Desertores surgem quando um grupo que possui um Desertor dominante é escolhido pela fêmea. Então, considerando um nascimento por passo de tempo, o número médio de novos Desertores pode ser escrito como

(36) $${⟨D^{+}⟩}_t={⟨\sum _{i:D\in n_i}{P}^s\left(n_i\right)⟩}_t,$$

sendo a soma realizada apenas sobre os grupos $i$, de tamanho $n_i$, que contêm um Desertor. Esse somatório descreve o valor esperado de novos Desertores no tempo $t+1$ dada uma população no tempo $t$, enquanto a média $⟨\cdot ⟩$ é realizada sobre todas as possíveis trajetórias/simulações do sistema no tempo $t$. Substituindo a (equação 32) na (equação 36), temos

(37) $${⟨D^{+}⟩}_t=\frac{1}{N}{⟨\sum _{i:D\in n_i}{n}_i⟩}_t.$$

Uma vez que cada grupo de Desertor tem apenas um indivíduo desse tipo (o dominante) e todos os demais são traidores, a soma de todos os $n_i$ é igual ao número de desertores mais o número de Traidores na população,

(38) $${⟨D^{+}⟩}_t=\frac{{⟨D⟩}_t+{⟨T⟩}_t}{N}.$$

Desertores irão desaparecer da população devido apenas à mortalidade. Considerando apenas um evento de morte por unidade de tempo ((equação 34)), chegamos ao número médio de mortes:

(39) $${⟨D^{-}⟩}_t={⟨D{P}^m⟩}_t=\frac{{⟨D⟩}_t}{N}.$$

Substituindo as (equações 38) e (39) na (equação 35), chega-se à expressão:

(40) $$\begin{array}{l}{⟨D⟩}_{t+1}={⟨D⟩}_t+\frac{{⟨D⟩}_t+{⟨T⟩}_t}{N}-\frac{{⟨D⟩}_t}{N},\\ {⟨D⟩}_{t+1}={⟨D⟩}_t+\frac{{⟨T⟩}_t}{N}.\end{array}$$

Note que não há nenhum termo que reduz a população de desertores nesta equação, uma vez que ${⟨T⟩}_t\ge 0$. Isso acontece porque as probabilidades de morte e de escolha da fêmea são uniformes, e é um indício do que irá acontecer em termos de estabilidade no sistema.

Similarmente, o número médio de Traidores em um tempo qualquer depende do número de Traidores no tempo anterior, do número de novos Traidores que surgem e do número de Traidores que desaparecemml:

(41) $${⟨T⟩}_{t+1}={⟨T⟩}_t+{⟨T^{+}⟩}_t-{⟨T^{-}⟩}_t.$$

Diferentemente dos Desertores, o surgimento de novos Traidores depende de dois eventos em conjunto: (i) um novo indivíduo com fenótipo de cooperar (1) precisa surgir; e (ii) esse novo indivíduo precisa escolher se associar a um grupo com Desertor dominante. O evento (i) depende da probabilidade da fêmea escolher um grupo dominado por um Cooperador, e o evento (ii) depende da probabilidade $P^c\left(n\right)$ ((equação 33)) desse novo indivíduo se associar a um grupo que gere Desertores. Portanto, dada uma população de machos no tempo $t$, a probabilidade dos dois eventos acontecerem é o produto das probabilidades individuais, e

(42) $${⟨T^{+}⟩}_t={⟨\left(\sum _{i:D\notin n_i}{P}^s\left(n_i\right)\right)\times \left(\sum _{i:D\in n_i}{P}^c\left(n_i\right)\right)⟩}_t,$$

sendo o primeiro somatório feito apenas nos grupos onde não há Desertores (portanto, são dominados por Cooperadores), e o segundo somatório entre os grupos em que há Desertores (dominados por Desertores). Substituindo as (equações 32) e (33) na (equação 42), temos um somatório no tamanho dos grupos

(43) $${⟨T^{+}⟩}_t={⟨\left(\frac{1}{N}\sum _{i:D\notin n_i}{n}_i\right)\times \left(\frac{1}{N}\sum _{i:D\in n_i}{n}_i\right)⟩}_t.$$

Como antes, sabemos que a soma de todos os grupos do primeiro somatório resulta na população total de Cooperadores, enquanto a do segundo somatório resulta na soma do número de Desertores e Traidores. Portanto,

(44) $${⟨T^{+}⟩}_t=\frac{{⟨C⟩}_t\left({⟨D⟩}_t+{⟨T⟩}_t\right)}{N^2}.$$

Note que aqui aproximamos a média do produto como o produto das médias, ${⟨C\left(D+T\right)⟩}_t\approx {⟨C⟩}_t\left({⟨D⟩}_t+{⟨D⟩}_t\right)$. Esta aproximação é exata para populações infinitas.

A remoção de Traidores da população pode ocorrer de duas formas independentes: um Traidor pode morrer, ou um conjunto de Traidores associados a um Desertor específico pode ir para a população de Cooperadores (no caso da morte de seu dominante, como exemplificado na Figura 4). Portanto,

(45) $$\begin{array}{l}{⟨T^{-}⟩}_t={⟨T{P}^m⟩}_t+{⟨T\to C⟩}_t,\\ {⟨T^{-}⟩}_t=\frac{{⟨T⟩}_t}{N}+{⟨T\to C⟩}_t.\end{array}$$

O termo ${⟨T\to C⟩}_t$ pode ser lido como “o número médio de Traidores do tempo $t$ que se tornarão Cooperadores no tempo $t+1$”. Em um grupo de tamanho $n_i$, existe no máximo um Desertor, sendo os demais $n_i-1$ indivíduos Traidores. Podemos então escrever

(46) $$\begin{array}{l}{⟨T\to C⟩}_t={⟨{\displaystyle \sum _{i:D\in n_i}}\left(n_i-1\right)P^m⟩}_t,\\ {⟨T\to C⟩}_t=\frac{1}{N}{⟨{\displaystyle \sum _{i:D\in n_i}}\left(n_i-1\right)⟩}_t,\\ {⟨T\to C⟩}_t=\frac{{⟨T⟩}_t}{N}.\end{array}$$

Utilizando as (equações 41), (44), (45) e (46), chegamos à equação para a evolução temporal da população média de Traidores

(47) $${⟨T⟩}_{t+1}={⟨T⟩}_t+\frac{{⟨C⟩}_t\left({⟨T⟩}_t+{⟨D⟩}_t\right)}{N^2}-2\frac{{⟨T⟩}_t}{N}.$$

Por fim, analogamente ao realizado para outras classes, o número médio de Cooperadores em um tempo $t+1$ é expresso como

(48) $${⟨C⟩}_{t+1}={⟨C⟩}_t+{⟨C^{+}⟩}_t-{⟨C^{-}⟩}_t.$$

O surgimento de novos Cooperadores depende da fêmea escolher um grupo de Cooperadores para que um novo indivíduo com fenótipo de cooperar (1) apareça ((equação 32)), e de que esse se associe a Cooperadores já existentes ((equação 33)). Além disso, essa classe pode aumentar por conta do recrutamento de Traidores ((equação 46)) na eventual morte de um Desertor dominante. Portanto, podemos escrever ${⟨C^{+}⟩}_t$ como:

(49) $$\begin{array}{l}{⟨C^{+}⟩}_t={⟨\left({\displaystyle \sum _{i:D\notin n_i}}{P}^s\left(n_i\right)\right)\times \left({\displaystyle \sum _{i:D\notin n_i}}{P}^c\left(n_i\right)\right)⟩}_t\\ +{⟨T\to C⟩}_t,\\ {⟨C^{+}⟩}_t=\frac{1}{N^2}{⟨\left({\displaystyle \sum _{i:D\notin n_i}}{n}_i\right)\times \left({\displaystyle \sum _{i:D\notin n_i}}{n}_i\right)⟩}_t+\frac{{⟨T⟩}_t}{N},\\ {⟨C^{+}⟩}_t={\left(\frac{{⟨C⟩}_t}{N}\right)}^2+\frac{{⟨T⟩}_t}{N}.\end{array}$$

Analogamente ao realizado na (equação 44), aproximamos o quadrado da média como sendo igual à média do quadrado, ${⟨C^2⟩}_t\approx {⟨C⟩}_t^2$.

Já o desaparecimento de Cooperadores, similarmente ao desaparecimento de Desertores, depende apenas da mortalidade (equação 34),

(50) $${⟨C^{-}⟩}_t={⟨C{P}^m⟩}_t=\frac{⟨C_t⟩}{N}.$$

Combinando as (equações 48), (49) e (50), temos:

(51) $${⟨C⟩}_{t+1}={⟨C⟩}_t+{\left(\frac{{⟨C⟩}_t}{N}\right)}^2+\frac{{⟨T⟩}_t-{⟨C⟩}_t}{N}.$$

Com isso, chegamos ao sistema de Equações Mestras de tempo discreto que descreve, de forma exata, o número médio de indivíduos de cada classe ao longo do tempo ((equações 40), (47) e (51)):

(52) $$\{\begin{array}{cc}{⟨D⟩}_{t+1}={⟨D⟩}_t+\frac{{⟨T⟩}_t}{N}\hfill & \\ {⟨T⟩}_{t+1}={⟨T⟩}_t+\frac{{⟨C⟩}_t\left({⟨D⟩}_t+{⟨T⟩}_t\right)}{N^2}-\frac{2{⟨T⟩}_t}{N}\hfill & \\ {⟨C⟩}_{t+1}={⟨C⟩}_t+{\left(\frac{{⟨C⟩}_t}{N}\right)}^2+\frac{{⟨T⟩}_t-{⟨C⟩}_t}{N}\hfill & \end{array}.$$

Vale ressaltar que, como a população é constante, a soma dos números (ou números médios) de indivíduos em cada classe em qualquer tempo deve resultar em $N$,

(53) $$N=D_t+T_t+C_t={⟨D⟩}_t+{⟨T⟩}_t+{⟨C⟩}_t.$$

Em nossas análises utilizaremos apenas equações que envolvem os termos da média das populações ((equações 52) e (53)). A fim de simplificar a notação, nas seções seguintes omitiremos o $⟨\cdot ⟩$, e os valores médios das populações de cada classe passam a ser escritos como $D_t$, $T_t$ e $C_t$.

3.4. Pontos fixos

Vamos agora calcular os pontos fixos do modelo partindo das (equações 52) e (53). Para além das possibilidades matemáticas, só temos interesse nos casos em que as grandezas $D_t$, $T_t$ e $C_t$ são maiores ou iguais a $0$, e $N\ge 2$. Isso porque uma população real, cuja dinâmica estamos tentando reproduzir, sempre será constituída por ao menos dois indivíduos que podem interagir entre si (a depender dos fenótipos), e as subpopulações sempre serão constituídas por uma quantidade não negativa.

Pela (equação 12), e supondo a existência de um ponto fixo $P^{\ast }=(D^{\ast },T^{\ast },C^{\ast })$, temos $D_{t+1}=D_t$, se $D_t=D^{\ast }$. O mesmo vale para $T^{\ast }$ e $C^{\ast }$. Aplicar essa propriedade nas (equações 52) resulta em

(54) $$\{\begin{array}{cc}0=\frac{T^{\ast }}{N}\hfill & \\ 0=\frac{C^{\ast }\left(D^{\ast }+T^{\ast }\right)}{N^2}-\frac{2T^{\ast }}{N}\hfill & \\ 0={\left(\frac{C^{\ast }}{N}\right)}^2+\frac{T^{\ast }-C^{\ast }}{N}\hfill & \end{array}.$$

Manipulando as expressões algebricamente para eliminar os denominadores, obtemos

(55) $$\{\begin{array}{cc}0=T^{\ast }\hfill & \\ 0=C^{\ast }\left(D^{\ast }+T^{\ast }\right)-2N{T}^{\ast }\hfill & \\ 0=C^{\ast 2}+N\left(T^{\ast }-C^{\ast }\right)\hfill & \end{array}.$$

Agora, aplicando $T^{\ast }=0$ na segunda expressão das (equações 55), encontramos

(56) $$0=C^{\ast }{D}^{\ast }.$$

Essa igualdade é satisfeita quando $C^{\ast }=0$ ou $D^{\ast }=0$. O caso $C^{\ast }=D^{\ast }=0$ implica em $N=0$ ((equação 53)), e não é de nosso interesse. Para $C^{\ast }=0$, temos $D^{\ast }=N$. Logo,

(57) $$P_1^{\ast }=\left(N\mathrm{,0,0}\right)$$

é um ponto fixo factível, isto é, que possui significado biológico. Ele representa uma situação em que todos os indivíduos da população são Desertores, e podemos descrever qualitativamente $P_1^{\ast }$ como “Ponto de dominação dos Desertores”. Alternativamente, para $D^{\ast }=0$, temos $C^{\ast }=N$. Logo,

(58) $$P_2^{\ast }=(\mathrm{0,0},N)$$

é um segundo ponto fixo factível. Nessa situação, todos os indivíduos da população são Cooperadores, tornando $P_2^{\ast }$ o “Ponto de dominação dos Cooperadores”. Como, para todo ponto fixo, $T^{\ast }=0$, e $C^{\ast }$ e $D^{\ast }$ devem satisfazer a (equação 56), $P_1^{\ast }$ e $P_2^{\ast }$ são os únicos pontos fixos do sistema.

3.5. Estabilidade

Uma vez identificados os pontos fixos, resta determinar a estabilidade deles. Como visto anteriormente, isso requer o cálculo da Matriz Jacobiana nesses pontos e de seus autovalores. Para isso, definimos as funções

(59) $$\{\begin{array}{cc}f(D_t,T_t)=D_{t+1}=D_t+\frac{T_t}{N}\hfill & \\ g(D_t,T_t,C_t)=T_{t+1}=T_t+\frac{C_t\left(D_t+T_t\right)}{N^2}-\frac{2T_t}{N}\hfill & \\ h(T_t,C_t)=C_{t+1}=C_t+{\left(\frac{C_t}{N}\right)}^2+\frac{T_t}{N}-\frac{C_t}{N}\hfill & \end{array},$$

e com o objetivo de simplificar os cálculos fazemos a substituição $D_t=N-C_t-T_t$ ((equação 53)), obtendo

(60) $$\{\begin{array}{cc}f(T_t,C_t)=D_{t+1}=N-C_t-T_t+\frac{T_t}{N}\hfill & \\ g(T_t,C_t)=T_{t+1}=T_t+\frac{C_t\left(N-C_t\right)}{N^2}-\frac{2T_t}{N}\hfill & \\ h(T_t,C_t)=C_{t+1}=C_t+{\left(\frac{C_t}{N}\right)}^2+\frac{T_t}{N}-\frac{C_t}{N}\hfill & \end{array}.$$

Agora, analisando a evolução do sistema em termos de $T_t$ e $C_t$, podemos suprimir a primeira expressão, terminando com

(61) $$\{\begin{array}{cc}g(T_t,C_t)=T_{t+1}=T_t+\frac{C_t\left(N-C_t\right)}{N^2}-\frac{2T_t}{N}\hfill & \\ h(T_t,C_t)=C_{t+1}=C_t+{\left(\frac{C_t}{N}\right)}^2+\frac{T_t}{N}-\frac{C_t}{N}\hfill & \end{array},$$

cujos pontos fixos ((equações 57) e (58)) $P_i^{\ast }=(T_i^{\ast },C_i^{\ast })$ são

(62) $$P_1^{\ast }=\left(\mathrm{0,0}\right)$$

(63) $$P_2^{\ast }=(0,N).$$

Esse procedimento simplifica o problema, o que será muito conveniente nos cálculos seguintes, ao custo de não obtermos informações explícitas sobre como o sistema evolui na dimensão $D$.

A Jacobiana ((equação 27)) das (equações 61) é dada por

(64) $${\mathbf{\text{J}}}_i=\left[\begin{array}{ccc}{\frac{\partial g}{\partial T_t}|}_{{\overrightarrow{P}}_i^{\ast }}& & {\frac{\partial g}{\partial C_t}|}_{{\overrightarrow{P}}_i^{\ast }}\\ {\frac{\partial h}{\partial T_t}|}_{{\overrightarrow{P}}_i^{\ast }}& & {\frac{\partial h}{\partial C_t}|}_{{\overrightarrow{P}}_i^{\ast }}\end{array}\right],$$

onde o índice $i$ indica o ponto fixo ao qual a Jacobiana está associada. Calculando as derivadas parciais obtemos

(65) $${\mathbf{\text{J}}}_i=\left[\begin{array}{ccc}1-\frac{2}{N}& & \frac{1}{N}-\frac{2C_i^{\ast }}{N^2}\\ \frac{1}{N}& & 1-\frac{1}{N}+\frac{2C_i^{\ast }}{N^2}\end{array}\right],$$

e assim, para $P_1^{\ast }$ e $P_2^{\ast }$ respectivamente, temos

(66) $${\mathbf{\text{J}}}_1=\left[\begin{array}{ccc}1-\frac{2}{N}& & \frac{1}{N}\\ \frac{1}{N}& & 1-\frac{1}{N}\end{array}\right],$$

(67) $${\mathbf{\text{J}}}_2=\left[\begin{array}{ccc}1-\frac{2}{N}& & -\frac{1}{N}\\ \frac{1}{N}& & 1+\frac{1}{N}\end{array}\right].$$

Começamos calculando os autovalores de ${\mathbf{\text{J}}}_1$ ((equação 66)) através da equação característica $det\left(\lambda I_2-\mathbf{\text{J}}\right)=0$, onde $I_2$ é a matriz identidade de ordem $2$, resultando em

(68) $$det\left(\left[\begin{array}{ccc}{\lambda }_1-1+\frac{2}{N}& & -\frac{1}{N}\\ -\frac{1}{N}& & {\lambda }_1-1+\frac{1}{N}\end{array}\right]\right)=0,$$

(69) $${\lambda }_1^{\text{2}}+\left(-2+\frac{3}{N}\right){\lambda }_1+1-\frac{3}{N}+\frac{1}{N^2}=0,$$

e pela fórmula quadrática encontramos as seguintes soluções

(70) $$\begin{array}{cc}{\lambda }_1^{+}=1+\frac{-3+\sqrt{5}}{2N}& \\ {\lambda }_1^{-}=1+\frac{-3-\sqrt{5}}{2N}& \end{array}.$$

Como $N\ge 2$, temos que $\left|{\lambda }_1^{+}\right|<1$ e $\left|{\lambda }_1^{-}\right|<1$, de modo que $P_1^{\ast }$ é um atrator. Assim, o Ponto de dominação dos Desertores é um ponto fixo estável (Figura 5).

Para ${\mathbf{\text{J}}}_2$ ((equação 67)) obtemos

(71) $$\begin{array}{cc}{\lambda }_2^{+}=1+\frac{-1+\sqrt{5}}{2N}& \\ {\lambda }_2^{-}=1+\frac{-1-\sqrt{5}}{2N}& \end{array}.$$

Nesse caso, temos $\left|{\lambda }_2^{+}\right|>1$ e $\left|{\lambda }_2^{-}\right|<1$, dado que $N\ge 2$. Descobrimos então que $P_2^{\ast }$ é um ponto de sela e, portanto, o Ponto de dominação dos Cooperadores é um ponto fixo instável (Figura 5). Por completude, no apêndice B calculamos os autovetores associados a cada autovalor.

Em suma, o modelo apresentado, cujo comportamento está sintetizado na Figura 5, mostra que interações cooperativas não se mantêm nas condições modeladas. Foi possível chegar a esse resultado através da aplicação dos conceitos discutidos da teoria de sistemas dinâmicos.

4. Considerações Finais

Neste trabalho, iniciamos revisitando a teoria de sistemas dinâmicos aplicada a equações determinísticas de tempo discreto, e mostramos como estudar a estabilidade local a partir de relações geométricas (Figuras 1, 2 e 3). Ao recorrermos a argumentos geométricos, permitimos com que o trabalho seja autocontido, facilitando a compreensão a respeito dos procedimentos analíticos necessários, como o cálculo da Matriz Jacobiana e de seus autovalores. Assim, esse trabalho pode servir tanto como introdução a sistemas dinâmicos de tempo discreto para estudantes de física e matemática em seus primeiros anos, quanto como material de apoio para estudantes de biologia, em disciplinas como “Introdução a sistemas dinâmicos” e “Matemática para a biologia”.

Em um segundo momento, apresentamos um problema recorrente em estudos de ecologia e evolução, referente à existência de comportamentos cooperativos em comunidades biológicas. No caso de aves das espécies Chiroxiphia caudata e Chiroxiphia linearis, esse comportamento é observado durante o cortejo das fêmeas, no qual um grupo de machos se une em uma série de exibições que aparentam beneficiar apenas um dos membros do grupo, que por fim copula com a fêmea. Propusemos um modelo em que uma população de machos pode ser classificada em três classes: os Desertores, que não possuem fenótipo de cooperar; os Traidores, que possuem o fenótipo de cooperar, mas cooperam com os Desertores; e os Cooperadores, que possuem fenótipo de cooperar e cooperam com Cooperadores. O modelo foi construído seguindo eventos probabilísticos associados a cada indivíduo, portanto um modelo baseado em indivíduos (IBM). Escrevemos então as equações exatas para o comportamento médio do sistema (número médio de Desertores, Traidores e Cooperadores), chegando em um conjunto de equações de tempo discreto determinísticas ((equações 52)). Aplicamos os conceitos de estabilidade abordados na primeira parte do trabalho e verificamos a existência de dois pontos fixos: um em que o sistema é dominado exclusivamente por Cooperadores, e outro exclusivamente dominado por Desertores. Entretanto, apenas o segundo ponto se mostrou estável. Dado tempo suficiente, e havendo ao menos um Desertor, a dinâmica sempre evolui para a dominação completa dos Desertores.

Como foi apresentado, comportamentos cooperativos não devem ser positivamente selecionados na ausência de mecanismos que os tornem vantajosos [⁹]. No modelo, há uma vantagem do comportamento cooperativo para indivíduos dominantes, uma vez que as fêmeas acasalam apenas com o dominante do grupo e a escolha delas é proporcional ao tamanho do grupo. Os Traidores, por fim, perdem a chance de se reproduzir, sem obter nenhum tipo de vantagem evolutiva. Nessa versão do modelo, a escolha do acasalamento feita pelas fêmeas e a escolha do grupo a cooperar, feita por um descendente com fenótipo de cooperar, dependem exclusivamente do tamanho do grupo ((equações 32) e (33)). Estes elementos configuram um cenário evolutivo neutro, o qual acaba por desfavorecer a manutenção da cooperação. Na natureza, a escolha das fêmeas não é homogênea e depende de uma série de características, como sons, cor do macho e consistência da dança [²⁰, ²¹]. A escolha feita pelo descendente com fenótipo de cooperar constitui outra simplificação do modelo, uma vez que não impomos um limite máximo ao tamanho dos grupos. Na natureza, os grupos são geralmente formados por dois ou três indivíduos, e no máximo por seis (no caso dos Chiroxiphia caudata) [¹², ¹³, ¹⁹]. Em estudos futuros, a inclusão destes elementos no modelo pode revelar possíveis mecanismos associados à manutenção da cooperação.

Agradecimentos

Os autores agradecem ao Laboratório de Interações Biológicas e ao Laboratório de Ecologia Comportamental e Ornitologia, ambos da UFPR, por contribuírem na discussão deste trabalho; à professora Lilian Tonelli Manica, pelo auxílio com a literatura sobre comportamentos cooperativos, no desenvolvimento do projeto de pesquisa, e também por permitir a utilização do vídeo com a dança das aves; e a João Lucas Burginski Rosa pela revisão do texto. GB, NCD e SBLA agradecem também ao CNPq pelo apoio financeiro (projeto número 304094/2025−0).

A. Notas Sobre a Estabilidade dos Pontos Fixos

A.1. Caso unidimensional

Em um sistema unidimensional, além de determinar se um ponto fixo é estável ou instável a partir do módulo de $J$, há ainda outras informações que podemos extrair da (equação 9):

• Se $0<J<1$, $x^{\ast }$ é um atrator, e $x_t$ aproxima-se de $x^{\ast }$ monotonicamente (Figura A.1a).

• Se $-1<J<0$, $x^{\ast }$ é um atrator, e $x_t$ aproxima-se de $x^{\ast }$ de forma oscilatória (Figura A.1b).

• Se $J>1$, $x^{\ast }$ é um repulsor, e $x_t$ afasta-se de $x^{\ast }$ monotonicamente (Figura A.1c).

• Se $J<-1$, $x^{\ast }$ é um repulsor, e $x_t$ afasta-se de $x^{\ast }$ de forma oscilatória (Figura A.1d).

Figura A.1
Atratores e repulsores em uma dimensão. Na linha superior, $x^{\ast }$ é um atrator ($\left|J\right|<1$) e, na linha inferior, $x^{\ast }$ é um repulsor ($\left|J\right|>1$). Na coluna da esquerda, o sistema evolui monotonicamente ($0<J\ne 1$) e, na coluna da direita, o sistema evolui de forma oscilatória ($0>J\ne -1$).

A.2. Caso bidimensional

Em um sistema bidimensional, a evolução temporal em torno do ponto fixo depende de dois autovalores:

• Se $\left|{\lambda }_{\mathrm{1,2}}\right|<1$, ${\overrightarrow{x}}^{\ast }$ é um atrator, e ${\overrightarrow{x}}_t$ aproxima-se monotonicamente (Figura A.2a).

• Se $\left|{\lambda }_{\mathrm{1,2}}\right|>1$, ${\overrightarrow{x}}^{\ast }$ é um repulsor, e ${\overrightarrow{x}}_t$ afasta-se monotonicamente (Figura A.2b).

• Se $\left|{\lambda }_1\right|<1$ e $\left|{\lambda }_2\right|>1$, ou $\left|{\lambda }_1\right|>1$ e $\left|{\lambda }_2\right|<1$, ${\overrightarrow{x}}^{\ast }$ é um ponto de sela, e ${\overrightarrow{x}}_t$ afasta-se em uma direção, e aproxima-se em outra (Figura A.2c).

Figura A.2
Pontos fixos e trajetórias no espaço de fase. As retas pretas indicam as direções dos autovetores no espaço de fase, enquanto as interseções entre elas, o ponto fixo estudado. As setas de cada reta indicam como se comporta (naquela direção) uma trajetória próxima do ponto fixo, e assim, a sua estabilidade. No item a) temos um atrator, b) um repulsor e c) um ponto de sela. Os pontos laranja, ligados por curvas pontilhadas, representam trajetórias com base na estabilidade dos pontos fixos, e os números associados a eles indicam a evolução temporal dessas trajetórias.

A.3. Caso n-dimensional

As expressões para $n$ dimensões podem ser obtidas por extrapolação do caso bidimensional. As (equações 10), (12), (13), (14) e (15) permanecem válidas, já que não restringem a quantidade de dimensões. A forma matricial da (equação 10) nesse caso é dada por

(A.1) $$\left[\begin{array}{c}{x}_{t+1}^1\\ ⋮\\ x_{t+1}^n\end{array}\right]=\left[\begin{array}{c}{F}^1(x_t^1,\dots ,x_t^n)\\ ⋮\\ F^n(x_t^1,\dots ,x_t^n)\end{array}\right].$$

na qual temos $n$ variáveis independentes e $n$ funções $F^i$ que definem a evolução do sistema, sendo $x_t^i$ a componente $i$ do vetor ${\overrightarrow{x}}_t$.

A forma geral da Matriz Jacobiana é dada por

(A.2) $$\mathbf{\text{J}}=\left[\begin{array}{ccccc}{\frac{\partial F^1}{\partial x_t^1}|}_{{\overrightarrow{x}}^{\ast }}& & \dots & & {\frac{\partial F^1}{\partial x_t^n}|}_{{\overrightarrow{x}}^{\ast }}\\ ⋮& & \ddots & & ⋮\\ {\frac{\partial F^n}{\partial x_t^1}|}_{{\overrightarrow{x}}^{\ast }}& & \dots & & {\frac{\partial F^n}{\partial x_t^n}|}_{{\overrightarrow{x}}^{\ast }}\end{array}\right],$$

e a interpretação dos autovetores e autovalores é a mesma: dado um ponto fixo do sistema ${\overrightarrow{x}}^{\ast }$, se o autovalor $\left|{\lambda }_i\right|<1$, então pequenas perturbações ao redor de ${\overrightarrow{x}}^{\ast }$, na direção do correspondente autovetor ${\overrightarrow{v}}_i$, tendem a decrescer em módulo, enquanto se $\left|{\lambda }_i\right|>1$, elas tendem a aumentar. Quanto aos critérios de classificação dos pontos fixos temos que, se $\left|{\lambda }_i\right|>1$ para todas as direções $i$, então o ponto fixo é instável, ou repulsor. Se $\left|{\lambda }_i\right|<1$ para todas as direções $i$, o ponto fixo é assintoticamente estável, ou atrator. Quando parte dos autovalores (em módulo) é maior que $1$, e a outra parte é menor que $1$, o ponto fixo é um ponto de sela hiperbólico (também instável).

B. Cálculo dos autovetores

Na Figura 5b, representamos os conjuntos de autovetores associados a cada um dos autovalores dos dois pontos fixos. Apesar de não serem utilizados na análise da estabilidade, os autovetores são outra informação que pode ser extraída de um modelo como o que foi apresentado, e que, como discutido anteriormente, tem um significado próprio: eles informam as direções, no espaço de fase, pelas quais o sistema se afasta ou se aproxima de um ponto fixo dada uma pequena perturbação. Por completude, calcularemos os autovetores do sistema.

Partindo da definição de autovalores e autovetores ((equação 28)), chegamos à expressão

(B.1) $$\left({\mathbf{\text{J}}}_i-{\lambda }_i^{\pm }I\right){\overrightarrow{v}}_i^{\pm }=\overrightarrow{0}.$$

Assim, dada uma matriz ${\mathbf{\text{J}}}_i$ e um autovalor ${\lambda }_i^{\pm }$, podemos calcular o autovetor ${\overrightarrow{v}}_i^{\pm }=\left[\begin{array}{c}{v}_{iT}^{\pm }\\ v_{iC}^{\pm }\end{array}\right]$ correspondente, cujos índices $T$ e $C$ (Traidores e Cooperadores respectivamente) indicam a quais eixos do espaço de fase as componentes pertencem.

Começando com $P_1^{\ast }$, temos a Matriz Jacobiana ${\mathbf{\text{J}}}_1$ ((equação 66)) e os autovalores ${\lambda }_1^{+}$ e ${\lambda }_1^{-}$ ((equações 70)). Para ${\lambda }_1^{+}$

(B.2) $$\left({\mathbf{\text{J}}}_1-{\lambda }_1^{+}I\right){\overrightarrow{v}}_1^{+}=\overrightarrow{0},$$

(B.3) $$\left[\begin{array}{ccc}1-\frac{2}{N}-{\lambda }_1^{+}& & \frac{1}{N}\\ \frac{1}{N}& & 1-\frac{1}{N}-{\lambda }_1^{+}\end{array}\right]\left[\begin{array}{c}{v}_{1T}^{+}\\ v_{1C}^{+}\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right],$$

(B.4) $$\{\begin{array}{cc}\left(\frac{-1-\sqrt{5}}{2}\right)v_{1T}^{+}+v_{1C}^{+}=0\hfill & \\ v_{1T}^{+}+\left(\frac{1-\sqrt{5}}{2}\right)v_{1C}^{+}=0\hfill & \end{array}.$$

Resolvendo esse sistema de equações, obtemos as componentes do autovetor ${\overrightarrow{v}}_1^{+}$ associado a ${\lambda }_1^{+}$ e, repetindo esse mesmo processo para ${\lambda }_1^{-}$, encontramos os seguintes autovetores associados a $P_1^{\ast }$

(B.5) $$\begin{array}{cc}{\overrightarrow{v}}_1^{+}=\left[\begin{array}{c}{r}_1\\ \left(\frac{1+\sqrt{5}}{2}\right)r_1\end{array}\right]& \\ & \\ {\overrightarrow{v}}_1^{-}=\left[\begin{array}{c}{r}_2\\ \left(\frac{1-\sqrt{5}}{2}\right)r_2\end{array}\right]& \end{array},$$

onde $r_1,r_2\in {ℝ}^{\ast }$. As soluções dos sistemas de equações resultam em um conjunto infinito de autovetores parametrizados por $r_1$ e $r_2$, e eles delimitam duas direções ao todo. Similarmente, para $P_2^{\ast }$, aplicamos as (equações 67) e (71) à equação (B.1), e obtemos os autovetores

(B.6) $$\begin{array}{cc}{\overrightarrow{v}}_2^{+}=\left[\begin{array}{c}{r}_3\\ \left(\frac{-3-\sqrt{5}}{2}\right)r_3\end{array}\right]& \\ & \\ {\overrightarrow{v}}_2^{-}=\left[\begin{array}{c}{r}_4\\ \left(\frac{-3+\sqrt{5}}{2}\right)r_4\end{array}\right]& \end{array},$$

onde $r_3,r_4\in {ℝ}^{\ast }$.

Disponibilidade de Dados

Todo o conjunto de dados que dá suporte aos resultados deste estudo foi disponibilizado em um Repositório de Dados:

• Nome do repositório: GitHub

• URL ou DOI: https://github.com/Natale57/Individual-Based-Cooperation-Model

• Data do depósito: 26/09/25

Referências

• ^[1] J.D. Murray, Mathematical Biology I. An Introduction (Springer, New York, 2007).
• ^[2] F. Brauer, C. Castillo-Chavez e Z. Feng, Mathematical Models in Epidemiology (Springer, New York, 2019).
• ^[3] S. Strogatz, Sync: How Order Emerges from Chaos in the Universe, Nature, and Daily Life (Hachette Books, New York, 2004).
• ^[4] M. Ridley, Evolução (Artmed, Porto Alegre, 2006).
• ^[5] L.H.A. Monteiro, Sistemas Dinâmicos (LF Editorial, São Paulo, 2023).
• ^[6] D.L. DeAngelis e V. Grimm, F1000Prime Reports 6, 39 (2014).
• ^[7] C.R. Silva, O. Mesquita, E. Araújo e A.S. Mata, Revista Brasileira de Ensino de Física 47, e20240464 (2025).
• ^[8] N.G. van Kampen, Stochastic Processes in Physics and Chemistry (North Holland, Amsterdam, 2007).
• ^[9] M.A. Nowak, Science 314, 1560 (2006).
• ^[10] S.P. Otto e T. Day, A Biologist’s Guide to Mathematical Modeling in Ecology and Evolution (Princeton University Press, Princeton, 2007).
• ^[11] S.L. Díaz-Muñoz, E.H. DuVal, A.H. Krakauer e E.A. Lacey, Animal Behaviour 88, 67 (2014).
• ^[12] M.S. Foster, Behavioral Ecology and Sociobiology 9, 167 (1981).
• ^[13] M.S. Foster, The American Naturalist 111, 845 (1977).
• ^[14] S.A. West, G.A. Cooper, M.B. Ghoul e A.S. Griffin, Nature Ecology and Evolution 5, 419 (2021).
• ^[15] D. Mitrea, T. Cioara, I. Anghel e L. Toderean, Energy and Buildings 327, 115057 (2025).
• ^[16] J.L. Molina, M.J. Lubbers, H. Valenzuela-Garcia e S. Gómez-Mestres, Anthropology Today 33, 11 (2017).
• ^[17] L. Miranda, A.J.F. Souza, F.F. Ferreira e P.R.A. Campos, PLoS ONE 7, e39188 (2012).
• ^[18] L.A. Dugatkin, Principles of Animal Behavior (University of Chicago Press, 2020).
• ^[19] P.H.L. Ribeiro, A.C. Guaraldo, R.H. Macedo e L.T. Manica, Journal of Ornithology 160, 485 (2019).
• ^[20] P.H.L. Ribeiro, S.B.L. Araujo, T.L. Prado, S.R. Lopes e L.T. Manica, Royal Society Open Science 12, 10 (2025).
• ^[21] L.M. Schaedler, P.H.L. Ribeiro e L.T. Manica, Animal Behavior 171, 29 (2021).

• 1
  Veja um exemplo da exibição em https://youtu.be/AkcGKiwE9fY.

Editado por

• Editor-Chefe:
  Marcello Ferreira https://orcid.org/0000-0003-4945-3169

Datas de Publicação

• Publicação nesta coleção
  06 Mar 2026
• Data do Fascículo
  2026

Histórico

• Recebido
  13 Out 2025
• Revisado
  17 Dez 2025
• Aceito
  24 Jan 2026