Equazioni differenziali applicate alla biologia: 3 modelli svolti

Quando osserviamo la diffusione di un virus, le oscillazioni tra predatori e prede oppure la crescita di una massa di cellule tumorali, può sembrare difficile individuare una regola che descriva l’evoluzione del fenomeno nel tempo. Le equazioni differenziali permettono di affrontare proprio questo problema, mettendo in relazione lo stato di un sistema biologico con la velocità con cui esso cambia e costruendo modelli quantitativi dai quali ricavare previsioni, condizioni di equilibrio e possibili scenari futuri.

Le applicazioni sono numerose e riguardano ambiti diversi, dall’epidemiologia all’ecologia, fino allo studio della crescita cellulare. Il modello SIR, per esempio, descrive come varia il numero di suscettibili, infetti e rimossi durante un’epidemia; il sistema di Lotka-Volterra rappresenta l’interazione tra due popolazioni, una costituita dalle prede e l’altra dai predatori; i modelli logistico e di Gompertz, infine, consentono di descrivere curve di crescita biologica che rallentano progressivamente nel tempo.

In questo articolo analizziamo tre esercizi svolti, accompagnando i calcoli con l’interpretazione dei risultati e con alcune considerazioni sui limiti dei modelli. L’obiettivo è comprendere come le equazioni differenziali traducano ipotesi biologiche in relazioni matematiche e quali informazioni sia davvero possibile ricavare dai dati disponibili.

Esercizio 1 — Modello SIR per un focolaio influenzale in un campus

Contesto e ipotesi

Consideriamo un campus chiuso, composto da [math]N=5000[/math] persone, nel quale si diffonde un’infezione simil-influenzale. All’inizio dell’osservazione sono presenti [math]I_0=10[/math] individui infetti, mentre le persone suscettibili sono [math]S_0=4990[/math]; supponiamo inoltre che non vi siano individui già guariti o rimossi, per cui [math]R(0)=0[/math].

Utilizziamo il modello SIR, nel quale la popolazione viene suddivisa in tre gruppi: suscettibili [math]S(t)[/math], infetti [math]I(t)[/math] e rimossi [math]R(t)[/math], categoria che comprende, secondo le ipotesi del modello, le persone guarite e immuni oppure non più in grado di trasmettere l’infezione.

Il sistema differenziale è

[math]\displaystyle \begin{aligned} \frac{dS}{dt}&=-\frac{\beta SI}{N}, \\ \frac{dI}{dt}&=\frac{\beta SI}{N}-\gamma I, \\ \frac{dR}{dt}&=\gamma I. \end{aligned}[/math]

Assumiamo un tasso di trasmissione [math]\beta=0{,}45\ \text{giorni}^{-1}[/math] e un tasso di rimozione [math]\gamma=0{,}20\ \text{giorni}^{-1}[/math]. Questi valori sono ipotetici e servono a illustrare il funzionamento del modello, non a descrivere una specifica epidemia reale.

Richieste

• Calcolare il numero di riproduzione di base e la soglia teorica di immunità di gregge.
• Determinare quanti suscettibili sono presenti al picco epidemico.
• Stimare il numero massimo di infetti simultanei sfruttando una quantità conservata del modello.
• Calcolare la dimensione finale dell’epidemia.
• Descrivere l’andamento della curva degli infetti nel tempo.

Soluzione

Il numero di riproduzione di base

Il numero di riproduzione di base, indicato con [math]R_0[/math], rappresenta il numero medio di infezioni secondarie prodotte da un individuo infetto in una popolazione inizialmente suscettibile, in assenza di interventi e nelle condizioni previste dal modello.

Nel sistema SIR considerato vale

[math]\displaystyle \begin{aligned} R_0=\frac{\beta}{\gamma}. \end{aligned}[/math]

Sostituendo i parametri assegnati otteniamo

[math]\displaystyle \begin{aligned} R_0=\frac{0{,}45}{0{,}20}=2{,}25. \end{aligned}[/math]

Poiché il valore è maggiore di uno, l’infezione può inizialmente diffondersi nella popolazione. La crescita effettiva dipende tuttavia anche dalla quota di persone ancora suscettibili, dalle condizioni iniziali e dalle ipotesi epidemiologiche adottate.

Nel modello omogeneo, la soglia teorica di immunità di gregge è

