Sobre vagalumes, pedestres e neurônios: a sincronização de osciladores de fase

About fireflies, pedestrians and neurons: the synchronization of phase oscillators

Resumos

Osciladores de fase são sistemas dinâmicos dissipativos com aporte externo de energia, e matematicamente caracterizados pela presença de ciclos-limite, ao longo dos quais a dinâmica pode ser descrita por um ângulo de fase. No presente artigo, consideramos fenômenos de sincronização num sistema de muitos osciladores de fase acoplados. Adotamos uma descrição qualitativa, baseada num modelo paradigmático proposto por Y. Kuramoto em 1975. Discutimos, à luz dessa descrição, algumas aplicações em sistemas físicos (metrônomos e pedestres numa ponte pênsil) e biológicos (neurônios e vagalumes piscando).

Palavras-chave
Kuramoto; osciladores; sincronização; fase

Phase oscillators are dissipative dynamic systems with external energy input, and mathematically characterized by the presence of limit cycles, along which the dynamics can be described by a phase angle. In the present article, we consider synchronization phenomena in a system of many phase-coupled oscillators. We adopt a qualitative description, based on a paradigmatic model proposed by Y. Kuramoto in 1975. We discuss, in the light of this description, some applications in physical systems (metronomes and pedestrians on a suspension bridge) and biological (neurons and fireflies blinking).

Keywords
Kuramoto; oscillators; synchronization; phase

1. Introdução

Fenômenos periódicos caracterizam-se por uma repetição sistemática após um determinado período de tempo. Dentre os primeiros fenômenos periódicos investigados pelo ser humano estão as efemérides astronômicas, como as mudanças de fases da Lua e de estações do ano, horários de nascer e do pôr do Sol, dentre outras [¹]. Uma das motivações para este estudo foi, sem dúvida, a necessidade de medição do tempo para a qual foram usados vários dispositivos ao longo da História [²].

Ao longo do século XVII houve as primeiras investigações sistemáticas de fenômenos periódicos no contexto físico. Galileo constatou o isocronismo das oscilações do pêndulo, segundo consta, usando a própria pulsação, que é também um fenômeno periódico. Posteriormente, Huygens empregou essa ideia na sua concepção original do relógio de pêndulo, usando um mecanismo de escape [³]. Desde então um dos principais ramos da Física é o estudo de fenômenos periódicos em fenômenos mecânicos, eletromagnéticos, óticos, etc. [⁴].

A sincronização de osciladores é um ajuste dos seus ritmos devido a uma interação fraca entre eles, ou seja, ela depende de algum tipo de acoplamento entre os osciladores [⁵]. Desde que Huygens descreveu a sincronização entre relógios de pêndulo, diversos exemplos de sincronização têm sido descritos em sistemas físicos [⁶, ⁷] e biológicos [⁸, ⁹, ¹⁰].

Um dos mais espetaculares fenômenos de sincronização na natureza é observado no piscar de algumas espécies de vagalumes do Sudeste Asiático [¹¹]. Milhares de vagalumes machos juntam-se em árvores à noite e piscam de forma sincronizada para atrair a atenção das fêmeas. Nesse caso, podemos identificar o efeito do acoplamento na percepção visual que os vagalumes têm ao sincronizar seu piscar com o de seus colegas. Observa-se, também, sincronização de aplausos de espectadores numa sala de concerto [¹²], bem como de metrônomos [⁶, ⁷]

Outro fenômeno notável de sincronização foi registrado na abertura da Ponte (pênsil) do Milênio, que cruza o Rio Tâmisa em Londres [¹³, ¹⁴]. No dia da sua abertura, em 2000, a ponte atraiu um enorme número de visitantes, e cerca de duas mil pessoas podiam ocupar a ponte ao mesmo tempo. O caminhar dos pedestres gerou um pequeno deslocamento lateral da ponte. Para evitar quedas, outros pedestres acabaram ajustando o próprio ritmo de sua caminhada, o que só piorava a situação causando oscilações laterais de grande amplitude devido à sincronização entre a ponte pênsil e o modo de andar dos pedestres [¹⁴]. A ponte foi, então, fechada para a instalação de amortecedores que diminuissem essa oscilação lateral.

Uma dificuldade comum à descrição da sincronização de osciladores de tipos tão diversos como vagalumes piscando e pedestres deslocando-se sobre uma ponte, está na grande diferença de modelamento desses eventos. No entanto, uma descrição qualitativa pode ser empregada, que leva em conta o fato de vagalumes, pedestres, metrônomos, neurônios, relógios de pêndulo e outros poderem – a despeito de suas diferenças – serem todos osciladores de fase. Isso é possível porque o movimento de osciladores auto-sustentáveis, dá-se ao longo de uma curva fechada no plano de fase denominada ciclo limite[¹⁵].

A sincronização de dois osciladores de fase, sob este ponto de vista, foi analisada em um trabalho prévio [¹⁶]. Mesmo sem conhecer em detalhes o tipo de acoplamento entre os osciladores, é possível obter conclusões gerais sobre suas propriedades de sincronização, usando equações médias nas vizinhanças de uma ressonância exata entre as frequências dos osciladores. Em ([¹⁶]) empregamos estes resultados para tratar da sincronização entre dois relógios de pêndulo e também de metrônomos.

No presente artigo temos como objetivo principal estender o tratamento abordado em [¹⁶] para o caso de um número arbitrário N de osciladores de fase acoplados. Esse tratamento levar-nos-á a um modelo proposto em 1975 por Y. Kuramoto, e que tem servido de paradigma para o estudo da sincronização de osciladores acoplados. Faremos uma abordagem do modelo de Kuramoto do ponto de vista descritivo, mostrando resultados de simulações numéricas e mencionando alguns resultados analíticos que podem ser encontrados na literatura especializada [¹⁷, ¹⁸].

