# Sistema dinamico della larva del pino

In biologia teorica, il sistema dinamico della larva del pino è un modello matematico rappresentante la dinamica della popolazione di un particolare genere di lepidotteri chiamato in inglese spruce budworm, quello dei Choristoneura (famiglia Tortricidae), che vive nelle foreste del Nord America ed attacca le foglie dell'Abies balsamea per ottenere nutrimento e protezione dagli uccelli predatori della specie Setophaga tigrina presenti nelle stesse foreste.
Ogni 40 anni circa si assiste ad un'esplosione della popolazione con conseguente devastazione dei pini presenti. Una volta consumate quasi completamente le risorse della foresta le larve tornano ad un livello di rifugio dando l'impressione di scomparire dalla foresta.
Il sistema è descritto, in forma adimensionalizzata, dall'equazione differenziale:
    ${\frac {\mathrm {d} N}{\mathrm {d} t}}=rN\left(1-{\frac {N}{K}}\right)-{\frac {N^{2}}{1+N^{2}}}$
dove $N$ rappresenta la popolazione di larve, $r$ il loro tasso di crescita e $K$ il carrying capacity, ovvero le risorse dell'ambiente.

## Descrizione del sistema

Nello studio del sistema si considerano alcune variabili come veloci ed altre come lente. Queste ultime, in una fase iniziale del modello, possono essere considerate come costanti.
La densità delle larve può aumentare di centinaia di volte nel giro di pochi anni. Quindi ha senso considerare la popolazione dei budworm come una variabile veloce. La stessa tempistica si ha per il cambiamento di comportamento degli uccelli predatori, mentre il numero dei volatili si può considerare costante nel tempo.
Gli alberi invece hanno una vita media di 100-150 anni ed il tempo di ripristino totale del fogliame si aggira sui 7-10 anni. Per tale ragione il numero di foglie può essere considerata una variabile lenta.

## Formulazione matematica

Considerando costante la quantità di fogliame con la quale si cibano le larve, un buon modello per il sistema può essere dato dall'equazione differenziale:
    ${\frac {\mathrm {d} N}{\mathrm {d} t}}=r_{B}N\left(1-{\frac {N}{K_{B}}}\right)-P(N)$
ovvero da un sistema che in assenza di uccelli segua un'equazione logistica di parametri $r=r_{B}$ e $K=K_{B}$, più un termine negativo di predazione $P(N)$ dato dalla presenza dei volatili.

## Termine di predazione

Poiché gli uccelli restano in numero pressoché costante e poiché la consumazione di prede per ogni singolo volatile raggiungerà un valore di saturazione (il predatore si sazia), la funzione $P(N)$, all'aumentare della popolazione di larve $N$ si saturerà stabilizzandosi ad un valore limite $B$. In formule:
    $\lim _{N\to \infty }P(N)=B$
Tali uccelli si cibano anche di altre prede, quindi, per basse popolazioni di larve, concentreranno i loro sforzi su altri tipi di cibi più facilmente reperibili. Solo quando le larve del pino cresceranno di numero, i volatili cominceranno a cercare selettivamente questi. Matematicamente possiamo supporre che $P(N)$ diminuisca quadraticamente per $N\to 0$.
Una formulazione semplice per cui valga il termine di saturazione in $B$ e che in un intorno di 0 sia dell'ordine di $N^{2}$ è:
    $P(N)=B{\frac {N^{2}}{A^{2}+N^{2}}}$
dove il termine $A$ rappresenta la scala della densità di larve per cui comincia il fenomeno di saturazione. Matematicamente quando $N=A$ la predazione è metà di quella massima:
    $P(A)=B{\frac {A^{2}}{A^{2}+A^{2}}}={\frac {B}{2}}$

## Adimensionalizzazione

Il modello spruce budworm diventa dunque:
    ${\frac {\mathrm {d} N}{\mathrm {d} t}}=r_{B}N\left(1-{\frac {N}{K_{B}}}\right)-{\frac {BN^{2}}{A^{2}+N^{2}}}$