[math]\displaystyle \begin{aligned} p_c=1-\frac{1}{R_0}. \end{aligned}[/math]

Pertanto,

[math]\displaystyle \begin{aligned} p_c=1-\frac{1}{2{,}25}\approx0{,}5556. \end{aligned}[/math]

La soglia è quindi pari a circa il [math]55{,}6\%[/math]. Questo risultato rappresenta una stima teorica valida sotto ipotesi semplificate, tra cui immunità efficace e contatti omogenei; non costituisce automaticamente una percentuale vaccinale sufficiente in una popolazione reale, nella quale possono intervenire eterogeneità dei contatti, efficacia imperfetta dei vaccini e variazioni della trasmissibilità.

Il numero di suscettibili al picco epidemico

Il picco degli infetti si verifica quando la derivata di [math]I(t)[/math] si annulla e la curva degli infetti passa dalla crescita alla diminuzione. Dalla seconda equazione del sistema otteniamo

[math]\displaystyle \begin{aligned} \frac{dI}{dt}=I\left(\frac{\beta S}{N}-\gamma\right). \end{aligned}[/math]

Quando [math]I>0[/math], la condizione [math]\frac{dI}{dt}=0[/math] implica

[math]\displaystyle \begin{aligned} \frac{\beta S}{N}=\gamma, \end{aligned}[/math]

da cui segue

[math]\displaystyle \begin{aligned} S^*=\frac{\gamma N}{\beta}=\frac{N}{R_0}. \end{aligned}[/math]

Nel nostro caso,

[math]\displaystyle \begin{aligned} S^*=\frac{5000}{2{,}25}\approx2222{,}2. \end{aligned}[/math]

Il picco si verifica dunque quando i suscettibili sono circa 2222. Da quel momento, nell’ambito del modello, la riduzione della popolazione suscettibile rende la trasmissione insufficiente a compensare le persone che escono dal gruppo degli infetti.

La stima del picco degli infetti

Per ricavare il numero massimo di infetti senza integrare numericamente tutto il sistema, possiamo utilizzare una relazione che rimane costante lungo le traiettorie del modello SIR.

Dividendo la prima equazione per la seconda e semplificando [math]I[/math], si ottiene

[math]\displaystyle \begin{aligned} \frac{dI}{dS}=-1+\frac{\gamma N}{\beta S}. \end{aligned}[/math]

Integrando rispetto a [math]S[/math],

[math]\displaystyle \begin{aligned} I+S-\frac{N}{R_0}\ln S=C, \end{aligned}[/math]

dove [math]C[/math] è una costante. Valutando la relazione all’inizio e al picco, quando [math]S=S^*=N/R_0[/math], ricaviamo

[math]\displaystyle \begin{aligned} I_{\text{picco}}=I_0+S_0-S^*+\frac{N}{R_0}\ln\left(\frac{S^*}{S_0}\right). \end{aligned}[/math]

Sostituendo i valori numerici,

[math]\displaystyle \begin{aligned} I_{\text{picco}}&=10+4990-2222{,}22 \\ &\quad+2222{,}22\ln\left(\frac{2222{,}22}{4990}\right). \end{aligned}[/math]

Poiché

[math]\displaystyle \begin{aligned} \ln\left(\frac{2222{,}22}{4990}\right)\approx-0{,}809, \end{aligned}[/math]

otteniamo

[math]\displaystyle \begin{aligned} I_{\text{picco}}\approx980. \end{aligned}[/math]

Il modello stima quindi un massimo di circa 980 persone infette contemporaneamente. Il risultato è una previsione matematica basata sui parametri ipotizzati e non va interpretato come un valore certo per un focolaio reale.

La dimensione finale dell’epidemia

Quando l’epidemia si esaurisce, il numero di infetti tende a zero. Nel modello chiuso, in assenza di reinfezioni e con [math]R(0)=0[/math], la relazione finale è

[math]\displaystyle \begin{aligned} \ln\left(\frac{S_\infty}{S_0}\right)=-R_0\left(1-\frac{S_\infty}{N}\right). \end{aligned}[/math]

Questa equazione non si risolve in generale con operazioni algebriche elementari e richiede un metodo numerico. Introducendo la variabile adimensionale [math]x=S_\infty/N[/math], possiamo scriverla come

[math]\displaystyle \begin{aligned} \ln x-\ln\left(\frac{S_0}{N}\right)+R_0(1-x)=0. \end{aligned}[/math]