À guisa de ilustração, mostraremos também como é possível tratar os disparos elétricos de um neurônio como um oscilador de fase; e de como sua sincronização pode ser estudada a partir do modelo de Kuramoto [¹⁹]. Dessa forma, sistemas tão diversos como vagalumes piscando, pedestres, metrônomos, relógios de pêndulo e neurônios podem ser descritos de forma unificada.

O restante deste artigo está organizado da seguinte forma: na seção 2 introduzimos a ideia de oscilador de fase, mostrando definições diferentes de fase de acordo com o tipo de oscilador estudado. Na seção 2 revemos a teoria para a sincronização de dois osciladores de fase, enfatizando o papel das equações médias na vizinhança da ressonância exata. A sincronização de três ou mais osciladores de fase é discutida , na seção 3, no contexto do modelo de Kuramoto, para o qual fazemos uma apresentação descritiva. A seção 4 é dedicada à análise dos resultados, usando para isso resultados teóricos desenvolvidos no caso de um número infinito de osciladores. Detalhes dos cálculos para esse caso são deixados para um Apêndice. A última seção é dedicada às Conclusões.

2. Osciladores de fase

Um pêndulo simples que é posto a oscilar tenderá, com o passar do tempo, ao repouso, pois a sua energia total é dissipada devido a uma série de fatores, como o atrito com o ar e a fricção do fio no ponto de suspensão. Para garantir a presença de oscilações auto-sustentadas é, portanto, necessário compensar essa perda com um aporte externo de energia. Num relógio de pêndulo, por exemplo, esse aporte é fornecido pela lenta queda de um corpo, transmitida ao pêndulo de forma intermitente por um mecanismo do tipo âncora, como proposto originalmente por Huygens [¹⁶].

Uma característica comum a oscilações auto-sustentadas é que suas trajetórias no plano de fase (ângulo versus velocidade angular para o pêndulo, por exemplo) estão limitadas a uma curva fechada e isolada, conhecida na dinâmica não-linear, como ciclo limite[¹⁵]. Esse ciclo limite deve ser atrativo, para que pequenos desvios não tirem o movimento dessa curva fechada.

O movimento ao longo do ciclo-limite (no plano de fase) pode ser descrito por meio de um ângulo de fase $\theta (t)$, que aumenta em $2\pi$ radianos se percorremos uma volta completa no ciclo limite. A equação de evolução temporal da fase é

(1) $$\frac{d\theta }{dt}=\omega ,$$

onde $\omega$ é a frequência do oscilador. Para muitos osciladores mecânicos (como metrônomos e relógios de pêndulo) e eletromagnéticos (como circuitos eletrônicos) podemos definir essa fase geometricamente [¹⁶]. No plano descrito pelas variáveis $(x,y)$ a fase é definida como

(2) $$\tan \theta (t)=\frac{y(t)-y_0}{x(t)-x_0}$$

onde $(x_0,y_0)$ é um ponto no interior do ciclo limite [Fig.1]. Essa não é a única maneira de definir uma fase, havendo outras mais apropriadas conforme o contexto no qual trabalhamos [⁵].

Figura 1
Ilustração da determinação de uma fase geométrica para o movimento sobre um ciclo limite atrativo descrito no plano de fase $x,y$, sendo $x_0,y_0$, os pontos considerados no interior do plano de fase.

Dentre os muitos exemplos de osciladores existentes no mundo real, um de extrema importância é o neurônio, que consiste em uma célula do sistema nervoso responsável pela transmissão de impulsos nervosos. Estes impulsos podem ser medidos através da membrana neuronal, e caracterizam-se por um potencial de membrana $V(t)$ (supomos que o meio extra-celular tenha potencial nulo) [²⁰]. Se o potencial de membrana é superior a um dado limiar, o neurônio “dispara” provocando um pico na série temporal de $V(t)$.

Matematicamente, o comportamento dos disparos neuronais pode ser obtido através de diferentes modelos computacionais, como por exemplo, Hodgkin-Huxley [²¹] e Izhikevich [²²]. Outro modelo amplamente utilizado para imitar o comportamento de um neurônio biológico, é o modelo fenomenológico de Rulkov [²³]. Este, é escrito na forma de um mapa bidimensional com um termo não-linear e três parâmetros; que descrevem de forma adequada, o comportamento dos disparos observados em sistemas neuronais. O mapa proposto por Rulkov é representado pelo seguinte conjunto de equações a tempos discretos $t=\mathrm{0,1,2},\mathrm{\dots }$:

(3) $$x_{t+1}={\displaystyle \frac{\alpha }{1+x_t^2}}+y_t$$

(4) $$y_{t+1}=y_t-\sigma (x_t-\rho ),$$

onde x_t e y_t representam, respectivamente, as variáveis rápida e lenta num tempo discreto t. A variável rápida está relacionada ao potencial de ação $V(t)$ do neurônio, enquanto a variável lenta produz uma modulação temporal. As constantes $\alpha$, $\sigma$ e $\rho$ têm seus valores fixados pelo tipo de comportamento neuronal que se deseja investigar. Diferentes combinações de $\sigma$, $\rho$ e $\alpha$ conferem distintos padrões de disparos neuronais, que correspondem aos tipos de sinais enviados ou recebidos pelos neurônios tais como repouso, disparos e rajadas de disparos caóticos. Tanto as variáveis x e y como os parâmetros do modelo de Rulkov são adimensionais por construção.

As Figuras2(a)–(b) representam, respectivamente, a evolução temporal das variáveis rápida e lenta do modelo de Rulkov (3)–(4), para $\alpha =\mathrm{4,1}$ e $\sigma =\rho =\mathrm{0,001}$. Observamos, na Fig.2(a), que a variável rápida apresenta rajadas de disparos (bursting), caracterizadas por uma rápida sucessão de disparos consecutivos do neurônio [²⁴]. Muitos neurônios exibem este tipo de comportamento, e o modelo de Rulkov consegue reproduzi-lo para valores adequadamente escolhidos das constantes, de sorte que é um modelo que pode ser empregado para simulações computacionais de diversos fenômenos de interesse em Neurociência [²⁵, ²⁶]. É visível o caráter periódico das rajadas de disparos, ainda que os disparos propriamente ditos ocorram de forma caótica, devido às propriedades matemáticas do modelo de Rulkov [²³].