Dei quattro parametri presenti è possibile, adimensionalizzando il sistema, farne scomparire due. Una scelta comoda è quella di lasciare come parametri quelli relativi al primo membro, in modo da ricondurre lo studio su questa parte alla nota equazione logistica.
Come primo passo si può dividere numeratore e denominatore del termine di predazione per $A^{2}$ ottenendo:
    ${\frac {\mathrm {d} N}{\mathrm {d} t}}=r_{B}N\left(1-{\frac {N}{K_{B}}}\right)-{\frac {B\left({\frac {N}{A}}\right)^{2}}{1+\left({\frac {N}{A}}\right)^{2}}}$
Poiché dimensionalmente si ha:
* $[A]=[N]$
si considera la variabile adimensionale $\,{\frac {N}{A}}\,$ che ci permette di riscrivere:
    ${\frac {\mathrm {d} \left({\frac {N}{A}}\right)}{\mathrm {d} t}}=r_{B}{\frac {N}{A}}\left(1-{\frac {N}{A}}{\dot {\frac {A}{K_{B}}}}\right)-{\frac {B}{A}}{\dot {\frac {\left({\frac {N}{A}}\right)^{2}}{1+\left({\frac {N}{A}}\right)^{2}}}}$
ovvero, ponendo $\,\mu ={\frac {N}{A}}\,$:
    ${\frac {\mathrm {d} \mu }{\mathrm {d} t}}=r_{B}\mu \left(1-A{\frac {\mu }{K_{B}}}\right)-{\frac {B}{A}}{\frac {\mu ^{2}}{1+\mu ^{2}}}$
Dividendo il tutto per $\,{\frac {B}{A}}$ otteniamo:
    ${\frac {\mathrm {d} \mu }{{\frac {B}{A}}\mathrm {d} t}}={\frac {r_{B}A}{B}}\mu \left(1-{\frac {\mu }{\left({\frac {K_{B}}{A}}\right)}}\right)-{\frac {\mu ^{2}}{1+\mu ^{2}}}$
Notando che dimensionalmente si ha:
* $[B]=[t^{-1}N]$
e quindi:
* $\left[{\frac {B}{A}}\right]=\left[{\frac {t^{-1}N}{N}}\right]=[t^{-1}]$
si può definire la variabile adimensionale $\tau ={\frac {B}{A}}\,t$, riscrivendo il sistema come:
    ${\frac {\mathrm {d} \mu }{\mathrm {d} \tau }}={\frac {r_{B}A}{B}}\mu \left(1-{\frac {\mu }{\left({\frac {K_{B}}{A}}\right)}}\right)-{\frac {\mu ^{2}}{1+\mu ^{2}}}$
Notiamo che:
* $[r_{B}]=[t^{-1}]$
* $[K_{B}]=[N]$
possiamo quindi definire i parametri adimensionali $r={\frac {r_{B}A}{B}}$ e $K={\frac {K_{B}}{A}}$, riscrivendo infine il sistema:
    ${\frac {\mathrm {d} \mu }{\mathrm {d} \tau }}=r\mu \left(1-{\frac {\mu }{K}}\right)-{\frac {\mu ^{2}}{1+\mu ^{2}}}$

## Punti d'equilibrio

Cercare i punti d'equilibrio del sistema vuol dire trovare le soluzioni di:
    ${\frac {\mathrm {d} \mu }{\mathrm {d} \tau }}=0$
Lo studio della funzione
$r\mu \left(1-{\frac {\mu }{K}}\right)-{\frac {\mu ^{2}}{1+\mu ^{2}}}$
è però complicato.
Il valore $\mu _{0}=0$, ovvero l'estinzione delle larve, è banalmente punto d'equilibrio. Dallo studio delle derivate si vede inoltre che è un punto d'equilibrio instabile.
Per ricavarsi gli altri eventuali punti fissi si possono studiare le intersezioni delle due parti del modello (dividendo per $\mu$ il cui caso è già stato discusso), ovvero trovando i punti per cui vale:
    $r\left(1-{\frac {\mu }{K}}\right)={\frac {\mu }{1+\mu ^{2}}}$