La soluzione numerica coerente con i dati è

[math]\displaystyle \begin{aligned} x\approx0{,}1461, \end{aligned}[/math]

da cui

[math]\displaystyle \begin{aligned} S_\infty\approx0{,}1461\cdot5000\approx731. \end{aligned}[/math]

Il numero finale di persone rimosse è pertanto

[math]\displaystyle \begin{aligned} R_\infty=N-S_\infty\approx4269. \end{aligned}[/math]

La quota della popolazione che è transitata attraverso lo stato infetto durante l’epidemia, nel modello, è quindi circa

[math]\displaystyle \begin{aligned} \frac{R_\infty}{N}\cdot100\approx85{,}4\%. \end{aligned}[/math]

È importante distinguere questa percentuale dalla soglia teorica di immunità di gregge: il fatto che l’epidemia si arresti quando la popolazione suscettibile è diminuita non significa che il numero totale di persone infettate coincida con la soglia necessaria a interrompere la crescita iniziale.

Come si presenta il grafico degli infetti

La curva [math]I(t)[/math] cresce inizialmente perché il numero di nuove infezioni supera quello delle persone che guariscono o vengono rimosse dalla trasmissione. Raggiunge poi un massimo quando [math]S=N/R_0[/math] e comincia a diminuire, fino a tendere a zero quando l’infezione si esaurisce.

La curva è quindi tipicamente unimodale, con una fase iniziale di crescita, un picco e una fase di decrescita. La posizione temporale del massimo non può essere ricavata dai soli calcoli precedenti senza risolvere il sistema differenziale o disporre di ulteriori informazioni.

Osservazioni e limiti del modello

Il modello SIR suppone che la popolazione sia chiusa, che gli individui abbiano contatti omogenei e che i tassi di trasmissione e rimozione rimangano costanti. Non considera esplicitamente il periodo di incubazione, le reinfezioni, la struttura dei contatti o l’introduzione di misure di contenimento.

Se [math]R_0<1[/math], in una popolazione inizialmente tutta suscettibile la crescita iniziale degli infetti non può sostenersi e [math]I(t)[/math] diminuisce. Se invece il numero iniziale di infetti è molto piccolo ma positivo e [math]R_0>1[/math], la crescita può comunque verificarsi, anche se il tempo necessario per raggiungere numeri elevati può essere maggiore. Se [math]I_0=0[/math], in assenza di introduzioni successive il modello rimane privo di infezione.

Domande di riflessione

• Perché il picco degli infetti si verifica quando [math]S=N/R_0[/math], e quale significato epidemiologico assume questa condizione?
• Quali differenze ci si deve attendere tra la previsione del modello SIR e l’andamento di un’epidemia reale, nella quale i contatti e la trasmissibilità possono variare nel tempo?

Esercizio 2 — Modello predatore-preda di Lotka-Volterra

Contesto e ipotesi

Consideriamo un sistema ecologico sperimentale nel quale una popolazione di prede [math]x(t)[/math] interagisce con una popolazione di predatori [math]y(t)[/math]. Il modello classico di Lotka-Volterra descrive la crescita delle prede in assenza di predatori e la diminuzione dei predatori in assenza di cibo, introducendo un termine di interazione che rappresenta gli incontri tra le due popolazioni.

Il sistema è

[math]\displaystyle \begin{aligned} \frac{dx}{dt}&=\alpha x-\beta xy, \\ \frac{dy}{dt}&=\delta xy-\gamma y, \end{aligned}[/math]

dove [math]\alpha=0{,}8\ \text{giorni}^{-1}[/math] rappresenta il tasso di crescita delle prede, [math]\beta=0{,}02[/math] misura l’intensità della predazione, [math]\delta=0{,}01[/math] descrive l’efficienza con cui l’interazione contribuisce alla crescita dei predatori e [math]\gamma=0{,}4\ \text{giorni}^{-1}[/math] è il tasso di mortalità dei predatori. Le popolazioni iniziali sono [math]x(0)=50[/math] e [math]y(0)=30[/math].

Richieste

• Determinare gli equilibri del sistema e classificare l’equilibrio interno.
• Calcolare il periodo delle piccole oscillazioni attorno all’equilibrio.
• Individuare una quantità conservata e interpretarne il significato.
• Descrivere il ritratto di fase e i limiti biologici del modello.