Figura 2
Evolução temporal da (a) variável rápida (relacionada ao potencial de membrana) e (b) variável lenta para um neurônio descrito pelo modelo de Rulkov (3)–(4).

Na Fig.2(b) mostramos a evolução temporal da variável lenta (para as mesmas escolhas de parâmetros), onde é mais aparente o caráter periódico das oscilações características deste fenômeno de rajadas de disparos. De fato, os máximos locais da variável lenta (indicados pelas barras verticais vermelhas) marcam o início de cada ciclo de rajada de disparos da variável rápida. Denotando por t_k o início do k-ésimo ciclo, podemos definir uma fase $\theta (t)$ da seguinte forma

(5) $$\theta (t)=2\pi \left\{k+\frac{t-t_k}{t_{k+1}-t_k}\right\},$$

de modo que a fase $\theta$ aumenta de $2\pi$ entre os instantes de tempo t_k e t_k+1, correspondentes ao k-ésimo ciclo. A duração de cada ciclo de rajadas de disparos $\mathrm{\Delta }{t}_k=t_k-t_{n-1}$ não é rigorosamente constante, devido à dinâmica caótica experimentada pela variável rápida. Portanto, a frequência correspondente à fase definida por (5) deve ser obtida por meio de uma média temporal:

(6) $$\omega =\underset{T\to \mathrm{\infty }}{lim}\frac{\theta (t=T)-\theta (t=0)}{T},$$

onde admitimos tacitamente que este limite sempre exista nos casos de interesse. O modelo de Rulkov trata, originalmente, de um neurônio isolado. Para compreender o efeito coletivo dos neurônios, como por exemplo em uma rede, devemos acoplar ao modelo um termo de modo que possamos estudar os efeitos dinâmicos causados pela adição de N mapas. Na próxima seção, estudaremos o caso de um sistema composto por dois osciladores (modelo mais simples) para depois tratarmos de modelo de N osciladores acoplados.

3. Acoplamento de dois osciladores de fase

Em um trabalho anterior, abordamos o acoplamento de dois osciladores de fase, descritos pelos ângulos ${\theta }_1$ e ${\theta }_2$, com $-\pi <{\theta }_{\mathrm{1,2}}\le \pi$[¹⁶]. Se desacoplados, eles têm frequências ${\omega }_1={\dot{\theta }}_1$ e ${\omega }_2={\dot{\theta }}_2$. Caso as fases sejam definidas como em (5), isto é, por meio de um tipo de mapa de Poincaré, as frequências devem ser entendidas no espírito da Eq. (6), ou seja, como médias temporais.

O tipo de acoplamento entre os osciladores depende da natureza da interação entre eles. No caso de pedestres na ponte pênsil, metrônomos e relógios de pêndulo, por exemplo, as oscilações produzidas provocam perturbações mecânicas transmitidas ao longo do suporte comum a ambos. Essas perturbações influenciam, por sua vez, as oscilações dos demais osciladores produzindo um efeito coletivo, que é a sincronização. Para o caso da bioluminescência dos vagalumes, por outro lado, o acoplamento é mediado pela propagação de ondas luminosas (eletromagnéticas) que se propagam pelo ar. Em ambos os casos, a velocidade de propagação das perturbações é tão grande que os osciladores praticamente não têm efeitos de retardo [²⁷]. Já o acoplamento entre neurônios pode ocorrer tanto devido à troca direta de íons entre eles, via contato das membranas (sinapses elétricas), como devido à emissão e absorção de substâncias químicas chamadas neurotransmissores (sinapses químicas) [²⁸].

É possível contornar as complicações físicas inerentes ao tipo de acoplamento usando argumentos gerais, válidos para osciladores com frequências muito próximas (nas vizinhanças de uma ressonância). Inicialmente escreveremos as equações para os osciladores acoplados na forma [⁵]

(7) $$\frac{d{\theta }_1}{dt}={\omega }_1+K{Q}_1({\theta }_1,{\theta }_2),$$

(8) $$\frac{d{\theta }_2}{dt}={\omega }_2+K{Q}_2({\theta }_2,{\theta }_1),$$

onde Q_1,2 são funções periódicas em θ_1,2, de tal sorte que podemos expandi-las em séries de Fourier duplas, cujos coeficientes dependem das particularidades do acoplamento existente. O parâmetro K indica a intensidade do acoplamento entre os osciladores. Suporemos que o acoplamento seja sempre suficientemente fraco, de modo a valer a condição $|K|\ll 1$.

Observa-se, das equações7 e8, que uma ressonância simples entre os osciladores ocorre quando ${\omega }_1={\omega }_2$. Como as frequências naturais dependem dos parâmetros característicos de cada oscilador, a condição de ressonância só ocorre, rigorosamente falando, se os osciladores forem idênticos. Na prática, porém, pequenas diferenças entre os osciladores geram frequências que são muito próximas, na qual podemos considerar ${\omega }_1-{\omega }_2\approx 0$. Um exemplo deste pequeno descasamento de frequências ocorre em neurônios que apresentam rajadas de disparos, como ilustrado pelo modelo de Rulkov (3)–(4), quando o parâmetro $\alpha$ de cada neurônio é ligeiramente diferente.

Nas séries de Fourier que representam as funções de acoplamento Q_1,2, os termos ressonantes oscilam lentamente com o tempo na vizinhança da ressonância exata. Já os termos não-ressonantes, variam rapidamente e anulam-se quando tomada uma média no tempo. Dessa forma, chega-se às seguintes equações médias para os osciladores de fase [¹⁶]

(9) $$\frac{d{\theta }_1}{dt}={\omega }_1+K{q}_1(\phi ),$$

(10) $$\frac{d{\theta }_2}{dt}={\omega }_2+K{q}_2(-\phi ),$$

