Aparición y desaparición de equilibrios en la bifurcación de nodo-silla

Tres fenómenos que hay que separar

La manera más simple en que una familia de sistemas cambia de organización al variar un parámetro es que dos equilibrios nazcan juntos al atravesar un valor crítico. El capítulo 43 la mostró con el ejemplo $\dot x=\mu-x^2$, al definir la bifurcación como un cambio cualitativo en la organización de las soluciones. Este capítulo estudia ese mecanismo, llamado bifurcación de nodo-silla o de pliegue, junto con las condiciones que permiten reconocerlo en familias generales (Kuznetsov, 2023).

Su estudio obliga a separar tres fenómenos que en una serie temporal pueden parecer similares. Un equilibrio puede dejar de existir. Puede también seguir existiendo y cambiar de estabilidad. Cabe, además, que el sistema abandone la cuenca de un equilibrio que sigue existiendo, porque una perturbación lo llevó al otro lado de su frontera. Solo el primero corresponde a la nodo-silla. Reconstruiremos la bifurcación desde la ecuación, el diagrama de ramas y la línea de fase, formularemos sus condiciones genéricas y las aplicaremos a un modelo de población con cosecha. Después extenderemos el resultado al plano y compararemos la desaparición de un estado estable con el cambio de destino que produce una perturbación con el parámetro fijo.

La forma normal y sus líneas de fase

Consideremos otra vez

con $x\in\mathbb R$ y un parámetro real $\mu$. Los equilibrios $x^*_\pm=\pm\sqrt\mu$ existen para $\mu\ge0$. Con $\mu>0$ la rama positiva es estable, porque $f_x=-2\sqrt\mu<0$, mientras que la negativa es inestable, porque allí $f_x=2\sqrt\mu>0$. En $\mu=0$ hay un único equilibrio doble y para $\mu<0$ no hay ninguno.

La función $f$ muestra la reorganización sin integrar trayectorias (figura 44.1a). Para $\mu>0$ la parábola corta el eje en dos puntos. Entre ellos $f>0$ y $x$ aumenta, mientras que fuera de ellos $f<0$ y $x$ disminuye. El equilibrio negativo repele por ambos lados y el positivo atrae por ambos lados dentro de su entorno. En $\mu=0$ la parábola toca el eje sin cruzarlo, de modo que el signo del campo no cambia. Para $\mu<0$ la curva queda entera por debajo del eje y la velocidad no se anula en ningún punto. La figura 44.1b superpone al diagrama de ramas las líneas de fase de varios sistemas de la familia.

La palabra colisión requiere una aclaración. Los dos equilibrios no se desplazan uno hacia otro en el tiempo. Sus posiciones dependen del parámetro. La distancia entre ellos, $2\sqrt\mu$, tiende a cero cuando $\mu\to0^+$, por lo que la colisión ocurre en el espacio de estado y parámetro. Por la misma razón, la pendiente de la rama,

diverge en el pliegue sin que ninguna variable adquiera velocidad infinita. Es una derivada respecto del parámetro y nada dice sobre la rapidez del movimiento en el tiempo.

El punto crítico

En $\mu=0$ la ecuación es $\dot x=-x^2$. El origen es un equilibrio con $f_x(0;0)=0$, de modo que la linealización, $\dot\xi=0$, no decide su estabilidad. El término cuadrático contiene toda la información. Como $-x^2<0$ para $x\neq0$, el campo apunta hacia la izquierda a ambos lados del origen. Una trayectoria que parte de $x_0>0$ se acerca a cero, mientras que una que parte de $x_0<0$ se aleja hacia valores más negativos. El equilibrio es semiestable.

Separando variables se obtiene la solución exacta

Con $x_0>0$ se tiene $x(t)\sim1/t$, una aproximación algebraica más lenta que cualquier exponencial. Con $x_0<0$ el denominador se anula en $t=-1/x_0$ y la solución diverge hacia $-\infty$ en tiempo finito.

Esa divergencia no es una predicción general de las bifurcaciones de nodo-silla. La forma normal describe la dinámica cerca del punto crítico. Cuando la trayectoria sale de ese entorno, los términos que la forma normal omite pueden acotarla o llevarla a otro atractor, como ocurrirá en el ejemplo de la cosecha y en el modelo con dos atractores de este capítulo.

Ralentización antes del pliegue

