Modelo de Gray-Scott: mudanças entre as edições

De Física Computacional
Ir para navegação Ir para pesquisar
Linha 33: Linha 33:
Somando essas duas equações, relacionamos as variáveis <math>u</math> e <math>v</math>:
Somando essas duas equações, relacionamos as variáveis <math>u</math> e <math>v</math>:


:<math> F(1-u) - (F+k)v = 0 \Rightarrow u = 1 - \gamma v </math>
:<math> F(1-u) - (F+k)v = 0 \Rightarrow u = 1 - \gamma v \quad (3) </math>


onde definiu-se o parâmetro auxiliar <math>\gamma = \frac{F+k}{F}</math>.
onde definiu-se o parâmetro auxiliar <math>\gamma = \frac{F+k}{F}</math>.
Linha 39: Linha 39:
Substituindo <math>u</math> na segunda equação do sistema (2) (e reescrevendo <math>F+k =\gamma F</math>), ficamos com:
Substituindo <math>u</math> na segunda equação do sistema (2) (e reescrevendo <math>F+k =\gamma F</math>), ficamos com:


:<math>\left(1 - \gamma v \right)v^2 -\gamma F v = 0 \Rightarrow -\gamma v^3 + v^2 - \gamma F v = 0</math>
:<math>\left(1 - \gamma v \right)v^2 -\gamma F v = 0 \Rightarrow -\gamma v^3 + v^2 - \gamma F v = 0 \quad (4)</math>


Evidentemente, <math>v = 0</math> é solução dessa equação, implicando em <math>u = 1</math>, como já havíamos inspecionado. Alternativamente, considerando <math>v \neq 0</math>, podemos dividir a expressão acima por <math>v</math>, ficando com <math>-\gamma v^2 + v - \gamma F = 0</math>. Resolvendo esta equação quadrática, obtemos duas novas soluções estacionárias para <math>v</math>:
Evidentemente, <math>v = 0</math> é solução dessa equação, implicando em <math>u = 1</math>, como já havíamos inspecionado. Alternativamente, considerando <math>v \neq 0</math>, podemos dividir (4) por <math>v</math>, ficando com <math>-\gamma v^2 + v - \gamma F = 0</math>. Resolvendo esta equação quadrática, obtemos duas novas soluções estacionárias para <math>v</math>:


:<math>v_{\pm} = \frac{1}{2\gamma}(1 \pm \sqrt{1 - 4\gamma^2 F})</math>
:<math>v_{\pm} = \frac{1}{2\gamma}(1 \pm \sqrt{1 - 4\gamma^2 F})</math>


Disso, pela relação <math>u = 1 - \gamma v </math>, temos que os valores correspondentes para <math>u</math> são:
Disso, pela relação (3), temos que os valores correspondentes para <math>u</math> são:


:<math>u_{\mp} = \frac{1}{2}(1 \mp \sqrt{1 - 4\gamma^2 F})</math>
:<math>u_{\mp} = \frac{1}{2}(1 \mp \sqrt{1 - 4\gamma^2 F})</math>
Linha 56: Linha 56:


:<math>\begin{align}
:<math>\begin{align}
  & u^{*}_{0} = 1    & v^{*}_{0} = 0 \\
  & u^{*}_{0} = 1    & v^{*}_{0} & = 0 & \\
  & u^{*}_{1} = \frac{1}{2}(1 + \sqrt{1 - 4\gamma^2 F})  & v^{*}_{1} = \frac{1}{2\gamma}(1 - \sqrt{1 - 4\gamma^2 F})  \\
  & u^{*}_{1} = \frac{1}{2}(1 + \sqrt{1 - 4\gamma^2 F})  & v^{*}_{1} & = \frac{1}{2\gamma}(1 - \sqrt{1 - 4\gamma^2 F}) &\quad (5)  \\
  & u^{*}_{2} = \frac{1}{2}(1 - \sqrt{1 - 4\gamma^2 F})  & v^{*}_{2} = \frac{1}{2\gamma}(1 + \sqrt{1 - 4\gamma^2 F}) \\
  & u^{*}_{2} = \frac{1}{2}(1 - \sqrt{1 - 4\gamma^2 F})  & v^{*}_{2} & = \frac{1}{2\gamma}(1 + \sqrt{1 - 4\gamma^2 F}) & \\
