Modelo de agentes activos de Vicsek
23 Feb 2024Hace algunos meses me enteré de la existencia del Núcleo Milenio de Materia Activa, el cual se especializa en el estudio de la dinámica de un gran número de miembros o "agentes", los cuales son unidades independientes que interaccionan con otros agentes vecinos a través de reglas simples y empíricas, por lo menos en sus fundamentos. De esta forma es posible describir el movimiento espontáneo y coherente de ciertos sistemas biológicos, como las bandadas de los pájaros. En el año 1995, el Biólogo físico Tamás Vicsek publicó un estudio llamado NOVEL TYPE OF PHASE TRANSITION IN A SYSTEM OF SELF-DRIVEN PARTICLES donde describe un modelo para describir el movimiento coherente de los cardúmenes y de las bandadas de pájaros en el vuelo.
wikipedia lo describe de la siguiente forma:
La materia activa está compuesta de un gran número de «agentes activos», cada uno de los cuales consume energía para moverse o para ejercer fuerzas mecánicas. Debido al consumo de energía, estos sistemas están intrínsecamente fuera de equilibrio térmico. Ejemplos de materia activa son los bancos de peces, bandadas de aves, bacterias, partículas autopropulsadas artificiales, y la auto-organización de los bio-polímeros tales como los microtúbulos y la actina, siendo parte ambos del citoesqueleto de las células vivas. La mayoría de los ejemplos de materia activa son de origen biológico; sin embargo, una gran cantidad del trabajo experimental está dedicado a los sistemas sintéticos.
Sin duda es una definición muy interesante, tomada o inspirada desde la mecánica estadística y aplicada elegántemente en los sistemas complejos biológicos fuera de equilibrio.
El modelo de Vicsek pretende cuantificar matemáticamente la interacción de un miembro del grupo, considerando su vecindad dentro de un rango o distancia de interacción. En este caso, para cada pájaro dentro de la bandada tanto como su posición, velocidad y ángulo del camino serán determinadas por la cantidad de vecinos dentro de una distancia particular, y con eso se creará la "materia activa", generando dichas estructuras coherentes o "materia blanda" que se observan en muchas configuraciones a nuestro alrededor.
Dado que el sistema describe la evolución temporal de cada pájaro miembro de la bandada, dada su posición inicial $\vec{r}(t=0)$ o en cualquier momento $\vec{r}(t)$, la posición cambiará en el siguiente instante $t + \Delta t$ siguiendo la ecuación itinerario (considerando la aceleración $\vec{a}=\vec{0}$):
$$ {\displaystyle \vec{r}(t+\Delta t)=\vec{r}(t)+\vec{v}(\theta) \Delta t}, $$donde la velocidad $\vec{v}(\theta) = v_0 \left( \cos(\theta)\ \hat{i} + \sin(\theta)\ \hat{j}\right)$ es constante para cada individuo, y la dirección dependerá completamente del ángulo $\theta(t)$, la que a su vez es definida por los agentes en la vecindad local, la que en instante siguiente $t + \Delta t$ es definido de la forma:
$$ {\displaystyle \theta(t+\Delta t)=\langle \theta \rangle _{|r_{i}-r_{j}|<R}+\eta\ (t)}, $$donde $\langle \theta \rangle$ representa el valor promedio de todas las direcciones del movimiento de cada vecino dentro del radio de acción ${|r_{i}-r_{j}|<R}$. En la siguiente imagen se grafica la configuración y la cuantificación matemática del modelo de materia activa. El agente azul se está moviendo hacia la izquierda, y su dirección se modificará teniendo en cuenta los agentes dentro de la circunferencia al rededor del agente azul.
El ángulo de cada ave dentro del círculo es considerado y el valor promedio del ángulo es obtenido de la forma:
$$ \displaystyle \langle \theta \rangle = \arctan{ \left( \frac{\sum_{k=1}^{N} \sin{(\theta_k)} }{\sum_{k=1}^{N} \cos{(\theta_k)}} \right) } $$Además, el término $\eta\ (t)$ agrega una perturbación aleatoria en la dirección de movimiento, para así incluir algo de realismo al modelo. Todo lo anterior debe ser calculado para cada uno de los pájaros dentro de la bandada, por lo que la complejidad de este modelo recae en el poder de cómputo disponible para calcular cada uno de los anteriores parámetros a través de todos los momentos (o épocas) necesarios para ver la evolución temporal de la bandada.
Todo el procedimiento anterior es explicado a continuación en forma de pseudo-algoritmo de la ejecución del modelo de Vicsek:
Para cada uno de los pájaros dentro de la bandada:
- Definir de forma aleatoria las posiciones $\vec{r}_i=(x,y)$ y el angulo inicial $\theta_i$.
- Contamos cuántos vecinos están alrededor, dentro del radio $R$, y sus ángulos son considerados.
- Calculamos el valor promedio de los ángulos $\langle \theta \rangle$.
- Actualizamos la dirección de la la velocidad $\vec{v}$ utilizando el nuevo ángulo $\langle \theta \rangle$.
- Agregamos el término de perturbación random $\eta$.
- Finalmente, actualizamos la posición de cada pájaro usando la ecuación itinerario en el instante $t + \Delta t$.
- Repetir por cada época deseada.
Todo lo anterior es implementado en el siguiente script de python, usando como base la implementación de Philip Mocz (ver referencias). Para la siguiente simulación, se considerará una bandada de 3000 pájaros, en 200 épocas que se diferencian entre ellas por un tiempo $\Delta t=0.01$ segundos. No se diga más.
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
import sys
plt.style.use('dark_background')
#Parámetros básico del modelo de Vicsek
#velocidad inicial:
v0 = 2.0
#Número de pájaros:
N = 3000
#Parámetro de ruido angular:
eta = 0.9
#Tamaño de la caja:
Lx = 4
Ly = 5
#Radio de interacción:
R = 0.1
#paso temporal Delta t:
dt = 0.01
Dt=dt
#Número de épocas:
Nt = 200
np.random.seed(N)
#Se eligen las posiciones y velocidades de forma random:
x = np.random.rand(N)*Lx
y = np.random.rand(N)*Ly
v_init = np.random.rand(N)*v0
theta = 2*np.pi*np.random.rand(N)
#Calculamos las componentes individuales
vx = v_init * np.cos(theta)
vy = v_init * np.sin(theta)
#Iniciamos el plot y método quiver para visualizar:
fig, ax = plt.subplots(figsize=(4, 5), dpi=100)
quiver = ax.quiver(x, y, vx, vy, color='yellow',linewidth=3.9, headwidth=10, headlength=10)
#Función para crear el video:
def dibujar_plot(frame, quiver):
#Iniciamos los valores para modificiarlos:
global x, y, Lx, Ly, theta, vx, vy, Dt
#Limpiamos el gráfico:
ax.clear()
ax.set_title('Simulación del modelo Vicsek')
ax.set_xlim(0, Lx)
ax.set_ylim(0, Ly)
#Actualizamos la posicion de cada integrante
x += vx * dt
y += vy * dt
#Condiciones de borde periódicas:
#(esto me voló la cabeza)
#https://pythoninchemistry.org/sim_and_scat/important_considerations/pbc.html
x = x % Lx
y = y % Ly
#Contamos cuántos agentes están dentro del radio R:
mean_theta = theta[:]
for b in range(N):
neighbors = np.sqrt((x - x[b])**2 + (y - y[b])**2) < R
sx = np.mean(np.cos(theta[neighbors]))
sy = np.mean(np.sin(theta[neighbors]))
mean_theta[b] = np.arctan2(sy, sx)
#Añadimos el ruido, hay que restar 0.5 porque los valores random van entre 0 y 1:
theta = mean_theta + eta*(np.random.rand(N)-0.5)
#Actualizamos la velocidad:
vx = v_init * np.cos(theta)
vy = v_init * np.sin(theta)
#Llevamos la cuenta del tiempo:
Dt += dt
#Actualizamos el plot con la nueva época:
quiver = ax.quiver(x, y, vx, vy, color='yellow',linewidth=3.9, headwidth=10, headlength=10)
ax.set_aspect('equal')
ax.tick_params(axis="x", labelsize=0)
ax.get_yaxis().set_visible(False)
ax.set_xlabel('tiempo t=%.1f [s]' %Dt)
#Guardamos la simulación:
ani = animation.FuncAnimation(fig, dibujar_plot, fargs=(quiver,), frames=Nt, interval=100)
ani.save('simulation_vicsek.mp4', writer='ffmpeg', fps=100)
plt.close()
Un video extendido de esta simulación puede ser encontrado en mi en mi canal de youtube o en mi cuenta de instagram @nicomediap
¡Gracias por la atención!
Referencias;
Simulate Active Matter: Viscek model for flocking birds: Una gran fuente de inspiración, implementado por Philip Mocz (2021). Básicamente tomé este código y le agregé cosas de mi cosecha.
Active Matter: Modelling swarming behaviour using the Vicsek model: Una exposición en Youtube sobre el tema.