Nell'espressione di sopra il primo membro rappresenta la crescita pro capita della variabile adimensionale $\mu$ rispetto al tempo $\tau$, mentre il secondo membro è il tasso di morte pro capita causato dalla predazione, sempre riferito alle variabili adimensionalizzate.
Poiché il numero delle foglie di pino di cui si cibano le larve è pressoché costante, si considera $K$ come costante e si vede come cambia il modello al variare del valore del parametro $r$. In particolare si parte da un $r$ piccolo e si vede cosa succede facendolo aumentare.
Oltre al già dato $\mu _{0}$ i punti di equilibrio ricavabili dall'equazione di sopra, ovvero i punti di intersezione delle due curve, variano in numero da 1 a 3. Infatti vi sono due diverse biforcazioni del tipo saddle node in cui si assiste alla comparsa/scomparsa di punti fissi.

### Un punto d'equilibrio: Rifugio

Per bassi valori di $r$ si ha una sola intersezione $\mu _{-}$ che lo studio geometrico del grafico mostra essere stabile.
Tale punto d'equilibrio, pur spostandosi verso destra, rimane relativamente vicino all'origine.
Biologicamente la popolazione aumenta ma il numero di individui rimane sempre basso. Per tale ragione il punto $\mu _{-}$ è detto rifugio o refuge.

### Prima biforcazione saddle-node

All'aumentare di $r$ si raggiunge un valore critico in cui compare un nuovo punto d'equilibrio semistabile (instabile a sinistra e stabile a destra).
Poiché il punto semistabile non varia numericamente il flusso, la popolazione resta prossima al punto d'equilibrio $\mu _{-}$

### Tre punti d'equilibrio

Aumentando ancora $r$ il nuovo punto di equilibrio si scinde in due diversi punti $\mu _{c}$ instabile ed $\mu _{+}$ stabile.
Pur con la presenza di un nuovo punto fisso nel sistema, la popolazione resta prossima al punto fisso $\mu _{-}$. Ciò è dovuto al fatto che, per raggiungere il bacino d'attrazione di $\mu _{+}$, il numero di individui dovrebbe oltrepassare il punto instabile $\mu _{c}$, ovvero essere sottoposto ad una grande perturbazione.

### Seconda biforcazione saddle-node

All'aumentare di $r$ si raggiunge un secondo valore critico in cui i punti $\mu _{-}$ e $\mu _{c}$ collidono in un unico punto semistabile (stabile a sinistra, instabile a destra).
Il punto $\mu _{-}$ è divenuto numericamente instabile. Biologicamente, quindi, a seguito di ogni seppur minima perturbazione, il sistema esplode a destra, arrestandosi solo in prossimità del punto $\mu _{+}$.
È dovuto a questa scomparsa improvvisa del punto d'equilibrio rifugio che ogni 40 anni circa si verifica l'invasione delle larve nelle foreste. Per tale ragione il punto $\mu _{+}$ è detto anche esplosione o outbreak.

### Un punto d'equilibrio: Esplosione

Per alti valori di $r$ si ha, nuovamente, una sola intersezione. Questa è stavolta $\mu _{+}$ che risulta essere stabile.
Tale punto d'equilibrio, seppure inferiore, è abbastanza vicino alla carrying capacity $K$.

### Ritorno al Rifugio

A questo punto, poiché la popolazione di larve e quindi il consumo di fogliame è di molto superiore rispetto a prima, non si può più considerare $K$ costante.
Il numero di foglie nella foresta, infatti, diminuisce drasticamente e seppur il parametro $r$ resta inalterato, col diminuire di K la retta interseca nuovamente la curva di predazione facendo rispuntare i punti $\mu _{c}$ e $\mu _{-}$ e, infine, facendo scomparire il punto fisso $\mu _{+}$.
A tal punto la popolazione di larve ritorna al valore iniziale di rifugio, dando nuovamente inizio al ciclo.
Da notare il fenomeno dell'isteresi per cui nella zona con due punti d'equilibrio stabili (più un terzo instabile), la popolazione tende al primo punto fisso o al secondo in base alla posizione nei passi precedenti, ovvero a seconda che il parametro $r$ sia in fase di crescita o di diminuzione.