onde $\phi ={\theta }_1-{\theta }_2$ é a diferença de fase e os termos q_1,2 são médias temporais de Q_1,2.

Subtraindo as equações9 e10 obtém-se

(11) $$\frac{d\phi }{dt}=-\mathrm{\Delta }\omega +Kq(\phi ),$$

onde $\mathrm{\Delta }\omega ={\omega }_2-{\omega }_1$ é a diferença entre frequências e definimos $q(\phi )=q_1(\phi )-q_2(-\phi )$, tal que $q(\phi )=-q(-\phi )$. Uma função simples que é periódica, suave e anti-simétrica é: $q(\phi )=\mathrm{sen}\phi$, de modo que (11) torna-se

(12) $$\frac{d\phi }{dt}=-\mathrm{\Delta }\omega +K\mathrm{sen}\phi ,$$

chamada equação de Adler, e introduzida no contexto de fenômenos de sincronização em circuitos eletrônicos [²⁹, ³⁰]. Portanto, as equações médias para os osciladores de fase quando ${\omega }_1\approx {\omega }_2$ podem ser reescritas como

(13) $$\frac{d{\theta }_1}{dt}={\omega }_1+\frac{K}{2}\mathrm{sen}({\theta }_1-{\theta }_2),$$

(14) $$\frac{d{\theta }_2}{dt}={\omega }_2+\frac{K}{2}\mathrm{sen}({\theta }_2-{\theta }_1).$$

As características gerais da dinâmica dos osciladores de fase acoplados podem ser extraídas da análise qualitativa da equação de Adler (12) [³⁰]. Supondo, sem perda de generalidade, que $\mathrm{\Delta }\omega >0$, vamos considerar inicialmente o caso $K<0$. Se $\mathrm{\Delta }\omega \le |K|$ haverá um ponto fixo de (12) dado pela relação

(15) $$\mathrm{sen}{\phi }^{*}=\frac{\mathrm{\Delta }\omega }{\mathrm{K}}.$$

que corresponde a uma diferença de fase constante

(16) $${\theta }_1-{\theta }_2={\phi }^{*}$$

caracterizando o chamado travamento de fase (phase locking). Na condição de ressonância exata $\mathrm{\Delta }\omega =0$, como ${\phi }^{*}=0$ temos a igualdade entre as fases (${\theta }_1={\theta }_2$), e os osciladores sincronizam em fase, com frequências iguais a

(17) $${\mathrm{\Omega }}_1={\mathrm{\Omega }}_2=\frac{{\omega }_1+{\omega }_2}{2},$$

o caso $K>0$ é similar mas traz uma diferença importante. Conquanto também tenhamos um travamento de fase como em (16), se a condição de ressonância for satisfeita ($\mathrm{\Delta }\omega =0$) teremos ${\phi }^{*}=\pi$, ou seja, ${\theta }_1=\pi +{\theta }_2$: os osciladores estão sincronizados em anti-fase.

Já, se a diferença entre as frequências for suficientemente grande ($\mathrm{\Delta }\omega >|K|$), não haverá um ponto fixo na equação de Adler (12) e, consequentemente, não teremos mais travamento de fases nem sincronização entre suas frequências. Nesse caso, a diferença de fase $\phi$ entre os osciladores varia com o tempo, e dizemos que os osciladores estão à deriva (“drifting”) [¹⁷, ³¹].

4. Acoplamento de N osciladores de fase

A discussão da seção precedente pode ser generalizada para o caso de N osciladores de fase, cada qual com fase ${\theta }_i$ e frequência natural ${\omega }_i$, com $i=\mathrm{1,2},\mathrm{\dots }N$. Quando consideramos o acoplamento entre os osciladores, podemos generalizar as expressões (13)–(14) e escrever

(18) $$\frac{d{\theta }_i}{dt}={\omega }_i+\frac{K}{N}\sum _{j=1}^N\mathrm{sen}({\theta }_{\mathrm{j}}-{\theta }_{\mathrm{i}}),$$

onde K é a intensidade do acoplamento e $i=\mathrm{1,2},\mathrm{\dots }N$.

As frequências naturais dos osciladores dependem das suas características individuais. Como não existem, na prática, osciladores exatamente idênticos, é de se esperar que suas frequências não sejam precisamente iguais, mas tenham uma certa distribuição de probabilidade em torno de um valor médio $\overline{\omega }$.

Vamos denominar essa distribuição de probabilidade $G(\omega )$, tal que $G(\omega )d\omega$ seja o número de osciladores com frequências naturais entre $\omega$ e $\omega +d\omega$. Se integrarmos entre todas as frequências teremos o número total de osciladores

