Física Computacional - Contribuições do usuário [pt-br]

Transição de fase em dinâmicas de avaliação

2021-06-04T03:48:49Z

Almeidajuliano:

== Falência Coletiva de Empresas e Transição de fase em dinâmicas de avaliação ==
=== Introdução ===

Existem muitos riscos na atividade bancária. Neste trabalho foi investigado uma Dinâmica usada para estudar estes riscos. Poder lidar com isso é importante para que a economia ao todo não seja afetada quando os piores cenários aconteçam. Um deles é o risco de Moratórias (Defaults). Uma moratória é definida como um estado em que um banco determina que um devedor não é mais capaz de pagar uma dívida sem com que o banco tome outras ações (p. e. estender o prazo de pagamento). Com isso em mente, um banco sempre precisa ter uma estimativa das suas perdas dado uma moratória para cada entidade capaz de gerar uma dívida <ref>Basel II: International Convergence of Capital Measurement and Capital Standards: a Revised Framework, The Basel Committee for Banking Supervision, Basel (2004), http://www.bis.org/publ/bcbs107.htm</ref>. Uma dessa entidades são empresas. Porém, empresas não são entidades independentes. Elas tem relações com outras empresas e isso pode torná-las dependentes uma das outras. Baseado em o quão forte for essa dependência, uma moratória em uma empresa pode levar a uma moratória em outra. Duas empresas que dependem do produto uma outra se veem em uma situação que se caso qualquer uma delas tenha dificuldade financeira, a outra sentirá também e, no pior dos casos, pode acabar com dívidas assim como a sua parceira. Isso é chamado de moratórias em conjunto ou falência coletiva.

Dependências podem surgir de fontes diferentes e nem sempre serão positivas. Empresas podem ter uma parceria de trocas, onde uma depende do que outra produz e a que produz depende da compra da outra. Essa dependência também pode ser indireta, na forma de ambas dependerem de um mesmo recurso. Como ambas usam a mesma matéria prima, caso algo aconteça com a mesma, ambas serão afetadas. Uma competição entre empresas também é uma dependência. A diferença é que quando uma empresa for afetada de forma negativa, a outra será afetada de forma positiva e vice-versa.

Nesse trabalho foi utilizada uma dinâmica simples para identificar essas moratórias. Isso foi feito simulando essa interação entre empresas e criando uma avaliação da situação de cada uma delas. As mudanças nessas avaliações vão vir de duas fontes. A primeira será uma Dinâmica Individual de uma firma, isto é, o estado econômico da mesma e como ela age em relação a isso. A segunda fonte será uma Dinâmica Coletiva das interações entre as empresas semelhante ao [[Modelo de Potts]]. Com uma dinâmica desse tipo, veremos duas fases bem definidas no número de moratórias.

===Modelo===

A condição financeira de uma empresa será descrita por uma variável R com valores discretos 0, 1, …, <math>R_{max}</math>. Essa medida de avaliação financeira vai nos dizer que quando R = 0, teremos uma moratória. R é modelado pela Eq. (1).

<math>
R(t) = R(t-1) + s(t) + \eta(R(t-1), s(t))
</math>

Essa avaliação dependerá do seu valor anterior e em cada passo, ela só vai mudar de um em um, pois essa variável <math>s(t)</math> vai ser -1, 0 ou 1. Ou seja, só é possível diminuir em um o R(t) de uma empresa, ou aumenta-lo em 1. O <math>\eta</math> estará modelando as condições de fronteira, já que o espaço dos R's é limitado. A fronteira R = 0 (Moratória) é absorvente, pois uma empresa não irá se recuperar neste modelo. Já o limite superior, ou seja, <math>R_{max+1}</math> é refletivo, pois não se pode ir além. Isto será modelado com <math>\eta(R,s) = -s</math> se R=0, ou <math>\eta(R,s) = -1</math> se <math>R = R_{max}</math> e <math>s = 1</math>, ou <math>\eta(R,s) = 0</math> para os outros casos.

<math>
P(s_i | s_1, ..., s_{i-1}, s_{i+1}, ..., s_N, R_i) = \frac{1}{Z} exp(\sum_{j \neq i} J_{ij} \delta(s_i, s_j) + f(R_i, s_i))
</math>

A função s(t) é semelhante aos spins de partículas para um Modelo de Potts com q = 3. Para N empresas, definiremos uma probabilidade de uma variável <math>s_i</math> de uma empresa <math>i</math> mudar de estado no instante t. Isto é, será uma probabilidade condicional, onde iremos supor uma transição e analisar a probabilidade disso ocorrer no tempo seguinte, baseado em como o sistema está agora e repetir isso para todas as possíveis transições, vide Eq. (2). Tudo antes da barra em (2) é em <math>t+1</math> e tudo depois da barra é em <math>t</math>. Ou seja, assumiremos que em <math>t+1</math>, <math>s_i = -1</math>. Calculando <math>P(s_i = -1 | s_1, s_2, ...)</math> teremos a probabilidade disso acontecer, analisando em especial o Delta de Kronecker. Ele só será não nulo com os <math>s_j</math> que estão no mesmo estado, ou seja, estamos contando quantos estados são iguais ao <math>s_i</math> que foi suposto. O termo <math>J_{ij}</math> será uma matriz que vai modelar a interação entre as empresas, ou a dependência entre i e j. Cada valor nesta matriz é Gaussiano (Eq. 3) com média <math>J_0</math> e desvio <math>\sigma_J</math>. Além disso Z será uma constante de normalização, para que <math>P(s_i = -1 | s_1, …) + P(s_i = 0 | s_1, …) + P(s_i = 1 | s_1, …) = 1</math> e <math>f(R_i, s_i)</math> será o termo que vai modelar a dinâmica interna da empresa (Veremos o caso <math>f=0</math> e o caso <math>f \neq 0</math>). Com isso feito, teremos a possibilidade do nodo <math>s_i</math> ir para -1 em <math>t+1</math>. Estenderemos isso para <math>P(s_i = 0| s_1, ...)</math> e <math>P(s_i = 1| s_1, ...)</math> e teremos todas as possibilidades de transição de <math>s_i</math> em <math>t+1</math>. Basta agora repetir isso para as N empresas e escolher aleatoriamente o destino de cada nodo baseado nessas probabilidades obtidas.

Os valores <math>J_{ij}</math> podem ser tanto positivos quanto negativos. Quando <math>J_{ij} > 0</math>, haverá uma interação positiva entre duas empresas, ou uma cooperação. Essa cooperação pode ser na forma de uma relação de compra e troca, ou na dependência de mesmos recursos. Se um dos nodos tiver algum problema financeiro, ele terá o seu capital afetado e logo diminuirá os lucros dos seus parceiros. Uma mudança em i causará uma mudança em j na mesma direção. Caso <math>J_{ij} < 0</math>, haverá uma concorrência entre i e j. Uma mudança em i causará uma mudança em j na direção oposta.

<math>
P(J_{ij}) = \frac{1}{\sqrt{2\pi\sigma_J^2}} exp(-\frac{(J_{ij} - J_0)^2}{2\sigma_J^2})
</math>

Definiremos também o número de moratórias <math>ND</math> (Number of Defaults). Sabemos que há uma moratória quando <math>R=0</math>, então <math>ND</math> será quantos <math>R_i</math> nulos teremos nas N empresas após toda uma simulação. Mas, como <math>J_{ij}</math> será um valor diferente de simulação para simulação e isso levará a valores diferentes de <math>ND</math>, iremos definir uma média de <math>ND</math> como <math><ND></math>. Na equação abaixo, K é o número de simulações feitas e <math>ND_k</math> é o <math>ND</math> de uma específica simulação.

<math>
<ND> = \frac{1}{K} \sum_{k = 1}^K ND_k
</math>

===Simulações===
== Referências ==

<references/>

Transição de fase em dinâmicas de avaliação

2021-06-04T03:46:02Z

Almeidajuliano:

== Falência Coletiva de Empresas e Transição de fase em dinâmicas de avaliação ==
=== Introdução ===

Existem muitos riscos na atividade bancária. Neste trabalho foi investigado uma Dinâmica usada para estudar estes riscos. Poder lidar com isso é importante para que a economia ao todo não seja afetada quando os piores cenários aconteçam. Um deles é o risco de Moratórias (Defaults). Uma moratória é definida como um estado em que um banco determina que um devedor não é mais capaz de pagar uma dívida sem com que o banco tome outras ações (p. e. estender o prazo de pagamento). Com isso em mente, um banco sempre precisa ter uma estimativa das suas perdas dado uma moratória para cada entidade capaz de gerar uma dívida <ref>Basel II: International Convergence of Capital Measurement and Capital Standards: a Revised Framework, The Basel Committee for Banking Supervision, Basel (2004), http://www.bis.org/publ/bcbs107.htm</ref>. Uma dessa entidades são empresas. Porém, empresas não são entidades independentes. Elas tem relações com outras empresas e isso pode torná-las dependentes uma das outras. Baseado em o quão forte for essa dependência, uma moratória em uma empresa pode levar a uma moratória em outra. Duas empresas que dependem do produto uma outra se veem em uma situação que se caso qualquer uma delas tenha dificuldade financeira, a outra sentirá também e, no pior dos casos, pode acabar com dívidas assim como a sua parceira. Isso é chamado de moratórias em conjunto ou falência coletiva.

Dependências podem surgir de fontes diferentes e nem sempre serão positivas. Empresas podem ter uma parceria de trocas, onde uma depende do que outra produz e a que produz depende da compra da outra. Essa dependência também pode ser indireta, na forma de ambas dependerem de um mesmo recurso. Como ambas usam a mesma matéria prima, caso algo aconteça com a mesma, ambas serão afetadas. Uma competição entre empresas também é uma dependência. A diferença é que quando uma empresa for afetada de forma negativa, a outra será afetada de forma positiva e vice-versa.

Nesse trabalho foi utilizada uma dinâmica simples para identificar essas moratórias. Isso foi feito simulando essa interação entre empresas e criando uma avaliação da situação de cada uma delas. As mudanças nessas avaliações vão vir de duas fontes. A primeira será uma Dinâmica Individual de uma firma, isto é, o estado econômico da mesma e como ela age em relação a isso. A segunda fonte será uma Dinâmica Coletiva das interações entre as empresas semelhante ao [[Modelo de Potts]]. Com uma dinâmica desse tipo, veremos duas fases bem definidas no número de moratórias.

===Modelo===

A condição financeira de uma empresa será descrita por uma variável R com valores discretos 0, 1, …, <math>R_{max}</math>. Essa medida de avaliação financeira vai nos dizer que quando R = 0, teremos uma moratória. R é modelado pela Eq. (1).

<math>
R(t) = R(t-1) + s(t) + \eta(R(t-1), s(t)) \tag{1}
</math>

Essa avaliação dependerá do seu valor anterior e em cada passo, ela só vai mudar de um em um, pois essa variável <math>s(t)</math> vai ser -1, 0 ou 1. Ou seja, só é possível diminuir em um o R(t) de uma empresa, ou aumenta-lo em 1. O <math>\eta</math> estará modelando as condições de fronteira, já que o espaço dos R's é limitado. A fronteira R = 0 (Moratória) é absorvente, pois uma empresa não irá se recuperar neste modelo. Já o limite superior, ou seja, <math>R_{max+1}</math> é refletivo, pois não se pode ir além. Isto será modelado com <math>\eta(R,s) = -s</math> se R=0, ou <math>\eta(R,s) = -1</math> se <math>R = R_{max}</math> e <math>s = 1</math>, ou <math>\eta(R,s) = 0</math> para os outros casos.

<math>
P(s_i | s_1, ..., s_{i-1}, s_{i+1}, ..., s_N, R_i) = \frac{1}{Z} exp(\sum_{j \neq i} J_{ij} \delta(s_i, s_j) + f(R_i, s_i)) \tag{2}
</math>

A função s(t) é semelhante aos spins de partículas para um Modelo de Potts com q = 3. Para N empresas, definiremos uma probabilidade de uma variável <math>s_i</math> de uma empresa <math>i</math> mudar de estado no instante t. Isto é, será uma probabilidade condicional, onde iremos supor uma transição e analisar a probabilidade disso ocorrer no tempo seguinte, baseado em como o sistema está agora e repetir isso para todas as possíveis transições, vide Eq. (2). Tudo antes da barra em (2) é em <math>t+1</math> e tudo depois da barra é em <math>t</math>. Ou seja, assumiremos que em <math>t+1</math>, <math>s_i = -1</math>. Calculando <math>P(s_i = -1 | s_1, s_2, ...)</math> teremos a probabilidade disso acontecer, analisando em especial o Delta de Kronecker. Ele só será não nulo com os <math>s_j</math> que estão no mesmo estado, ou seja, estamos contando quantos estados são iguais ao <math>s_i</math> que foi suposto. O termo <math>J_{ij}</math> será uma matriz que vai modelar a interação entre as empresas, ou a dependência entre i e j. Cada valor nesta matriz é Gaussiano (Eq. 3) com média <math>J_0</math> e desvio <math>\sigma_J</math>. Além disso Z será uma constante de normalização, para que <math>P(s_i = -1 | s_1, …) + P(s_i = 0 | s_1, …) + P(s_i = 1 | s_1, …) = 1</math> e <math>f(R_i, s_i)</math> será o termo que vai modelar a dinâmica interna da empresa (Veremos o caso <math>f=0</math> e o caso <math>f \neq 0</math>). Com isso feito, teremos a possibilidade do nodo <math>s_i</math> ir para -1 em <math>t+1</math>. Estenderemos isso para <math>P(s_i = 0| s_1, ...)</math> e <math>P(s_i = 1| s_1, ...)</math> e teremos todas as possibilidades de transição de <math>s_i</math> em <math>t+1</math>. Basta agora repetir isso para as N empresas e escolher aleatoriamente o destino de cada nodo baseado nessas probabilidades obtidas.

Os valores <math>J_{ij}</math> podem ser tanto positivos quanto negativos. Quando <math>J_{ij} > 0</math>, haverá uma interação positiva entre duas empresas, ou uma cooperação. Essa cooperação pode ser na forma de uma relação de compra e troca, ou na dependência de mesmos recursos. Se um dos nodos tiver algum problema financeiro, ele terá o seu capital afetado e logo diminuirá os lucros dos seus parceiros. Uma mudança em i causará uma mudança em j na mesma direção. Caso <math>J_{ij} < 0</math>, haverá uma concorrência entre i e j. Uma mudança em i causará uma mudança em j na direção oposta.

<math>
P(J_{ij}) = \frac{1}{\sqrt{2\pi\sigma_J^2}} exp(-\frac{(J_{ij} - J_0)^2}{2\sigma_J^2}) \tag{3}
</math>

Definiremos também o número de moratórias <math>ND</math> (Number of Defaults). Sabemos que há uma moratória quando <math>R=0</math>, então <math>ND</math> será quantos <math>R_i</math> nulos teremos nas N empresas após toda uma simulação. Mas, como <math>J_{ij}</math> será um valor diferente de simulação para simulação e isso levará a valores diferentes de <math>ND</math>, iremos definir uma média de <math>ND</math> como <math><ND></math>. Na equação abaixo, K é o número de simulações feitas e <math>ND_k</math> é o <math>ND</math> de uma específica simulação.

<math>
<ND> = \frac{1}{K} \sum_{k = 1}^K ND_k \tag{4}
</math>

===Simulações===
== Referências ==

<references/>

Transição de fase em dinâmicas de avaliação

2021-05-29T03:51:52Z

Almeidajuliano:

Transição de fase em dinâmicas de avaliação

2021-05-29T03:51:08Z

Almeidajuliano:

Transição de fase em dinâmicas de avaliação

2021-05-29T03:46:08Z

Almeidajuliano: /* Falência Coletiva de Empresas e Transição de fase em dinâmicas de avaliação */

Transição de fase em dinâmicas de avaliação

2021-05-29T03:43:50Z

Almeidajuliano: /* Falência Coletiva de Empresas e Transição de fase em dinâmicas de avaliação */

Transição de fase em dinâmicas de avaliação

2021-05-29T03:42:52Z

Almeidajuliano: /* Introdução */

Transição de fase em dinâmicas de avaliação

2021-05-29T03:40:44Z

Almeidajuliano: /* Introdução */

Transição de fase em dinâmicas de avaliação

2021-05-29T03:36:47Z

Almeidajuliano:

Transição de fase em dinâmicas de avaliação

2021-05-29T03:04:06Z

Almeidajuliano: /* Introdução */

Transição de fase em dinâmicas de avaliação

2021-05-29T02:57:36Z

Almeidajuliano: Criou página com '== Falência Coletiva de Empresas e Transição de fase em dinâmicas de avaliação == === Introdução ==='

== Falência Coletiva de Empresas e Transição de fase em dinâmicas de avaliação ==
=== Introdução ===

Trabalhos 2020-2

2021-05-29T02:08:28Z

Almeidajuliano:

=== [[Modelo de Potts 2D]] ===
=== [[Misturas binárias na rede 2D]] ===
=== [[Clusterização]] ===
=== [[Transição de fase em dinâmicas de avaliação]] ===

Grupo - Solução Exata para Movimento Anisotrópico de Ornstein-Uhlenbeck

2020-01-10T02:24:38Z

Almeidajuliano:

Aluno: Juliano Almeida Machado - Física
==Introdução==

Movimento coletivo ordenado de diversos organismos é um fenômeno que pode ser observado em diversas situações, como em movimento de cardumes, migração de pássaros, migração celular, etc..

Aquilo que é mais curioso destes tipos de movimento é que o mesmo surge através de movimentos aleatórios dos indivíduos que o compõem, ou seja, o movimento individual e aleatório de um grupo de células quando postos a interagir entre si leva a um movimento ordenado em curtas escalas de tempo.