En la rama estable, $x^*=\sqrt\mu$, la linealización es $\dot\xi=-2\sqrt\mu\,\xi$. Una perturbación pequeña decae como $e^{-2\sqrt\mu\,t}$, con tiempo característico

Cuando $\mu\to0^+$ el valor propio $-2\sqrt\mu$ se acerca a cero y $\tau$ crece sin límite. Las perturbaciones tardan cada vez más en disiparse, un fenómeno llamado ralentización crítica.

La forma normal permite comprobarlo de manera exacta. Con $\mu=a^2>0$ y $x(0)=0$, la solución es

como se verifica al sustituir, porque $\dot x=a^2\operatorname{sech}^2(at)=a^2-x^2$. Para tiempos grandes, $a-x(t)\approx2a\,e^{-2at}$, que decae con la tasa $2a=2\sqrt\mu$ predicha por la linealización.

La expresión de $\tau$ vale para la linealización de esta forma normal y para perturbaciones pequeñas alrededor de la rama estable. De ella no se sigue que un sistema empírico que se recupera lentamente esté cerca de una nodo-silla. Volveremos sobre este punto al final del capítulo.

El cuello después del pliegue

La desaparición de los equilibrios no borra de inmediato su influencia. Tomemos $\mu=-\delta^2<0$. Entonces $\dot x=-(\delta^2+x^2)$. No hay equilibrio, aunque cerca de $x=0$ la velocidad es pequeña si $\delta$ lo es. Una trayectoria que atraviesa esa región avanza despacio. La región recibe el nombre de cuello o de fantasma de la nodo-silla (Strogatz, 2015, sec. 4.3).

El tiempo necesario para ir desde $x=L$ hasta $x=-L$, con $L>0$, se calcula sin simular,

Cuando $\delta\to0^+$ el arcotangente tiende a $\pi/2$ y

La escala temporal también diverge al acercarse al valor crítico desde el lado sin equilibrios, con el mismo exponente $-1/2$ que el tiempo de recuperación del otro lado (figura 44.2a). Casi todo el tiempo de paso transcurre cerca de $x=0$. Con $\mu=-0{,}01$, la trayectoria que parte de $x=1$ tarda unas 29 unidades de tiempo en llegar a $x=-1$ y permanece la mayor parte de ese intervalo muy cerca del origen (figura 44.2b).

Esta demora puede producir una serie temporal engañosa. Después de cruzar el valor crítico, el sistema puede permanecer mucho tiempo cerca de donde estaba el equilibrio y solo después alejarse. En un registro corto parecería que el estado estacionario persiste. El retraso proviene de la dinámica determinista, sin intervención del ruido ni de errores numéricos. Tampoco es el retraso dinámico del capítulo 43, que aparece cuando el parámetro cambia durante la evolución y crece como la potencia $2/3$ de la velocidad del barrido. Aquí el parámetro está fijo y lo que crece es el tiempo de tránsito por la región donde el campo casi se anula.

Condiciones genéricas de una nodo-silla

La forma normal resume el comportamiento de una clase amplia de sistemas cerca de su punto crítico. Consideremos $\dot x=F(x,\mu)$ con $F$ suficientemente suave y un punto $(x_c,\mu_c)$ con $F(x_c,\mu_c)=0$. Para que un equilibrio hiperbólico desaparezca localmente debe perder la hiperbolicidad, lo que en una dimensión exige $F_x(x_c,\mu_c)=0$. Esa condición no basta. Una nodo-silla genérica requiere además dos condiciones de no degeneración,

La condición sobre $F_{xx}$ asegura que el término cuadrático en el estado domina la no linealidad. La condición sobre $F_\mu$ exige que la variación del parámetro desplace el campo de manera transversal a través de la condición crítica (Kuznetsov, 2023).

Con $u=x-x_c$ y $\nu=\mu-\mu_c$, el desarrollo de Taylor alrededor del punto crítico es

porque los términos constante y lineal en $u$ se anulan. Si $a=F_\mu\neq0$ y $b=\tfrac12F_{xx}\neq0$, los cambios de escala y de orientación reducen la dinámica local a $\dot z=\eta-z^2$ o a una de sus variantes con signos invertidos. Los signos de $a$ y $b$ deciden de qué lado del valor crítico existen las ramas y en qué sentido fluye el campo, sin alterar el mecanismo de creación o destrucción del par.