\end{align}</math>   
\end{align}</math>   


Linha 65: Linha 65:
Para avaliar a estabilidade das soluções acima, faz-se necessário obter a matriz Jacobiana dos termos de reação, <math>R_{i}(u,v)</math>. Explicitamente, analisando o sistema (1) de equações, temos que <math>R_{1}(u,v) = -uv^2 + F(1-u)</math> e <math>R_{2}(u,v) = uv^2 - (F+k)v</math>. A matriz Jacobiana do sistema é então dada por:
Para avaliar a estabilidade das soluções acima, faz-se necessário obter a matriz Jacobiana dos termos de reação, <math>R_{i}(u,v)</math>. Explicitamente, analisando o sistema (1) de equações, temos que <math>R_{1}(u,v) = -uv^2 + F(1-u)</math> e <math>R_{2}(u,v) = uv^2 - (F+k)v</math>. A matriz Jacobiana do sistema é então dada por:


:<math>J_{R}(u,v) = \begin{bmatrix}\partial_{u} R_{1}& \partial_{v} R_{1}\\ \partial_{u} R_{2}& \partial_{v} R_{2}\end{bmatrix} = \begin{bmatrix} -v^2 -F& -2uv\\ v^2 & 2uv - (F+k)\end{bmatrix}</math>
:<math>J_{R}(u,v) = \begin{bmatrix}\partial_{u} R_{1}& \partial_{v} R_{1}\\ \partial_{u} R_{2}& \partial_{v} R_{2}\end{bmatrix} = \begin{bmatrix} -v^2 -F& -2uv\\ v^2 & 2uv - (F+k)\end{bmatrix} \quad (6)</math>


Analisemos a estabilidade para os três pares <math>(u^{*}_{i},v^{*}_{i})</math> de soluções estacionárias:
Analisemos a estabilidade para os três pares <math>(u^{*}_{i},v^{*}_{i})</math> de soluções estacionárias:
Linha 71: Linha 71:
* Para <math>(u^{*}_{0},v^{*}_{0}) = (1, 0)</math>:
* Para <math>(u^{*}_{0},v^{*}_{0}) = (1, 0)</math>:


:<math>J_{R}(u^{*}_{0},v^{*}_{0}) = \begin{bmatrix}-F& 0\\ 0 & -(F+k)\end{bmatrix}</math>
:<math>J_{R}(u^{*}_{0},v^{*}_{0}) = \begin{bmatrix}-F& 0\\ 0 & -(F+k)\end{bmatrix} \quad (7)</math>


:Por essa ser uma matriz diagonal, os autovalores <math>\lambda_{i}</math> são justamente as entradas das diagonais; ou seja, <math>\lambda_{1} = -F</math> e <math>\lambda_{2} = -(F+k)</math>. Uma vez que <math>F</math> e <math>k</math> são parâmetros positivos, os dois autovalores são reais e negativos, e portanto o ponto <math>(u^{*}_{0},v^{*}_{0})</math> é '''sempre estável'''.
:Por essa ser uma matriz diagonal, os autovalores <math>\lambda_{i}</math> são justamente as entradas das diagonais; ou seja, <math>\lambda_{1} = -F</math> e <math>\lambda_{2} = -(F+k)</math>. Uma vez que <math>F</math> e <math>k</math> são parâmetros positivos, os dois autovalores são reais e negativos, e portanto o ponto <math>(u^{*}_{0},v^{*}_{0})</math> é '''sempre estável'''.
Linha 87: Linha 87:
:Assim, a matriz jacobiana desses pontos fica:
:Assim, a matriz jacobiana desses pontos fica:


:<math>J_{R}(u^{*}_{i},v^{*}_{i}) = \begin{bmatrix} -(v^{*}_{i})^2 - F & -2(F+k)\\ (v^{*}_{i})^2 & (F+k)\end{bmatrix}, \quad i = 1,2.</math>
:<math>J_{R}(u^{*}_{i},v^{*}_{i}) = \begin{bmatrix} -(v^{*}_{i})^2 - F & -2(F+k)\\ (v^{*}_{i})^2 & (F+k)\end{bmatrix}, \quad i = 1,2 \quad (8)</math>


:Sabemos que o produto dos autovalores dessa matriz é igual ao seu determinante. Calculando-o, obtém-se:
:Sabemos que o produto dos autovalores dessa matriz é igual ao seu determinante. Calculando-o, obtém-se:


:<math>\Delta_{i} := \operatorname{det}\left(J_{R}(u^{*}_{i},v^{*}_{i})\right) = (F+k)\left[(v^{*}_{i})^2 -F \right]</math>
:<math>\Delta_{i} := \operatorname{det}\left(J_{R}(u^{*}_{i},v^{*}_{i})\right) = (F+k)\left[(v^{*}_{i})^2 -F \right] \quad (9)</math>


:Dividindo por <math>F(F+k)</math>:
:Dividindo por <math>F(F+k)</math>:
Linha 99: Linha 99:
:onde se definiu <math>a = \gamma \sqrt{F}</math> ('''observação:''' este é o <math>\gamma</math> definido no Gros). Nota-se que a condição de existência <math>a^2 \leq 1/4</math> para os dois pontos não-triviais é equivalente a <math>F \geq 4(F+k)^2</math>. Expandindo os termos, é possível mostrar que a expressão acima pode ser reescrita como:
:onde se definiu <math>a = \gamma \sqrt{F}</math> ('''observação:''' este é o <math>\gamma</math> definido no Gros). Nota-se que a condição de existência <math>a^2 \leq 1/4</math> para os dois pontos não-triviais é equivalente a <math>F \geq 4(F+k)^2</math>. Expandindo os termos, é possível mostrar que a expressão acima pode ser reescrita como:


:<math>\frac{\Delta_{i}}{F(F+k)} = \frac{\sqrt{1-4a^2}}{2a^2} \left(\sqrt{1-4a^2} \mp 1\right)</math>
:<math>\frac{\Delta_{i}}{F(F+k)} = \frac{\sqrt{1-4a^2}}{2a^2} \left(\sqrt{1-4a^2} \mp 1\right) \quad (10)</math>


:* Para o caso <math>i = 1</math> (sinal negativo), temos a cota superior <math>1-4a^2 < 1</math>. Portanto, <math>\Delta_{1} < 0</math> para todo <math>a</math> que satisfaça a condição de existência. Como o determinante é negativo, sabemos que os autovalores são reais ('''comentário:''' como as entradas da matriz são reais, se os autovalores fossem complexos, seriam também conjugados, de modo que o produto deles fosse igual ao módulo ao quadrado de qualquer um, que seria um valor positivo). Ademais, como seu produto é negativo, eles têm sinais opostos; isto é, um deles é positivo, de modo que o ponto <math>(u^{*}_{1}, v^{*}_{1})</math> '''nunca seja estável'''. Depreendemos desse raciocínio que o determinante da matriz jacobiana de entradas reais ser positivo é uma condição necessária para que haja estabilidade do ponto.
:* Para o caso <math>i = 1</math> (sinal negativo em (10)), temos a cota superior <math>1-4a^2 < 1</math>. Portanto, <math>\Delta_{1} < 0</math> para todo <math>a</math> que satisfaça a condição de existência. Como o determinante é negativo, sabemos que os autovalores são reais ('''comentário:''' como as entradas da matriz são reais, se os autovalores fossem complexos, seriam também conjugados, de modo que o produto deles fosse igual ao módulo ao quadrado de qualquer um, que seria um valor positivo). Ademais, como seu produto é negativo, eles têm sinais opostos; isto é, um deles é positivo, de modo que o ponto <math>(u^{*}_{1}, v^{*}_{1})</math> '''nunca seja estável'''. Depreendemos desse raciocínio que o determinante da matriz jacobiana de entradas reais ser positivo é uma condição necessária para que haja estabilidade do ponto.