(19) $${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\mathit{d}\omega G(\omega )=N.$$

Se definirmos a densidade de probabilidade

(20) $$g(\omega )=\frac{1}{N}G(\omega ),$$

a eq. (19) leva à condição de normalização

(21) $${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\mathit{d}\omega g(\omega )=1,$$

de forma que o valor médio da frequência é

(22) $$\overline{\omega }={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\mathit{d}\omega \omega g(\omega ).$$

É conveniente redefinir as frequências como

(23) $$\tilde{\omega }=\omega -\overline{\omega },$$

tal que o valor médio é transladado para a origem. Frequências positivas e negativas são interpretadas como resultado de rotações nos sentido anti-horário e horário, respectivamente. Essa redefinição implica em fazermos a seguinte transformação de fase

(24) $${\theta }_i\to {\tilde{\theta }}_i={\theta }_i-\overline{\omega }t,$$

que equivale a mudar a descrição das fases para um referencial que gira com velocidade angular $-\overline{\omega }$.

Vamos supor, também, que a densidade de probabilidade tenha duas características: (i) seja simétrica

(25) $$g(\tilde{\omega })=g(-\tilde{\omega });$$

e (ii) unimodal, ou seja, $g(\tilde{\omega })$ tem um máximo em $\tilde{\omega }=0$. Um exemplo de ambas é a distribuição Lorentziana [Fig.3]:

Figura 3
Histograma para as frequências naturais de N=1000 osciladores de fase. A curva vermelha representa uma distribuição Lorentziana (26) com $\gamma =\mathrm{0,5}$.

(26) $$g(\tilde{\omega })=\frac{\gamma }{\pi ({\gamma }^2+{\tilde{\omega }}^2)},$$

onde $\gamma$ é a largura característica da distribuição, ou seja,

(27) $$g(\tilde{\omega }=\pm \gamma )=\frac{g_{max}}{2}=\frac{1}{2\pi \gamma }.$$

Nestas novas variáveis (fases e frequências) o modelo de Kuramoto (18) é escrito como

(28) $$\frac{d{\tilde{\theta }}_i}{dt}={\tilde{\omega }}_i+\frac{K}{N}\sum _{j=1}^N\mathrm{sen}({\tilde{\theta }}_{\mathrm{j}}-{\tilde{\theta }}_{\mathrm{i}}),(\mathrm{i}=\mathrm{1,2},\mathrm{\dots }\mathrm{N})$$

Doravante, por simplicidade, iremos substituir a notação $\frac{d{\tilde{\theta }}_i}{dt}$ por $\frac{d{\theta }_i}{dt}$ sobre as fases, e ${\tilde{\omega }}_i$ por ${\omega }_i$ sobre as frequências.

5. Sincronização no modelo de Kuramoto

O modelo de Kuramoto

(29) $$\frac{d{\theta }_i}{dt}={\omega }_i+\frac{K}{N}\sum _{j=1}^N\mathrm{sen}({\theta }_{\mathrm{j}}-{\theta }_{\mathrm{i}}),(\mathrm{i}=\mathrm{1,2},\mathrm{\dots }\mathrm{N})$$

(18) é um sistema de N equações diferenciais ordinárias e acopladas de primeira ordem em relação ao tempo [³², ³³] e que, neste caso, são escritas considerando um referencial girante. Portanto, para resolver esse sistema precisamos empregar métodos computacionais, no nosso caso, o método de Runge-Kutta de quarta ordem, com passo fixo $\mathrm{0,01}$.

Como discutimos anteriomente no caso de dois osciladores, para um conjunto de N osciladores é plausível que as frequências ${\omega }_i$ sejam suficientemente próximas e escolhidas aleatoriamente num dado intervalo, segundo uma distribuição de probabilidade. É comum considerarmos frequências que satisfaçam uma distribuição de probabilidade Lorenztiana do tipo (26) com média zero e largura $\gamma =\mathrm{0,5}$ [Fig.3]. Usualmente escolhemos as condições iniciais ${\theta }_i(t=0)$ aleatoriamente com distribuição uniforme. Assim esperamos que os resultados obtidos não dependam, em média, da forma específica que supomos para as condições iniciais.

Resolvendo numericamente as equações (29) obtemos ${\theta }_i(t)$ para cada oscilador $i=\mathrm{1,2},\mathrm{\dots }N$ em um tempo qualquer $t>0$. Nesse caso, as frequências perturbadas serão dadas por

(30) $${\mathrm{\Omega }}_i={\omega }_i+\frac{K}{N}\sum _{j=1}^N\mathrm{sen}({\theta }_{\mathrm{j}}-{\theta }_{\mathrm{i}}),(\mathrm{i}=\mathrm{1,2},\mathrm{\dots }\mathrm{N})$$

e diferem, em geral, das frequências naturais ${\omega }_i$, para $K\ne 0$.

Vamos considerar, inicialmente, o caso $K=\mathrm{0,5}$. Após um intervalo de tempo longo (t = 100) observamos que as fases dos osciladores distribuem-se de maneira relativamente uniforme ao longo do intervalo $[\mathrm{0,2}\pi )$ conforme visto na parte superior da Fig.4(a). Já as frequências perturbadas (30) estão distribuídas de forma simétrica em torno de ${\mathrm{\Omega }}_i=0$, de forma semelhante àquela dos osciladores não-perturbados, ilustrado abaixo na Fig.4(a).

Figura 4
Fases (cor verde, acima) e frequências perturbadas (cor vermelha, abaixo) e não-perturbadas (cor azul, abaixo) para o modelo de Kuramoto com N=100 osciladores nos casos (a) $K=\mathrm{0,5}$; (b) $\mathrm{1,5}$; (c) $\mathrm{3,0}$. Nos painéis acima indicamos os parâmetros de ordem correspondentes aos três casos.

Como generalizamos a diferença de fase para N osciladores acoplados? Dadas as fases dos osciladores num instante t, ${\theta }_i(t)$, Kuramoto introduziu o parâmetro de ordem complexo [³², ³³]

(31) $$z(t)=r(t)e^{i\psi (t)}=\frac{1}{N}\sum _{j=1}^N{e}^{i{\theta }_j(t)}$$

onde r e $\psi$ são o módulo e o argumento, respectivamente. Podemos interpretar z como um vetor no interior do círculo unitário, que é a resultante dos vetores correspondentes a cada oscilador. Se os osciladores estão distribuidos uniformemente ao longo do círculo unitário esse vetor resultante terá um módulo próximo a zero. De fato, no caso ilustrado pela Fig.4(a), obtemos $r=\mathrm{0,14}$ e $\psi =\mathrm{1,56}$ radianos.

Aumentando o valor de K para $\mathrm{1,5}$ observamos uma concentração das fases dos osciladores na parte superior da Fig.4(b), o que aumenta o módulo do parâmetro de ordem para $r=\mathrm{0,69}$ e $\psi =\mathrm{2,19}$. Há uma sincronização entre as frequências perturbadas de um grande número de osciladores em torno de ${\mathrm{\Omega }}_i=0$ mostrado na parte inferior da Fig.4(b). Para valores ainda maiores da intensidade de acoplamento ($K=\mathrm{3,0}$) aumenta a concentração das fases na parte superior da Fig.4(c), tal que $r=\mathrm{0,92}$ e $\psi =\mathrm{5,48}$. Haverá a sincronização de quase todas as frequências conforme mostrado no quadro inferior da Fig.4(c). No limite, se todas as fases fossem exatamente iguais, teríamos r = 1.

Para o caso em que $K=\mathrm{0,5}$, a frequência dos osciladores acoplados Ω_i (vermelho) segue a mesma curva da frequência natural ω_i (azul). Quando $K=\mathrm{1,5}$, nota-se que a maioria dos osciladores apresenta, praticamente a mesma frequência, de modo que, ${\omega }_i-{\omega }_j\approx 0$, ou seja, quase todos os osciladores apresentam frequências próximas a zero. Quando $K=\mathrm{3,0}$ a diferença na frequência para todos os osciladores é praticamente inexistente, o que corrobora os valores do parâmetro de ordem na parte superior da Figura4.

Para apreciarmos de forma mais completa o efeito do aumento da intensidade do acoplamento, calculamos o módulo do parâmetro de ordem para diferentes valores de K em um sistema de $N=\mathrm{100,200,500,1000}$ e 2000 osciladores de fase. É importante mencionar que os valores de r dependem de forma bastante sensível tanto do número de osciladores N como das condições iniciais (${\theta }_i(0)$). Para cada valor de K, escolhemos 10 conjuntos diferentes de condições iniciais aleatórias (alterando a “semente” usada no gerador de números pseudo-randômicos).

Na Fig.5 mostramos o valor da magnitude do parâmetro de ordem r em função da intensidade do acoplamento K, para conjuntos com diferentes números de osciladores acoplados. Para intensidades de acoplamento relativamente pequenas o parâmetro de ordem assume valores entre zero e $\mathrm{0,2}$, indicando ausência de sincronização de fases (pode haver sincronizações por acaso entre alguns osciladores, mas espera-se esse efeito torne-se irrelevante à medida que N cresce). Quando a constante de acoplamento aproxima-se de K=1, o módulo do parâmetro de ordem começa a crescer e atinge valores da ordem de $\mathrm{0,85}$ para $K>3$. De fato, se $K\gg 1$ temos que $r\to 1$. Conforme o número de osciladores aumenta, a curva do parâmetro de ordem se aproxima da previsão teórica para o limite do contínuo $(N\to \mathrm{\infty })$, representada pela linha tracejada preta na Fig.5. A derivação matemática da curva teórica será apresentada no apêndice A. Os resultados numéricos aproximam-se desta curva teórica à medida que o número de osciladores aumenta, como pode ser visto para os conjuntos com N = 1000 e 2000 osciladores de fase, representados pelas curvas azul e amarela, respectivamente, na Figura5

Figura 5
Módulo do parâmetro de ordem em função da intensidade de acoplamento em um sistema de $N=\mathrm{100,200,500,1000}$ e 2000 osciladores de fase acoplados segundo o modelo de Kuramoto. A linha tracejada preta representa a previsão teórica (ver material suplementar) no limite de N grande.

6. Discussão dos resultados

Os resultados mostrados na seção anterior indicam que para um sistema de N osciladores acoplados de acordo com a Equação (29), o aumento da intensidade de acoplamento tem duas consequências: (a) há uma aproximação das fases dos osciladores, até a situação onde elas são iguais; (b) a progressiva sincronização das frequências dos osciladores no seu valor médio. Há, portanto, dois fenômenos ocorrendo: uma sincronização parcial das fases e uma sincronização quase completa das frequências respectivas.

Para uma melhor compreensão desses fenômenos, vamos multiplicar o parâmetro de ordem (31) por $e^{-i{\theta }_i}$:

(32) $$r(t)e^{i(\psi -{\theta }_i)}=\frac{1}{N}\sum _{j=1}^N{e}^{i({\theta }_j-{\theta }_i)},$$

cuja parte imaginária é

(33) $$r(t)\text{sen}(\psi -{\theta }_i)=\frac{1}{N}\sum _{j=1}^N\text{sen}({\theta }_j-{\theta }_i).$$

Substituindo (33) em (29) e considerando um referencial não girante, podemos reescrever o modelo de Kuramoto na forma

(34) $$\frac{d{\theta }_i}{dt}={\omega }_i+Kr\text{sen}(\psi -{\theta }_i),(i=\mathrm{1,2},\mathrm{\dots }N).$$

A despeito dessa forma, as equações para os osciladores de fase não estão desacopladas. Tudo o que fizemos foi embutir o termo de acoplamento em r e $\psi$ (para implementar computacionalmente o modelo de Kuramoto é recomendável usar (34)). Observe que, na Eq.34, estamos utilizando a mesma notação do referencial girante por conveniência.

Uma simplificação considerável no nosso estudo virá se levarmos em consideração apenas os estados estacionários, ou seja, aqueles para os quais o parâmetro de ordem tem módulo r constante e gira com velocidade angular constante $\mathrm{\Omega }$[¹⁷]. Na prática, os estados estacionários podem ser obtidos para tempos suficientemente grandes.

Nessa situação podemos fazer uma nova mudança para um referencial girante com a mesma velocidade angular $\mathrm{\Omega }$, o que equivale à um valor constante $\psi$, que podemos escolher como $\psi =0$. No referencial girante, as equações (34) tornam-se

(35) $$\frac{d{\theta }_i}{dt}={\omega }_i-Kr\text{sen}{\theta }_i,(i=\mathrm{1,2},\mathrm{\dots }N).$$

Impondo, agora, a condição de fases constantes (no referencial girante) $\dot{{\theta }_i^{*}}=0$, obtemos a relação

(36) $$\text{sen}{\theta }_i^{*}=\frac{{\omega }_i}{Kr}.$$

Como $|\text{sen}{\theta }_i^{*}|\le 1$, concluímos que será possível satisfazer a condição de fases constantes desde que

(37) $$-Kr\le {\omega }_i\le Kr,$$

para valores dados de K e r.

Os osciladores que satisfazem (37) têm fases que giram com velocidade angular $\mathrm{\Omega }$ no referencial fixo. Como $\mathrm{\Omega }$ é constante, haverá travamento de fases entre esses osciladores, de modo que as diferenças de fase mútuas ${\phi }_{ij}={\theta }_i^{*}-{\theta }_j^{*}$ são independentes do tempo, como no caso de dois osciladores visto anteriormente. Para saber quanto é essa diferença de fase nós usamos (36) e o valor da frequência natural dos osciladores:

(38) $${\phi }_{ij}=\text{arcsen}\left(\frac{{\omega }_i}{Kr}\right)-\text{arcsen}\left(\frac{{\omega }_j}{Kr}\right).$$

As frequências perturbadas dos osciladores que exibem travamento de fase são dadas, no referencial girante, por

(39) $${\mathrm{\Omega }}_i={\omega }_i-Kr\text{sen}{\theta }_i^{*}=0,(i=\mathrm{1,2},\mathrm{\dots }N)$$

ou seja, esses osciladores sincronizam à frequência zero. É o que vemos nos painéis exibidos na Fig.4(b): uma parte dos osciladores que se concentram em torno de $\theta =\pi$ têm frequência nula. Quanto maior o valor de K maior será o valor de r e, portanto, mais osciladores passam a satisfazer a condição (37). Já os osciladores que não satisfazem (37), por terem $|{\omega }_i|>Kr$, não exibem travamento de fase nem sincronização entre suas frequências, o que explica os osciladores “rebeldes” que vemos nas Figs.4(b) e (c). Em analogia com o caso de dois osciladores diremos que estes osciladores estão ‘à deriva’; pois suas diferenças de fase variam com o tempo. Evidentemente, para K e r suficientemente grandes, a condição (37) será satisfeita para todos os osciladores.

A Fig.5, contendo a variação do módulo do parâmetro de ordem r com a intensidade do acoplamento K, apresenta uma transição entre um estado onde não há travamento de fase nem sincronização de frequência (para K baixo) e um estado com travamento de fase e sincronização de frequência parciais, à medida em que K aumenta. As flutuações observadas na Fig.5 são devidas tanto ao efeito das condições iniciais como ao tamanho finito da rede N. Se a figura for obtida para um número maior de osciladores, essas flutuações serão significativamente menores.

Desse modo, podemos perguntar o que ocorre no limite “termodinâmico”, ou seja, quando $N\to \mathrm{\infty }$. Nesse limite é possível usar a teoria do campo médio para determinar as propriedades gerais do modelo de Kuramoto em termos da intensidade de acoplamento K. Essa análise emprega um ferramental emprestado da mecânica estatística, e os resultados principais são deduzidos no Apêndice A. Vamos resumir os resultados principais da análise do caso $N\to \mathrm{\infty }$.

Quando $0<K<K_c$, onde K_c é um valor crítico dado por

(40) $$K_c=\frac{2}{\pi g(0)},$$

os osciladores estão completamente dessincronizados, ou seja, temos r = 0[¹⁷]. Se $K>K_c$ os osciladores começam a apresentar travamento parcial de fases e sincronização de frequências, de forma que r cresça com a raiz quadrada de K – K_c.

No caso de uma distribuição Lorentziana de frequências com largura $\gamma$, dada por (26), onde $g(0)=1/\pi \gamma$, de modo que $K_c=2\gamma$. Se $\gamma =\mathrm{0,5}$ (como na Fig.5), então K_{c = 1}. Além disso, para $K>K_C$ temos a seguinte relação (exata) entre o parâmetro de ordem e a intensidade do acoplamento [¹⁷]

(41) $$r=\sqrt{1-\frac{2\gamma }{K}}=\sqrt{1-\frac{K_c}{K}}.$$

indicada na Fig.5 por uma curva tracejada. Se $K\to \mathrm{\infty }$, então, $r\to 1$, que corresponde ao perfeito travamento de fases e sincronização de frequência de todos os osciladores. Para outras distribuições de probabilidade $g(\omega )$, é possível obter expressões gerais nas vizinhanças do ponto crítico. O modelo de Kuramoto exibe, assim, um tipo de transição de fase de segunda ordem (contínua) no ponto K_c, semelhante à transição paramagnética-ferromagnética que ocorre em certos materiais.

7. Conclusões

A sincronização de um conjunto de osciladores acoplados é um fenômeno bastante comum em aplicações físicas e biológicas. Como há vários tipos de osciladores, assim como diferentes formas de acoplamento, é interessante abordar esse assunto sob um ponto de vista geral, de modo a salientar os seus aspectos essenciais, sem a necessidade de levar em conta detalhes dos modelos. Esse foi o ponto de vista prevalente no presente trabalho.

A aproximação de osciladores de fase é importante pois reduz a dinâmica a uma única dimensão, simplificando a descrição matemática. Já o tipo de acoplamento pode ser, também, simplificado usando equações médias, donde é possível motivar o modelo de Kuramoto, no qual cada oscilador é influenciado pela média de todos os demais. Esse tipo de sistema tem sido usado como paradigma para o estudo de osciladores acoplados.

Na ausência de acoplamento supomos que as frequências correspondentes são ligeiramente diferentes, satisfazendo uma distribuição de probabilidade simétrica e unimodal. Se a intensidade de acoplamento for suficientemente pequena, os osciladores apresentam um estado incoerente. Em contrapartida, se a intensidade for superior a um dado valor crítico, aparecem duas populações de osciladores: (i) aqueles que apresentam travamento de fase (ou seja, suas diferenças de fase permanecem constantes com o tempo) e sincronização de frequências; (ii) aqueles cujas fases estão à deriva, ou seja, variam com o tempo de forma independente.

Observou-se que o crescimento da fração de osciladores travados/sincronizados cresce monotonicamente com a intensidade do acoplamento. Essa constatação também surge a partir da definição de um parâmetro de ordem, que cresce com a intensidade do acoplamento a partir do seu valor crítico. No limite de infinitos osciladores acoplados, uma teoria de campo médio prevê que essa dependência é da forma ${(K-K_c)}^{1/2}$, onde K_c é um valor crítico que depende do valor máximo da distribuição de frequências dos osciladores. Esse resultado faz com que a transição para a sincronização no modelo de Kuramoto seja comparada à transição de fase de segunda ordem entre sistemas paramagnéticos e ferromagnéticos.

As características observadas no modelo de Kuramoto, tanto o travamento de fase como a sincronização de frequências, podem ser empregadas para explicar qualitativamente sistemas tão diversos como vagalumes piscando, neurônios em rajadas de disparos, e pedestres em uma ponte pênsil. Para cada um destes exemplos, os fenômenos de sincronização resultam em consequências observáveis.

Agradecimentos

Os autores agradecem ao CNPq (processos: 301019/ 2019-3; 161949/2020-7; 302665/2017-0 e 407299/ 2018-1), à CAPES e à Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) – Auxílio No. 2018/ 03211-6 pelo apoio financeiro para a realização deste trabalho.

Apêndice A:

O limite do contínuo

Referências

• ^1.
  S.O. Kepler e M.F.O. Saraiva, Astronomia e Astrofísica (Livraria da Física, São Paulo, 2004), 2^a ed.
• ^2.
  M. Bennett, M.F. Schatz, H. Rockwood e K. Wiesenfeld, Proc. R. Soc. Lond. A 458, 563 (2002).
• ^3.
  M. Denny, European Journal of Physics 23, 449 (2002).
• ^4.
  H.M. Nussenzveig, Curso de Física Básica (Edgard Blücher, São Paulo, 2014), v. 2.
• ^5.
  A. Pikowsky, M. Rosenblum e J. Kurths, Synchronization: a Universal Concept in Nonlinear Sciences (Cambridge University Press, Cambridge, 2001).
• ^6.
  J. Pantaleone, American Journal of Physics 70, 992 (2002).
• ^7.
  M. Kapitaniak, K. Czolczynski, P. Perlikowski, A. Stefanski e T. Kapitaniak, Physics Reports 517, 1 (2012).
• ^8.
  J.L. Cloudsley-Thompson, em: The Study of Time III, editado por J.T. Fraser, N. Lawrence e D.A. Park, (Springer, New York, 1978).
• ^9.
  J. Yan e A. Goldbeter, J. R. Soc. Interface. 16, 20190376 (2019).
• ^10.
  R.C. Elson, A.I. Selverston, R. Huerta, N.F. Rulkov, M.I. Rabinovich e H.D.I. Abarbanel, Phys. Rev. Lett. 81, 5692 (1998).
• ^11.
  S. Strogatz, Sync: How Order Emerges from Chaos in the Universe, Nature, and Daily Life (Hachette Book, Paris, 2012).
• ^12.
  Z. Néda, E. Ravasz, Y. Brechet, T. Viczek e A.L. Barabási, Nature 403, 849 (2000).
• ^13.
  D. Eroglu, J.S.W. Lamb e T. Pereira, Contemporary Physics 58, 207 (2017).
• ^14.
  S.H. Strogatz, D.M. Abrams, A. McRobie, B. Eckhardt e E. Ott, Nature 438, 43 (2005).
• ^15.
  L.H.A. Monteiro, Sistemas Dinâmicos (Editora da Física, São Paulo, 2007), 2^a ed.
• ^16.
  R.L. Viana, Rev. Bras. Ens. Fís. 42, e20200272 (2020).
• ^17.
  J.A. Acebrón, L.L. Bonilla, C.J. Pérez Vicente, F. Ritort e R. Spigler, Reviews of Modern Physics 77, 137 (2005).
• ^18.
  F.A. Rodrigues, T.K. Peron, P. Jie e J. Kurths, Physics Reports 610, 1 (2016).
• ^19.
  F.A.S. Ferrari, R.L. Viana, S.R. Lopes e R. Stoop, Neural Networks 66, 107 (2005).
• ^20.
  C. Koch e I. Segev, Methods in neuronal modeling: from ions to networks (The MIT Press, Cambridge, 1999), 2^a ed.
• ^21.
  A.L. Hodgkin e A.F. Huxley, J Physiol 117, 500 (1952).
• ^22.
  E.M. Izhikevich, IEEE Transactions on Neural Networks 14, 1569 (2003).
• ^23.
  N.F. Rulkov, Phys Rev E 65, 041922 (2002).
• ^24.
  T.L. Prado, S.R. Lopes, C.A.S. Batista, J. Kurths e R.L. Viana, Phys Rev E 90, 032818 (2014).
• ^25.
  F.A.S. Ferrari, R.L. Viana, F. Gomez, T. Lorimer e R. Stoop, New J. Phys. 17, 055024 (2015).
• ^26.
  J.A.P. Silveira, P.R. Protachevicz, R.L. Viana e A.M. Batista, Neurocomputing 436, 126 (2021).
• ^27.
  R.P. Aristides e R.L. Viana, Int. J. Nonlinear Dynamics and Control 1, 376 (2019).
• ^28.
  M.F. Bear, B.W. Connors e M.A. Paradiso, Neurociências: Desvendando o sistema nervoso (Artmed, Porto Alegre, 2002), 2^a ed.
• ^29.
  R. Adler, Proceedings of the IRE34, 351 (1946).
• ^30.
  R.L. Viana e F.F. de Carvalho, Rev. Bras. Ens. Fís. 39, e3306 (2017).
• ^31.
  S.H. Strogatz, Physica D 143, 1 (2000).
• ^32.
  Y. Kuramoto, em Lecture Notes in Physics, International Symposium on Mathematical Problems in Theoretical Physics, editado por H. Araki (Springer-Verlag, New York, 1975), v. 39, p. 420.
• ^33.
  Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence (Springer-Verlag, New York, 1984).

Datas de Publicação

• Publicação nesta coleção
  11 Fev 2022
• Data do Fascículo
  2022

Histórico

• Recebido
  15 Out 2021
• Revisado
  05 Jan 2022
• Aceito
  07 Jan 2022