Dos contraejemplos muestran por qué importan esas condiciones. En $\dot x=\mu-x^3$ se cumplen $F=0$ y $F_x=0$ en el origen, pero también $F_{xx}=0$. El equilibrio $x^*=\sqrt[3]\mu$ existe para todo $\mu$ y ningún par aparece ni desaparece. En $\dot x=\mu^2-x^2$ falla la transversalidad, porque $F_\mu(0,0)=0$. Para todo $\mu\neq0$ hay dos equilibrios, $\pm|\mu|$, que coinciden en $\mu=0$. El retrato de fase en $\mu=0$ difiere del de los valores vecinos, aunque la familia vuelve al mismo régimen a ambos lados y el número de equilibrios no cambia al atravesar el valor crítico. Un valor propio nulo señala que la linealización deja de decidir. El tipo de bifurcación lo fijan los términos y las condiciones adicionales.

La nodo-silla genérica es una bifurcación de codimensión uno. En el punto crítico deben cumplirse $F=0$ y $F_x=0$, dos ecuaciones que, para una familia escalar con un parámetro, pueden resolverse para las dos incógnitas $(x_c,\mu_c)$. Las condiciones $F_{xx}\neq0$ y $F_\mu\neq0$ garantizan que no haga falta ningún ajuste adicional. Por eso la nodo-silla aparece en modelos muy distintos sin necesidad de una simetría especial. Si se perturban un poco los coeficientes de una familia que cumple las condiciones, el punto crítico se desplaza y el mecanismo persiste. Esta robustez es una propiedad de la familia. En el valor crítico mismo el equilibrio no es hiperbólico y el sistema individual no es estructuralmente estable.

Un ejemplo aplicado con cosecha constante

Una población con crecimiento logístico sometida a una extracción constante $h\ge0$ cumple, en variables adimensionales,

Los equilibrios satisfacen $x^2-x+h=0$, de donde $x_\pm=\tfrac12\bigl(1\pm\sqrt{1-4h}\bigr)$, que existen para $h\le1/4$. Las condiciones del punto crítico, $F=x(1-x)-h=0$ y $F_x=1-2x=0$, dan $x_c=1/2$ y $h_c=1/4$. Como $F_{xx}=-2\neq0$ y $F_h=-1\neq0$, el valor $h=1/4$ es una nodo-silla genérica.

En este caso la reducción a la forma normal es exacta. Con $u=x-1/2$ se tiene $x(1-x)=1/4-u^2$, de modo que

con parámetro efectivo $1/4-h$. Por debajo del umbral, la rama estable $x_+$ tiene valor propio $F_x(x_+)=-\sqrt{1-4h}$. Con $h=0{,}24$ el tiempo de recuperación vale 5 unidades y con $h=0{,}2475$ vale 10. Por encima del umbral no hay equilibrio positivo y la población decrece hasta extinguirse. Con $h=0{,}26$ tarda unas 27,5 unidades en caer desde $x=1$ hasta cero. De ese tiempo, $5\pi\approx15{,}7$ unidades transcurren entre $x=0{,}6$ y $x=0{,}4$, cerca del antiguo punto de pliegue.

El ejemplo muestra las dos caras del mecanismo. Antes del umbral, la población se recupera cada vez más despacio de las perturbaciones. Después, el colapso es inevitable dentro del modelo, aunque su comienzo puede parecer lento. El modelo supone una extracción fija e independiente del tamaño de la población. Su utilidad para una pesquería concreta depende de que esa hipótesis sea razonable.

Por qué se llama nodo-silla

La forma escalar puede insertarse en un sistema plano muy simple,

Para $\mu>0$ los equilibrios son $E_\pm=(\pm\sqrt\mu,0)$ y el jacobiano es diagonal, con elementos $-2x$ y $-1$. En $E_+$ los valores propios son $-2\sqrt\mu$ y $-1$, un nodo estable. En $E_-$ son $2\sqrt\mu$ y $-1$, una silla. En $\mu=0$ ambos equilibrios coinciden en el origen, con valores propios $0$ y $-1$. La dirección $y$ sigue contrayéndose mientras la dirección $x$ deja de ser hiperbólica, por lo que la dinámica decisiva se concentra en una dimensión central. El nombre de la bifurcación proviene de este choque entre un nodo y una silla.