:* Já para <math>i = 2</math> (sinal positivo), temos sempre que <math>\Delta_{2} > 0</math>. Para verificar a estabilidade, temos que agora calcular o traço da matriz jacobiana, pois o traço é a soma dos autovalores: se os autovalores são reais, eles têm o mesmo sinal por seu determinante ser positivo, de modo que o traço compartilhe o sinal com os dois autovalores; se os autovalores <math>\lambda_{i}</math> são complexos, eles serão conjugados e o traço será <math>\operatorname{Tr}(J_{R}) = 2 \operatorname{Re}(\lambda_{i})</math>, de modo que a parte real dos autovalores tenha o mesmo sinal do traço. Assim, ''basta que o traço seja negativo para que o ponto seja estável, e que seja positivo para que seja instável''.
:* Já para <math>i = 2</math> (sinal positivo em (10)), temos sempre que <math>\Delta_{2} > 0</math>. Para verificar a estabilidade, temos que agora calcular o traço da matriz jacobiana, pois o traço é a soma dos autovalores: se os autovalores são reais, eles têm o mesmo sinal por seu determinante ser positivo, de modo que o traço compartilhe o sinal com os dois autovalores; se os autovalores <math>\lambda_{i}</math> são complexos, eles serão conjugados e o traço será <math>\operatorname{Tr}(J_{R}) = 2 \operatorname{Re}(\lambda_{i})</math>, de modo que a parte real dos autovalores tenha o mesmo sinal do traço. Assim, ''basta que o traço seja negativo para que o ponto seja estável, e que seja positivo para que seja instável''.


::No caso, temos que <math>\operatorname{Tr}(J_{R}(u^{*}_{2},v^{*}_{2})) = k - (v^{*}_{2})^2</math>. Esse traço é negativo quando <math>(v^{*}_{2})^2 > k</math> e positivo quando <math>(v^{*}_{2})^2 < k</math>; ou seja, <math>(u^{*}_{2}, v^{*}_{2})</math> '''é estável quando''' <math>v^{*}_{2} > \sqrt{k}</math> '''e instável quando''' <math>v^{*}_{2} < \sqrt{k}</math> (lembrando que <math>v^{*}_{2} > 0</math> para todo <math>F</math> e <math>k</math>). Desse modo, pode-se caracterizar uma transição de estabilidade quando <math>v^{*}_{2} = \sqrt{k}</math>.
::No caso, temos que <math>\operatorname{Tr}(J_{R}(u^{*}_{2},v^{*}_{2})) = k - (v^{*}_{2})^2</math>. Esse traço é negativo quando <math>(v^{*}_{2})^2 > k</math> e positivo quando <math>(v^{*}_{2})^2 < k</math>; ou seja, <math>(u^{*}_{2}, v^{*}_{2})</math> '''é estável quando''' <math>v^{*}_{2} > \sqrt{k}</math> '''e instável quando''' <math>v^{*}_{2} < \sqrt{k}</math> (lembrando que <math>v^{*}_{2} > 0</math> para todo <math>F</math> e <math>k</math>). Desse modo, pode-se caracterizar uma transição de estabilidade quando <math>v^{*}_{2} = \sqrt{k}</math>.

Edição das 21h26min de 25 de fevereiro de 2022

Introdução

Descrição do Modelo

O modelo de Gray-Scott descreve uma reação autocatalítica. Sejam duas substâncias químicas cujas concentrações em um dado ponto do espaço são dadas pelas variáveis e , a reação pode ser representada comoFalhou ao verificar gramática (erro de sintaxe): {\displaystyle Inserir fórmula aqui}

Isso significa que uma molécula da substância é transformada em uma molécula da substância por meio da ação de outras duas moléculas da substância , ou seja, é um catalisador de sua própria produção (daí o termo autocatálise). Além dessa reação, ambas substâncias se difundem pelo meio (por isso esse modelo pertence à classe mais geral de modelos reativos-difusivos) e, portanto, as concentrações e mudam com o tempo e diferem em cada ponto. Por simplicidade, assume-se que a reação reversa (i.e., ) não ocorre.

O comportamento geral do sistema pode ser descrito pelas equações abaixo:

A primeira das equações acima pode ser interpretada da seguinte forma. Em um dado ponto, a variação na concentração aumenta: 1) proporcionalmente ao laplaciano de naquele ponto (termo de difusão), i.e., quando a concentração na vizinhança desse ponto é alta; 2) proporcionalmente à taxa de reposição de (taxa de alimentação, ou feed rate). A concentração diminui com o termo reativo , que representa a reação .

De outro lado, a concentração de aumenta com o termo reativo e com o laplaciano de naquele ponto (termo de difusão), mas diminui com a remoção de a uma taxa , mais rápida do que a reposição de .

Análise de estabilidade

Nota: A análise em toda esta seção pressupõe sempre que os parâmetros e coeficientes de difusão são positivos.

Soluções estacionárias sem difusão