Soluzione

Gli equilibri del sistema

Un equilibrio è una configurazione nella quale entrambe le popolazioni rimangono costanti, almeno finché il sistema conserva le condizioni previste dal modello. Per individuarlo imponiamo che le due derivate siano nulle:

[math]\displaystyle \begin{aligned} x(\alpha-\beta y)=0,\qquad y(\delta x-\gamma)=0. \end{aligned}[/math]

Il primo equilibrio è l’origine,

[math]\displaystyle \begin{aligned} (x^*,y^*)=(0,0), \end{aligned}[/math]

che rappresenta l’assenza di entrambe le popolazioni. Quando entrambe le popolazioni sono positive, le condizioni di equilibrio richiedono invece

[math]\displaystyle \begin{aligned} y^*=\frac{\alpha}{\beta},\qquad x^*=\frac{\gamma}{\delta}. \end{aligned}[/math]

Sostituendo i parametri assegnati otteniamo

[math]\displaystyle \begin{aligned} x^*=\frac{0{,}4}{0{,}01}=40,\qquad y^*=\frac{0{,}8}{0{,}02}=40. \end{aligned}[/math]

L’equilibrio interno è dunque [math](40,40)[/math], nel quale entrambe le popolazioni rimangono costanti perché la crescita naturale delle prede è bilanciata dalla predazione e la mortalità dei predatori è compensata dalle risorse ottenute.

La linearizzazione attorno all’equilibrio

Per studiare il comportamento locale del sistema calcoliamo la matrice Jacobiana, che contiene le derivate parziali delle funzioni che descrivono le due dinamiche:

[math]\displaystyle \begin{aligned} J(x,y)=\begin{pmatrix} \alpha-\beta y & -\beta x \\ \delta y & \delta x-\gamma \end{pmatrix}. \end{aligned}[/math]

Valutandola in [math](40,40)[/math] otteniamo

[math]\displaystyle \begin{aligned} J(40,40)=\begin{pmatrix} 0 & -0{,}8 \\ 0{,}4 & 0 \end{pmatrix}. \end{aligned}[/math]

Gli autovalori si ricavano dall’equazione caratteristica

[math]\displaystyle \begin{aligned} \lambda^2+0{,}32=0, \end{aligned}[/math]

da cui

[math]\displaystyle \begin{aligned} \lambda=\pm i\sqrt{0{,}32}\approx\pm0{,}5657i. \end{aligned}[/math]

Gli autovalori immaginari puri indicano che il sistema linearizzato presenta oscillazioni attorno all’equilibrio, senza una crescita o un decadimento esponenziale dell’ampiezza. Nel sistema non lineare completo, tuttavia, questa informazione da sola non basta a dimostrare che le orbite siano chiuse; per il modello classico di Lotka-Volterra tale proprietà può essere verificata attraverso una quantità conservata.

Il periodo delle piccole oscillazioni

La pulsazione delle oscillazioni linearizzate è

[math]\displaystyle \begin{aligned} \omega=\sqrt{\alpha\gamma}. \end{aligned}[/math]

Il periodo corrispondente è

[math]\displaystyle \begin{aligned} T=\frac{2\pi}{\sqrt{\alpha\gamma}}. \end{aligned}[/math]

Inserendo i parametri,

[math]\displaystyle \begin{aligned} T&=\frac{2\pi}{\sqrt{0{,}8\cdot0{,}4}} \\ &=\frac{2\pi}{\sqrt{0{,}32}}\approx11{,}1\ \text{giorni}. \end{aligned}[/math]

Questo valore descrive il periodo delle oscillazioni nella linearizzazione e approssima quello delle traiettorie vicine all’equilibrio. Nel modello classico, il periodo delle orbite non lineari può variare con la traiettoria, perciò non va interpretato come un periodo universale valido per ogni condizione iniziale.

La quantità conservata

Una funzione che rimane costante lungo le soluzioni del sistema è

[math]\displaystyle \begin{aligned} V(x,y)=\delta x-\gamma\ln x+\beta y-\alpha\ln y, \end{aligned}[/math]

definita per [math]x>0[/math] e [math]y>0[/math]. Per verificarne la conservazione, deriviamo rispetto al tempo:

[math]\displaystyle \begin{aligned} \frac{dV}{dt}=\left(\delta-\frac{\gamma}{x}\right)\frac{dx}{dt}+\left(\beta-\frac{\alpha}{y}\right)\frac{dy}{dt}. \end{aligned}[/math]

Sostituendo le equazioni del modello,

[math]\displaystyle \begin{aligned} \frac{dV}{dt}&=(\delta x-\gamma)(\alpha-\beta y) \\ &\quad+(\beta y-\alpha)(\delta x-\gamma) \\ &=0. \end{aligned}[/math]

La derivata è identicamente nulla, quindi [math]V(x(t),y(t))[/math] conserva il proprio valore lungo ogni traiettoria positiva.

Le curve di livello di questa funzione, per il sistema classico, descrivono orbite chiuse attorno all’equilibrio interno. Le oscillazioni non vengono smorzate e non si amplificano progressivamente: il sistema è neutro rispetto alle perturbazioni dell’equilibrio, nel senso che una traiettoria vicina non converge all’equilibrio, ma continua a percorrere la propria orbita.

Il ritratto di fase e il significato biologico

Nel piano con [math]x[/math] sull’asse orizzontale e [math]y[/math] su quello verticale, le traiettorie positive circondano il punto [math](40,40)[/math]. Quando le prede sono numerose, i predatori dispongono di maggiori risorse e la loro popolazione tende ad aumentare; l’aumento dei predatori intensifica però la predazione, facendo diminuire le prede. La scarsità di prede riduce successivamente le risorse disponibili per i predatori, la cui popolazione diminuisce, consentendo infine alle prede di riprendersi.

Il ciclo si ripete nel modello ideale, ma il comportamento reale delle popolazioni può essere molto diverso. Il modello classico non include una capacità portante per le prede, una risposta funzionale dei predatori alla disponibilità di cibo, ritardi riproduttivi, competizione o variazioni ambientali.

Se si introduce una crescita logistica delle prede, sostituendo il primo termine con

[math]\displaystyle \begin{aligned} \frac{dx}{dt}=\alpha x\left(1-\frac{x}{K}\right)-\beta xy, \end{aligned}[/math]

la conservazione appena dimostrata non vale più, in generale. La stabilità del nuovo equilibrio dipende dai parametri: possono comparire equilibri stabili o instabili e, in determinate condizioni, cicli limite. Non è corretto concludere che l’aggiunta della capacità portante renda sempre le oscillazioni smorzate, perché il risultato dipende dalla struttura del sistema modificato.

Domande di riflessione

• Perché gli autovalori immaginari puri descrivono oscillazioni nel sistema linearizzato, ma non sono sufficienti da soli a dimostrare la chiusura delle orbite del sistema non lineare?
• In che modo l’introduzione di una capacità portante per le prede può modificare gli equilibri e il ritratto di fase?

Esercizio 3 — Crescita di uno sferoide tumorale: confronto tra modello logistico e modello di Gompertz

Contesto e dati sperimentali

In un esperimento in vitro si misura il volume [math]V(t)[/math] di uno sferoide tumorale durante venti giorni di osservazione. Il volume è espresso in millimetri cubi e le misure sono soggette a variabilità sperimentale. L’obiettivo è confrontare due modelli di crescita comunemente utilizzati per descrivere processi biologici che aumentano inizialmente in modo rapido e rallentano progressivamente.

Il modello logistico è

[math]\displaystyle \begin{aligned} V_{\mathrm{log}}(t)=\frac{K}{1+\frac{K-V_0}{V_0}e^{-rt}}, \end{aligned}[/math]

dove [math]V_0[/math] è il volume iniziale, [math]r[/math] è il tasso di crescita e [math]K[/math] è il livello di saturazione previsto dal modello.

Il modello di Gompertz è

[math]\displaystyle \begin{aligned} V_{\mathrm{G}}(t)=V_0\exp\left[\frac{A}{B}\left(1-e^{-Bt}\right)\right], \end{aligned}[/math]

dove [math]A[/math] rappresenta il tasso di crescita relativo iniziale e [math]B[/math] controlla il rallentamento progressivo della crescita.

I dati sperimentali sono riportati nella tabella.

Richieste

• Stimare i parametri dei due modelli mediante minimi quadrati non lineari.
• Confrontare gli adattamenti usando l’AIC.
• Calcolare il tempo previsto per raggiungere un volume di [math]2{,}0\ \text{mm}^3[/math].
• Discutere l’identificabilità dei parametri e i limiti delle previsioni nella fase iniziale della crescita.