En dimensión mayor, una nodo-silla genérica con un único valor propio nulo simple puede reducirse localmente a una ecuación escalar sobre una variedad central (Guckenheimer y Holmes, 1983). Las cantidades de no degeneración se calculan con los vectores nulos derecho e izquierdo del jacobiano crítico $J=D_{\mathbf x}\mathbf F$, que cumplen

Con ellos se definen

Si $a\neq0$ y $b\neq0$, la dinámica reducida sobre la dirección crítica tiene la forma $\dot u=a\nu+bu^2+\cdots$ (Kuznetsov, 2023). La fórmula permite reconocer el mecanismo cuando el equilibrio no puede aislarse sobre una línea de fase, como en el sistema acoplado del ejercicio 44.3. No reemplaza la comprobación de las hipótesis ni informa sobre las cuencas una vez que la trayectoria abandona el entorno local.

Desaparición de un equilibrio y cruce de una cuenca

Para mostrar la diferencia hace falta un sistema con más de un atractor. Consideremos

Cerca de $x=0$ el factor $x+2$ es positivo y próximo a 2, así que la estructura local de la nodo-silla se conserva salvo por una reescala temporal suave. El factor adicional introduce un equilibrio persistente en $x=-2$ (figura 44.3a). Con $\mu=0{,}25$ los equilibrios son $-2$, $-0{,}5$ y $0{,}5$. Las derivadas $G_x$ valen $-3{,}75$, $1{,}5$ y $-2{,}5$, por lo que $-2$ y $0{,}5$ son estables y $-0{,}5$ es inestable. En una dimensión ese equilibrio inestable separa las cuencas, $(-\infty,-0{,}5)$ para $-2$ y $(-0{,}5,\infty)$ para $0{,}5$.

Comparemos dos intervenciones. Con $\mu=0{,}25$ fijo, un sistema que estaba en $x=0{,}5$ y recibe una perturbación que lo lleva a $x=-0{,}4$ permanece en la cuenca derecha y vuelve a $0{,}5$. Si la perturbación lo lleva a $x=-0{,}6$, cruza la frontera y termina en $-2$ (figura 44.3b). La trayectoria cambia de destino sin que la familia atraviese ninguna bifurcación, porque el parámetro no cambió y los tres equilibrios siguen existiendo. Si en cambio el sistema parte de la rama estable derecha y el parámetro se reduce a $\mu=-0{,}25$, el par $\pm\sqrt\mu$ deja de existir y solo queda el equilibrio $-2$. Desde $x=0{,}5$ la trayectoria llega a $-1{,}99$ en unas 3,5 unidades de tiempo. Ese cambio sí resulta de atravesar la nodo-silla.

La distinción importa al interpretar datos. Una transición brusca no identifica por sí sola el mecanismo. Puede corresponder a la desaparición de un atractor, al cruce de una frontera de cuenca o a una perturbación externa, entre otras posibilidades. El diagrama de bifurcación es una propiedad de un modelo parametrizado. Atribuirlo a un fenómeno real exige evidencia independiente de que el modelo y su parámetro describen adecuadamente ese fenómeno.

Qué permite inferir la ralentización

La nodo-silla produce dos escalas temporales que divergen cerca del punto crítico. Desde el lado con equilibrio estable, $\tau\propto\mu^{-1/2}$. Desde el lado sin equilibrios, el tiempo de paso por el cuello crece como $(-\mu)^{-1/2}$. Estas leyes son consecuencias de la forma normal. Por eso en un modelo gobernado por una nodo-silla es razonable esperar tiempos de recuperación crecientes antes del umbral y permanencias prolongadas después.

El razonamiento inverso es más débil. Una recuperación lenta no demuestra que haya una nodo-silla cerca. Otras bifurcaciones también producen ralentización, como la transcrítica, en la que un valor propio cruza el cero sin que desaparezca el equilibrio, o la de Hopf, en la que la parte real de un par complejo se acerca a cero. En datos con ruido, la estimación del tiempo de recuperación depende además de la frecuencia de muestreo y de la longitud de la ventana de estimación.

La idea de señal temprana debe mantenerse, por eso, separada de la demostración de una bifurcación. Scheffer y colaboradores (2009) propusieron el aumento de la autocorrelación y de la varianza como indicadores genéricos de la cercanía de una transición crítica. Van de Leemput y colaboradores (2014) aplicaron esa idea a registros de emociones tomados en la vida diaria y encontraron valores más altos de esos indicadores en personas que después cambiaron de estado anímico. Boettiger y Hastings (2012) mostraron que la fiabilidad de estos indicadores varía mucho entre conjuntos de datos y que sus tasas de falsas alarmas y de detecciones fallidas deben cuantificarse. En un modelo conocido podemos derivar qué magnitud se acerca a cero y con qué ley. En un sistema empírico hay que establecer además que el parámetro pertinente está cambiando, que la dinámica local es aproximadamente estacionaria mientras se estima y que otros mecanismos no explican el patrón observado.