O modelo de Gray-Scott depende dos parâmetros e dos coeficientes de difusão das espécies químicas. Ignorando em um primeiro momento os termos de difusão, percebe-se que, por inspeção, o sistema possui uma solução estacionária em para quaisquer valores dos parâmetros. Esse ponto, no entanto, não é a única solução estacionária do sistema; para encontrar as outras, é necessário impor nas equações do sistema. Fazendo isso e dispensando os termos de difusão (), obtém-se o seguinte sistema de equações:

Somando essas duas equações, relacionamos as variáveis e :

onde definiu-se o parâmetro auxiliar .

Substituindo na segunda equação do sistema (2) (e reescrevendo ), ficamos com:

Evidentemente, é solução dessa equação, implicando em , como já havíamos inspecionado. Alternativamente, considerando , podemos dividir (4) por , ficando com . Resolvendo esta equação quadrática, obtemos duas novas soluções estacionárias para :

Disso, pela relação (3), temos que os valores correspondentes para são:

É necessário apontar que, para que as duas últimas soluções (não-triviais) existam — isto é, sejam números reais — o fator dentro da raiz quadrada tem de ser positivo ( ). Por consequência:

, para que existam as soluções não-triviais.

Portanto, há três soluções estacionárias do sistema:[1]

Estabilidade dos estados estacionários (sem difusão)

Para avaliar a estabilidade das soluções acima, faz-se necessário obter a matriz Jacobiana dos termos de reação, . Explicitamente, analisando o sistema (1) de equações, temos que e . A matriz Jacobiana do sistema é então dada por:

Analisemos a estabilidade para os três pares de soluções estacionárias:

  • Para :
Por essa ser uma matriz diagonal, os autovalores são justamente as entradas das diagonais; ou seja, e . Uma vez que e são parâmetros positivos, os dois autovalores são reais e negativos, e portanto o ponto é sempre estável.
  • Para , podemos utilizar uma estratégia que simplifica as contas. Em particular, nota-se que os dois pontos obedecem à segunda equação do sistema (2) com . Desse modo, se dividirmos tal equação por , percebemos que ambos os pontos obedecem a:
Dessa equação, podemos calcular as entradas da segunda coluna da matriz jacobiana com facilidade:
Assim, a matriz jacobiana desses pontos fica:
Sabemos que o produto dos autovalores dessa matriz é igual ao seu determinante. Calculando-o, obtém-se:
Dividindo por :
onde se definiu (observação: este é o definido no Gros). Nota-se que a condição de existência para os dois pontos não-triviais é equivalente a . Expandindo os termos, é possível mostrar que a expressão acima pode ser reescrita como:
  • Para o caso (sinal negativo em (10)), temos a cota superior . Portanto, para todo que satisfaça a condição de existência. Como o determinante é negativo, sabemos que os autovalores são reais (comentário: como as entradas da matriz são reais, se os autovalores fossem complexos, seriam também conjugados, de modo que o produto deles fosse igual ao módulo ao quadrado de qualquer um, que seria um valor positivo). Ademais, como seu produto é negativo, eles têm sinais opostos; isto é, um deles é positivo, de modo que o ponto nunca seja estável. Depreendemos desse raciocínio que o determinante da matriz jacobiana de entradas reais ser positivo é uma condição necessária para que haja estabilidade do ponto.
  • Já para (sinal positivo em (10)), temos sempre que . Para verificar a estabilidade, temos que agora calcular o traço da matriz jacobiana, pois o traço é a soma dos autovalores: se os autovalores são reais, eles têm o mesmo sinal por seu determinante ser positivo, de modo que o traço compartilhe o sinal com os dois autovalores; se os autovalores são complexos, eles serão conjugados e o traço será , de modo que a parte real dos autovalores tenha o mesmo sinal do traço. Assim, basta que o traço seja negativo para que o ponto seja estável, e que seja positivo para que seja instável.
No caso, temos que . Esse traço é negativo quando e positivo quando ; ou seja, é estável quando e instável quando (lembrando que para todo e ). Desse modo, pode-se caracterizar uma transição de estabilidade quando .

(...)

Estabilidade dos estados estacionários (com difusão)

Precisamos agora analisar a estabilidade dos pontos estacionários na presença de difusão, como prescreve o sistema de equações (1), que descreve o modelo. Para isso, é necessário levar em consideração, para cada um dos estados de equilíbrio, os autovalores da matriz , em que é a matriz diagonal cujas entradas são e :[2]