Um dos modelos anteriores a este, o modelo de ''Vicsek'' [https://en.wikipedia.org/wiki/Vicsek_model], visa descrever o movimento coletivo ordenado de partículas “inteligentes”, isto é, partículas que podem observar as outras partículas na vizinhança e mover-se de acordo com a média do movimento do grupo. Um bom exemplo disso é o movimento de pássaros ou peixes, da onde movimentos ordenados podem emergir de movimentos aleatórios dos indivíduos autopropulsionados. Curiosamente, mesmo que estes sistemas estejam longe de um estado de equilíbrio, os mesmos apresentam características típicas de sistemas no equilíbrio, como transições de fase bem definidas.

Contudo, partículas “não inteligentes” já haviam demonstrado o mesmo tipo de movimento coletivo experimentalmente, com as mesmas mudanças de fase bem definidas. É para explicar isso, que o modelo de ''Szabó'' ([[Grupo - Modelo de Szabó]]) foi criado. Nele, o mesmo fenômeno é descrito utilizando somente as forças de interação intercelulares já que isso é algo que independe da capacidade da célula de poder observar a média do movimento da vizinhança, como em ''Vicsek''.

Essas forças são repulsivas para distâncias (<math>d_{ij} = |r_{i}-r_{j}|</math>) menores que <math>R_{eq}</math>, sendo <math>F_{rep}</math> o seu máximo para esse caso, atrativa para <math>R_{eq} \le d_{ij} \le R_{o}</math>, sendo <math>F_{adh}</math> o máximo para esse caso, e zero para <math>d_{ij}</math> maior que um certo <math>R_{o}</math>.

<math>
\textbf{F}(\textbf{r}_{i}, \textbf{r}_{j} ) = \textbf{e}_{ij}\begin{cases}
F_{rep}\frac{d_{ij} - R_{eq}}{R_{eq}},\quad d_{ij} < R_{eq} \\
F_{adh}\frac{d_{ij} - R_{eq}}{R_{0} - Re_{eq}}, \; \; R_{eq} \le d_{ij} \le R_{0}\\
0, \quad \quad \quad \quad \quad \; d_{ij} > R_{0}\\
\end{cases}
</math>

[[Arquivo:Forcaszabo.png|frame|center|Figura 1: Esquematização da força intercelular de Szabó.]]

Essas forças levam às seguintes equações de movimento:

<math>\frac{dr_{i}(t)}{dt} = v_{o}n_{i}(t) + \mu \sum_{j=1}^N F(r_{i}, r{j})</math>

Onde <math>N</math> é o número de partículas e <math>\sum_{j=1}^N F(r_{i}, r{j})</math> é o somatório da todas as interações sobre uma partícula ''i''.

<math> \frac{d_{i}^{n}(t)}{dt} = \frac{1}{\tau}arcsin[(n_{i}(t)\frac{v_{i}(t)}{|v_{i}(t)|})e_{z}] + \xi </math>

Onde <math>\tau</math> é o tempo de relaxamento, <math>n_{i}(t)</math> é o vetor na direção de polarização, tal que <math>n_{i}(t) = (cos(\theta_{i}^{n}(t)), sen(\theta_{i}^{n}(t)))</math>, <math>e_{z}</math> é o vetor normal, e <math>\xi</math> é um ruído branco.

Neste modelo, o movimento bidimensional das células é descrito por uma dinâmica superamortecida. Na mesma, cada célula tenta manter uma velocidade de magnitude <math>v_{o}</math> na direção do vetor <math>n_{i}(t)</math>, sendo que a direção desse vetor é descrita pelo ângulo <math>\theta_{i}^{n}(t)</math> que tenta relaxar na direção da velocidade enquanto sofre ruído gaussiano.

Embora descreva bem este movimento, o modelo de Szabó fornece somente soluções numéricas para o problema. Para obtermos uma solução analítica, o problema precisa ser descrito de outra forma. É com esse objetivo que o modelo de Rita et al. foi feito.

[[Arquivo:Szaboilus.png|frame|center|Figura 2: Ilustração do resultado obtido por Szabó.]]

==Modelo==
Neste modelo, em vez de ter a relaxação da velocidade na direção da polarização, a mesma é projetada sobre a polarização, enquanto ela muda com ruído. Essa direção de polarização interna pode representar, p.e., a polarização dos filamentos de actina na direção preferencial em migração celular.

Com isso pode-se observar duas transições, ou seja, uma transição de um estado desordenado para um estado ordenado e para um estado desordenado de novo. Essas duas transições representam a aleatoriedade do movimento inicial das células, depois a escolha de uma direção preferencial e por último a variação desta direção preferencial.

Contudo, neste modelo deve-se assumir que, para a partícula, seu movimento é anisotrópico, ou seja, apresenta assimetria nas dinâmicas de movimento em direções diferentes, uma na direção paralela à direção de polarização e outra perpendicular.

===Movimento Paralelo===
Na direção paralela à polarização, temos o sistema agindo numa dinâmica semelhante a de ''Langevin'' (<math>m\mathbf{a} = m\frac{d\mathbf{v}}{dt} = F(\mathbf{x}) - \beta \mathbf{v} + \eta(t)</math>)[https://en.wikipedia.org/wiki/Langevin_equation], mais especificamente, um processo de ''Orstein-Uhlenbeck'' amortecido. Nesta direção a velocidade (paralela) é bem definida. Tanto para o movimento paralelo, quanto para o movimento perpendicular, <math>p_{i}(t) = (cos(\theta_{i}(t)), sen(\theta_{i}(t)))</math> é o vetor que indica a direção de polarização e <math>n_{i}(t) = (sen(\theta_{i}(t)), -cos(\theta_{i}(t)))</math> é o vetor perpendicular a <math>p_{i}(t)</math>.

<math>\frac{dv_{||}}{dt} = -\frac{\gamma}{m}v_{||} + \frac{\xi_{||}(t)}{m} </math>
<math>v_{||}(t+\Delta t)p(t+\Delta t) = [(1-\gamma\Delta t)v_{||}(t) + \int_{t}^{t+\Delta t} \xi_{||}(t)dt](p_{i}(t)p_{i}(t+\Delta t))p_{i}(t+\Delta t) </math>

Onde <math>v_{||}</math> é a velocidade paralela, <math>\gamma</math> é o coeficiente de atrito e <math>\xi_{||}</math> é o ruído gaussiano paralelo.

A velocidade diminui por causa do atrito com o substrato, já que pode-se observar que <math>(1-\gamma\Delta t)</math> tende a <math>e^{-\gamma\Delta t}</math> se <math>\Delta t</math> for pequeno, e diminui ou aumenta por causa do termo estocástico <math>\int_{t}^{t+\Delta t} \xi_{||}(t)dt = W(t)</math>. A primeira equação assume que durante o intervalo <math>\Delta t</math>, a direção de polarização é constante. Somente em <math>t+\Delta t</math>, a velocidade na direção <math>p_{i}(t)</math> é projetada sobre a nova direção de polarização <math>p_{i}(t+\Delta t)</math>.

===Movimento Perpendicular===
Na direção perpendicular à polarização, o sistema segue um processo de ''Wiener'' (movimento aleatório) por causa das flutuações na direção de polarização.

<math>\frac{dr_{\perp}}{dt} = n_{i}(t)\xi_{\perp}(t) </math>
<math>[r_{\perp}(t+\Delta t) - r_{\perp}(t)]n_{i}(t) = n_{i}(t)\int_{t}^{t+\Delta t} \xi_{\perp}(t)dt </math>

Onde <math>\xi_{\perp}</math> é o ruído gaussiano na direção perpendicular. Pode-se observar que <math>\int_{t}^{t+\Delta t} \xi_{\perp}(t)dt = W(t) </math>. Este movimento, por ser puramente um processo de ''Wiener'', é melhor descrito pelo deslocamento do que pela velocidade.

===Ângulo e Ruídos===
A direção de polarização é definida pelo ângulo teta que também apresenta ruído branco. Esse ruído <math>\beta(t)</math> é relacionado à <math> \xi_{\perp}(t)</math> através de <math>\xi_{perp}(t) = \sqrt q\beta(t)</math>. É essa mudança na direção de polarização que causa o deslocamento aleatório perpendicular.

Embora o ruído no ângulo e o ruído na direção perpendicular estejam relacionados, os mesmos são independentes do ruído na direção paralela. Abaixo pode-se ver a correlação dos mesmos. Os três são ruídos gaussianos, então a média deles é zero.

<math> <\xi_{||}(t)\xi_{||}(t')> = g\delta(t-t') </math>
<math> <\beta(t)\beta(t')> = 2k\delta(t-t') </math>
<math> <\xi_{\perp}(t)\xi_{\perp}(t')> = 2qk\delta(t-t') </math>

Onde <math>k</math> tem unidades de 1/tempo, <math>\sqrt q</math> tem unidades de comprimento, e tanto <math>g</math> quanto <math>qk</math> tem unidades de [comprimento]^2 / tempo.

Em resumo, este modelo considera uma partícula com dois graus de liberdades espaciais e um grau de liberdade interna, que considera a polarização interna. Esse grau de liberdade interna quebra a simetria espacial de tal forma que a dinâmica na direção da polarização é semelhante a ''Langevin'' e na direção perpendicular a dinâmica segue um processo de ''Wiener''. Existem duas fontes independentes de ruído: Uma age na dinâmica na direção de polarização e a outra age na polarização em si. Além da força de dissipação, a mudança na polarização age como mais um termo na perda de energia cinética.

[[Arquivo:Wforce.gif|frame|350x350px|center|Figura 3: Simulação do modelo de Rita et al. A seta azul representa a direção de polarização e a seta verde a velocidade.]]

===Mean Square Displacement(MSD)===
O MSD é uma quantificação do movimento líquido das partículas. O mesmo foi utilizado para medir a aglomeração e a transição de fase das mesmas. Para isso foi feita uma média do MSD de 100 trajetórias sendo que o Método das Janelas foi feito para cada uma das mesmas. Através do resultados analíticos para as equações de movimento pode-se obter um resultado analítico para o MSD também:

<math><r^{2}> = \frac{4DP}{1-S}[\tau - (1-S)(1-e^{-\tau})]</math>

Onde <math>P</math> é o tempo de persistência do movimento em um estado cinético, <math>S</math> é uma fração deste tempo, relacionado a quanto tempo o movimento persiste no primeiro estado cinético e <math>D =\frac{g}{2m^{2}\tau\tau_{b}}</math>, sendo que <math>\tau = \frac{\gamma}{m} + k </math> e <math>\tau_{b} = \frac{\gamma}{m} + 2k </math>.

Para testar este modelo resolveu-se numericamente as equações diferenciais de movimento com um algoritmo simples de Euler-Maruyama e comparadas com os resultados analíticos do MSD para diferentes valores do parâmetro <math>k</math>.

O resultado obtido consta abaixo. Nele os pontos são os valores obtidos pela simulação e as linhas são as previsões da solução analítica para o MSD. No mesmo pode-se claramente ver as duas transições de estado cinético, especialmente para valores pequenos de <math>k</math>.

[[Arquivo:Msd_v2.png|frame|center|Figura 4: MSD analítico para 100 trajetórias]]

===O Programa===
O programa criado para gerar simulações do modelo descrito e calcular o MSD do mesmo foi feito em C. O mesmo começa com uma inialização de todas as variáveis necessárias e a definição das mesmas de acordo com aquilo que foi descrito acima sendo os parâmetros iniciais:<math>k, q, g</math> e <math> \gamma</math> e os parâmetros termodinâmicos <math> S,P </math> e <math>D</math> que descrevem o comportamento macroscópico da célula.

 
<source lang="c">
//Escrito por Mendeli Vainstein;

#define N 100

#define GAMMA 1.0 //3.0
#define M 1.0 //5.0
#define G 1.0 //7.0
#define K 0.00 // 0.01 //0.001

#define T_MAX (1e2) // 1e2 //5e1 //1e2 //5e2
#define DT (1e-4) // 1e-5//1e-5 //1e-6
#define TRANSIENT (2e1)

#define NUM_CORR_MEASURES 200
#define NUM_VCORR_MEASURES 200
#define T_MAX_VCORR 20

unsigned long int rseed;

double *rx,*dr_par,*vpar,*dx,*dy;
double *ry,*dr_perp,*vperp;
double *v2, *corr_v, *r2, *v2_ave,*v_ave, *corr_v_ave;

double K;
int N, MAX_STEPS, TRANSIENT_STEPS, STEP_SIZE_VCORR;
long corr_time[NUM_CORR_MEASURES];//, vcorr_time[NUM_VCORR_MEASURES];

double TAUB, GAMMA_M, TAU2,GAMMA_M2,V02;

double DT_DIV_2,M2, TAU, Q, SQRT_Q;
double EXP_VPAR, SQRT_DT;
double XI_PAR,XI_PAR_SQRT_DT,XI_PAR_SQRT_DT_DIV_M;
double XI_THETA,XI_THETA_SQRT_DT; // Used if langevin theta dynamics

double *theta, *d_theta;
double two_PI = 2*M_PI;
const gsl_rng_type * T;
gsl_rng * rand_vec;

char filename[200];
FILE *fcorr, *ftraj, *fr2;

/************************************************************************************************
* *
************************************************************************************************/
void allocate_pointers(int n);
void free_pointers(int n);
void measure(double *rx, double *ry,double *rx0, double *ry0, double *v2, unsigned long int t);
void set_initial_conditions(int N);
void set_correlation_vecs();
void calculate_correlations(double *rx, double *ry, double *v);
void print_correlations(void);
void reset_initial_condition(int t_counter);
void set_gsl_rng(void);
double wrap_angle(double theta);
void create_time_table(long *t1, int total_time, int measures);

/************************************************************************************************
* *
************************************************************************************************/
</source> 

Inicialização:
 
<source lang="c">
void set_constants(void)
{
MAX_STEPS = (int)(T_MAX/DT) + 1;
TRANSIENT_STEPS = TRANSIENT/DT;

DT_DIV_2 = 0.5*DT;
SQRT_DT = sqrt(DT);

M2 = M*M;
TAU = GAMMA/M + K;
TAU2 = TAU*TAU;
TAUB = GAMMA/M + 2*K;
GAMMA_M = GAMMA/M;
GAMMA_M2 = GAMMA_M*GAMMA_M;

XI_THETA = sqrt(2.0*K);
XI_THETA_SQRT_DT = XI_THETA * SQRT_DT; // Used if langevin theta dynamics

XI_PAR = sqrt(G); //sqrt(2.0*(GAMMA+K*M)*kT);
XI_PAR_SQRT_DT = XI_PAR * SQRT_DT;
XI_PAR_SQRT_DT_DIV_M = XI_PAR * SQRT_DT/M;

EXP_VPAR = exp(-DT_DIV_2*GAMMA/M); //For langevin algorithm

V02 = G/(2*M2*TAU); //G/(2*M2*(GAMMA/M + K))

Q = V02/TAU2; //(kT/M)/TAU2;//v2/(TAU*TAU);
SQRT_Q = sqrt(Q);

create_time_table(corr_time, MAX_STEPS, NUM_CORR_MEASURES);
//create_time_table(vcorr_time, T_MAX_VCORR/DT, NUM_VCORR_MEASURES);
STEP_SIZE_VCORR = T_MAX_VCORR/(NUM_VCORR_MEASURES*DT);

return;
}
</source> 

Em seguida foram definidas as funções que iriam calcular a dinâmica descrita acima, sendo que cada parte foi feita numa função separada (p.e. a dinâmica paralela e a perpendicular). Com isso polarizações e velocidades iniciais foram atribuidas, tal que as mesmas começassem em equilibrio.

 
<source lang="c">
double msd_analitico (double t)
{
return (G/(M2*TAUB*TAU)) * (t - (1.0-exp(-TAUB*t))/TAUB) + 2*Q*K*t;
}

/************************************************************************************************
* *
************************************************************************************************/
void set_gsl_rng(void)
{
#ifdef DEBUG
rseed=0;
#else
rseed=time(NULL);
#endif

gsl_rng_env_setup();
T = gsl_rng_default;
rand_vec = gsl_rng_alloc (T);
gsl_rng_set (rand_vec, rseed);

return;
}

/************************************************************************************************
* *
************************************************************************************************/
void open_output_files(void)
{

sprintf(filename,"corr_K_%.2e_N_%d.txt",K,N);
fcorr = fopen(filename,"w");
fprintf(fcorr,"#Delta t\t vcorr\n");

sprintf(filename,"traj_K_%.2e.txt",K);
ftraj = fopen(filename,"w");
fprintf(ftraj,"#t\t x\t y\n");

sprintf(filename,"r2_K_%.2e_N_%d.txt",K,N);
fr2 = fopen(filename,"w");
fprintf(fr2,"#Delta t\t r2\t r2_exact\t vcorr\n");

return;
}
</source> 

Dinâmicas (Euler-Maruyama[https://en.wikipedia.org/wiki/Euler%E2%80%93Maruyama_method]):
 
<source lang="c">
/************************************************************************************************
* *
************************************************************************************************/
void dynamics_r_parallel(int tstep)
{
/**********************************************************************
* dv_par/dt = -gamma/m v_par + (1/m) xi(t)
*
* < xi(t2)x(t1) > = G delta (t2-t1)
*
**********************************************************************/
double dr_par_1k,vpar_1k,vpar_2k,eta_par;

// Wiener process propto sqrt(dt)
eta_par = XI_PAR_SQRT_DT_DIV_M * gsl_ran_gaussian_ziggurat(rand_vec,1.0);

// Langevin eq. in par direction // EXP_VPAR = exp(-DT_DIV_2*GAMMA/M);
vpar_1k = vpar[tstep-1]*EXP_VPAR; dr_par_1k = DT_DIV_2*vpar_1k; //rpar_1k = rpar[i-1] + DT_DIV_2*vpar_1k;
vpar_2k = vpar_1k + eta_par;
vpar[tstep] = vpar_2k*EXP_VPAR; dr_par[tstep] = dr_par_1k + DT_DIV_2*vpar_2k; //rpar[i] = rpar_1k + DT_DIV_2*vpar_2k;

return;
}

/************************************************************************************************
* *
************************************************************************************************/
void dynamics_r_perpendicular(int tstep)
{
/**
// Overdamped langevin eq. in perp direction

dr_perp[i] = eta_perp; //rperp[i] += eta_perp;
rB[i] += eta_perp;
vperp[i] = INV_DT*eta_perp;
*/

dr_perp[tstep] = SQRT_Q*sin(d_theta[tstep]);
//dr_perp[tstep] = SQRT_Q*d_theta[tstep]; // Does not change aparently

return;
}

/************************************************************************************************
* *
************************************************************************************************/
double wrap_angle(double theta)
{
return theta - two_PI*floor(theta/two_PI);
}

/************************************************************************************************
* *
************************************************************************************************/
void dynamics_theta(int tstep)
{
// Overdamped langevin eq. in THETA

double eta_theta;

eta_theta = XI_THETA_SQRT_DT * gsl_ran_gaussian_ziggurat(rand_vec,1.0);

d_theta[tstep] = eta_theta;
theta[tstep] = theta[tstep-1] + eta_theta;
theta[tstep] = wrap_angle(theta[tstep]);

return;
}

/************************************************************************************************
* *
************************************************************************************************/
void dynamics(int tstep)
{
double sin_alpha,cos_alpha;
double dr,dr2;
double cos_theta_old,sin_theta_old,cos_theta_dr,sin_theta_dr;

dynamics_r_parallel(tstep);

cos_theta_old = cos(theta[tstep-1]);
sin_theta_old = sin(theta[tstep-1]);

dynamics_theta(tstep);

vpar[tstep] *= cos(d_theta[tstep]);

dynamics_r_perpendicular(tstep);

//Find total displacement
dr2 = dr_par[tstep]*dr_par[tstep] + dr_perp[tstep]*dr_perp[tstep];
dr = sqrt(dr2);

// Find displacement direction relative to original (old) direction (orientation)
if (dr != 0.0)
{
sin_alpha = dr_perp[tstep]/dr;
cos_alpha = dr_par[tstep] /dr;
}
else // Maintain direction (alpha = 0.0)
{
//fprintf(stderr,"=====>>> dr == 0.0\n");fflush(stderr);
sin_alpha = 0.0;
cos_alpha = 1.0;
}

// Find displacement angle relative to x-axis
// + is the counterclockwise direction: (theta_old + alpha)
cos_theta_dr = cos_theta_old*cos_alpha - sin_alpha*sin_theta_old;
sin_theta_dr = sin_theta_old*cos_alpha + sin_alpha*cos_theta_old;

dx[tstep] = dr*cos_theta_dr;
dy[tstep] = dr*sin_theta_dr;

rx[tstep] = rx[tstep-1] + dx[tstep];
ry[tstep] = ry[tstep-1] + dy[tstep];

// Find velocity
v2[tstep] = vpar[tstep]*vpar[tstep];

return;

}
</source> 

Por fim, colocou-se a ordem de execução das funções na ordem correta e calculou-se os parâmetros termodinâmicos relevantes.
 
<source lang="c">

int main(int argc, char *argv[])
{
int npart,k;
unsigned long int t_counter;

if (TRANSIENT_STEPS > MAX_STEPS)
{
printf("Number of transient steps must be smaller than MAX_STEPS\n");
exit(1);
}
if (argc != 3)
{
printf("Program needs 2 arguments:\n 1) Value of N,\n 2) Value of K.\n");
exit(1);
}
N=atoi(argv[1]);
K=atof(argv[2]);

fprintf(stderr,"Setting constants... "); fflush(stderr);
set_constants();
fprintf(stderr,"Done\n"); fflush(stderr);
printf("Allocating pointers...\n");fflush(stdout);
allocate_pointers(MAX_STEPS);
set_correlation_vecs();
printf("Done: Allocating pointers\n");fflush(stdout);
set_gsl_rng();
printf("Done: Seting initial conditions.\n");fflush(stdout);
open_output_files();
printf("Done: Opening files.\n");fflush(stdout);

for(npart=0; npart<N; ++npart)
{
set_initial_conditions(MAX_STEPS);
for(t_counter=1;t_counter<TRANSIENT/DT;++t_counter)
dynamics(t_counter);

--t_counter; // == TRANSIENT_STEPS - 1 // fprintf(stderr,"Done transient.\n");fflush(stderr);
reset_initial_condition(t_counter);

//fprintf(stderr,"Done transient\n"); fflush(stderr);

for(t_counter=1;t_counter<MAX_STEPS;++t_counter)
dynamics(t_counter);

if (npart==0)
{
for(k=0;k<MAX_STEPS;++k)
fprintf(ftraj,"%f %e %e\n", k*DT, rx[k], ry[k]);
fclose(ftraj);
}
if (npart%10==0)
{
printf("Done particle %d dynamics. Calculating correlations... ", npart);
fflush(stdout);
}
calculate_correlations(rx,ry,vpar);
if (npart%10==0) printf("Done.\n");
}

print_correlations();
fclose(fcorr);
fclose(fr2);

//for(i=0;i<NUM_CORR_MEASURES;++i) printf("%ld\n",corr_time[i]);

printf("Freeing pointers\n");fflush(stdout);
free_pointers(MAX_STEPS);

return 0;
}

/************************************************************************************************
* *
************************************************************************************************/
void set_correlation_vecs()
{
int i;

for(i=0;i<NUM_CORR_MEASURES;++i)
{
r2[i] = 0.0;
v2_ave[i] = 0.0;
v_ave[i] = 0.0;
corr_v_ave[i] = 0.0;
}

for(i=0;i<NUM_VCORR_MEASURES;++i)
corr_v[i] = 0.0;

return;

}
/************************************************************************************************
* *
************************************************************************************************/
void calculate_correlations(double *rx, double *ry, double *v)
{
double c_v, dr2, dx, dy, v_ave_temp,temp;
int i,j,k,num_points;

for(k=0; k<NUM_VCORR_MEASURES;++k)
{
i=k*STEP_SIZE_VCORR; //vcorr_time[k];
c_v = 0.0;
num_points = 0;
for (j=0; j < MAX_STEPS - i; ++j)
{
c_v += v[j+i] * v[j];
++num_points;
}
corr_v[k] += c_v/num_points;
//fprintf(stderr," Done corr delta t %f\n",i*DT);fflush(stderr);
}

//printf("Done 1\n");fflush(stdout);
//for (i=0; i < MAX_STEPS; ++i)
for(k=0; k<NUM_CORR_MEASURES;++k)
{
i=corr_time[k];

dr2 = 0.0;
v_ave_temp = 0.0;
v2_ave[k] += v[i]*v[i];

num_points = 0;
for (j=0; j < MAX_STEPS - i; ++j)
{
dx = rx[j+i]-rx[j];
dy = ry[j+i]-ry[j];
temp = dx*dx + dy*dy;
dr2 += temp; // dx*dx + dy*dy;
v_ave_temp += sqrt(temp)/(i*DT);

++num_points;
}
r2[k] += dr2/num_points;
v_ave[k] += v_ave_temp/num_points;

//fprintf(stderr," Done corr delta t %f\n",i*DT);fflush(stderr);
}

//printf("Done 2\n");fflush(stdout);

/*
* for(k=0; k<NUM_VCORR_MEASURES;++k)
{
i=corr_time[k];

c_v = 0.0;

num_points = 0;
for (j=0; j < MAX_STEPS - i; ++j)
{
c_v += v_ave[j+i] * v_ave[j];
++num_points;
}
corr_v_ave[k] += c_v/num_points;

//fprintf(stderr," Done corr delta t %f\n",i*DT);fflush(stderr);
}
*/
//printf("Done 3\n");fflush(stdout);

return;
}
/************************************************************************************************
* *
************************************************************************************************/
void print_correlations(void)
{
int i,k;
double t;

//for(i=0;i<MAX_STEPS;++i)
for(k=0;k<NUM_VCORR_MEASURES;++k)
{
i=k*STEP_SIZE_VCORR;//vcorr_time[k];
t=i*DT;
corr_v[k] /= N;
fprintf(fcorr,"%f %e\n", t, corr_v[k]);
}

for(k=0;k<NUM_CORR_MEASURES;++k)
{
i=corr_time[k];
t=i*DT;
r2[k] /= N;
v2_ave[k] /= N;
v_ave[k] /= N;
corr_v_ave[k] /= N;
//fprintf(fr2,"%f %e %e %e %e %e\n", t, r2[k], msd_analitico (t), v2_ave[k], v_ave[k],corr_v_ave[k]);
fprintf(fr2,"%f %e %e %e %e\n", t, r2[k], msd_analitico (t), v2_ave[k], v_ave[k]);
}

return;
}
/************************************************************************************************
* *
************************************************************************************************/
void set_initial_conditions(int N)
{
int i;

for (i=0;i<N;++i)
{
rx[i] = 0.0;
ry[i] = 0.0;
dr_par[i] = 0.0;
dr_perp[i] = 0.0;
vpar[i] = 0.0; //(G/(2*M2*(GAMMA_M+K))) * gsl_ran_gaussian_ziggurat(rand_vec,1.0);
v2[i] = 0.0; //vpar[i]*vpar[i]+vperp[i]*vperp[i];
theta[i] = 0.0; //two_PI*gsl_rng_uniform(rand_vec);
}

vpar[0] = (G/(2*M2*(GAMMA_M+K))) * gsl_ran_gaussian_ziggurat(rand_vec,1.0); //raiz, nao?!?!

//Isotropic Initial condition
//Choose a rand_vec with direction uniformly distributed in [0,2pi)
theta[0] = two_PI*gsl_rng_uniform(rand_vec);

return;
}
/************************************************************************************************
* *
************************************************************************************************/
void reset_initial_condition(int t_counter)
{
rx[0] = 0.0; //rx[t_counter];
ry[0] = 0.0; //ry[t_counter];
dr_par[0] = dr_par[t_counter];
dr_perp[0] = dr_perp[t_counter];
theta[0] = theta[t_counter];
vpar[0] = vpar[t_counter];
v2[0] = v2[t_counter];

return;
}

/************************************************************************************************
* *
************************************************************************************************/
void allocate_pointers(int n)
{
rx = create_double_pointer(n);
ry = create_double_pointer(n);
dx = create_double_pointer(n);
dy = create_double_pointer(n);
dr_par = create_double_pointer(n);
dr_perp = create_double_pointer(n);
vpar = create_double_pointer(n);
v2 = create_double_pointer(n);
theta = create_double_pointer(n);
d_theta = create_double_pointer(n);

r2 = create_double_pointer(NUM_CORR_MEASURES);
v2_ave = create_double_pointer(NUM_CORR_MEASURES);
corr_v = create_double_pointer(NUM_VCORR_MEASURES);

v_ave = create_double_pointer(NUM_CORR_MEASURES);
corr_v_ave = create_double_pointer(NUM_CORR_MEASURES);

return;
}

/************************************************************************************************
* *
************************************************************************************************/
void free_pointers(int n)
{

free(rx);
free(ry);
free(dx);
free(dy);
free(dr_par);
free(dr_perp);
free(vpar);
free(corr_v);
free(theta);
free(d_theta);

free(r2);
free(v2);
free(v2_ave);

free(v_ave);
free(corr_v_ave);

gsl_rng_free (rand_vec);

return;

}
/*********************************************************************
*** Time Table ***
*** The total number of time steps and the number of measures are ***
*** specified. ***
*********************************************************************/
void create_time_table(long *t1, int total_time, int measures)
{
unsigned long i,k;
double temp;

t1[0] = 0;
temp = pow((double) total_time,1.0/(measures-1));
k=0;
for (i=1; i<measures; ++i)
{
t1[i] = (int) pow(temp, (double) i);
if (t1[i]<=k) t1[i]=k+1;
k = t1[i];
}
t1[--i] = total_time-2;

return;
}
</source>

Grupo - Solução Exata para Movimento Anisotrópico de Ornstein-Uhlenbeck

2020-01-10T02:21:08Z

Almeidajuliano:

Aluno: Juliano Almeida Machado - Física
==Introdução==

Movimento coletivo ordenado de diversos organismos é um fenômeno que pode ser observado em diversas situações, como em movimento de cardumes, migração de pássaros, migração celular, etc..

Aquilo que é mais curioso destes tipos de movimento é que o mesmo surge através de movimentos aleatórios dos indivíduos que o compõem, ou seja, o movimento individual e aleatório de um grupo de células quando postos a interagir entre si leva a um movimento ordenado em curtas escalas de tempo.

Um dos modelos anteriores a este, o modelo de ''Vicsek'' [https://en.wikipedia.org/wiki/Vicsek_model], visa descrever o movimento coletivo ordenado de partículas “inteligentes”, isto é, partículas que podem observar as outras partículas na vizinhança e mover-se de acordo com a média do movimento do grupo. Um bom exemplo disso é o movimento de pássaros ou peixes, da onde movimentos ordenados podem emergir de movimentos aleatórios dos indivíduos autopropulsionados. Curiosamente, mesmo que estes sistemas estejam longe de um estado de equilíbrio, os mesmos apresentam características típicas de sistemas no equilíbrio, como transições de fase bem definidas.

Contudo, partículas “não inteligentes” já haviam demonstrado o mesmo tipo de movimento coletivo experimentalmente, com as mesmas mudanças de fase bem definidas. É para explicar isso, que o modelo de ''Szabó'' ([[Grupo - Modelo de Szabó]]) foi criado. Nele, o mesmo fenômeno é descrito utilizando somente as forças de interação intercelulares já que isso é algo que independe da capacidade da célula de poder observar a média do movimento da vizinhança, como em ''Vicsek''.

Essas forças são repulsivas para distâncias (<math>d_{ij} = |r_{i}-r_{j}|</math>) menores que <math>R_{eq}</math>, sendo <math>F_{rep}</math> o seu máximo para esse caso, atrativa para <math>R_{eq} \le d_{ij} \le R_{o}</math>, sendo <math>F_{adh}</math> o máximo para esse caso, e zero para <math>d_{ij}</math> maior que um certo <math>R_{o}</math>.

<math>
\textbf{F}(\textbf{r}_{i}, \textbf{r}_{j} ) = \textbf{e}_{ij}\begin{cases}
F_{rep}\frac{d_{ij} - R_{eq}}{R_{eq}},\quad d_{ij} < R_{eq} \\
F_{adh}\frac{d_{ij} - R_{eq}}{R_{0} - Re_{eq}}, \; \; R_{eq} \le d_{ij} \le R_{0}\\
0, \quad \quad \quad \quad \quad \; d_{ij} > R_{0}\\
\end{cases}
</math>

[[Arquivo:Forcaszabo.png|frame|center|Figura 1: Esquematização da força intercelular de Szabó.]]

Essas forças levam às seguintes equações de movimento:

<math>\frac{dr_{i}(t)}{dt} = v_{o}n_{i}(t) + \mu \sum_{j=1}^N F(r_{i}, r{j})</math>

Onde <math>N</math> é o número de partículas e <math>\sum_{j=1}^N F(r_{i}, r{j})</math> é o somatório da todas as interações sobre uma partícula ''i''.

<math> \frac{d_{i}^{n}(t)}{dt} = \frac{1}{\tau}arcsin[(n_{i}(t)\frac{v_{i}(t)}{|v_{i}(t)|})e_{z}] + \xi </math>

Onde <math>\tau</math> é o tempo de relaxamento, <math>n_{i}(t)</math> é o vetor na direção de polarização, tal que <math>n_{i}(t) = (cos(\theta_{i}^{n}(t)), sen(\theta_{i}^{n}(t)))</math>, <math>e_{z}</math> é o vetor normal, e <math>\xi</math> é um ruído branco.

Neste modelo, o movimento bidimensional das células é descrito por uma dinâmica superamortecida. Na mesma, cada célula tenta manter uma velocidade de magnitude <math>v_{o}</math> na direção do vetor <math>n_{i}(t)</math>, sendo que a direção desse vetor é descrita pelo ângulo <math>\theta_{i}^{n}(t)</math> que tenta relaxar na direção da velocidade enquanto sofre ruído gaussiano.

Embora descreva bem este movimento, o modelo de Szabó fornece somente soluções numéricas para o problema. Para obtermos uma solução analítica, o problema precisa ser descrito de outra forma. É com esse objetivo que o modelo de Rita et al. foi feito.

[[Arquivo:Szaboilus.png|frame|center|Figura 2: Ilustração do resultado obtido por Szabó.]]

==Modelo==
Neste modelo, em vez de ter a relaxação da velocidade na direção da polarização, a mesma é projetada sobre a polarização, enquanto ela muda com ruído. Essa direção de polarização interna pode representar, p.e., a polarização dos filamentos de actina na direção preferencial em migração celular.

Com isso pode-se observar duas transições, ou seja, uma transição de um estado desordenado para um estado ordenado e para um estado desordenado de novo. Essas duas transições representam a aleatoriedade do movimento inicial das células, depois a escolha de uma direção preferencial e por último a variação desta direção preferencial.

Contudo, neste modelo deve-se assumir que, para a partícula, seu movimento é anisotrópico, ou seja, apresenta assimetria nas dinâmicas de movimento em direções diferentes, uma na direção paralela à direção de polarização e outra perpendicular.

===Movimento Paralelo===
Na direção paralela à polarização, temos o sistema agindo numa dinâmica semelhante a de ''Langevin'' (<math>m\mathbf{a} = m\frac{d\mathbf{v}}{dt} = F(\mathbf{x}) - \beta \mathbf{v} + \eta(t)</math>)[https://en.wikipedia.org/wiki/Langevin_equation], mais especificamente, um processo de ''Orstein-Uhlenbeck'' amortecido. Nesta direção a velocidade (paralela) é bem definida. Tanto para o movimento paralelo, quanto para o movimento perpendicular, <math>p_{i}(t) = (cos(\theta_{i}(t)), sen(\theta_{i}(t)))</math> é o vetor que indica a direção de polarização e <math>n_{i}(t) = (sen(\theta_{i}(t)), -cos(\theta_{i}(t)))</math> é o vetor perpendicular a <math>p_{i}(t)</math>.

<math>\frac{dv_{||}}{dt} = -\frac{\gamma}{m}v_{||} + \frac{\xi_{||}(t)}{m} </math>
<math>v_{||}(t+\Delta t)p(t+\Delta t) = [(1-\gamma\Delta t)v_{||}(t) + \int_{t}^{t+\Delta t} \xi_{||}(t)dt](p_{i}(t)p_{i}(t+\Delta t))p_{i}(t+\Delta t) </math>

Onde <math>v_{||}</math> é a velocidade paralela, <math>\gamma</math> é o coeficiente de atrito e <math>\xi_{||}</math> é o ruído gaussiano paralelo.

A velocidade diminui por causa do atrito com o substrato, já que pode-se observar que <math>(1-\gamma\Delta t)</math> tende a <math>e^{-\gamma\Delta t}</math> se <math>\Delta t</math> for pequeno, e diminui ou aumenta por causa do termo estocástico <math>\int_{t}^{t+\Delta t} \xi_{||}(t)dt = W(t)</math>. A primeira equação assume que durante o intervalo <math>\Delta t</math>, a direção de polarização é constante. Somente em <math>t+\Delta t</math>, a velocidade na direção <math>p_{i}(t)</math> é projetada sobre a nova direção de polarização <math>p_{i}(t+\Delta t)</math>.

===Movimento Perpendicular===
Na direção perpendicular à polarização, o sistema segue um processo de ''Wiener'' (movimento aleatório) por causa das flutuações na direção de polarização.

<math>\frac{dr_{\perp}}{dt} = n_{i}(t)\xi_{\perp}(t) </math>
<math>[r_{\perp}(t+\Delta t) - r_{\perp}(t)]n_{i}(t) = n_{i}(t)\int_{t}^{t+\Delta t} \xi_{\perp}(t)dt </math>

Onde <math>\xi_{\perp}</math> é o ruído gaussiano na direção perpendicular. Pode-se observar que <math>\int_{t}^{t+\Delta t} \xi_{\perp}(t)dt = W(t) </math>. Este movimento, por ser puramente um processo de ''Wiener'', é melhor descrito pelo deslocamento do que pela velocidade.

===Ângulo e Ruídos===
A direção de polarização é definida pelo ângulo teta que também apresenta ruído branco. Esse ruído <math>\beta(t)</math> é relacionado à <math> \xi_{\perp}(t)</math> através de <math>\xi_{perp}(t) = \sqrt q\beta(t)</math>. É essa mudança na direção de polarização que causa o deslocamento aleatório perpendicular.

Embora o ruído no ângulo e o ruído na direção perpendicular estejam relacionados, os mesmos são independentes do ruído na direção paralela. Abaixo pode-se ver a correlação dos mesmos. Os três são ruídos gaussianos, então a média deles é zero.

<math> <\xi_{||}(t)\xi_{||}(t')> = g\delta(t-t') </math>
<math> <\beta(t)\beta(t')> = 2k\delta(t-t') </math>
<math> <\xi_{\perp}(t)\xi_{\perp}(t')> = 2qk\delta(t-t') </math>

Onde <math>k</math> tem unidades de 1/tempo, <math>\sqrt q</math> tem unidades de comprimento, e tanto <math>g</math> quanto <math>qk</math> tem unidades de [comprimento]^2 / tempo.

Em resumo, este modelo considera uma partícula com dois graus de liberdades espaciais e um grau de liberdade interna, que considera a polarização interna. Esse grau de liberdade interna quebra a simetria espacial de tal forma que a dinâmica na direção da polarização é semelhante a ''Langevin'' e na direção perpendicular a dinâmica segue um processo de ''Wiener''. Existem duas fontes independentes de ruído: Uma age na dinâmica na direção de polarização e a outra age na polarização em si. Além da força de dissipação, a mudança na polarização age como mais um termo na perda de energia cinética.

[[Arquivo:Wforce.gif|frame|350x350px|center|Figura 3: Simulação do modelo de Rita et al. A seta azul representa a direção de polarização e a seta verde a velocidade.]]

===Mean Square Displacement(MSD)===
O MSD é uma quantificação do movimento líquido das partículas. O mesmo foi utilizado para medir a aglomeração e a transição de fase das mesmas. Para isso foi feita uma média do MSD de 100 trajetórias sendo que o Método das Janelas foi feito para cada uma das mesmas. Através do resultados analíticos para as equações de movimento pode-se obter um resultado analítico para o MSD também:

<math><r^{2}> = \frac{4DP}{1-S}[\tau - (1-S)(1-e^{-\tau})]</math>

Onde <math>P</math> é o tempo de persistência do movimento em um estado cinético, <math>S</math> é uma fração deste tempo, relacionado a quanto tempo o movimento persiste no primeiro estado cinético e <math>D =\frac{g}{2m^{2}\tau\tau_{b}}</math>, sendo que <math>\tau = \frac{\gamma}{m} + k </math> e <math>\tau_{b} = \frac{\gamma}{m} + 2k </math>.

Para testar este modelo resolveu-se numericamente as equações diferenciais de movimento com um algoritmo simples de Euler-Maruyama e comparadas com os resultados analíticos do MSD para diferentes valores do parâmetro <math>k</math>.

O resultado obtido consta abaixo. Nele os pontos são os valores obtidos pela simulação e as linhas são as previsões da solução analítica para o MSD.

[[Arquivo:Msd_v2.png|frame|center|Figura 4: MSD analítico para 100 trajetórias]]

===O Programa===
O programa criado para gerar simulações do modelo descrito e calcular o MSD do mesmo foi feito em C. O mesmo começa com uma inialização de todas as variáveis necessárias e a definição das mesmas de acordo com aquilo que foi descrito acima sendo os parâmetros iniciais:<math>k, q, g</math> e <math> \gamma</math> e os parâmetros termodinâmicos <math> S,P </math> e <math>D</math> que descrevem o comportamento macroscópico da célula.

 
<source lang="c">
#define N 100

#define GAMMA 1.0 //3.0
#define M 1.0 //5.0
#define G 1.0 //7.0
#define K 0.00 // 0.01 //0.001

#define T_MAX (1e2) // 1e2 //5e1 //1e2 //5e2
#define DT (1e-4) // 1e-5//1e-5 //1e-6
#define TRANSIENT (2e1)

#define NUM_CORR_MEASURES 200
#define NUM_VCORR_MEASURES 200
#define T_MAX_VCORR 20

unsigned long int rseed;

double *rx,*dr_par,*vpar,*dx,*dy;
double *ry,*dr_perp,*vperp;
double *v2, *corr_v, *r2, *v2_ave,*v_ave, *corr_v_ave;

double K;
int N, MAX_STEPS, TRANSIENT_STEPS, STEP_SIZE_VCORR;
long corr_time[NUM_CORR_MEASURES];//, vcorr_time[NUM_VCORR_MEASURES];

double TAUB, GAMMA_M, TAU2,GAMMA_M2,V02;

double DT_DIV_2,M2, TAU, Q, SQRT_Q;
double EXP_VPAR, SQRT_DT;
double XI_PAR,XI_PAR_SQRT_DT,XI_PAR_SQRT_DT_DIV_M;
double XI_THETA,XI_THETA_SQRT_DT; // Used if langevin theta dynamics

double *theta, *d_theta;
double two_PI = 2*M_PI;
const gsl_rng_type * T;
gsl_rng * rand_vec;

char filename[200];
FILE *fcorr, *ftraj, *fr2;

/************************************************************************************************
* *
************************************************************************************************/
void allocate_pointers(int n);
void free_pointers(int n);
void measure(double *rx, double *ry,double *rx0, double *ry0, double *v2, unsigned long int t);
void set_initial_conditions(int N);
void set_correlation_vecs();
void calculate_correlations(double *rx, double *ry, double *v);
void print_correlations(void);
void reset_initial_condition(int t_counter);
void set_gsl_rng(void);
double wrap_angle(double theta);
void create_time_table(long *t1, int total_time, int measures);

/************************************************************************************************
* *
************************************************************************************************/
</source> 

Inicialização:
 
<source lang="c">
void set_constants(void)
{
MAX_STEPS = (int)(T_MAX/DT) + 1;
TRANSIENT_STEPS = TRANSIENT/DT;

DT_DIV_2 = 0.5*DT;
SQRT_DT = sqrt(DT);

M2 = M*M;
TAU = GAMMA/M + K;
TAU2 = TAU*TAU;
TAUB = GAMMA/M + 2*K;
GAMMA_M = GAMMA/M;
GAMMA_M2 = GAMMA_M*GAMMA_M;

XI_THETA = sqrt(2.0*K);
XI_THETA_SQRT_DT = XI_THETA * SQRT_DT; // Used if langevin theta dynamics

XI_PAR = sqrt(G); //sqrt(2.0*(GAMMA+K*M)*kT);
XI_PAR_SQRT_DT = XI_PAR * SQRT_DT;
XI_PAR_SQRT_DT_DIV_M = XI_PAR * SQRT_DT/M;

EXP_VPAR = exp(-DT_DIV_2*GAMMA/M); //For langevin algorithm

V02 = G/(2*M2*TAU); //G/(2*M2*(GAMMA/M + K))

Q = V02/TAU2; //(kT/M)/TAU2;//v2/(TAU*TAU);
SQRT_Q = sqrt(Q);

create_time_table(corr_time, MAX_STEPS, NUM_CORR_MEASURES);
//create_time_table(vcorr_time, T_MAX_VCORR/DT, NUM_VCORR_MEASURES);
STEP_SIZE_VCORR = T_MAX_VCORR/(NUM_VCORR_MEASURES*DT);

return;
}
</source> 

Em seguida foram definidas as funções que iriam calcular a dinâmica descrita acima, sendo que cada parte foi feita numa função separada (p.e. a dinâmica paralela e a perpendicular). Com isso polarizações e velocidades iniciais foram atribuidas, tal que as mesmas começassem em equilibrio.

 
<source lang="c">
double msd_analitico (double t)
{
return (G/(M2*TAUB*TAU)) * (t - (1.0-exp(-TAUB*t))/TAUB) + 2*Q*K*t;
}

/************************************************************************************************
* *
************************************************************************************************/
void set_gsl_rng(void)
{
#ifdef DEBUG
rseed=0;
#else
rseed=time(NULL);
#endif

gsl_rng_env_setup();
T = gsl_rng_default;
rand_vec = gsl_rng_alloc (T);
gsl_rng_set (rand_vec, rseed);

return;
}

/************************************************************************************************
* *
************************************************************************************************/
void open_output_files(void)
{

sprintf(filename,"corr_K_%.2e_N_%d.txt",K,N);
fcorr = fopen(filename,"w");
fprintf(fcorr,"#Delta t\t vcorr\n");

sprintf(filename,"traj_K_%.2e.txt",K);
ftraj = fopen(filename,"w");
fprintf(ftraj,"#t\t x\t y\n");

sprintf(filename,"r2_K_%.2e_N_%d.txt",K,N);
fr2 = fopen(filename,"w");
fprintf(fr2,"#Delta t\t r2\t r2_exact\t vcorr\n");

return;
}
</source> 

Dinâmicas (Euler-Maruyama[https://en.wikipedia.org/wiki/Euler%E2%80%93Maruyama_method]):
 
<source lang="c">
/************************************************************************************************
* *
************************************************************************************************/
void dynamics_r_parallel(int tstep)
{
/**********************************************************************
* dv_par/dt = -gamma/m v_par + (1/m) xi(t)
*
* < xi(t2)x(t1) > = G delta (t2-t1)
*
**********************************************************************/
double dr_par_1k,vpar_1k,vpar_2k,eta_par;

// Wiener process propto sqrt(dt)
eta_par = XI_PAR_SQRT_DT_DIV_M * gsl_ran_gaussian_ziggurat(rand_vec,1.0);

// Langevin eq. in par direction // EXP_VPAR = exp(-DT_DIV_2*GAMMA/M);
vpar_1k = vpar[tstep-1]*EXP_VPAR; dr_par_1k = DT_DIV_2*vpar_1k; //rpar_1k = rpar[i-1] + DT_DIV_2*vpar_1k;
vpar_2k = vpar_1k + eta_par;
vpar[tstep] = vpar_2k*EXP_VPAR; dr_par[tstep] = dr_par_1k + DT_DIV_2*vpar_2k; //rpar[i] = rpar_1k + DT_DIV_2*vpar_2k;

return;
}

/************************************************************************************************
* *
************************************************************************************************/
void dynamics_r_perpendicular(int tstep)
{
/**
// Overdamped langevin eq. in perp direction

dr_perp[i] = eta_perp; //rperp[i] += eta_perp;
rB[i] += eta_perp;
vperp[i] = INV_DT*eta_perp;
*/

dr_perp[tstep] = SQRT_Q*sin(d_theta[tstep]);
//dr_perp[tstep] = SQRT_Q*d_theta[tstep]; // Does not change aparently

return;
}

/************************************************************************************************
* *
************************************************************************************************/
double wrap_angle(double theta)
{
return theta - two_PI*floor(theta/two_PI);
}

/************************************************************************************************
* *
************************************************************************************************/
void dynamics_theta(int tstep)
{
// Overdamped langevin eq. in THETA

double eta_theta;

eta_theta = XI_THETA_SQRT_DT * gsl_ran_gaussian_ziggurat(rand_vec,1.0);

d_theta[tstep] = eta_theta;
theta[tstep] = theta[tstep-1] + eta_theta;
theta[tstep] = wrap_angle(theta[tstep]);

return;
}

/************************************************************************************************
* *
************************************************************************************************/
void dynamics(int tstep)
{
double sin_alpha,cos_alpha;
double dr,dr2;
double cos_theta_old,sin_theta_old,cos_theta_dr,sin_theta_dr;

dynamics_r_parallel(tstep);

cos_theta_old = cos(theta[tstep-1]);
sin_theta_old = sin(theta[tstep-1]);

dynamics_theta(tstep);

vpar[tstep] *= cos(d_theta[tstep]);

dynamics_r_perpendicular(tstep);

//Find total displacement
dr2 = dr_par[tstep]*dr_par[tstep] + dr_perp[tstep]*dr_perp[tstep];
dr = sqrt(dr2);

// Find displacement direction relative to original (old) direction (orientation)
if (dr != 0.0)
{
sin_alpha = dr_perp[tstep]/dr;
cos_alpha = dr_par[tstep] /dr;
}
else // Maintain direction (alpha = 0.0)
{
//fprintf(stderr,"=====>>> dr == 0.0\n");fflush(stderr);
sin_alpha = 0.0;
cos_alpha = 1.0;
}

// Find displacement angle relative to x-axis
// + is the counterclockwise direction: (theta_old + alpha)
cos_theta_dr = cos_theta_old*cos_alpha - sin_alpha*sin_theta_old;
sin_theta_dr = sin_theta_old*cos_alpha + sin_alpha*cos_theta_old;

dx[tstep] = dr*cos_theta_dr;
dy[tstep] = dr*sin_theta_dr;

rx[tstep] = rx[tstep-1] + dx[tstep];
ry[tstep] = ry[tstep-1] + dy[tstep];

// Find velocity
v2[tstep] = vpar[tstep]*vpar[tstep];

return;

}
</source> 

Por fim, colocou-se a ordem de execução das funções na ordem correta e calculou-se os parâmetros termodinâmicos relevantes.
 
<source lang="c">

int main(int argc, char *argv[])
{
int npart,k;
unsigned long int t_counter;

if (TRANSIENT_STEPS > MAX_STEPS)
{
printf("Number of transient steps must be smaller than MAX_STEPS\n");
exit(1);
}
if (argc != 3)
{
printf("Program needs 2 arguments:\n 1) Value of N,\n 2) Value of K.\n");
exit(1);
}
N=atoi(argv[1]);
K=atof(argv[2]);

fprintf(stderr,"Setting constants... "); fflush(stderr);
set_constants();
fprintf(stderr,"Done\n"); fflush(stderr);
printf("Allocating pointers...\n");fflush(stdout);
allocate_pointers(MAX_STEPS);
set_correlation_vecs();
printf("Done: Allocating pointers\n");fflush(stdout);
set_gsl_rng();
printf("Done: Seting initial conditions.\n");fflush(stdout);
open_output_files();
printf("Done: Opening files.\n");fflush(stdout);

for(npart=0; npart<N; ++npart)
{
set_initial_conditions(MAX_STEPS);
for(t_counter=1;t_counter<TRANSIENT/DT;++t_counter)
dynamics(t_counter);

--t_counter; // == TRANSIENT_STEPS - 1 // fprintf(stderr,"Done transient.\n");fflush(stderr);
reset_initial_condition(t_counter);

//fprintf(stderr,"Done transient\n"); fflush(stderr);

for(t_counter=1;t_counter<MAX_STEPS;++t_counter)
dynamics(t_counter);

if (npart==0)
{
for(k=0;k<MAX_STEPS;++k)
fprintf(ftraj,"%f %e %e\n", k*DT, rx[k], ry[k]);
fclose(ftraj);
}
if (npart%10==0)
{
printf("Done particle %d dynamics. Calculating correlations... ", npart);
fflush(stdout);
}
calculate_correlations(rx,ry,vpar);
if (npart%10==0) printf("Done.\n");
}

print_correlations();
fclose(fcorr);
fclose(fr2);

//for(i=0;i<NUM_CORR_MEASURES;++i) printf("%ld\n",corr_time[i]);

printf("Freeing pointers\n");fflush(stdout);
free_pointers(MAX_STEPS);

return 0;
}

/************************************************************************************************
* *
************************************************************************************************/
void set_correlation_vecs()
{
int i;

for(i=0;i<NUM_CORR_MEASURES;++i)
{
r2[i] = 0.0;
v2_ave[i] = 0.0;
v_ave[i] = 0.0;
corr_v_ave[i] = 0.0;
}

for(i=0;i<NUM_VCORR_MEASURES;++i)
corr_v[i] = 0.0;

return;

}
/************************************************************************************************
* *
************************************************************************************************/
void calculate_correlations(double *rx, double *ry, double *v)
{
double c_v, dr2, dx, dy, v_ave_temp,temp;
int i,j,k,num_points;

for(k=0; k<NUM_VCORR_MEASURES;++k)
{
i=k*STEP_SIZE_VCORR; //vcorr_time[k];
c_v = 0.0;
num_points = 0;
for (j=0; j < MAX_STEPS - i; ++j)
{
c_v += v[j+i] * v[j];
++num_points;
}
corr_v[k] += c_v/num_points;
//fprintf(stderr," Done corr delta t %f\n",i*DT);fflush(stderr);
}

//printf("Done 1\n");fflush(stdout);
//for (i=0; i < MAX_STEPS; ++i)
for(k=0; k<NUM_CORR_MEASURES;++k)
{
i=corr_time[k];

dr2 = 0.0;
v_ave_temp = 0.0;
v2_ave[k] += v[i]*v[i];

num_points = 0;
for (j=0; j < MAX_STEPS - i; ++j)
{
dx = rx[j+i]-rx[j];
dy = ry[j+i]-ry[j];
temp = dx*dx + dy*dy;
dr2 += temp; // dx*dx + dy*dy;
v_ave_temp += sqrt(temp)/(i*DT);

++num_points;
}
r2[k] += dr2/num_points;
v_ave[k] += v_ave_temp/num_points;

//fprintf(stderr," Done corr delta t %f\n",i*DT);fflush(stderr);
}

//printf("Done 2\n");fflush(stdout);

/*
* for(k=0; k<NUM_VCORR_MEASURES;++k)
{
i=corr_time[k];

c_v = 0.0;

num_points = 0;
for (j=0; j < MAX_STEPS - i; ++j)
{
c_v += v_ave[j+i] * v_ave[j];
++num_points;
}
corr_v_ave[k] += c_v/num_points;

//fprintf(stderr," Done corr delta t %f\n",i*DT);fflush(stderr);
}
*/
//printf("Done 3\n");fflush(stdout);

return;
}
/************************************************************************************************
* *
************************************************************************************************/
void print_correlations(void)
{
int i,k;
double t;

//for(i=0;i<MAX_STEPS;++i)
for(k=0;k<NUM_VCORR_MEASURES;++k)
{
i=k*STEP_SIZE_VCORR;//vcorr_time[k];
t=i*DT;
corr_v[k] /= N;
fprintf(fcorr,"%f %e\n", t, corr_v[k]);
}

for(k=0;k<NUM_CORR_MEASURES;++k)
{
i=corr_time[k];
t=i*DT;
r2[k] /= N;
v2_ave[k] /= N;
v_ave[k] /= N;
corr_v_ave[k] /= N;
//fprintf(fr2,"%f %e %e %e %e %e\n", t, r2[k], msd_analitico (t), v2_ave[k], v_ave[k],corr_v_ave[k]);
fprintf(fr2,"%f %e %e %e %e\n", t, r2[k], msd_analitico (t), v2_ave[k], v_ave[k]);
}

return;
}
/************************************************************************************************
* *
************************************************************************************************/
void set_initial_conditions(int N)
{
int i;

for (i=0;i<N;++i)
{
rx[i] = 0.0;
ry[i] = 0.0;
dr_par[i] = 0.0;
dr_perp[i] = 0.0;
vpar[i] = 0.0; //(G/(2*M2*(GAMMA_M+K))) * gsl_ran_gaussian_ziggurat(rand_vec,1.0);
v2[i] = 0.0; //vpar[i]*vpar[i]+vperp[i]*vperp[i];
theta[i] = 0.0; //two_PI*gsl_rng_uniform(rand_vec);
}

vpar[0] = (G/(2*M2*(GAMMA_M+K))) * gsl_ran_gaussian_ziggurat(rand_vec,1.0); //raiz, nao?!?!

//Isotropic Initial condition
//Choose a rand_vec with direction uniformly distributed in [0,2pi)
theta[0] = two_PI*gsl_rng_uniform(rand_vec);

return;
}
/************************************************************************************************
* *
************************************************************************************************/
void reset_initial_condition(int t_counter)
{
rx[0] = 0.0; //rx[t_counter];
ry[0] = 0.0; //ry[t_counter];
dr_par[0] = dr_par[t_counter];
dr_perp[0] = dr_perp[t_counter];
theta[0] = theta[t_counter];
vpar[0] = vpar[t_counter];
v2[0] = v2[t_counter];

return;
}

/************************************************************************************************
* *
************************************************************************************************/
void allocate_pointers(int n)
{
rx = create_double_pointer(n);
ry = create_double_pointer(n);
dx = create_double_pointer(n);
dy = create_double_pointer(n);
dr_par = create_double_pointer(n);
dr_perp = create_double_pointer(n);
vpar = create_double_pointer(n);
v2 = create_double_pointer(n);
theta = create_double_pointer(n);
d_theta = create_double_pointer(n);

r2 = create_double_pointer(NUM_CORR_MEASURES);
v2_ave = create_double_pointer(NUM_CORR_MEASURES);
corr_v = create_double_pointer(NUM_VCORR_MEASURES);

v_ave = create_double_pointer(NUM_CORR_MEASURES);
corr_v_ave = create_double_pointer(NUM_CORR_MEASURES);

return;
}

/************************************************************************************************
* *
************************************************************************************************/
void free_pointers(int n)
{

free(rx);
free(ry);
free(dx);
free(dy);
free(dr_par);
free(dr_perp);
free(vpar);
free(corr_v);
free(theta);
free(d_theta);

free(r2);
free(v2);
free(v2_ave);

free(v_ave);
free(corr_v_ave);

gsl_rng_free (rand_vec);

return;

}
/*********************************************************************
*** Time Table ***
*** The total number of time steps and the number of measures are ***
*** specified. ***
*********************************************************************/
void create_time_table(long *t1, int total_time, int measures)
{
unsigned long i,k;
double temp;

t1[0] = 0;
temp = pow((double) total_time,1.0/(measures-1));
k=0;
for (i=1; i<measures; ++i)
{
t1[i] = (int) pow(temp, (double) i);
if (t1[i]<=k) t1[i]=k+1;
k = t1[i];
}
t1[--i] = total_time-2;

return;
}
</source> 
Código escrito pelo Prof. Mendeli Vainstein.

Grupo - Solução Exata para Movimento Anisotrópico de Ornstein-Uhlenbeck

2020-01-10T02:11:27Z

Almeidajuliano:

Aluno: Juliano Almeida Machado - Física
==Introdução==

Movimento coletivo ordenado de diversos organismos é um fenômeno que pode ser observado em diversas situações, como em movimento de cardumes, migração de pássaros, migração celular, etc..

Aquilo que é mais curioso destes tipos de movimento é que o mesmo surge através de movimentos aleatórios dos indivíduos que o compõem, ou seja, o movimento individual e aleatório de um grupo de células quando postos a interagir entre si leva a um movimento ordenado em curtas escalas de tempo.

Um dos modelos anteriores a este, o modelo de ''Vicsek'' [https://en.wikipedia.org/wiki/Vicsek_model], visa descrever o movimento coletivo ordenado de partículas “inteligentes”, isto é, partículas que podem observar as outras partículas na vizinhança e mover-se de acordo com a média do movimento do grupo. Um bom exemplo disso é o movimento de pássaros ou peixes, da onde movimentos ordenados podem emergir de movimentos aleatórios dos indivíduos autopropulsionados. Curiosamente, mesmo que estes sistemas estejam longe de um estado de equilíbrio, os mesmos apresentam características típicas de sistemas no equilíbrio, como transições de fase bem definidas.

Contudo, partículas “não inteligentes” já haviam demonstrado o mesmo tipo de movimento coletivo experimentalmente, com as mesmas mudanças de fase bem definidas. É para explicar isso, que o modelo de ''Szabó'' ([[Grupo - Modelo de Szabó]]) foi criado. Nele, o mesmo fenômeno é descrito utilizando somente as forças de interação intercelulares já que isso é algo que independe da capacidade da célula de poder observar a média do movimento da vizinhança, como em ''Vicsek''.

Essas forças são repulsivas para distâncias (<math>d_{ij} = |r_{i}-r_{j}|</math>) menores que <math>R_{eq}</math>, sendo <math>F_{rep}</math> o seu máximo para esse caso, atrativa para <math>R_{eq} \le d_{ij} \le R_{o}</math>, sendo <math>F_{adh}</math> o máximo para esse caso, e zero para <math>d_{ij}</math> maior que um certo <math>R_{o}</math>.

<math>
\textbf{F}(\textbf{r}_{i}, \textbf{r}_{j} ) = \textbf{e}_{ij}\begin{cases}
F_{rep}\frac{d_{ij} - R_{eq}}{R_{eq}},\quad d_{ij} < R_{eq} \\
F_{adh}\frac{d_{ij} - R_{eq}}{R_{0} - Re_{eq}}, \; \; R_{eq} \le d_{ij} \le R_{0}\\
0, \quad \quad \quad \quad \quad \; d_{ij} > R_{0}\\
\end{cases}
</math>

[[Arquivo:Forcaszabo.png|frame|center|Figura 1: Esquematização da força intercelular de Szabó.]]

Essas forças levam às seguintes equações de movimento:

<math>\frac{dr_{i}(t)}{dt} = v_{o}n_{i}(t) + \mu \sum_{j=1}^N F(r_{i}, r{j})</math>

Onde <math>N</math> é o número de partículas e <math>\sum_{j=1}^N F(r_{i}, r{j})</math> é o somatório da todas as interações sobre uma partícula ''i''.

<math> \frac{d_{i}^{n}(t)}{dt} = \frac{1}{\tau}arcsin[(n_{i}(t)\frac{v_{i}(t)}{|v_{i}(t)|})e_{z}] + \xi </math>

Onde <math>\tau</math> é o tempo de relaxamento, <math>n_{i}(t)</math> é o vetor na direção de polarização, tal que <math>n_{i}(t) = (cos(\theta_{i}^{n}(t)), sen(\theta_{i}^{n}(t)))</math>, <math>e_{z}</math> é o vetor normal, e <math>\xi</math> é um ruído branco.

Neste modelo, o movimento bidimensional das células é descrito por uma dinâmica superamortecida. Na mesma, cada célula tenta manter uma velocidade de magnitude <math>v_{o}</math> na direção do vetor <math>n_{i}(t)</math>, sendo que a direção desse vetor é descrita pelo ângulo <math>\theta_{i}^{n}(t)</math> que tenta relaxar na direção da velocidade enquanto sofre ruído gaussiano.

Embora descreva bem este movimento, o modelo de Szabó fornece somente soluções numéricas para o problema. Para obtermos uma solução analítica, o problema precisa ser descrito de outra forma. É com esse objetivo que o modelo de Rita et al. foi feito.

[[Arquivo:Szaboilus.png|frame|center|Figura 2: Ilustração do resultado obtido por Szabó.]]

==Modelo==
Neste modelo, em vez de ter a relaxação da velocidade na direção da polarização, a mesma é projetada sobre a polarização, enquanto ela muda com ruído. Essa direção de polarização interna pode representar, p.e., a polarização dos filamentos de actina na direção preferencial em migração celular.

Com isso pode-se observar duas transições, ou seja, uma transição de um estado desordenado para um estado ordenado e para um estado desordenado de novo. Essas duas transições representam a aleatoriedade do movimento inicial das células, depois a escolha de uma direção preferencial e por último a variação desta direção preferencial.

Contudo, neste modelo deve-se assumir que, para a partícula, seu movimento é anisotrópico, ou seja, apresenta assimetria nas dinâmicas de movimento em direções diferentes, uma na direção paralela à direção de polarização e outra perpendicular.

===Movimento Paralelo===
Na direção paralela à polarização, temos o sistema agindo numa dinâmica semelhante a de ''Langevin'' (<math>m\mathbf{a} = m\frac{d\mathbf{v}}{dt} = F(\mathbf{x}) - \beta \mathbf{v} + \eta(t)</math>), mais especificamente, um processo de ''Orstein-Uhlenbeck'' amortecido. Nesta direção a velocidade (paralela) é bem definida. Tanto para o movimento paralelo, quanto para o movimento perpendicular, o vetor <math>p_{i}(t) = (cos(\theta_{i}(t)), sen(\theta_{i}(t)))</math> é o vetor que indica a direção de polarização e <math>n_{i}(t) = (sen(\theta_{i}(t)), -cos(\theta_{i}(t)))</math> é o vetor perpendicular a <math>p_{i}(t)</math>.

<math>\frac{dv_{||}}{dt} = -\frac{\gamma}{m}v_{||} + \frac{\xi_{||}(t)}{m} </math>
<math>v_{||}(t+\Delta t)p(t+\Delta t) = [(1-\gamma\Delta t)v_{||}(t) + \int_{t}^{t+\Delta t} \xi_{||}(t)dt](p_{i}(t)p_{i}(t+\Delta t))p_{i}(t+\Delta t) </math>

Onde <math>v_{||}</math> é a velocidade paralela, <math>\gamma</math> é o coeficiente de atrito e <math>\xi_{||}</math> é o ruído gaussiano paralelo.

A velocidade diminui por causa do atrito com o substrato, já que pode-se observar que <math>(1-\gamma\Delta t)</math> tende a <math>e^{-\gamma\Delta t}</math> se <math>\gamma</math> for pequeno, e diminui ou aumenta por causa do termo estocástico <math>\int_{t}^{t+\Delta t} \xi_{||}(t)dt = W(t)</math>. A primeira equação assume que durante o intervalo <math>\Delta t</math>, a direção de polarização é constante. Somente em <math>t+\Delta t</math>, a velocidade na direção <math>p_{i}(t)</math> é projetada sobre a nova direção de polarização <math>p_{i}(t+\Delta t)</math>.

===Movimento Perpendicular===
Na direção perpendicular à polarização, o sistema segue um processo de ''Wiener'' (movimento aleatório) por causa das flutuações na direção de polarização.

<math>\frac{dr_{\perp}}{dt} = n_{i}(t)\xi_{\perp}(t) </math>
<math>[r_{\perp}(t+\Delta t) - r_{\perp}(t)]n_{i}(t) = n_{i}(t)\int_{t}^{t+\Delta t} \xi_{\perp}(t)dt </math>

Onde <math>\xi_{\perp}</math> é o ruído gaussiano na direção perpendicular. Pode-se observar que <math>\int_{t}^{t+\Delta t} \xi_{\perp}(t)dt = W(t) </math>. Este movimento, por ser um processo de ''Wiener'', é melhor descrito pelo deslocamento do que pela velocidade.

===Ângulo e Ruídos===
A direção de polarização é definida pelo ângulo teta que também apresenta ruído branco. Esse ruído <math>\beta(t)</math> é relacionado à <math> \xi_{\perp}(t)</math> através de <math>\xi_{perp}(t) = \sqrt(q)\beta(t)</math>. É essa mudança na direção de polarização que causa o deslocamento aleatório perpendicular.

Embora o ruído no ângulo e o ruído na direção perpendicular estejam relacionados, os mesmos são independentes do ruído na direção paralela. Abaixo pode-se ver a correlação dos mesmos. Os três são ruídos gaussianos, então a média deles é zero.

<math> <\xi_{||}(t)\xi_{||}(t')> = g\delta(t-t') </math>
<math> <\beta(t)\beta(t')> = 2k\delta(t-t') </math>
<math> <\xi_{\perp}(t)\xi_{\perp}(t')> = 2qk\delta(t-t') </math>

Onde <math>k</math> tem unidades de 1/tempo, <math>\sqrt q</math> tem unidades de comprimento, e tanto <math>g</math> quanto <math>qk</math> tem unidades de [comprimento]^2 / tempo.

Em resumo, este modelo considera uma partícula com dois graus de liberdades espaciais e um grau de liberdade interna, que considera a polarização interna. Esse grau de liberdade interna quebra a simetria espacial de tal forma que a dinâmica na direção da polarização é semelhante a ''Langevin'' e na direção perpendicular a dinâmica segue um processo de ''Wiener''. Existem duas fontes independentes de ruído: Uma age na dinâmica na direção de polarização e a outra age na polarização em si. Além da força de dissipação, a mudança na polarização age como mais um termo na perda de energia cinética.

[[Arquivo:Wforce.gif|frame|350x350px|center|Figura 3: Simulação do modelo de Rita et al. A seta azul representa a direção de polarização e a seta verde a velocidade.]]

===Mean Square Displacement(MSD)===
O MSD é uma quantificação do movimento líquido das partículas. O mesmo foi utilizado para medir a aglomeração e a transição de fase das partículas. Para isso foi feita uma média do MSD de 100 trajetórias sendo que o Método das Janelas foi feito para cada uma das mesmas. Através do resultados analíticos para as equações de movimento pode-se obter um resultado analítico para o MSD também:

<math><r^{2}> = \frac{4DP}{1-S}[\tau - (1-S)(1-e^{-\tau})]</math>

Onde <math>P</math> é o tempo de persistência do movimento em um estado cinético, <math>S</math> é uma fração deste tempo, relacionado a quanto tempo o movimento persiste no primeiro estado cinético e <math>D =\frac{g}{2m^{2}\tau\tau_{b}}</math>, sendo que <math>\tau = \frac{\gamma}{m} + k </math> e <math>\tau_{b} = \frac{\gamma}{m} + 2k </math>.

Para testar este modelo resolveu-se numericamente as equações diferenciais de movimento com um algoritmo de Euler-Maruyama e comparadas com os resultados analíticos do MSD para diferentes valores do parâmetro <math>k</math>.

O resultado obtido consta abaixo. Nele os pontos são os valores obtidos pela simulação e as linhas são as previsões da solução analítica para o MSD.

[[Arquivo:Msd_v2.png|frame|center|Figura 4: MSD analítico para 100 trajetórias]]

===O Programa===
O programa criado para gerar simulações do modelo descrito e calcular o MSD do mesmo foi feito em C. O mesmo começa com uma inialização de todas as variáveis necessárias e a definição das mesmas de acordo com aquilo que foi descrito acima sendo os parâmetros iniciais:<math>k, q, g</math> e <math> \gamma</math> e os parâmetros termodinâmicos <math> S,P </math> e <math>D</math> que descrevem o comportamento macroscópico da célula.

 
<source lang="c">
#define N 100

#define GAMMA 1.0 //3.0
#define M 1.0 //5.0
#define G 1.0 //7.0
#define K 0.00 // 0.01 //0.001

#define T_MAX (1e2) // 1e2 //5e1 //1e2 //5e2
#define DT (1e-4) // 1e-5//1e-5 //1e-6
#define TRANSIENT (2e1)

#define NUM_CORR_MEASURES 200
#define NUM_VCORR_MEASURES 200
#define T_MAX_VCORR 20

unsigned long int rseed;

double *rx,*dr_par,*vpar,*dx,*dy;
double *ry,*dr_perp,*vperp;
double *v2, *corr_v, *r2, *v2_ave,*v_ave, *corr_v_ave;

double K;
int N, MAX_STEPS, TRANSIENT_STEPS, STEP_SIZE_VCORR;
long corr_time[NUM_CORR_MEASURES];//, vcorr_time[NUM_VCORR_MEASURES];

double TAUB, GAMMA_M, TAU2,GAMMA_M2,V02;

double DT_DIV_2,M2, TAU, Q, SQRT_Q;
double EXP_VPAR, SQRT_DT;
double XI_PAR,XI_PAR_SQRT_DT,XI_PAR_SQRT_DT_DIV_M;
double XI_THETA,XI_THETA_SQRT_DT; // Used if langevin theta dynamics

double *theta, *d_theta;
double two_PI = 2*M_PI;
const gsl_rng_type * T;
gsl_rng * rand_vec;

char filename[200];
FILE *fcorr, *ftraj, *fr2;

/************************************************************************************************
* *
************************************************************************************************/
void allocate_pointers(int n);
void free_pointers(int n);
void measure(double *rx, double *ry,double *rx0, double *ry0, double *v2, unsigned long int t);
void set_initial_conditions(int N);
void set_correlation_vecs();
void calculate_correlations(double *rx, double *ry, double *v);
void print_correlations(void);
void reset_initial_condition(int t_counter);
void set_gsl_rng(void);
double wrap_angle(double theta);
void create_time_table(long *t1, int total_time, int measures);

/************************************************************************************************
* *
************************************************************************************************/
</source> 

Inicialização:
 
<source lang="c">
void set_constants(void)
{
MAX_STEPS = (int)(T_MAX/DT) + 1;
TRANSIENT_STEPS = TRANSIENT/DT;

DT_DIV_2 = 0.5*DT;
SQRT_DT = sqrt(DT);

M2 = M*M;
TAU = GAMMA/M + K;
TAU2 = TAU*TAU;
TAUB = GAMMA/M + 2*K;
GAMMA_M = GAMMA/M;
GAMMA_M2 = GAMMA_M*GAMMA_M;

XI_THETA = sqrt(2.0*K);
XI_THETA_SQRT_DT = XI_THETA * SQRT_DT; // Used if langevin theta dynamics

XI_PAR = sqrt(G); //sqrt(2.0*(GAMMA+K*M)*kT);
XI_PAR_SQRT_DT = XI_PAR * SQRT_DT;
XI_PAR_SQRT_DT_DIV_M = XI_PAR * SQRT_DT/M;

EXP_VPAR = exp(-DT_DIV_2*GAMMA/M); //For langevin algorithm

V02 = G/(2*M2*TAU); //G/(2*M2*(GAMMA/M + K))

Q = V02/TAU2; //(kT/M)/TAU2;//v2/(TAU*TAU);
SQRT_Q = sqrt(Q);

create_time_table(corr_time, MAX_STEPS, NUM_CORR_MEASURES);
//create_time_table(vcorr_time, T_MAX_VCORR/DT, NUM_VCORR_MEASURES);
STEP_SIZE_VCORR = T_MAX_VCORR/(NUM_VCORR_MEASURES*DT);

return;
}
</source> 

Em seguida foram definidas as funções que iriam calcular a dinâmica descrita acima, sendo que cada parte foi feita numa função separada (p.e. a dinâmica paralela e a perpendicular). Com isso polarizações e velocidades iniciais foram atribuidas, tal que as mesmas começassem em equilibrio.

 
<source lang="c">
double msd_analitico (double t)
{
return (G/(M2*TAUB*TAU)) * (t - (1.0-exp(-TAUB*t))/TAUB) + 2*Q*K*t;
}

/************************************************************************************************
* *
************************************************************************************************/
void set_gsl_rng(void)
{
#ifdef DEBUG
rseed=0;
#else
rseed=time(NULL);
#endif

gsl_rng_env_setup();
T = gsl_rng_default;
rand_vec = gsl_rng_alloc (T);
gsl_rng_set (rand_vec, rseed);

return;
}

/************************************************************************************************
* *
************************************************************************************************/
void open_output_files(void)
{

sprintf(filename,"corr_K_%.2e_N_%d.txt",K,N);
fcorr = fopen(filename,"w");
fprintf(fcorr,"#Delta t\t vcorr\n");

sprintf(filename,"traj_K_%.2e.txt",K);
ftraj = fopen(filename,"w");
fprintf(ftraj,"#t\t x\t y\n");

sprintf(filename,"r2_K_%.2e_N_%d.txt",K,N);
fr2 = fopen(filename,"w");
fprintf(fr2,"#Delta t\t r2\t r2_exact\t vcorr\n");

return;
}
</source> 

Dinâmicas:
 
<source lang="c">
/************************************************************************************************
* *
************************************************************************************************/
void dynamics_r_parallel(int tstep)
{
/**********************************************************************
* dv_par/dt = -gamma/m v_par + (1/m) xi(t)
*
* < xi(t2)x(t1) > = G delta (t2-t1)
*
**********************************************************************/
double dr_par_1k,vpar_1k,vpar_2k,eta_par;

// Wiener process propto sqrt(dt)
eta_par = XI_PAR_SQRT_DT_DIV_M * gsl_ran_gaussian_ziggurat(rand_vec,1.0);

// Langevin eq. in par direction // EXP_VPAR = exp(-DT_DIV_2*GAMMA/M);
vpar_1k = vpar[tstep-1]*EXP_VPAR; dr_par_1k = DT_DIV_2*vpar_1k; //rpar_1k = rpar[i-1] + DT_DIV_2*vpar_1k;
vpar_2k = vpar_1k + eta_par;
vpar[tstep] = vpar_2k*EXP_VPAR; dr_par[tstep] = dr_par_1k + DT_DIV_2*vpar_2k; //rpar[i] = rpar_1k + DT_DIV_2*vpar_2k;

return;
}

/************************************************************************************************
* *
************************************************************************************************/
void dynamics_r_perpendicular(int tstep)
{
/**
// Overdamped langevin eq. in perp direction

dr_perp[i] = eta_perp; //rperp[i] += eta_perp;
rB[i] += eta_perp;
vperp[i] = INV_DT*eta_perp;
*/

dr_perp[tstep] = SQRT_Q*sin(d_theta[tstep]);
//dr_perp[tstep] = SQRT_Q*d_theta[tstep]; // Does not change aparently

return;
}

/************************************************************************************************
* *
************************************************************************************************/
double wrap_angle(double theta)
{
return theta - two_PI*floor(theta/two_PI);
}

/************************************************************************************************
* *
************************************************************************************************/
void dynamics_theta(int tstep)
{
// Overdamped langevin eq. in THETA

double eta_theta;

eta_theta = XI_THETA_SQRT_DT * gsl_ran_gaussian_ziggurat(rand_vec,1.0);

d_theta[tstep] = eta_theta;
theta[tstep] = theta[tstep-1] + eta_theta;
theta[tstep] = wrap_angle(theta[tstep]);

return;
}

/************************************************************************************************
* *
************************************************************************************************/
void dynamics(int tstep)
{
double sin_alpha,cos_alpha;
double dr,dr2;
double cos_theta_old,sin_theta_old,cos_theta_dr,sin_theta_dr;

dynamics_r_parallel(tstep);

cos_theta_old = cos(theta[tstep-1]);
sin_theta_old = sin(theta[tstep-1]);

dynamics_theta(tstep);

vpar[tstep] *= cos(d_theta[tstep]);

dynamics_r_perpendicular(tstep);

//Find total displacement
dr2 = dr_par[tstep]*dr_par[tstep] + dr_perp[tstep]*dr_perp[tstep];
dr = sqrt(dr2);

// Find displacement direction relative to original (old) direction (orientation)
if (dr != 0.0)
{
sin_alpha = dr_perp[tstep]/dr;
cos_alpha = dr_par[tstep] /dr;
}
else // Maintain direction (alpha = 0.0)
{
//fprintf(stderr,"=====>>> dr == 0.0\n");fflush(stderr);
sin_alpha = 0.0;
cos_alpha = 1.0;
}

// Find displacement angle relative to x-axis
// + is the counterclockwise direction: (theta_old + alpha)
cos_theta_dr = cos_theta_old*cos_alpha - sin_alpha*sin_theta_old;
sin_theta_dr = sin_theta_old*cos_alpha + sin_alpha*cos_theta_old;

dx[tstep] = dr*cos_theta_dr;
dy[tstep] = dr*sin_theta_dr;

rx[tstep] = rx[tstep-1] + dx[tstep];
ry[tstep] = ry[tstep-1] + dy[tstep];

// Find velocity
v2[tstep] = vpar[tstep]*vpar[tstep];

return;

}
</source> 

Por fim, colocou-se a ordem de execução das funções na ordem correta e calculou-se os parâmetros termodinâmicos relevantes.
 
<source lang="c">

int main(int argc, char *argv[])
{
int npart,k;
unsigned long int t_counter;

if (TRANSIENT_STEPS > MAX_STEPS)
{
printf("Number of transient steps must be smaller than MAX_STEPS\n");
exit(1);
}
if (argc != 3)
{
printf("Program needs 2 arguments:\n 1) Value of N,\n 2) Value of K.\n");
exit(1);
}
N=atoi(argv[1]);
K=atof(argv[2]);

fprintf(stderr,"Setting constants... "); fflush(stderr);
set_constants();
fprintf(stderr,"Done\n"); fflush(stderr);
printf("Allocating pointers...\n");fflush(stdout);
allocate_pointers(MAX_STEPS);
set_correlation_vecs();
printf("Done: Allocating pointers\n");fflush(stdout);
set_gsl_rng();
printf("Done: Seting initial conditions.\n");fflush(stdout);
open_output_files();
printf("Done: Opening files.\n");fflush(stdout);

for(npart=0; npart<N; ++npart)
{
set_initial_conditions(MAX_STEPS);
for(t_counter=1;t_counter<TRANSIENT/DT;++t_counter)
dynamics(t_counter);

--t_counter; // == TRANSIENT_STEPS - 1 // fprintf(stderr,"Done transient.\n");fflush(stderr);
reset_initial_condition(t_counter);

//fprintf(stderr,"Done transient\n"); fflush(stderr);

for(t_counter=1;t_counter<MAX_STEPS;++t_counter)
dynamics(t_counter);

if (npart==0)
{
for(k=0;k<MAX_STEPS;++k)
fprintf(ftraj,"%f %e %e\n", k*DT, rx[k], ry[k]);
fclose(ftraj);
}
if (npart%10==0)
{
printf("Done particle %d dynamics. Calculating correlations... ", npart);
fflush(stdout);
}
calculate_correlations(rx,ry,vpar);
if (npart%10==0) printf("Done.\n");
}

print_correlations();
fclose(fcorr);
fclose(fr2);

//for(i=0;i<NUM_CORR_MEASURES;++i) printf("%ld\n",corr_time[i]);

printf("Freeing pointers\n");fflush(stdout);
free_pointers(MAX_STEPS);

return 0;
}

/************************************************************************************************
* *
************************************************************************************************/
void set_correlation_vecs()
{
int i;

for(i=0;i<NUM_CORR_MEASURES;++i)
{
r2[i] = 0.0;
v2_ave[i] = 0.0;
v_ave[i] = 0.0;
corr_v_ave[i] = 0.0;
}

for(i=0;i<NUM_VCORR_MEASURES;++i)
corr_v[i] = 0.0;

return;

}
/************************************************************************************************
* *
************************************************************************************************/
void calculate_correlations(double *rx, double *ry, double *v)
{
double c_v, dr2, dx, dy, v_ave_temp,temp;
int i,j,k,num_points;

for(k=0; k<NUM_VCORR_MEASURES;++k)
{
i=k*STEP_SIZE_VCORR; //vcorr_time[k];
c_v = 0.0;
num_points = 0;
for (j=0; j < MAX_STEPS - i; ++j)
{
c_v += v[j+i] * v[j];
++num_points;
}
corr_v[k] += c_v/num_points;
//fprintf(stderr," Done corr delta t %f\n",i*DT);fflush(stderr);
}

//printf("Done 1\n");fflush(stdout);
//for (i=0; i < MAX_STEPS; ++i)
for(k=0; k<NUM_CORR_MEASURES;++k)
{
i=corr_time[k];

dr2 = 0.0;
v_ave_temp = 0.0;
v2_ave[k] += v[i]*v[i];

num_points = 0;
for (j=0; j < MAX_STEPS - i; ++j)
{
dx = rx[j+i]-rx[j];
dy = ry[j+i]-ry[j];
temp = dx*dx + dy*dy;
dr2 += temp; // dx*dx + dy*dy;
v_ave_temp += sqrt(temp)/(i*DT);

++num_points;
}
r2[k] += dr2/num_points;
v_ave[k] += v_ave_temp/num_points;

//fprintf(stderr," Done corr delta t %f\n",i*DT);fflush(stderr);
}

//printf("Done 2\n");fflush(stdout);

/*
* for(k=0; k<NUM_VCORR_MEASURES;++k)
{
i=corr_time[k];

c_v = 0.0;

num_points = 0;
for (j=0; j < MAX_STEPS - i; ++j)
{
c_v += v_ave[j+i] * v_ave[j];
++num_points;
}
corr_v_ave[k] += c_v/num_points;

//fprintf(stderr," Done corr delta t %f\n",i*DT);fflush(stderr);
}
*/
//printf("Done 3\n");fflush(stdout);

return;
}
/************************************************************************************************
* *
************************************************************************************************/
void print_correlations(void)
{
int i,k;
double t;

//for(i=0;i<MAX_STEPS;++i)
for(k=0;k<NUM_VCORR_MEASURES;++k)
{
i=k*STEP_SIZE_VCORR;//vcorr_time[k];
t=i*DT;
corr_v[k] /= N;
fprintf(fcorr,"%f %e\n", t, corr_v[k]);
}

for(k=0;k<NUM_CORR_MEASURES;++k)
{
i=corr_time[k];
t=i*DT;
r2[k] /= N;
v2_ave[k] /= N;
v_ave[k] /= N;
corr_v_ave[k] /= N;
//fprintf(fr2,"%f %e %e %e %e %e\n", t, r2[k], msd_analitico (t), v2_ave[k], v_ave[k],corr_v_ave[k]);
fprintf(fr2,"%f %e %e %e %e\n", t, r2[k], msd_analitico (t), v2_ave[k], v_ave[k]);
}

return;
}
/************************************************************************************************
* *
************************************************************************************************/
void set_initial_conditions(int N)
{
int i;

for (i=0;i<N;++i)
{
rx[i] = 0.0;
ry[i] = 0.0;
dr_par[i] = 0.0;
dr_perp[i] = 0.0;
vpar[i] = 0.0; //(G/(2*M2*(GAMMA_M+K))) * gsl_ran_gaussian_ziggurat(rand_vec,1.0);
v2[i] = 0.0; //vpar[i]*vpar[i]+vperp[i]*vperp[i];
theta[i] = 0.0; //two_PI*gsl_rng_uniform(rand_vec);
}

vpar[0] = (G/(2*M2*(GAMMA_M+K))) * gsl_ran_gaussian_ziggurat(rand_vec,1.0); //raiz, nao?!?!

//Isotropic Initial condition
//Choose a rand_vec with direction uniformly distributed in [0,2pi)
theta[0] = two_PI*gsl_rng_uniform(rand_vec);

return;
}
/************************************************************************************************
* *
************************************************************************************************/
void reset_initial_condition(int t_counter)
{
rx[0] = 0.0; //rx[t_counter];
ry[0] = 0.0; //ry[t_counter];
dr_par[0] = dr_par[t_counter];
dr_perp[0] = dr_perp[t_counter];
theta[0] = theta[t_counter];
vpar[0] = vpar[t_counter];
v2[0] = v2[t_counter];

return;
}

/************************************************************************************************
* *
************************************************************************************************/
void allocate_pointers(int n)
{
rx = create_double_pointer(n);
ry = create_double_pointer(n);
dx = create_double_pointer(n);
dy = create_double_pointer(n);
dr_par = create_double_pointer(n);
dr_perp = create_double_pointer(n);
vpar = create_double_pointer(n);
v2 = create_double_pointer(n);
theta = create_double_pointer(n);
d_theta = create_double_pointer(n);

r2 = create_double_pointer(NUM_CORR_MEASURES);
v2_ave = create_double_pointer(NUM_CORR_MEASURES);
corr_v = create_double_pointer(NUM_VCORR_MEASURES);

v_ave = create_double_pointer(NUM_CORR_MEASURES);
corr_v_ave = create_double_pointer(NUM_CORR_MEASURES);

return;
}

/************************************************************************************************
* *
************************************************************************************************/
void free_pointers(int n)
{

free(rx);
free(ry);
free(dx);
free(dy);
free(dr_par);
free(dr_perp);
free(vpar);
free(corr_v);
free(theta);
free(d_theta);

free(r2);
free(v2);
free(v2_ave);

free(v_ave);
free(corr_v_ave);

gsl_rng_free (rand_vec);

return;

}
/*********************************************************************
*** Time Table ***
*** The total number of time steps and the number of measures are ***
*** specified. ***
*********************************************************************/
void create_time_table(long *t1, int total_time, int measures)
{
unsigned long i,k;
double temp;

t1[0] = 0;
temp = pow((double) total_time,1.0/(measures-1));
k=0;
for (i=1; i<measures; ++i)
{
t1[i] = (int) pow(temp, (double) i);
if (t1[i]<=k) t1[i]=k+1;
k = t1[i];
}
t1[--i] = total_time-2;

return;
}
</source> 
Código escrito pelo Prof. Mendeli Vainstein.

Grupo - Solução Exata para Movimento Anisotrópico de Ornstein-Uhlenbeck

2020-01-10T01:54:53Z

Almeidajuliano:

Arquivo:Szaboilus.png

2020-01-10T01:45:23Z

Almeidajuliano:

Grupo - Solução Exata para Movimento Anisotrópico de Ornstein-Uhlenbeck

2020-01-10T01:38:21Z

Almeidajuliano:

==Introdução==

Movimento coletivo ordenado de diversos organismos é um fenômeno que pode ser observado em diversas situações, como em movimento de cardumes, migração de pássaros, migração celular, etc..

Aquilo que é mais curioso destes tipos de movimento é que o mesmo surge através de movimentos aleatórios dos indivíduos que o compõem, ou seja, o movimento individual e a aleatório de um grupo de células quando postos a interagir entre si leva a um movimento ordenado em curtas escalas de tempo.

Um dos modelos anteriores a este, o modelo de ''Vicsek'', visa descrever o movimento coletivo ordenado de partículas “inteligentes”, isto é, partículas que podem observar as outras partículas na vizinhança e mover-se de acordo com a média do movimento do grupo. Um bom exemplo deste tipo de movimento é o movimento de pássaros ou peixes, da onde movimentos ordenados podem emergir de movimentos aleatórios dos indivíduos autopropulsionados. Curiosamente, mesmo que estes sistemas estejam longe de um estado de equilíbrio, os mesmos apresentam características típicas de sistemas no equilíbrio, como transições de fase bem definidas.

Contudo, partículas “não inteligentes” já haviam demonstrado o mesmo tipo de movimento coletivo experimentalmente, com as mesmas mudanças de fase bem definidas. É para explicar isso, que o modelo de ''Szabó'' foi criado. Nele, o mesmo fenômeno é descrito utilizando somente as forças de interação intercelulares já que isso é algo que independe da capacidade da célula de poder observar a média do movimento da vizinhança, como era em ''Vicsek''.

Essas forças são repulsivas para distâncias (<math>d_{ij} = |r_{i}-r_{j}|</math>) menores que <math>R_{eq}</math>, sendo <math>F_{rep}</math> o seu máximo para esse caso, atrativa para <math>R_{eq} \le d_{ij} \le R_{o}</math>, sendo <math>F_{adh}</math> o máximo para esse caso, e zero para <math>d_{ij}</math> maior que um certo <math>R_{o}</math>.

<math>
\textbf{F}(\textbf{r}_{i}, \textbf{r}_{j} ) = \textbf{e}_{ij}\begin{cases}
F_{rep}\frac{d_{ij} - R_{eq}}{R_{eq}},\quad d_{ij} < R_{eq} \\
F_{adh}\frac{d_{ij} - R_{eq}}{R_{0} - Re_{eq}}, \; \; R_{eq} \le d_{ij} \le R_{0}\\
0, \quad \quad \quad \quad \quad \; d_{ij} > R_{0}\\
\end{cases}
</math>

[[Arquivo:Forcaszabo.png|frame|center|Figura 1: Esquematização da força intercelular de Szabó.]]

Essas forças levam às seguintes equações de movimento:

<math>\frac{dr_{i}(t)}{dt} = v_{o}n_{i}(t) + \mu \sum_{j=1}^N F(r_{i}, r{j})</math>

Onde <math>N</math> é o número de partículas e <math>\sum_{j=1}^N F(r_{i}, r{j})</math> é o somatório da todas as interações sobre uma partícula ''i''.

<math> \frac{d_{i}^{n}(t)}{dt} = \frac{1}{\tau}arcsin[(n_{i}(t)\frac{v_{i}(t)}{|v_{i}(t)|})e_{z}] + \xi </math>

Onde <math>\tau</math> é o tempo de relaxamento, <math>n_{i}(t)</math> é o vetor na direção de polarização, tal que <math>n_{i}(t) = (cos(\theta_{i}^{n}(t)), sen(\theta_{i}^{n}(t)))</math>, <math>e_{z}</math> é o vetor normal, e <math>\xi</math> é um ruído branco.

Neste modelo, o movimento bidimensional das células é descrito por uma dinâmica superamortecida. Na mesma, cada célula tenta manter uma velocidade de magnitude <math>v_{o}</math> na direção do vetor <math>n_{i}(t)</math>, sendo que a direção desse vetor é descrita pelo ângulo <math>\theta_{i}^{n}(t)</math> que tenta relaxar na direção da velocidade enquanto sofre ruído gaussiano.

Embora descreva bem este movimento, o modelo de Szabó fornece somente soluções numéricas para o problema. Para obtermos uma solução analítica, o problema precisa ser descrito de outra forma. É com esse objetivo que o modelo de Rita et al. foi feito.

==Modelo==
Neste modelo, em vez de ter a relaxação da velocidade na direção da polarização, a mesma é projetada sobre a polarização, enquanto ela muda com ruído. Essa direção de polarização interna pode representar, p.e., a polarização dos filamentos de actina na direção preferencial em migração celular.

Com isso pode-se observar duas transições, ou seja, uma transição de um estado desordenado para um estado ordenado e para um estado desordenado de novo. Essas duas transições representam a aleatoriedade do movimento inicial das células, depois a escolha de uma direção preferencial e por último a variação desta direção preferencial.

Contudo, neste modelo deve-se assumir que, para a partícula, seu movimento é anisotrópico, ou seja, apresenta assimetria nas dinâmicas de movimento em direções diferentes, uma na direção paralela à direção de polarização e outra perpendicular.

===Movimento Paralelo===
Na direção paralela à polarização, temos o sistema agindo numa dinâmica semelhante a de ''Langevin'' (<math>m\mathbf{a} = m\frac{d\mathbf{v}}{dt} = F(\mathbf{x}) - \beta \mathbf{v} + \eta(t)</math>), mais especificamente, um processo de ''Orstein-Uhlenbeck'' amortecido. Nesta direção a velocidade (paralela) é bem definida. Tanto para o movimento paralelo, quanto para o movimento perpendicular, o vetor <math>p_{i}(t) = (cos(\theta_{i}(t)), sen(\theta_{i}(t)))</math> é o vetor que indica a direção de polarização e <math>n_{i}(t) = (sen(\theta_{i}(t)), -cos(\theta_{i}(t)))</math> é o vetor perpendicular a <math>p_{i}(t)</math>.

<math>\frac{dv_{||}}{dt} = -\frac{\gamma}{m}v_{||} + \frac{\xi_{||}(t)}{m} </math>
<math>v_{||}(t+\Delta t)p(t+\Delta t) = [(1-\gamma\Delta t)v_{||}(t) + \int_{t}^{t+\Delta t} \xi_{||}(t)dt](p_{i}(t)p_{i}(t+\Delta t))p_{i}(t+\Delta t) </math>

Onde <math>v_{||}</math> é a velocidade paralela, <math>\gamma</math> é o coeficiente de atrito e <math>\xi_{||}</math> é o ruído gaussiano paralelo.

A velocidade diminui por causa do atrito com o substrato, já que pode-se observar que <math>(1-\gamma\Delta t)</math> tende a <math>e^{-\gamma\Delta t}</math> se <math>\gamma</math> for pequeno, e diminui ou aumenta por causa do termo estocástico <math>\int_{t}^{t+\Delta t} \xi_{||}(t)dt = W(t)</math>. A primeira equação assume que durante o intervalo <math>\Delta t</math>, a direção de polarização é constante. Somente em <math>t+\Delta t</math>, a velocidade na direção <math>p_{i}(t)</math> é projetada sobre a nova direção de polarização <math>p_{i}(t+\Delta t)</math>.

===Movimento Perpendicular===
Na direção perpendicular à polarização, o sistema segue um processo de ''Wiener'' (movimento aleatório) por causa das flutuações na direção de polarização.

<math>\frac{dr_{\perp}}{dt} = n_{i}(t)\xi_{\perp}(t) </math>
<math>[r_{\perp}(t+\Delta t) - r_{\perp}(t)]n_{i}(t) = n_{i}(t)\int_{t}^{t+\Delta t} \xi_{\perp}(t)dt </math>

Onde <math>\xi_{\perp}</math> é o ruído gaussiano na direção perpendicular. Pode-se observar que <math>\int_{t}^{t+\Delta t} \xi_{\perp}(t)dt = W(t) </math>. Este movimento, por ser um processo de ''Wiener'', é melhor descrito pelo deslocamento do que pela velocidade.

===Ângulo e Ruídos===
A direção de polarização é definida pelo ângulo teta que também apresenta ruído branco. Esse ruído <math>\beta(t)</math> é relacionado à <math> \xi_{\perp}(t)</math> através de <math>\xi_{perp}(t) = \sqrt(q)\beta(t)</math>. É essa mudança na direção de polarização que causa o deslocamento aleatório perpendicular.

Embora o ruído no ângulo e o ruído na direção perpendicular estejam relacionados, os mesmos são independentes do ruído na direção paralela. Abaixo pode-se ver a correlação dos mesmos. Os três são ruídos gaussianos, então a média deles é zero.

<math> <\xi_{||}(t)\xi_{||}(t')> = g\delta(t-t') </math>
<math> <\beta(t)\beta(t')> = 2k\delta(t-t') </math>
<math> <\xi_{\perp}(t)\xi_{\perp}(t')> = 2qk\delta(t-t') </math>

Onde <math>k</math> tem unidades de 1/tempo, <math>\sqrt q</math> tem unidades de comprimento, e tanto <math>g</math> quanto <math>qk</math> tem unidades de [comprimento]^2 / tempo.

Em resumo, este modelo considera uma partícula com dois graus de liberdades espaciais e um grau de liberdade interna, que considera a polarização interna. Esse grau de liberdade interna quebra a simetria espacial de tal forma que a dinâmica na direção da polarização é semelhante a ''Langevin'' e na direção perpendicular a dinâmica segue um processo de ''Wiener''. Existem duas fontes independentes de ruído: Uma age na dinâmica na direção de polarização e a outra age na polarização em si. Além da força de dissipação, a mudança na polarização age como mais um termo na perda de energia cinética.

[[Arquivo:Wforce.gif|frame|350x350px|center|Figura 2: Simulação do modelo de Rita et al. A seta azul representa a direção de polarização e a seta verde a velocidade.]]

===Mean Square Displacement(MSD)===
O MSD é uma quantificação do movimento líquido das partículas. O mesmo foi utilizado para medir a aglomeração e a transição de fase das partículas. Para isso foi feita uma média do MSD de 100 trajetórias sendo que o Método das Janelas foi feito para cada uma das mesmas. Através do resultados analíticos para as equações de movimento pode-se obter um resultado analítico para o MSD também:

<math><r^{2}> = \frac{4DP}{1-S}[\tau - (1-S)(1-e^{-\tau})]</math>

Onde <math>P</math> é o tempo de persistência do movimento em um estado cinético, <math>S</math> é uma fração deste tempo, relacionado a quanto tempo o movimento persiste no primeiro estado cinético e <math>D =\frac{g}{2m^{2}\tau\tau_{b}}</math>, sendo que <math>\tau = \frac{\gamma}{m} + k </math> e <math>\tau_{b} = \frac{\gamma}{m} + 2k </math>.

Para testar este modelo resolveu-se numericamente as equações diferenciais de movimento com um algoritmo de Euler-Maruyama e comparadas com os resultados analíticos do MSD para diferentes valores do parâmetro <math>k</math>.

O resultado obtido consta abaixo. Nele os pontos são os valores obtidos pela simulação e as linhas são as previsões da solução analítica para o MSD.

[[Arquivo:Msd_v2.png|frame|center|Figura 3: MSD analítico para 100 trajetórias]]

===O Programa===
O programa criado para gerar simulações do modelo descrito e calcular o MSD do mesmo foi feito em C. O mesmo começa com uma inialização de todas as variáveis necessárias e a definição das mesmas de acordo com aquilo que foi descrito acima sendo os parâmetros iniciais:<math>k, q, g</math> e <math> \gamma</math> e os parâmetros termodinâmicos <math> S,P </math> e <math>D</math> que descrevem o comportamento macroscópico da célula.

Em seguida foram definidas as funções que iriam calcular a dinâmica descrita acima, sendo que cada parte foi feita numa função separada (p.e. a dinâmica paralela e a perpendicular). Com isso polarizações e velocidades iniciais foram atribuidas, tal que as mesmas começassem em equilibrio.

Por fim, colocou-se a ordem de execução das funções na ordem correta e calculou-se os parâmetros termodinâmicos relevantes.

Grupo - Solução Exata para Movimento Anisotrópico de Ornstein-Uhlenbeck

2020-01-10T01:26:30Z

Almeidajuliano:

Arquivo:Forcaszabo.png

2020-01-10T01:23:44Z

Almeidajuliano:

Grupo - Solução Exata para Movimento Anisotrópico de Ornstein-Uhlenbeck

2020-01-10T01:20:24Z

Almeidajuliano:

==Introdução==

Movimento coletivo ordenado de diversos organismos é um fenômeno que pode ser observado em diversas situações, como em movimento de cardumes, migração de pássaros, migração celular, etc..

Aquilo que é mais curioso destes tipos de movimento é que o mesmo surge através de movimentos aleatórios dos indivíduos que o compõem, ou seja, o movimento individual e a aleatório de um grupo de células quando postos a interagir entre si leva a um movimento ordenado em curtas escalas de tempo.

Um dos modelos anteriores a este, o modelo de ''Vicsek'', visa descrever o movimento coletivo ordenado de partículas “inteligentes”, isto é, partículas que podem observar as outras partículas na vizinhança e mover-se de acordo com a média do movimento do grupo. Um bom exemplo deste tipo de movimento é o movimento de pássaros ou peixes, da onde movimentos ordenados podem emergir de movimentos aleatórios dos indivíduos autopropulsionados. Curiosamente, mesmo que estes sistemas estejam longe de um estado de equilíbrio, os mesmos apresentam características típicas de sistemas no equilíbrio, como transições de fase bem definidas.

Contudo, partículas “não inteligentes” já haviam demonstrado o mesmo tipo de movimento coletivo experimentalmente, com as mesmas mudanças de fase bem definidas. É para explicar isso, que o modelo de ''Szabó'' foi criado. Nele, o mesmo fenômeno é descrito utilizando somente as forças de interação intracelulares já que isso é algo que independe da capacidade da célula de poder observar a média do movimento da vizinhança, como era em ''Vicsek''.

Essas forças são repulsivas para distâncias (<math>d_{ij} = |r_{i}-r_{j}|</math>) menores que <math>R_{eq}</math>, sendo <math>F_{rep}</math> o seu máximo para esse caso, atrativa para <math>R_{eq} \le d_{ij} \le R_{o}</math>, sendo <math>F_{adh}</math> o máximo para esse caso, e zero para <math>d_{ij}</math> maior que um certo <math>R_{o}</math>.

<math>
\textbf{F}(\textbf{r}_{i}, \textbf{r}_{j} ) = \textbf{e}_{ij}\begin{cases}
F_{rep}\frac{d_{ij} - R_{eq}}{R_{eq}},\quad d_{ij} < R_{eq} \\
F_{adh}\frac{d_{ij} - R_{eq}}{R_{0} - Re_{eq}}, \; \; R_{eq} \le d_{ij} \le R_{0}\\
0, \quad \quad \quad \quad \quad \; d_{ij} > R_{0}\\
\end{cases}
</math>

Essas forças levam às seguintes equações de movimento:

<math>\frac{dr_{i}(t)}{dt} = v_{o}n_{i}(t) + \mu \sum_{j=1}^N F(r_{i}, r{j})</math>

Onde <math>N</math> é o número de partículas e <math>\sum_{j=1}^N F(r_{i}, r{j})</math> é o somatório da todas as interações sobre uma partícula ''i''.

<math> \frac{d_{i}^{n}(t)}{dt} = \frac{1}{\tau}arcsin[(n_{i}(t)\frac{v_{i}(t)}{|v_{i}(t)|})e_{z}] + \xi </math>

Onde <math>\tau</math> é o tempo de relaxamento, <math>n_{i}(t)</math> é o vetor na direção de polarização, tal que <math>n_{i}(t) = (cos(\theta_{i}^{n}(t)), sen(\theta_{i}^{n}(t)))</math>, <math>e_{z}</math> é o vetor normal, e <math>\xi</math> é um ruído branco.

Neste modelo, o movimento bidimensional das células é descrito por uma dinâmica superamortecida. Na mesma, cada célula tenta manter uma velocidade de magnitude <math>v_{o}</math> na direção do vetor <math>n_{i}(t)</math>, sendo que a direção desse vetor é descrita pelo ângulo <math>\theta_{i}^{n}(t)</math> que tenta relaxar na direção da velocidade enquanto sofre ruído gaussiano.

Embora descreva bem este movimento, o modelo de Szabó fornece somente soluções numéricas para o problema. Para obtermos uma solução analítica, o problema precisa ser descrito de outra forma. É com esse objetivo que o modelo de Rita et al. foi feito.

==Modelo==
Neste modelo, em vez de ter a relaxação da velocidade na direção da polarização, a mesma é projetada sobre a polarização, enquanto ela muda com ruído. Essa direção de polarização interna pode representar, p.e., a polarização dos filamentos de actina na direção preferencial em migração celular.

Com isso pode-se observar duas transições, ou seja, uma transição de um estado desordenado para um estado ordenado e para um estado desordenado de novo. Essas duas transições representam a aleatoriedade do movimento inicial das células, depois a escolha de uma direção preferencial e por último a variação desta direção preferencial.

Contudo, neste modelo deve-se assumir que, para a partícula, seu movimento é anisotrópico, ou seja, apresenta assimetria nas dinâmicas de movimento em direções diferentes, uma na direção paralela à direção de polarização e outra perpendicular.

===Movimento Paralelo===
Na direção paralela à polarização, temos o sistema agindo numa dinâmica semelhante a de ''Langevin'' (<math>m\mathbf{a} = m\frac{d\mathbf{v}}{dt} = F(\mathbf{x}) - \beta \mathbf{v} + \eta(t)</math>), mais especificamente, um processo de ''Orstein-Uhlenbeck'' amortecido. Nesta direção a velocidade (paralela) é bem definida. Tanto para o movimento paralelo, quanto para o movimento perpendicular, o vetor <math>p_{i}(t) = (cos(\theta_{i}(t)), sen(\theta_{i}(t)))</math> é o vetor que indica a direção de polarização e <math>n_{i}(t) = (sen(\theta_{i}(t)), -cos(\theta_{i}(t)))</math> é o vetor perpendicular a <math>p_{i}(t)</math>.

<math>\frac{dv_{||}}{dt} = -\frac{\gamma}{m}v_{||} + \frac{\xi_{||}(t)}{m} </math>
<math>v_{||}(t+\Delta t)p(t+\Delta t) = [(1-\gamma\Delta t)v_{||}(t) + \int_{t}^{t+\Delta t} \xi_{||}(t)dt](p_{i}(t)p_{i}(t+\Delta t))p_{i}(t+\Delta t) </math>

Onde <math>v_{||}</math> é a velocidade paralela, <math>\gamma</math> é o coeficiente de atrito e <math>\xi_{||}</math> é o ruído gaussiano paralelo.

A velocidade diminui por causa do atrito com o substrato, já que pode-se observar que <math>(1-\gamma\Delta t)</math> tende a <math>e^{-\gamma\Delta t}</math> se <math>\gamma</math> for pequeno, e diminui ou aumenta por causa do termo estocástico <math>\int_{t}^{t+\Delta t} \xi_{||}(t)dt = W(t)</math>. A primeira equação assume que durante o intervalo <math>\Delta t</math>, a direção de polarização é constante. Somente em <math>t+\Delta t</math>, a velocidade na direção <math>p_{i}(t)</math> é projetada sobre a nova direção de polarização <math>p_{i}(t+\Delta t)</math>.

===Movimento Perpendicular===
Na direção perpendicular à polarização, o sistema segue um processo de ''Wiener'' (movimento aleatório) por causa das flutuações na direção de polarização.

<math>\frac{dr_{\perp}}{dt} = n_{i}(t)\xi_{\perp}(t) </math>
<math>[r_{\perp}(t+\Delta t) - r_{\perp}(t)]n_{i}(t) = n_{i}(t)\int_{t}^{t+\Delta t} \xi_{\perp}(t)dt </math>

Onde <math>\xi_{\perp}</math> é o ruído gaussiano na direção perpendicular. Pode-se observar que <math>\int_{t}^{t+\Delta t} \xi_{\perp}(t)dt = W(t) </math>. Este movimento, por ser um processo de ''Wiener'', é melhor descrito pelo deslocamento do que pela velocidade.

===Ângulo e Ruídos===
A direção de polarização é definida pelo ângulo teta que também apresenta ruído branco. Esse ruído <math>\beta(t)</math> é relacionado à <math> \xi_{\perp}(t)</math> através de <math>\xi_{perp}(t) = \sqrt(q)\beta(t)</math>. É essa mudança na direção de polarização que causa o deslocamento aleatório perpendicular.

Embora o ruído no ângulo e o ruído na direção perpendicular estejam relacionados, os mesmos são independentes do ruído na direção paralela. Abaixo pode-se ver a correlação dos mesmos. Os três são ruídos gaussianos, então a média deles é zero.

<math> <\xi_{||}(t)\xi_{||}(t')> = g\delta(t-t') </math>
<math> <\beta(t)\beta(t')> = 2k\delta(t-t') </math>
<math> <\xi_{\perp}(t)\xi_{\perp}(t')> = 2qk\delta(t-t') </math>

Onde <math>k</math> tem unidades de 1/tempo, <math>\sqrt q</math> tem unidades de comprimento, e tanto <math>g</math> quanto <math>qk</math> tem unidades de [comprimento]^2 / tempo.

Em resumo, este modelo considera uma partícula com dois graus de liberdades espaciais e um grau de liberdade interna, que considera a polarização interna. Esse grau de liberdade interna quebra a simetria espacial de tal forma que a dinâmica na direção da polarização é semelhante a ''Langevin'' e na direção perpendicular a dinâmica segue um processo de ''Wiener''. Existem duas fontes independentes de ruído: Uma age na dinâmica na direção de polarização e a outra age na polarização em si. Além da força de dissipação, a mudança na polarização age como mais um termo na perda de energia cinética.

[[Arquivo:Wforce.gif|frame|350x350px|center|Figura 1: Simulação do modelo de Rita et al. A seta azul representa a direção de polarização e a seta verde a velocidade.]]

===Mean Square Displacement(MSD)===
O MSD é uma quantificação do movimento líquido das partículas. O mesmo foi utilizado para medir a aglomeração e a transição de fase das partículas. Para isso foi feita uma média do MSD de 100 trajetórias sendo que o Método das Janelas foi feito para cada uma das mesmas. Através do resultados analíticos para as equações de movimento pode-se obter um resultado analítico para o MSD também:

<math><r^{2}> = \frac{4DP}{1-S}[\tau - (1-S)(1-e^{-\tau})]</math>

Onde <math>P</math> é o tempo de persistência do movimento em um estado cinético, <math>S</math> é uma fração deste tempo, relacionado a quanto tempo o movimento persiste no primeiro estado cinético e <math>D =\frac{g}{2m^{2}\tau\tau_{b}}</math>, sendo que <math>\tau = \frac{\gamma}{m} + k </math> e <math>\tau_{b} = \frac{\gamma}{m} + 2k </math>.

Para testar este modelo resolveu-se numericamente as equações diferenciais de movimento com um algoritmo de Euler-Maruyama e comparadas com os resultados analíticos do MSD para diferentes valores do parâmetro <math>k</math>.

O resultado obtido consta abaixo. Nele os pontos são os valores obtidos pela simulação e as linhas são as previsões da solução analítica para o MSD.

[[Arquivo:Msd_v2.png|frame|center|Figura 2: MSD analítico para 100 trajetórias]]

Arquivo:Wforce.gif

2020-01-10T01:13:09Z

Almeidajuliano:

Grupo - Solução Exata para Movimento Anisotrópico de Ornstein-Uhlenbeck

2020-01-10T01:12:02Z

Almeidajuliano:

==Introdução==

Movimento coletivo ordenado de diversos organismos é um fenômeno que pode ser observado em diversas situações, como em movimento de cardumes, migração de pássaros, migração celular, etc..

Aquilo que é mais curioso destes tipos de movimento é que o mesmo surge através de movimentos aleatórios dos indivíduos que o compõem, ou seja, o movimento individual e a aleatório de um grupo de células quando postos a interagir entre si leva a um movimento ordenado em curtas escalas de tempo.

Um dos modelos anteriores a este, o modelo de ''Vicsek'', visa descrever o movimento coletivo ordenado de partículas “inteligentes”, isto é, partículas que podem observar as outras partículas na vizinhança e mover-se de acordo com a média do movimento do grupo. Um bom exemplo deste tipo de movimento é o movimento de pássaros ou peixes, da onde movimentos ordenados podem emergir de movimentos aleatórios dos indivíduos autopropulsionados. Curiosamente, mesmo que estes sistemas estejam longe de um estado de equilíbrio, os mesmos apresentam características típicas de sistemas no equilíbrio, como transições de fase bem definidas.

Contudo, partículas “não inteligentes” já haviam demonstrado o mesmo tipo de movimento coletivo experimentalmente, com as mesmas mudanças de fase bem definidas. É para explicar isso, que o modelo de ''Szabó'' foi criado. Nele, o mesmo fenômeno é descrito utilizando somente as forças de interação intracelulares já que isso é algo que independe da capacidade da célula de poder observar a média do movimento da vizinhança, como era em ''Vicsek''.

Essas forças são repulsivas para distâncias (<math>d_{ij} = |r_{i}-r_{j}|</math>) menores que <math>R_{eq}</math>, sendo <math>F_{rep}</math> o seu máximo para esse caso, atrativa para <math>R_{eq} \le d_{ij} \le R_{o}</math>, sendo <math>F_{adh}</math> o máximo para esse caso, e zero para <math>d_{ij}</math> maior que um certo <math>R_{o}</math>.

<math>
\textbf{F}(\textbf{r}_{i}, \textbf{r}_{j} ) = \textbf{e}_{ij}\begin{cases}
F_{rep}\frac{d_{ij} - R_{eq}}{R_{eq}},\quad d_{ij} < R_{eq} \\
F_{adh}\frac{d_{ij} - R_{eq}}{R_{0} - Re_{eq}}, \; \; R_{eq} \le d_{ij} \le R_{0}\\
0, \quad \quad \quad \quad \quad \; d_{ij} > R_{0}\\
\end{cases}
</math>

Essas forças levam às seguintes equações de movimento:

<math>\frac{dr_{i}(t)}{dt} = v_{o}n_{i}(t) + \mu \sum_{j=1}^N F(r_{i}, r{j})</math>

Onde <math>N</math> é o número de partículas e <math>\sum_{j=1}^N F(r_{i}, r{j})</math> é o somatório da todas as interações sobre uma partícula ''i''.

<math> \frac{d_{i}^{n}(t)}{dt} = \frac{1}{\tau}arcsin[(n_{i}(t)\frac{v_{i}(t)}{|v_{i}(t)|})e_{z}] + \xi </math>

Onde <math>\tau</math> é o tempo de relaxamento, <math>n_{i}(t)</math> é o vetor na direção de polarização, tal que <math>n_{i}(t) = (cos(\theta_{i}^{n}(t)), sen(\theta_{i}^{n}(t)))</math>, <math>e_{z}</math> é o vetor normal, e <math>\xi</math> é um ruído branco.

Neste modelo, o movimento bidimensional das células é descrito por uma dinâmica superamortecida. Na mesma, cada célula tenta manter uma velocidade de magnitude <math>v_{o}</math> na direção do vetor <math>n_{i}(t)</math>, sendo que a direção desse vetor é descrita pelo ângulo <math>\theta_{i}^{n}(t)</math> que tenta relaxar na direção da velocidade enquanto sofre ruído gaussiano.

Embora descreva bem este movimento, o modelo de Szabó fornece somente soluções numéricas para o problema. Para obtermos uma solução analítica, o problema precisa ser descrito de outra forma. É com esse objetivo que o modelo de Rita et al. foi feito.

==Modelo==
Neste modelo, em vez de ter a relaxação da velocidade na direção da polarização, a mesma é projetada sobre a polarização, enquanto ela muda com ruído. Essa direção de polarização interna pode representar, p.e., a polarização dos filamentos de actina na direção preferencial em migração celular.

Com isso pode-se observar duas transições, ou seja, uma transição de um estado desordenado para um estado ordenado e para um estado desordenado de novo. Essas duas transições representam a aleatoriedade do movimento inicial das células, depois a escolha de uma direção preferencial e por último a variação desta direção preferencial.

Contudo, neste modelo deve-se assumir que, para a partícula, seu movimento é anisotrópico, ou seja, apresenta assimetria nas dinâmicas de movimento em direções diferentes, uma na direção paralela à direção de polarização e outra perpendicular.

===Movimento Paralelo===
Na direção paralela à polarização, temos o sistema agindo numa dinâmica semelhante a de ''Langevin'' (<math>m\mathbf{a} = m\frac{d\mathbf{v}}{dt} = F(\mathbf{x}) - \beta \mathbf{v} + \eta(t)</math>), mais especificamente, um processo de ''Orstein-Uhlenbeck'' amortecido. Nesta direção a velocidade (paralela) é bem definida. Tanto para o movimento paralelo, quanto para o movimento perpendicular, o vetor <math>p_{i}(t) = (cos(\theta_{i}(t)), sen(\theta_{i}(t)))</math> é o vetor que indica a direção de polarização e <math>n_{i}(t) = (sen(\theta_{i}(t)), -cos(\theta_{i}(t)))</math> é o vetor perpendicular a <math>p_{i}(t)</math>.

<math>\frac{dv_{||}}{dt} = -\frac{\gamma}{m}v_{||} + \frac{\xi_{||}(t)}{m} </math>
<math>v_{||}(t+\Delta t)p(t+\Delta t) = [(1-\gamma\Delta t)v_{||}(t) + \int_{t}^{t+\Delta t} \xi_{||}(t)dt](p_{i}(t)p_{i}(t+\Delta t))p_{i}(t+\Delta t) </math>

Onde <math>v_{||}</math> é a velocidade paralela, <math>\gamma</math> é o coeficiente de atrito e <math>\xi_{||}</math> é o ruído gaussiano paralelo.

A velocidade diminui por causa do atrito com o substrato, já que pode-se observar que <math>(1-\gamma\Delta t)</math> tende a <math>e^{-\gamma\Delta t}</math> se <math>\gamma</math> for pequeno, e diminui ou aumenta por causa do termo estocástico <math>\int_{t}^{t+\Delta t} \xi_{||}(t)dt = W(t)</math>. A primeira equação assume que durante o intervalo <math>\Delta t</math>, a direção de polarização é constante. Somente em <math>t+\Delta t</math>, a velocidade na direção <math>p_{i}(t)</math> é projetada sobre a nova direção de polarização <math>p_{i}(t+\Delta t)</math>.

===Movimento Perpendicular===
Na direção perpendicular à polarização, o sistema segue um processo de ''Wiener'' (movimento aleatório) por causa das flutuações na direção de polarização.

<math>\frac{dr_{\perp}}{dt} = n_{i}(t)\xi_{\perp}(t) </math>
<math>[r_{\perp}(t+\Delta t) - r_{\perp}(t)]n_{i}(t) = n_{i}(t)\int_{t}^{t+\Delta t} \xi_{\perp}(t)dt </math>

Onde <math>\xi_{\perp}</math> é o ruído gaussiano na direção perpendicular. Pode-se observar que <math>\int_{t}^{t+\Delta t} \xi_{\perp}(t)dt = W(t) </math>. Este movimento, por ser um processo de ''Wiener'', é melhor descrito pelo deslocamento do que pela velocidade.

===Ângulo e Ruídos===
A direção de polarização é definida pelo ângulo teta que também apresenta ruído branco. Esse ruído <math>\beta(t)</math> é relacionado à <math> \xi_{\perp}(t)</math> através de <math>\xi_{perp}(t) = \sqrt(q)\beta(t)</math>. É essa mudança na direção de polarização que causa o deslocamento aleatório perpendicular.

Embora o ruído no ângulo e o ruído na direção perpendicular estejam relacionados, os mesmos são independentes do ruído na direção paralela. Abaixo pode-se ver a correlação dos mesmos. Os três são ruídos gaussianos, então a média deles é zero.

<math> <\xi_{||}(t)\xi_{||}(t')> = g\delta(t-t') </math>
<math> <\beta(t)\beta(t')> = 2k\delta(t-t') </math>
<math> <\xi_{\perp}(t)\xi_{\perp}(t')> = 2qk\delta(t-t') </math>

Onde <math>k</math> tem unidades de 1/tempo, <math>\sqrt q</math> tem unidades de comprimento, e tanto <math>g</math> quanto <math>qk</math> tem unidades de [comprimento]^2 / tempo.

Em resumo, este modelo considera uma partícula com dois graus de liberdades espaciais e um grau de liberdade interna, que considera a polarização interna. Esse grau de liberdade interna quebra a simetria espacial de tal forma que a dinâmica na direção da polarização é semelhante a ''Langevin'' e na direção perpendicular a dinâmica segue um processo de ''Wiener''. Existem duas fontes independentes de ruído: Uma age na dinâmica na direção de polarização e a outra age na polarização em si. Além da força de dissipação, a mudança na polarização age como mais um termo na perda de energia cinética.

===Mean Square Displacement(MSD)===
O MSD é uma quantificação do movimento líquido das partículas. O mesmo foi utilizado para medir a aglomeração e a transição de fase das partículas. Para isso foi feita uma média do MSD de 100 trajetórias sendo que o Método das Janelas foi feito para cada uma das mesmas. Através do resultados analíticos para as equações de movimento pode-se obter um resultado analítico para o MSD também:

<math><r^{2}> = \frac{4DP}{1-S}[\tau - (1-S)(1-e^{-\tau})]</math>

Onde <math>P</math> é o tempo de persistência do movimento em um estado cinético, <math>S</math> é uma fração deste tempo, relacionado a quanto tempo o movimento persiste no primeiro estado cinético e <math>D =\frac{g}{2m^{2}\tau\tau_{b}}</math>, sendo que <math>\tau = \frac{\gamma}{m} + k </math> e <math>\tau_{b} = \frac{\gamma}{m} + 2k </math>.

Para testar este modelo resolveu-se numericamente as equações diferenciais de movimento com um algoritmo de Euler-Maruyama e comparadas com os resultados analíticos do MSD para diferentes valores do parâmetro <math>k</math>.

O resultado obtido consta abaixo. Nele os pontos são os valores obtidos pela simulação e as linhas são as previsões da solução analítica para o MSD.

[[Arquivo:Msd_v2.png|frame|center|Figura 2: MSD analítico para 100 trajetórias]]

Arquivo:Msd v2.png

2020-01-10T01:04:06Z

Almeidajuliano:

Grupo - Solução Exata para Movimento Anisotrópico de Ornstein-Uhlenbeck

2020-01-10T00:57:09Z

Almeidajuliano:

==Introdução==

Movimento coletivo ordenado de diversos organismos é um fenômeno que pode ser observado em diversas situações, como em movimento de cardumes, migração de pássaros, migração celular, etc..

Aquilo que é mais curioso destes tipos de movimento é que o mesmo surge através de movimentos aleatórios dos indivíduos que o compõem, ou seja, o movimento individual e a aleatório de um grupo de células quando postos a interagir entre si leva a um movimento ordenado em curtas escalas de tempo.

Um dos modelos anteriores a este, o modelo de ''Vicsek'', visa descrever o movimento coletivo ordenado de partículas “inteligentes”, isto é, partículas que podem observar as outras partículas na vizinhança e mover-se de acordo com a média do movimento do grupo. Um bom exemplo deste tipo de movimento é o movimento de pássaros ou peixes, da onde movimentos ordenados podem emergir de movimentos aleatórios dos indivíduos autopropulsionados. Curiosamente, mesmo que estes sistemas estejam longe de um estado de equilíbrio, os mesmos apresentam características típicas de sistemas no equilíbrio, como transições de fase bem definidas.

Contudo, partículas “não inteligentes” já haviam demonstrado o mesmo tipo de movimento coletivo experimentalmente, com as mesmas mudanças de fase bem definidas. É para explicar isso, que o modelo de ''Szabó'' foi criado. Nele, o mesmo fenômeno é descrito utilizando somente as forças de interação intracelulares já que isso é algo que independe da capacidade da célula de poder observar a média do movimento da vizinhança, como era em ''Vicsek''.

Essas forças são repulsivas para distâncias (<math>d_{ij} = |r_{i}-r_{j}|</math>) menores que <math>R_{eq}</math>, sendo <math>F_{rep}</math> o seu máximo para esse caso, atrativa para <math>R_{eq} \le d_{ij} \le R_{o}</math>, sendo <math>F_{adh}</math> o máximo para esse caso, e zero para <math>d_{ij}</math> maior que um certo <math>R_{o}</math>.

<math>
\textbf{F}(\textbf{r}_{i}, \textbf{r}_{j} ) = \textbf{e}_{ij}\begin{cases}
F_{rep}\frac{d_{ij} - R_{eq}}{R_{eq}},\quad d_{ij} < R_{eq} \\
F_{adh}\frac{d_{ij} - R_{eq}}{R_{0} - Re_{eq}}, \; \; R_{eq} \le d_{ij} \le R_{0}\\
0, \quad \quad \quad \quad \quad \; d_{ij} > R_{0}\\
\end{cases}
</math>

Essas forças levam às seguintes equações de movimento:

<math>\frac{dr_{i}(t)}{dt} = v_{o}n_{i}(t) + \mu \sum_{j=1}^N F(r_{i}, r{j})</math>

Onde <math>N</math> é o número de partículas e <math>\sum_{j=1}^N F(r_{i}, r{j})</math> é o somatório da todas as interações sobre uma partícula ''i''.

<math> \frac{d_{i}^{n}(t)}{dt} = \frac{1}{\tau}arcsin[(n_{i}(t)\frac{v_{i}(t)}{|v_{i}(t)|})e_{z}] + \xi </math>

Onde <math>\tau</math> é o tempo de relaxamento, <math>n_{i}(t)</math> é o vetor na direção de polarização, tal que <math>n_{i}(t) = (cos(\theta_{i}^{n}(t)), sen(\theta_{i}^{n}(t)))</math>, <math>e_{z}</math> é o vetor normal, e <math>\xi</math> é um ruído branco.

Neste modelo, o movimento bidimensional das células é descrito por uma dinâmica superamortecida. Na mesma, cada célula tenta manter uma velocidade de magnitude <math>v_{o}</math> na direção do vetor <math>n_{i}(t)</math>, sendo que a direção desse vetor é descrita pelo ângulo <math>\theta_{i}^{n}(t)</math> que tenta relaxar na direção da velocidade enquanto sofre ruído gaussiano.

Embora descreva bem este movimento, o modelo de Szabó fornece somente soluções numéricas para o problema. Para obtermos uma solução analítica, o problema precisa ser descrito de outra forma. É com esse objetivo que o modelo de Rita et al. foi feito.

==Modelo==
Neste modelo, em vez de ter a relaxação da velocidade na direção da polarização, a mesma é projetada sobre a polarização, enquanto ela muda com ruído. Essa direção de polarização interna pode representar, p.e., a polarização dos filamentos de actina na direção preferencial em migração celular.

Com isso pode-se observar duas transições, ou seja, uma transição de um estado desordenado para um estado ordenado e para um estado desordenado de novo. Essas duas transições representam a aleatoriedade do movimento inicial das células, depois a escolha de uma direção preferencial e por último a variação desta direção preferencial.

Contudo, neste modelo deve-se assumir que, para a partícula, seu movimento é anisotrópico, ou seja, apresenta assimetria nas dinâmicas de movimento em direções diferentes, uma na direção paralela à direção de polarização e outra perpendicular.

===Movimento Paralelo===
Na direção paralela à polarização, temos o sistema agindo numa dinâmica semelhante a de ''Langevin'' (<math>m\mathbf{a} = m\frac{d\mathbf{v}}{dt} = F(\mathbf{x}) - \beta \mathbf{v} + \eta(t)</math>), mais especificamente, um processo de ''Orstein-Uhlenbeck'' amortecido. Nesta direção a velocidade (paralela) é bem definida. Tanto para o movimento paralelo, quanto para o movimento perpendicular, o vetor <math>p_{i}(t) = (cos(\theta_{i}(t)), sen(\theta_{i}(t)))</math> é o vetor que indica a direção de polarização e <math>n_{i}(t) = (sen(\theta_{i}(t)), -cos(\theta_{i}(t)))</math> é o vetor perpendicular a <math>p_{i}(t)</math>.

<math>\frac{dv_{||}}{dt} = -\frac{\gamma}{m}v_{||} + \frac{\xi_{||}(t)}{m} </math>
<math>v_{||}(t+\Delta t)p(t+\Delta t) = [(1-\gamma\Delta t)v_{||}(t) + \int_{t}^{t+\Delta t} \xi_{||}(t)dt](p_{i}(t)p_{i}(t+\Delta t))p_{i}(t+\Delta t) </math>

Onde <math>v_{||}</math> é a velocidade paralela, <math>\gamma</math> é o coeficiente de atrito e <math>\xi_{||}</math> é o ruído gaussiano paralelo.

A velocidade diminui por causa do atrito com o substrato, já que pode-se observar que <math>(1-\gamma\Delta t)</math> tende a <math>e^{-\gamma\Delta t}</math> se <math>\gamma</math> for pequeno, e diminui ou aumenta por causa do termo estocástico <math>\int_{t}^{t+\Delta t} \xi_{||}(t)dt = W(t)</math>. A primeira equação assume que durante o intervalo <math>\Delta t</math>, a direção de polarização é constante. Somente em <math>t+\Delta t</math>, a velocidade na direção <math>p_{i}(t)</math> é projetada sobre a nova direção de polarização <math>p_{i}(t+\Delta t)</math>.

===Movimento Perpendicular===
Na direção perpendicular à polarização, o sistema segue um processo de ''Wiener'' (movimento aleatório) por causa das flutuações na direção de polarização.

<math>\frac{dr_{\perp}}{dt} = n_{i}(t)\xi_{\perp}(t) </math>
<math>[r_{\perp}(t+\Delta t) - r_{\perp}(t)]n_{i}(t) = n_{i}(t)\int_{t}^{t+\Delta t} \xi_{\perp}(t)dt </math>

Onde <math>\xi_{\perp}</math> é o ruído gaussiano na direção perpendicular. Pode-se observar que <math>\int_{t}^{t+\Delta t} \xi_{\perp}(t)dt = W(t) </math>. Este movimento, por ser um processo de ''Wiener'', é melhor descrito pelo deslocamento do que pela velocidade.

===Ângulo e Ruídos===
A direção de polarização é definida pelo ângulo teta que também apresenta ruído branco. Esse ruído <math>\beta(t)</math> é relacionado à <math> \xi_{\perp}(t)</math> através de <math>\xi_{perp}(t) = \sqrt(q)\beta(t)</math>. É essa mudança na direção de polarização que causa o deslocamento aleatório perpendicular.

Embora o ruído no ângulo e o ruído na direção perpendicular estejam relacionados, os mesmos são independentes do ruído na direção paralela. Abaixo pode-se ver a correlação dos mesmos. Os três são ruídos gaussianos, então a média deles é zero.

<math> <\xi_{||}(t)\xi_{||}(t')> = g\delta(t-t') </math>
<math> <\beta(t)\beta(t')> = 2k\delta(t-t') </math>
<math> <\xi_{\perp}(t)\xi_{\perp}(t')> = 2qk\delta(t-t') </math>

Onde <math>k</math> tem unidades de 1/tempo, <math>\sqrt q</math> tem unidades de comprimento, e tanto <math>g</math> quanto <math>qk</math> tem unidades de [comprimento]^2 / tempo.

Em resumo, este modelo considera uma partícula com dois graus de liberdades espaciais e um grau de liberdade interna, que considera a polarização interna. Esse grau de liberdade interna quebra a simetria espacial de tal forma que a dinâmica na direção da polarização é semelhante a ''Langevin'' e na direção perpendicular a dinâmica segue um processo de ''Wiener''. Existem duas fontes independentes de ruído: Uma age na dinâmica na direção de polarização e a outra age na polarização em si. Além da força de dissipação, a mudança na polarização age como mais um termo na perda de energia cinética.

===Mean Square Displacement(MSD)===
O MSD é uma quantificação do movimento líquido das partículas. O mesmo foi utilizado para medir a aglomeração e a transição de fase das partículas. Para isso foi feita uma média do MSD de 100 trajetórias sendo que o Método das Janelas foi feito para cada uma das mesmas. Através do resultados analíticos para as equações de movimento pode-se obter um resultado analítico para o MSD também:

<math><r^{2}> = \frac{4DP}{1-S}[\tau - (1-S)(1-e^{-\tau})]</math>

Onde <math>P</math> é o tempo de persistência do movimento em um estado cinético, <math>S</math> é uma fração deste tempo, relacionado a quanto tempo o movimento persiste no primeiro estado cinético e <math>D =\frac{g}{2m^{2}\tau\tau_{b}}</math>, sendo que <math>\tau = \frac{\gamma}{m} + k </math> e <math>\tau_{b} = \frac{\gamma}{m} + 2k </math>.

Grupo - Solução Exata para Movimento Anisotrópico de Ornstein-Uhlenbeck

2020-01-10T00:52:00Z

Almeidajuliano:

==Introdução==

Movimento coletivo ordenado de diversos organismos é um fenômeno que pode ser observado em diversas situações, como em movimento de cardumes, migração de pássaros, migração celular, etc..

Aquilo que é mais curioso destes tipos de movimento é que o mesmo surge através de movimentos aleatórios dos indivíduos que o compõem, ou seja, o movimento individual e a aleatório de um grupo de células quando postos a interagir entre si leva a um movimento ordenado em curtas escalas de tempo.

Um dos modelos anteriores a este, o modelo de ''Vicsek'', visa descrever o movimento coletivo ordenado de partículas “inteligentes”, isto é, partículas que podem observar as outras partículas na vizinhança e mover-se de acordo com a média do movimento do grupo. Um bom exemplo deste tipo de movimento é o movimento de pássaros ou peixes, da onde movimentos ordenados podem emergir de movimentos aleatórios dos indivíduos autopropulsionados. Curiosamente, mesmo que estes sistemas estejam longe de um estado de equilíbrio, os mesmos apresentam características típicas de sistemas no equilíbrio, como transições de fase bem definidas.

Contudo, partículas “não inteligentes” já haviam demonstrado o mesmo tipo de movimento coletivo experimentalmente, com as mesmas mudanças de fase bem definidas. É para explicar isso, que o modelo de ''Szabó'' foi criado. Nele, o mesmo fenômeno é descrito utilizando somente as forças de interação intracelulares já que isso é algo que independe da capacidade da célula de poder observar a média do movimento da vizinhança, como era em ''Vicsek''.

Essas forças são repulsivas para distâncias (<math>d_{ij} = |r_{i}-r_{j}|</math>) menores que <math>R_{eq}</math>, sendo <math>F_{rep}</math> o seu máximo para esse caso, atrativa para <math>R_{eq} \le d_{ij} \le R_{o}</math>, sendo <math>F_{adh}</math> o máximo para esse caso, e zero para <math>d_{ij}</math> maior que um certo <math>R_{o}</math>.

<math>
\textbf{F}(\textbf{r}_{i}, \textbf{r}_{j} ) = \textbf{e}_{ij}\begin{cases}
F_{rep}\frac{d_{ij} - R_{eq}}{R_{eq}},\quad d_{ij} < R_{eq} \\
F_{adh}\frac{d_{ij} - R_{eq}}{R_{0} - Re_{eq}}, \; \; R_{eq} \le d_{ij} \le R_{0}\\
0, \quad \quad \quad \quad \quad \; d_{ij} > R_{0}\\
\end{cases}
</math>

Essas forças levam às seguintes equações de movimento:

<math>\frac{dr_{i}(t)}{dt} = v_{o}n_{i}(t) + \mu \sum_{j=1}^N F(r_{i}, r{j})</math>

Onde <math>N</math> é o número de partículas e <math>\sum_{j=1}^N F(r_{i}, r{j})</math> é o somatório da todas as interações sobre uma partícula ''i''.

<math> \frac{d_{i}^{n}(t)}{dt} = \frac{1}{\tau}arcsin[(n_{i}(t)\frac{v_{i}(t)}{|v_{i}(t)|})e_{z}] + \xi </math>

Onde <math>\tau</math> é o tempo de relaxamento, <math>n_{i}(t)</math> é o vetor na direção de polarização, tal que <math>n_{i}(t) = (cos(\theta_{i}^{n}(t)), sen(\theta_{i}^{n}(t)))</math>, <math>e_{z}</math> é o vetor normal, e <math>\xi</math> é um ruído branco.

Neste modelo, o movimento bidimensional das células é descrito por uma dinâmica superamortecida. Na mesma, cada célula tenta manter uma velocidade de magnitude <math>v_{o}</math> na direção do vetor <math>n_{i}(t)</math>, sendo que a direção desse vetor é descrita pelo ângulo <math>\theta_{i}^{n}(t)</math> que tenta relaxar na direção da velocidade enquanto sofre ruído gaussiano.

Embora descreva bem este movimento, o modelo de Szabó fornece somente soluções numéricas para o problema. Para obtermos uma solução analítica, o problema precisa ser descrito de outra forma. É com esse objetivo que o modelo de Rita et al. foi feito.

==Modelo==
Neste modelo, em vez de ter a relaxação da velocidade na direção da polarização, a mesma é projetada sobre a polarização, enquanto ela muda com ruído. Essa direção de polarização interna pode representar, p.e., a polarização dos filamentos de actina na direção preferencial em migração celular.

Com isso pode-se observar duas transições, ou seja, uma transição de um estado desordenado para um estado ordenado e para um estado desordenado de novo. Essas duas transições representam a aleatoriedade do movimento inicial das células, depois a escolha de uma direção preferencial e por último a variação desta direção preferencial.

Contudo, neste modelo deve-se assumir que, para a partícula, seu movimento é anisotrópico, ou seja, apresenta assimetria nas dinâmicas de movimento em direções diferentes, uma na direção paralela à direção de polarização e outra perpendicular.

===Movimento Paralelo===
Na direção paralela à polarização, temos o sistema agindo numa dinâmica semelhante a de ''Langevin'' (<math>m\mathbf{a} = m\frac{d\mathbf{v}}{dt} = F(\mathbf{x}) - \beta \mathbf{v} + \eta(t)</math>), mais especificamente, um processo de ''Orstein-Uhlenbeck'' amortecido. Nesta direção a velocidade (paralela) é bem definida. Tanto para o movimento paralelo, quanto para o movimento perpendicular, o vetor <math>p_{i}(t) = (cos(\theta_{i}(t)), sen(\theta_{i}(t)))</math> é o vetor que indica a direção de polarização e <math>n_{i}(t) = (sen(\theta_{i}(t)), -cos(\theta_{i}(t)))</math> é o vetor perpendicular a <math>p_{i}(t)</math>.

<math>\frac{dv_{||}}{dt} = -\frac{\gamma}{m}v_{||} + \frac{\xi_{||}(t)}{m} </math>
<math>v_{||}(t+\Delta t)p(t+\Delta t) = [(1-\gamma\Delta t)v_{||}(t) + \int_{t}^{t+\Delta t} \xi_{||}(t)dt](p_{i}(t)p_{i}(t+\Delta t))p_{i}(t+\Delta t) </math>

Onde <math>v_{||}</math> é a velocidade paralela, <math>\gamma</math> é o coeficiente de atrito e <math>\xi_{||}</math> é o ruído gaussiano paralelo.

A velocidade diminui por causa do atrito com o substrato, já que pode-se observar que <math>(1-\gamma\Delta t)</math> tende a <math>e^{-\gamma\Delta t}</math> se <math>\gamma</math> for pequeno, e diminui ou aumenta por causa do termo estocástico <math>\int_{t}^{t+\Delta t} \xi_{||}(t)dt = W(t)</math>. A primeira equação assume que durante o intervalo <math>\Delta t</math>, a direção de polarização é constante. Somente em <math>t+\Delta t</math>, a velocidade na direção <math>p_{i}(t)</math> é projetada sobre a nova direção de polarização <math>p_{i}(t+\Delta t)</math>.

===Movimento Perpendicular===
Na direção perpendicular à polarização, o sistema segue um processo de ''Wiener'' (movimento aleatório) por causa das flutuações na direção de polarização.

<math>\frac{dr_{\perp}}{dt} = n_{i}(t)\xi_{\perp}(t) </math>
<math>[r_{\perp}(t+\Delta t) - r_{\perp}(t)]n_{i}(t) = n_{i}(t)\int_{t}^{t+\Delta t} \xi_{\perp}(t)dt </math>

Onde <math>\xi_{\perp}</math> é o ruído gaussiano na direção perpendicular. Pode-se observar que <math>\int_{t}^{t+\Delta t} \xi_{\perp}(t)dt = W(t) </math>. Este movimento, por ser um processo de ''Wiener'', é melhor descrito pelo deslocamento do que pela velocidade.

===Ângulo e Ruídos===
A direção de polarização é definida pelo ângulo teta que também apresenta ruído branco. Esse ruído <math>\beta(t)</math> é relacionado à <math> \xi_{\perp}(t)</math> através de <math>\xi_{perp}(t) = \sqrt(q)\beta(t)</math>. É essa mudança na direção de polarização que causa o deslocamento aleatório perpendicular.

Embora o ruído no ângulo e o ruído na direção perpendicular estejam relacionados, os mesmos são independentes do ruído na direção paralela. Abaixo pode-se ver a correlação dos mesmos. Os três são ruídos gaussianos, então a média deles é zero.

<math> <\xi_{||}(t)\xi_{||}(t')> = g\delta(t-t') </math>
<math> <\beta(t)\beta(t')> = 2k\delta(t-t') </math>
<math> <\xi_{\perp}(t)\xi_{\perp}(t')> = 2qk\delta(t-t') </math>

Onde <math>k</math> tem unidades de 1/tempo, <math>\sqrt q</math> tem unidades de comprimento, e tanto <math>g</math> quanto <math>qk</math> tem unidades de [comprimento]^2 / tempo.

Em resumo, este modelo considera uma partícula com dois graus de liberdades espaciais e um grau de liberdade interna, que considera a polarização interna. Esse grau de liberdade interna quebra a simetria espacial de tal forma que a dinâmica na direção da polarização é semelhante a ''Langevin'' e na direção perpendicular a dinâmica segue um processo de ''Wiener''. Existem duas fontes independentes de ruído: Uma age na dinâmica na direção de polarização e a outra age na polarização em si. Além da força de dissipação, a mudança na polarização age como mais um termo na perda de energia cinética.

===Mean Square Displacement(MSD)===
O MSD é uma quantificação do movimento líquido das partículas. O mesmo foi utilizado para medir a aglomeração e a transição de fase das partículas. Para isso foi feita uma média do MSD de 100 trajetórias sendo que o Método das Janelas foi feito para cada uma das mesmas. Através do resultados analíticos para as equações de movimento pode-se obter um resultado analítico para o MSD também:

<math><r^{2}> = \frac{4DP}{1-S}[\tau - (1-S)(1-e^{-\tau})]</math>

Grupo - Solução Exata para Movimento Anisotrópico de Ornstein-Uhlenbeck

2020-01-10T00:43:13Z

Almeidajuliano:

==Introdução==

Movimento coletivo ordenado de diversos organismos é um fenômeno que pode ser observado em diversas situações, como em movimento de cardumes, migração de pássaros, migração celular, etc..

Aquilo que é mais curioso destes tipos de movimento é que o mesmo surge através de movimentos aleatórios dos indivíduos que o compõem, ou seja, o movimento individual e a aleatório de um grupo de células quando postos a interagir entre si leva a um movimento ordenado em curtas escalas de tempo.

Um dos modelos anteriores a este, o modelo de ''Vicsek'', visa descrever o movimento coletivo ordenado de partículas “inteligentes”, isto é, partículas que podem observar as outras partículas na vizinhança e mover-se de acordo com a média do movimento do grupo. Um bom exemplo deste tipo de movimento é o movimento de pássaros ou peixes, da onde movimentos ordenados podem emergir de movimentos aleatórios dos indivíduos autopropulsionados. Curiosamente, mesmo que estes sistemas estejam longe de um estado de equilíbrio, os mesmos apresentam características típicas de sistemas no equilíbrio, como transições de fase bem definidas.

Contudo, partículas “não inteligentes” já haviam demonstrado o mesmo tipo de movimento coletivo experimentalmente, com as mesmas mudanças de fase bem definidas. É para explicar isso, que o modelo de ''Szabó'' foi criado. Nele, o mesmo fenômeno é descrito utilizando somente as forças de interação intracelulares já que isso é algo que independe da capacidade da célula de poder observar a média do movimento da vizinhança, como era em ''Vicsek''.

Essas forças são repulsivas para distâncias (<math>d_{ij} = |r_{i}-r_{j}|</math>) menores que <math>R_{eq}</math>, sendo <math>F_{rep}</math> o seu máximo para esse caso, atrativa para <math>R_{eq} \le d_{ij} \le R_{o}</math>, sendo <math>F_{adh}</math> o máximo para esse caso, e zero para <math>d_{ij}</math> maior que um certo <math>R_{o}</math>.

<math>
\textbf{F}(\textbf{r}_{i}, \textbf{r}_{j} ) = \textbf{e}_{ij}\begin{cases}
F_{rep}\frac{d_{ij} - R_{eq}}{R_{eq}},\quad d_{ij} < R_{eq} \\
F_{adh}\frac{d_{ij} - R_{eq}}{R_{0} - Re_{eq}}, \; \; R_{eq} \le d_{ij} \le R_{0}\\
0, \quad \quad \quad \quad \quad \; d_{ij} > R_{0}\\
\end{cases}
</math>

Essas forças levam às seguintes equações de movimento:

<math>\frac{dr_{i}(t)}{dt} = v_{o}n_{i}(t) + \mu \sum_{j=1}^N F(r_{i}, r{j})</math>

Onde <math>N</math> é o número de partículas e <math>\sum_{j=1}^N F(r_{i}, r{j})</math> é o somatório da todas as interações sobre uma partícula ''i''.

<math> \frac{d_{i}^{n}(t)}{dt} = \frac{1}{\tau}arcsin[(n_{i}(t)\frac{v_{i}(t)}{|v_{i}(t)|})e_{z}] + \xi </math>

Onde <math>\tau</math> é o tempo de relaxamento, <math>n_{i}(t)</math> é o vetor na direção de polarização, tal que <math>n_{i}(t) = (cos(\theta_{i}^{n}(t)), sen(\theta_{i}^{n}(t)))</math>, <math>e_{z}</math> é o vetor normal, e <math>\xi</math> é um ruído branco.

Neste modelo, o movimento bidimensional das células é descrito por uma dinâmica superamortecida. Na mesma, cada célula tenta manter uma velocidade de magnitude <math>v_{o}</math> na direção do vetor <math>n_{i}(t)</math>, sendo que a direção desse vetor é descrita pelo ângulo <math>\theta_{i}^{n}(t)</math> que tenta relaxar na direção da velocidade enquanto sofre ruído gaussiano.

Embora descreva bem este movimento, o modelo de Szabó fornece somente soluções numéricas para o problema. Para obtermos uma solução analítica, o problema precisa ser descrito de outra forma. É com esse objetivo que o modelo de Rita et al. foi feito.

==Modelo==
Neste modelo, em vez de ter a relaxação da velocidade na direção da polarização, a mesma é projetada sobre a polarização, enquanto ela muda com ruído. Essa direção de polarização interna pode representar, p.e., a polarização dos filamentos de actina na direção preferencial em migração celular.

Com isso pode-se observar duas transições, ou seja, uma transição de um estado desordenado para um estado ordenado e para um estado desordenado de novo. Essas duas transições representam a aleatoriedade do movimento inicial das células, depois a escolha de uma direção preferencial e por último a variação desta direção preferencial.

Contudo, neste modelo deve-se assumir que, para a partícula, seu movimento é anisotrópico, ou seja, apresenta assimetria nas dinâmicas de movimento em direções diferentes, uma na direção paralela à direção de polarização e outra perpendicular.

===Movimento Paralelo===
Na direção paralela à polarização, temos o sistema agindo numa dinâmica semelhante a de ''Langevin'' (<math>m\mathbf{a} = m\frac{d\mathbf{v}}{dt} = F(\mathbf{x}) - \beta \mathbf{v} + \eta(t)</math>), mais especificamente, um processo de ''Orstein-Uhlenbeck'' amortecido. Nesta direção a velocidade (paralela) é bem definida. Tanto para o movimento paralelo, quanto para o movimento perpendicular, o vetor <math>p_{i}(t) = (cos(\theta_{i}(t)), sen(\theta_{i}(t)))</math> é o vetor que indica a direção de polarização e <math>n_{i}(t) = (sen(\theta_{i}(t)), -cos(\theta_{i}(t)))</math> é o vetor perpendicular a <math>p_{i}(t)</math>.

<math>\frac{dv_{||}}{dt} = -\frac{\gamma}{m}v_{||} + \frac{\xi_{||}(t)}{m} </math>
<math>v_{||}(t+\Delta t)p(t+\Delta t) = [(1-\gamma\Delta t)v_{||}(t) + \int_{t}^{t+\Delta t} \xi_{||}(t)dt](p_{i}(t)p_{i}(t+\Delta t))p_{i}(t+\Delta t) </math>

Onde <math>v_{||}</math> é a velocidade paralela, <math>\gamma</math> é o coeficiente de atrito e <math>\xi_{||}</math> é o ruído gaussiano paralelo.

A velocidade diminui por causa do atrito com o substrato, já que pode-se observar que <math>(1-\gamma\Delta t)</math> tende a <math>e^{-\gamma\Delta t}</math> se <math>\gamma</math> for pequeno, e diminui ou aumenta por causa do termo estocástico <math>\int_{t}^{t+\Delta t} \xi_{||}(t)dt = W(t)</math>. A primeira equação assume que durante o intervalo <math>\Delta t</math>, a direção de polarização é constante. Somente em <math>t+\Delta t</math>, a velocidade na direção <math>p_{i}(t)</math> é projetada sobre a nova direção de polarização <math>p_{i}(t+\Delta t)</math>.

===Movimento Perpendicular===
Na direção perpendicular à polarização, o sistema segue um processo de ''Wiener'' (movimento aleatório) por causa das flutuações na direção de polarização.

<math>[r_{\perp}(t+\Delta t) - r_{\perp}(t)]n_{i}(t) = n_{i}(t)\int_{t}^{t+\Delta t} \xi_{\perp}(t)dt </math>
<math>\frac{dr_{\perp}}{dt} = n_{i}(t)\xi_{\perp}(t) </math>

Onde <math>\xi_{\perp}</math> é o ruído gaussiano na direção perpendicular. Pode-se observar que <math>\int_{t}^{t+\Delta t} \xi_{\perp}(t)dt = W(t) </math>. Este movimento, por ser um processo de ''Wiener'', é melhor descrito pelo deslocamento do que pela velocidade.

===Ângulo e Ruídos===
A direção de polarização é definida pelo ângulo teta que também apresenta ruído branco. Esse ruído <math>\beta(t)</math> é relacionado à <math> \xi_{\perp}(t)</math> através de <math>\xi_{perp}(t) = \sqrt(q)\beta(t)</math>. É essa mudança na direção de polarização que causa o deslocamento aleatório perpendicular.

Embora o ruído no ângulo e o ruído na direção perpendicular estejam relacionados, os mesmos são independentes do ruído na direção paralela. Abaixo pode-se ver a correlação dos mesmos. Os três são ruídos gaussianos, então a média deles é zero.

<math> <\xi_{||}(t)\xi_{||}(t')> = g\delta(t-t') </math>
<math> <\beta(t)\beta(t')> = 2k\delta(t-t') </math>
<math> <\xi_{\perp}(t)\xi_{\perp}(t')> = 2qk\delta(t-t') </math>

Onde <math>k</math> tem unidades de 1/tempo, <math>\sqrt q</math> tem unidades de comprimento, e tanto <math>g</math> quanto <math>qk</math> tem unidades de [comprimento]^2 / tempo.

Em resumo, este modelo considera uma partícula com dois graus de liberdades espaciais e um grau de liberdade interna, que considera a polarização interna. Esse grau de liberdade interna quebra a simetria espacial de tal forma que a dinâmica na direção da polarização é semelhante a ''Langevin'' e na direção perpendicular a dinâmica segue um processo de ''Wiener''. Existem duas fontes independentes de ruído: Uma age na dinâmica na direção de polarização e a outra age na polarização em si. Além da força de dissipação, a mudança na polarização age como mais um termo na perda de energia cinética.

===Mean Square Displacement(MSD)===
O MSD é uma quantificação do movimento líquido das partículas. O mesmo foi utilizado para medir a aglomeração e a transição de fase das partículas. Para isso foi feita uma média do MSD de 100 trajetórias sendo que o Método das Janelas foi feito para cada uma das mesmas.

Grupo - Solução Exata para Movimento Anisotrópico de Ornstein-Uhlenbeck

2020-01-10T00:35:30Z

Almeidajuliano:

==Introdução==

Movimento coletivo ordenado de diversos organismos é um fenômeno que pode ser observado em diversas situações, como em movimento de cardumes, migração de pássaros, migração celular, etc..

Aquilo que é mais curioso destes tipos de movimento é que o mesmo surge através de movimentos aleatórios dos indivíduos que o compõem, ou seja, o movimento individual e a aleatório de um grupo de células quando postos a interagir entre si leva a um movimento ordenado em curtas escalas de tempo.

Um dos modelos anteriores a este, o modelo de ''Vicsek'', visa descrever o movimento coletivo ordenado de partículas “inteligentes”, isto é, partículas que podem observar as outras partículas na vizinhança e mover-se de acordo com a média do movimento do grupo. Um bom exemplo deste tipo de movimento é o movimento de pássaros ou peixes, da onde movimentos ordenados podem emergir de movimentos aleatórios dos indivíduos autopropulsionados. Curiosamente, mesmo que estes sistemas estejam longe de um estado de equilíbrio, os mesmos apresentam características típicas de sistemas no equilíbrio, como transições de fase bem definidas.

Contudo, partículas “não inteligentes” já haviam demonstrado o mesmo tipo de movimento coletivo experimentalmente, com as mesmas mudanças de fase bem definidas. É para explicar isso, que o modelo de ''Szabó'' foi criado. Nele, o mesmo fenômeno é descrito utilizando somente as forças de interação intracelulares já que isso é algo que independe da capacidade da célula de poder observar a média do movimento da vizinhança, como era em ''Vicsek''.

Essas forças são repulsivas para distâncias (<math>d_{ij} = |r_{i}-r_{j}|</math>) menores que <math>R_{eq}</math>, sendo <math>F_{rep}</math> o seu máximo para esse caso, atrativa para <math>R_{eq} \le d_{ij} \le R_{o}</math>, sendo <math>F_{adh}</math> o máximo para esse caso, e zero para <math>d_{ij}</math> maior que um certo <math>R_{o}</math>.

<math>
\textbf{F}(\textbf{r}_{i}, \textbf{r}_{j} ) = \textbf{e}_{ij}\begin{cases}
F_{rep}\frac{d_{ij} - R_{eq}}{R_{eq}},\quad d_{ij} < R_{eq} \\
F_{adh}\frac{d_{ij} - R_{eq}}{R_{0} - Re_{eq}}, \; \; R_{eq} \le d_{ij} \le R_{0}\\
0, \quad \quad \quad \quad \quad \; d_{ij} > R_{0}\\
\end{cases}
</math>

Essas forças levam às seguintes equações de movimento:

<math>\frac{dr_{i}(t)}{dt} = v_{o}n_{i}(t) + \mu \sum_{j=1}^N F(r_{i}, r{j})</math>

Onde <math>N</math> é o número de partículas e <math>\sum_{j=1}^N F(r_{i}, r{j})</math> é o somatório da todas as interações sobre uma partícula ''i''.

<math> \frac{d_{i}^{n}(t)}{dt} = \frac{1}{\tau}arcsin[(n_{i}(t)\frac{v_{i}(t)}{|v_{i}(t)|})e_{z}] + \xi </math>

Onde <math>\tau</math> é o tempo de relaxamento, <math>n_{i}(t)</math> é o vetor na direção de polarização, tal que <math>n_{i}(t) = (cos(\theta_{i}^{n}(t)), sen(\theta_{i}^{n}(t)))</math>, <math>e_{z}</math> é o vetor normal, e <math>\xi</math> é um ruído branco.

Neste modelo, o movimento bidimensional das células é descrito por uma dinâmica superamortecida. Na mesma, cada célula tenta manter uma velocidade de magnitude <math>v_{o}</math> na direção do vetor <math>n_{i}(t)</math>, sendo que a direção desse vetor é descrita pelo ângulo <math>\theta_{i}^{n}(t)</math> que tenta relaxar na direção da velocidade enquanto sofre ruído gaussiano.

Embora descreva bem este movimento, o modelo de Szabó fornece somente soluções numéricas para o problema. Para obtermos uma solução analítica, o problema precisa ser descrito de outra forma. É com esse objetivo que o modelo de Rita et al. foi feito.

==Modelo==
Neste modelo, em vez de ter a relaxação da velocidade na direção da polarização, a mesma é projetada sobre a polarização, enquanto ela muda com ruído. Essa direção de polarização interna pode representar, p.e., a polarização dos filamentos de actina na direção preferencial em migração celular.

Com isso pode-se observar duas transições, ou seja, uma transição de um estado desordenado para um estado ordenado e para um estado desordenado de novo. Essas duas transições representam a aleatoriedade do movimento inicial das células, depois a escolha de uma direção preferencial e por último a variação desta direção preferencial.

Contudo, neste modelo deve-se assumir que, para a partícula, seu movimento é anisotrópico, ou seja, apresenta assimetria nas dinâmicas de movimento em direções diferentes, uma na direção paralela à direção de polarização e outra perpendicular.

===Movimento Paralelo===
Na direção paralela à polarização, temos o sistema agindo numa dinâmica semelhante a de ''Langevin'' (<math>m\mathbf{a} = m\frac{d\mathbf{v}}{dt} = F(\mathbf{x}) - \beta \mathbf{v} + \eta(t)</math>), mais especificamente, um processo de ''Orstein-Uhlenbeck'' amortecido. Nesta direção a velocidade (paralela) é bem definida. Tanto para o movimento paralelo, quanto para o movimento perpendicular, o vetor <math>p_{i}(t) = (cos(\theta_{i}(t)), sen(\theta_{i}(t)))</math> é o vetor que indica a direção de polarização e <math>n_{i}(t) = (sen(\theta_{i}(t)), -cos(\theta_{i}(t)))</math> é o vetor perpendicular a <math>p_{i}(t)</math>.

<math>\frac{dv_{||}}{dt} = -\frac{\gamma}{m}v_{||} + \frac{\xi_{||}(t)}{m} </math>
<math>v_{||}(t+\Delta t)p(t+\Delta t) = [(1-\gamma\Delta t)v_{||}(t) + \int_{t}^{t+\Delta t} \xi_{||}(t)dt](p_{i}(t)p_{i}(t+\Delta t))p_{i}(t+\Delta t) </math>

Onde <math>v_{||}</math> é a velocidade paralela, <math>\gamma</math> é o coeficiente de atrito e <math>\xi_{||}</math> é o ruído gaussiano paralelo.

A velocidade diminui por causa do atrito com o substrato, já que pode-se observar que <math>(1-\gamma\Delta t)</math> tende a <math>e^{-\gamma\Delta t}</math> se <math>\gamma</math> for pequeno, e diminui ou aumenta por causa do termo estocástico <math>\int_{t}^{t+\Delta t} \xi_{||}(t)dt = W(t)</math>. A primeira equação assume que durante o intervalo <math>\Delta t</math>, a direção de polarização é constante. Somente em <math>t+\Delta t</math>, a velocidade na direção <math>p_{i}(t)</math> é projetada sobre a nova direção de polarização <math>p_{i}(t+\Delta t)</math>.

===Movimento Perpendicular===
Na direção perpendicular à polarização, o sistema segue um processo de ''Wiener'' (movimento aleatório) por causa das flutuações na direção de polarização.

<math>[r_{\perp}(t+\Delta t) - r_{\perp}(t)]n_{i}(t)\int_{t}^{t+\Delta t} \xi_{\perp}(t)dt </math>
<math>\frac{dr_{\perp}}{dt} = n_{i}(t)\xi_{\perp}(t) </math>

Pode-se observar que <math>\int_{t}^{t+\Delta t} \xi_{\perp}(t)dt = W(t) </math>. Este movimento, por ser um processo de ''Wiener'', é melhor descrito pelo deslocamento do que pela velocidade.

===Ângulo e Ruídos===
A direção de polarização é definida pelo ângulo teta que também apresenta ruído branco. Esse ruído <math>\beta(t)</math> é relacionado à <math> \xi_{\perp}(t)</math> através de <math>\xi_{perp}(t) = \sqrt(q)\beta(t)</math>. É essa mudança na direção de polarização que causa o deslocamento aleatório perpendicular.

Embora o ruído no ângulo e o ruído na direção perpendicular estejam relacionados, os mesmos são independentes do ruído na direção paralela. Abaixo pode-se ver a correlação dos mesmos. Os três são ruídos gaussianos, então a média deles é zero.

<math> <\xi_{||}(t)\xi_{||}(t')> = g\delta(t-t') </math>
<math> <\beta(t)\beta(t')> = 2k\delta(t-t') </math>
<math> <\xi_{\perp}(t)\xi_{\perp}(t')> = 2qk\delta(t-t') </math>

Onde <math>k</math> tem unidades de 1/tempo, <math>\sqrt q</math> tem unidades de comprimento, e tanto <math>g</math> quanto <math>qk</math> tem unidades de [comprimento]^2 / tempo.

Em resumo, este modelo considera uma partícula com dois graus de liberdades espaciais e um grau de liberdade interna, que considera a polarização interna. Esse grau de liberdade interna quebra a simetria espacial de tal forma que a dinâmica na direção da polarização é semelhante a ''Langevin'' e na direção perpendicular a dinâmica segue um processo de ''Wiener''. Existem duas fontes independentes de ruído: Uma age na dinâmica na direção de polarização e a outra age na polarização em si. Além da força de dissipação, a mudança na polarização age como mais um termo na perda de energia cinética.

Grupo - Solução Exata para Movimento Anisotrópico de Ornstein-Uhlenbeck

2020-01-10T00:09:49Z

Almeidajuliano:

Grupo - Solução Exata para Movimento Anisotrópico de Ornstein-Uhlenbeck

2020-01-09T23:41:05Z

Almeidajuliano:

Grupo - Solução Exata para Movimento Anisotrópico de Ornstein-Uhlenbeck

2020-01-09T23:37:19Z

Almeidajuliano:

Grupo - Solução Exata para Movimento Anisotrópico de Ornstein-Uhlenbeck

2020-01-09T23:30:56Z

Almeidajuliano:

Grupo - Solução Exata para Movimento Anisotrópico de Ornstein-Uhlenbeck

2020-01-09T23:30:09Z

Almeidajuliano:

Grupo - Solução Exata para Movimento Anisotrópico de Ornstein-Uhlenbeck

2020-01-09T23:28:04Z

Almeidajuliano:

Grupo - Solução Exata para Movimento Anisotrópico de Ornstein-Uhlenbeck

2020-01-09T23:05:00Z

Almeidajuliano: Criou página com '==Introdução== Movimento coletivo ordenado de diversos organismos é um fenômeno que pode ser observado em diversas situações, como em movimento de cardumes, migração...'

Métodos computacionais

2020-01-09T22:46:43Z

Almeidajuliano: /* Métodos Computacionais C */

Física computacional é uma abordagem da '''física teórica''' com o auxílio do computador essencialmente
quando a complexidade do problema impossibilita o avanço pela via analítica e/ou porque os cálculos
numéricos são longos demais para serem feitos sem automação.
Alguns consideram a '''física computacional''' um terceiro (e mais recente) vértice do triângulo da maneira de se fazer física, onde os outros dois vértices são a '''física teórica''' e a '''física experimental'''.

'''Métodos computacionais''' é a disciplina onde estudamos ferramentas, métodos e algoritmos numéricos para a resolução de problemas de física onde uma abordagem analítica é extremamente complexa ou impossível.

Alguns exemplos de aplicação são: a solução numérica de equações diferenciais ordinárias, integração numérica via métodos de aproximação ou estatísticos como método de Monte Carlo, equações diferencias parciais como as equações de Maxwell e de Schroedinger, métodos matriciais para a solução de problemas de autovalor e autovetor como os encontrados na Mecânica Quântica.

==Breve Historia da Computação==

De Conrad Zuse (1941) ao IBM Blue/Gene (2006)

==Arquitectura==

[http://en.wikipedia.org/wiki/Image:Personal_computer%2C_exploded_5.svg Diagrama de PC]

== Ferramentas ==
===== [[Comandos Unix/Linux]] =====
===== [[Gnuplot]] e [[xmgrace]] =====
===== [[LaTex]] =====
===== [[FORTRAN]] =====
===== [[C]] =====
===== [[Julia]] =====

== Métodos Computacionais A ==
===== [[Derivada Numérica]] =====
===== [[Integração Numérica]] =====
===== [[Interpolação e extrapolação]] =====
====== [[Fórmula de Lagrange]]======
====== [[Spline cúbico]]======

===== [[Zeros de Funções]]=====
===== [[Mínimos Quadrados]] =====
===== [[Listas de exercícios]] =====
====== [[Área 1]] ======
====== [[Área 2]] ======
====== [[Área 3]] ======

== Métodos Computacionais B ==

===== [[Integração numérica de equações diferenciais ordinárias]] =====

===== [[Métodos multipassos]]=====

===== [[Métodos de passo variável]]=====

===== [[Aplicações (Mapas)]]=====

===== [[Números Aleatórios]]=====

===== [[Histogramas e Densidade de Probabilidade ]] =====

===== [[Método de Monte Carlo (integração numérica)]]=====

== Métodos Computacionais C ==

===== [[Equações Diferenciais Parciais]] =====

=====[[Grupo1 - Dif em 2D]] =====

=====[[Grupo2 - Ondas1]]=====

=====[[Grupo3 - Ondas2]]=====

=====[[Grupo4 - FFT]]=====

=====[[Grupo5 - Eq. Schroedinger]]=====

=====[[Grupo - Ising 2D]]=====

===== [[Monte Carlo]] =====

=====[[Grupo - Lennard Jones]]=====

=====[[Grupo - BOIDS]]=====

=====[[Grupo - Tráfego]]=====

===== [[Teste_conv]] =====

===== [[Teste2]] =====

===== [[Grupo - Dilema Do Prisioneiro]] =====

===== [[Grupo - Modelo Sznajd]] =====

=====[[Grupo - Modelo de Potts]]=====

=====[[Grupo - Conservação do Parâmetro de Ordem]]=====

=====[[ Movimento Coletivo ]] =====

=====[[Grupo - Modelo de Szabó]]=====

=====[[ Ressonância Estocástica ]] =====

=====[[Grupo - Equações de Schrödinger não-lineares acopladas]]=====

=====[[Grupo - Correlações no Movimento de Átomos em Argônio]]=====

=====[[Grupo - Solução Exata para Movimento Anisotrópico de Ornstein-Uhlenbeck]] =====