Comprobación numérica

Las integraciones del capítulo usaron el método DOP853 (Hairer et al., 1993) de SciPy (Virtanen et al., 2020), con tolerancia relativa $10^{-12}$ y absoluta $10^{-14}$. Con $\mu=0{,}25$ y $x(0)=0$, la solución numérica difiere de $0{,}5\tanh(0{,}5t)$ en menos de $3{,}1\times10^{-12}$ entre $t=0$ y $t=20$. La solución crítica desde $x_0=0{,}8$ difiere de la exacta en menos de $4\times10^{-13}$. En el modelo de cosecha con $h=0{,}26$, los tiempos de paso numéricos coinciden con la fórmula del cuello en todas las cifras mostradas en el texto.

El bloque siguiente estima la tasa de recuperación a partir de una trayectoria simulada, mediante un ajuste lineal del logaritmo de la distancia al equilibrio. También mide el tiempo de paso por el cuello con un evento que detiene la integración en $x=-1$.

import numpy as np
from scipy.integrate import solve_ivp

opciones = dict(method="DOP853", rtol=1e-12, atol=1e-14)

print("Recuperación en la rama estable, mu > 0")
print("    mu    tasa estimada   2*sqrt(mu)")
for mu in (1.0, 0.25, 0.04, 0.01):
    xe = np.sqrt(mu)
    sol = solve_ivp(lambda t, x: mu - x**2, (0, 10 / xe), [xe * 1.001],
                    dense_output=True, **opciones)
    t = np.linspace(2 / xe, 6 / xe, 200)     # ventana alejada del inicio
    d = sol.sol(t)[0] - xe                   # distancia al equilibrio
    tasa = -np.polyfit(t, np.log(d), 1)[0]   # pendiente de log(d)
    print(f"{mu:6.2f}   {tasa:12.6f}   {2 * xe:10.6f}")

print("Paso por el cuello, mu < 0, desde x = 1 hasta x = -1")
print("     mu     T numérico      fórmula   pi/sqrt(-mu)")
for mu in (-0.1, -0.01, -0.001):
    e = np.sqrt(-mu)
    llegada = lambda t, x: x[0] + 1.0
    llegada.terminal = True
    sol = solve_ivp(lambda t, x: mu - x**2, (0, 1e4), [1.0],
                    events=llegada, **opciones)
    T = sol.t_events[0][0]
    print(f"{mu:7.3f}  {T:12.6f}  {2 / e * np.arctan(1 / e):11.6f}  {np.pi / e:12.6f}")

La ejecución produce la salida siguiente.

Recuperación en la rama estable, mu > 0
    mu    tasa estimada   2*sqrt(mu)
  1.00       2.000002     2.000000
  0.25       1.000001     1.000000
  0.04       0.400000     0.400000
  0.01       0.200000     0.200000
Paso por el cuello, mu < 0, desde x = 1 hasta x = -1
     mu     T numérico      fórmula   pi/sqrt(-mu)
 -0.100      7.997520     7.997520      9.934588
 -0.010     29.422553    29.422553     31.415927
 -0.001     97.346549    97.346549     99.345883

Las tasas estimadas coinciden con $2\sqrt\mu$ hasta el quinto decimal y los tiempos de paso coinciden con la fórmula exacta en todas las cifras impresas. La aproximación $\pi/\sqrt{-\mu}$ mejora a medida que $\mu$ se acerca a cero, porque el cuello ocupa una fracción creciente del recorrido.

Ejercicios

Ejercicio 44.1. Considere la familia $\dot x=\mu-2x+x^2$. (a) Encuentre el punto crítico y compruebe las condiciones de no degeneración. (b) ¿De qué lado del valor crítico existen los equilibrios y cuál de ellos es estable? (c) Obtenga la forma normal con el cambio $u=x-1$.