Soluzione

La stima dei parametri mediante minimi quadrati

Per ciascun modello calcoliamo i volumi previsti [math]\widehat V(t_i)[/math] e confrontiamo tali valori con le osservazioni [math]V_i[/math]. La funzione obiettivo è la somma dei quadrati dei residui:

[math]\displaystyle \begin{aligned} RSS=\sum_{i=1}^{n}\left[V_i-\widehat V(t_i)\right]^2. \end{aligned}[/math]

La stima ai minimi quadrati consiste nel trovare i parametri che rendono minima questa quantità. Poiché i due modelli sono non lineari rispetto ai parametri, è opportuno utilizzare un algoritmo numerico, specificando le condizioni iniziali dell’ottimizzazione e gli eventuali vincoli imposti ai parametri.

Un esempio riproducibile in Python è il seguente:

import numpy as np
from scipy.optimize import curve_fit

t = np.arange(0, 22, 2, dtype=float)

V = np.array([
    0.05, 0.08, 0.13, 0.21, 0.32, 0.48,
    0.69, 0.95, 1.25, 1.58, 1.92
])

def logistico(t, K, r, V0):
    return K / (
        1 + ((K - V0) / V0) * np.exp(-r * t)
    )

def gompertz(t, V0, A, B):
    return V0 * np.exp(
        (A / B) * (1 - np.exp(-B * t))
    )

p_log, _ = curve_fit(
    logistico, t, V,
    p0=[3.0, 0.2, 0.05],
    bounds=([1.92, 0.001, 0.001],
            [100, 3, 1]),
    maxfev=100000
)

p_gomp, _ = curve_fit(
    gompertz, t, V,
    p0=[0.05, 0.3, 0.05],
    bounds=([0.001, 0.001, 0.001],
            [100, 3, 2]),
    maxfev=100000
)

def metriche(modello, parametri):
    residui = V - modello(t, *parametri)
    rss = np.sum(residui**2)
    n = len(V)
    p = len(parametri)
    aic = n * np.log(rss / n) + 2 * p
    return rss, aic

print("Logistico:", p_log)
print("RSS e AIC:", metriche(logistico, p_log))

print("Gompertz:", p_gomp)
print("RSS e AIC:", metriche(gompertz, p_gomp))

L’esempio utilizza tre parametri liberi per ciascun modello e vincoli numerici positivi. I risultati possono dipendere dall’algoritmo, dalle inizializzazioni e dai vincoli; per questo motivo è importante riportare le impostazioni utilizzate insieme alle stime.

Con questa procedura, i risultati indicativi sono:

I valori di [math]V_0[/math] sono parametri stimati liberamente e non coincidono necessariamente con la prima misura osservata, pari a [math]0{,}05\ \text{mm}^3[/math]. Se si decide di fissare [math]V_0[/math] a tale valore, bisogna eseguire nuovamente l’ottimizzazione per entrambi i modelli e confrontare i risultati ottenuti con la stessa impostazione.

Il confronto mediante AIC

Il criterio di informazione di Akaike, nella forma utilizzata per errori gaussiani indipendenti con varianza stimata, è

[math]\displaystyle \begin{aligned} AIC=n\ln\left(\frac{RSS}{n}\right)+2p, \end{aligned}[/math]

dove [math]n[/math] è il numero di osservazioni e [math]p[/math] è il numero di parametri stimati secondo la convenzione adottata. A parità di dati e di struttura statistica, un AIC più basso indica un compromesso migliore tra adattamento e complessità.

Nel nostro esempio entrambi i modelli hanno tre parametri liberi. Il modello logistico presenta un RSS inferiore e un AIC di circa 5,7 unità più basso rispetto al modello di Gompertz, perciò risulta preferibile tra i due modelli e con questa specifica procedura di stima.

La differenza non dimostra che il modello logistico rappresenti necessariamente il processo biologico reale, né garantisce che le sue previsioni siano sempre più accurate. Il confronto dipende dai dati disponibili, dalla funzione di errore, dai vincoli di stima e dalle ipotesi statistiche; inoltre, entrambi i modelli possono risultare inadeguati se la crescita osservata presenta caratteristiche non rappresentate dalle rispettive equazioni.

Il tempo necessario per raggiungere [math]2{,}0\ \text{mm}^3[/math]