Se escrevermos, genericamente, que , teremos a seguinte matriz jacobiana de reação-difusão:

Como já detalhado acima, para que o ponto seja estável, tal matriz tem que ter a parte real de todos os seus autovalores negativa, de modo que seu determinante seja positivo () e seu traço negativo ().[3] Impondo tais condições à matriz acima, obteremos, após manipulações:[2]

Se o traço é negativo



Se agora incluímos os termos de difusão , deve-se levar em consideração a matriz .

Aqui, é a matriz jacobiana dos termos de reação, é a matriz diagonal dos termos de difusão e é o parâmetro que determina a frequência espacial das perturbações. A demonstração da validade desse método pode ser encontrada na referência[4]. Aplicando ao modelo de Gray-Scott em :



Para que o estado de equilíbrio seja estável é necessário que o determinante da matriz acima seja positivo e o traço seja negativo. Obtém-se então



Ambas desigualdades são imediatamente satisfeitas para quaisquer valores de , e . Portanto, o estado de equilíbrio permanece estável no modelo de Gray-Scott mesmo após a inclusão dos coeficientes de difusão, sejam quais forem os valores desses coeficientes (lembrando que estamos nos restringindo a valores positivos dos parâmetros e coeficientes).

Esse é um resultado à primeira vista surpreendente. Em geral, o surgimento de padrões complexos e não homogêneos em sistemas reativos-difusivos está relacionado à desestabilização de um ou mais estados de equilíbrio homogêneo causada pela introdução dos coeficientes de difusão (conhecida como instabilidade de Turing)[5].

Entretanto, no caso do modelo de Gray-Scott, o surgimento de padrões complexos e não homogêneos não decorre da instabilidade de Turing, uma vez que o surgimento de padrões não triviais nesse modelo ocorre mesmo quando apenas o estado de equilíbrio trivial está presente [1].

Implementação

Será usado o método FTCS (Foward Time Central Space) para integrar as equações do modelo. Como existem explicações do método em toda literatura e em outras entradas da Wiki (ver, por exemplo, Modelo de Turing), a explicação aqui será sucinta.

O método consiste em discretizar a derivada parcial em relação ao tempo para frente e discretizar as derivadas parciais de segunda ordem em relação ao espaço centralmente. Para uma função :





A partir das duas últimas equações acima é fácil mostrar que o laplaciano em duas dimensões, como será usado no presente trabalho, pode ser escrito como



Fazendo , pode-se simplificar a discretização do laplaciano para



Usando a notação é possível então escrever as equações do modelo de forma discretizada:



Utilizou-se uma rede quadrada de tamanho . O estado do inicial do sistema é aquele em que todos os pontos estão no estado de equilíbrio estável trivial , exceto o ponto central, em que é introduzida uma perturbação com . Foram usadas condições de fronteira conforme a Figura 1.

  • Figura 1 - Grid para exemplificar as condições de fronteira usadas na simulação.
  • Cada elemento na matriz tem quatro vizinhos que são denominados por U (Up), D (Down), L (Left), R (Right). Na Figura 1, o elemento , possui os vizinhos , ; o elemento possui como vizinhos , e ; o elemento tem vizinhos , e ; e, finalmente, os vizinhos de são , , e .

    Essas condições de fronteira e a condição inicial explicada acima buscam reproduzir as mesmas condições usadas na simulação de Sayama[6].

    Resultados e discussão

    Modelo de Gray-Scott com
    Alt text
    Concentração de para , de t=0 até t=2000.
    Alt text
    Concentração de para , de t=0 até t=2000.

    Programa

    Simulação do Modelo de Gray-Scott


    Referências

    1. 1,0 1,1 Gros, p
    2. 2,0 2,1 Sayama, pp. 287-289
    3. Sayama, p. 124
    4. Erro de citação: Marca <ref> inválida; não foi fornecido texto para as refs chamadas Sayama260
    5. Week 13, MCB111: Mathematics in Biology (Fall 2021)
    6. Sayama, p. 268

    Bibliografia

    • C. Gros, "Complex and Adaptive Dynamical Systems". Springer-Verlag, Berlim, 2015.
    • H. Sayama, "Introduction to the Modeling and Analysis of Complex Systems". Open SUNY Textbooks, Geneseo, NY, 2015.