Ejercicio 44.2. Decida si cada una de las familias siguientes tiene una nodo-silla genérica en $(x,\mu)=(0,0)$ y describa el número de equilibrios cercanos al origen a ambos lados de $\mu=0$. (a) $\dot x=\mu-x^2-x^4$. (b) $\dot x=\mu-x^3$. (c) $\dot x=\mu^2-x^2$.

Ejercicio 44.3. Considere el sistema acoplado $\dot x=y-x^2$, $\dot y=\mu-y$. (a) Encuentre los equilibrios y el valor crítico. (b) Calcule el jacobiano crítico y los vectores nulos $v$ y $w$ con $w^{\mathsf T}v=1$. (c) Calcule $a$ y $b$ y escriba la ecuación reducida. (d) Clasifique los equilibrios para $\mu>0$.

Ejercicio 44.4. En el modelo de cosecha $\dot x=x(1-x)-h$, (a) calcule el tiempo de recuperación de la rama estable para $h=0{,}2$ y para $h=0{,}249$. (b) Con $h=0{,}26$, calcule el tiempo de paso entre $x=0{,}6$ y $x=0{,}4$. (c) Con $h=0{,}2501$, estime el tiempo que la población pasa en el cuello y compárelo con el caso $h=0{,}26$.

Ejercicio 44.5. Considere el tiempo de recuperación $\tau$ de la forma normal. (a) ¿Por qué factor aumenta $\tau$ cuando $\mu$ se divide por 100? (b) Si en un sistema que se supone gobernado por la forma normal se mide primero $\tau=2$ y después $\tau=4$, ¿qué razón entre los valores de $\mu$ se deduce? (c) Dé dos razones por las que un aumento observado de $\tau$ no basta para afirmar que un sistema se acerca a una nodo-silla.

Ejercicio 44.6. En $\dot x=(x+2)(\mu-x^2)$ con $\mu=0{,}25$, (a) verifique la estabilidad de los tres equilibrios. (b) ¿Qué perturbación negativa mínima lleva el sistema desde $0{,}5$ hasta el atractor $-2$? (c) Explique por qué, con $\mu=-0{,}25$, ninguna perturbación puede devolver el sistema a un equilibrio cercano a $0{,}5$.

Referencias

Boettiger, C. y Hastings, A. (2012). Quantifying limits to detection of early warning for critical transitions. Journal of the Royal Society Interface, 9(75), 2527–2539. https://doi.org/10.1098/rsif.2012.0125

Guckenheimer, J. y Holmes, P. (1983). Nonlinear oscillations, dynamical systems, and bifurcations of vector fields. Springer. https://doi.org/10.1007/978-1-4612-1140-2

Hairer, E., Nørsett, S. P. y Wanner, G. (1993). Solving ordinary differential equations I: Nonstiff problems (2.ª ed. rev.). Springer. https://doi.org/10.1007/978-3-540-78862-1

Kuznetsov, Y. A. (2023). Elements of applied bifurcation theory (4.ª ed.). Springer. https://doi.org/10.1007/978-3-031-22007-4

Scheffer, M., Bascompte, J., Brock, W. A., Brovkin, V., Carpenter, S. R., Dakos, V., Held, H., van Nes, E. H., Rietkerk, M. y Sugihara, G. (2009). Early-warning signals for critical transitions. Nature, 461(7260), 53–59. https://doi.org/10.1038/nature08227

Strogatz, S. H. (2015). Nonlinear dynamics and chaos: With applications to physics, biology, chemistry, and engineering (2.ª ed.). Westview Press.

van de Leemput, I. A., Wichers, M., Cramer, A. O. J., Borsboom, D., Tuerlinckx, F., Kuppens, P., van Nes, E. H., Viechtbauer, W., Giltay, E. J., Aggen, S. H., Derom, C., Jacobs, N., Kendler, K. S., van der Maas, H. L. J., Neale, M. C., Peeters, F., Thiery, E., Zachar, P. y Scheffer, M. (2014). Critical slowing down as early warning for the onset and termination of depression. Proceedings of the National Academy of Sciences, 111(1), 87–92. https://doi.org/10.1073/pnas.1312114110

Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S. J., Brett, M., Wilson, J., Millman, K. J., Mayorov, N., Nelson, A. R. J., Jones, E., Kern, R., Larson, E., … SciPy 1.0 Contributors. (2020). SciPy 1.0: Fundamental algorithms for scientific computing in Python. Nature Methods, 17(3), 261–272. https://doi.org/10.1038/s41592-019-0686-2