Per il modello logistico stimato, utilizziamo

[math]\displaystyle \begin{aligned} V(t)=\frac{K}{1+\frac{K-V_0}{V_0}e^{-rt}}. \end{aligned}[/math]

Ponendo [math]V(t)=2[/math], otteniamo

[math]\displaystyle \begin{aligned} 2=\frac{K}{1+\frac{K-V_0}{V_0}e^{-rt}}. \end{aligned}[/math]

Isolando il termine esponenziale e prendendo il logaritmo,

[math]\displaystyle \begin{aligned} t=\frac{1}{r}\ln\left[\frac{K/V_0-1}{K/2-1}\right]. \end{aligned}[/math]

Inserendo i parametri stimati [math]K\approx2{,}966[/math], [math]r\approx0{,}2258[/math] e [math]V_0\approx0{,}0577[/math], si ricava

[math]\displaystyle \begin{aligned} t\approx20{,}6\ \text{giorni}. \end{aligned}[/math]

Il modello logistico prevede quindi il raggiungimento di [math]2{,}0\ \text{mm}^3[/math] poco dopo il ventesimo giorno. Si tratta di un’estrapolazione breve oltre l’ultima osservazione, che vale [math]1{,}92\ \text{mm}^3[/math], ma la previsione resta subordinata all’adeguatezza del modello e all’incertezza dei parametri stimati.

Identificabilità dei parametri e limiti della previsione

Quando i dati riguardano soltanto la fase iniziale della crescita, la curva può apparire quasi esponenziale e fornire poche informazioni sul livello di saturazione [math]K[/math]. In tale situazione, diversi valori di [math]K[/math], accompagnati da differenti valori di [math]r[/math], possono produrre andamenti molto simili nell’intervallo osservato. La stima della capacità portante diventa così incerta e può risultare fortemente correlata agli altri parametri.

Per stimare [math]K[/math] con maggiore affidabilità è utile disporre di osservazioni che comprendano anche la fase di rallentamento e, possibilmente, l’avvicinamento alla saturazione. Se si utilizzassero soltanto i primi quattro punti, la stima del livello massimo sarebbe generalmente molto più dipendente dalle ipotesi del modello, perché nessuna misura mostrerebbe ancora in modo convincente l’eventuale appiattimento della curva.

Occorre infine ricordare che uno sferoide tumorale coltivato in vitro è un sistema sperimentale semplificato. I modelli logistico e di Gompertz descrivono l’andamento del volume, ma non rappresentano esplicitamente fenomeni come vascolarizzazione, necrosi, eterogeneità cellulare, disponibilità di nutrienti o interazione con il sistema immunitario. Le stime non devono quindi essere trasferite direttamente alla crescita di un tumore in un organismo.

Domande di riflessione

• Perché il criterio AIC penalizza i modelli più complessi e quali condizioni devono essere rispettate per confrontare correttamente due modelli?
• Per quale motivo la capacità portante [math]K[/math] risulta difficile da stimare quando le osservazioni coprono soltanto la fase iniziale della crescita?

Che cosa ci insegnano questi tre modelli

I tre casi di studio mostrano come le equazioni differenziali possano descrivere fenomeni biologici differenti, purché siano formulate ipotesi coerenti con il sistema osservato e i parametri siano stimati o scelti con attenzione.

Nel modello SIR, la relazione tra trasmissione e rimozione permette di individuare la condizione del picco epidemico e di stimare la dimensione finale dell’epidemia. Nel sistema di Lotka-Volterra, lo studio degli equilibri e della quantità conservata chiarisce come possano emergere oscillazioni periodiche tra due popolazioni. Nei modelli logistico e di Gompertz, infine, il confronto con i dati mostra che la scelta di una curva non dipende soltanto dall’aspetto grafico, ma anche dalla qualità dell’adattamento, dalla stabilità delle stime e dalla capacità di descrivere il processo biologico.

Un modello matematico non elimina l’incertezza e non restituisce una descrizione completa della realtà: rende esplicite le ipotesi, permette di valutare le conseguenze quantitative di tali ipotesi e offre un metodo per confrontare previsioni e osservazioni. La sua utilità dipende quindi non soltanto dalle equazioni utilizzate, ma anche dalla qualità dei dati e dalla consapevolezza dei limiti entro i quali i risultati possono essere interpretati.

---