Un modello matematico per la ME/CFS

Un modello matematico per la ME/CFS

La versione in inglese di questo articolo è disponibile qui.

Introduzione

Molti dei miei lettori sono probabilmente a conoscenza dei tentativi attualmente fatti per simulare matematicamente il metabolismo energetico dei pazienti ME/CFS, integrando i dati metabolici con i dati genetici. In particolare, il dr. Robert Phair ha sviluppato un modello matematico delle principali vie metaboliche coinvolte nella conversione dell’energia, dall’energia immagazzinata nei legami chimici di grandi molecole come glucosio, acidi grassi e amminoacidi, all’energia immagazzinata nell’adenosina trifosfato (ATP), pronta per l’uso. Phair, che è un ingegnere, ha determinato le equazioni differenziali che regolano questa enorme quantità di reazioni chimiche e le ha adattate al profilo genetico trovato nei pazienti ME/CFS. Ma già alcuni anni fa due fisici pubblicarono un interessante modello matematico del metabolismo energetico durante e dopo l’esercizio, nei pazienti ME/CFS (Lengert N. et Drossel B. 2015). In quanto segue descriverò questo modello e le sue previsioni e vedremo da vicino queste equazioni differenziali.

Le vie metaboliche che sono state analizzate

Il modello di Lengert e Drossel estende due sistemi di equazioni differenziali precedentemente pubblicati che descrivono il comportamento della glicolisi, del ciclo di Krebs (enormemente semplificato come una singola reazione!), della catena di trasporto degli elettroni mitocondriale (descritta in dettaglio), del sistema della creatina chinasi e della conversione di adenosina difosfato (ADP) in ATP, nei muscoli scheletrici (Korzeniewski B. et Zoladz JA. 2001), (Korzeniewski B. et Liguzinski P. 2004). Gli autori hanno aggiunto equazioni per l’accumulo di lattato e il suo efflusso fuori dalla cellula, per la sintesi de novo di inosina monofosfato (IMP) durante il recupero, per la degradazione dell’adenosina monofosfato (AMP) in IMP, per la degradazione di IMP in inosina e ipoxantina. Tutte le vie coinvolte sono raccolte nella figura 1. Queste reazioni sono descritte da 15 equazioni differenziali e la soluzione è un insieme di 15 funzioni del tempo che rappresentano la concentrazione dei principali metaboliti coinvolti (come il lattato, il piruvato, l’ATP, ecc.). Diamo ora uno sguardo più da vicino a una di queste equazioni e alla struttura generale dell’intero sistema di equazioni.

1-s2.0-S0301462215000630-fx1.jpg
Figura 1. Questa è una rappresentazione schematica dei percorsi metabolici descritti dal modello matematico sviluppato da Lengert e Drossel. In dettaglio: sintesi citosolica e degradazione di ADP, AMP e IMP (a sinistra), via della protein chinasi e glicolisi (centro), catena di trasporto degli elettroni e ciclo TCA (a destra). Da Lengert N. et Drossel B. 2015.
lactate dehydrogenase.PNG
Figure 2. La lattato deidrogenasi è l’enzima coinvolto nella catalisi della conversione del lattato in piruvato. Questa reazione procede in entrambe le direzioni.

Equazioni differenziali per reazioni chimiche

Consideriamo l’equazione utilizzata dagli autori per la reazione catalizzata dalla lattato deidrogenasi (la trasformazione del piruvato in lattato, figura 2) dove si è anche tenuto conto dell’efflusso di lattato dal citosol. L’equazione differenziale è la seguente:

equation.PNG

dove i tre parametri sono determinati sperimentalmente e i loro valori sono

equations.PNG

Il primo descrive l’attività dell’enzima lattato deidrogenasi: più questo parametro è elevato, più l’enzima è attivo. Il secondo descrive la reazione inversa (dal lattato al piruvato). Il terzo è una misura di quanto lattato la cellula è in grado di trasportare al di fuori della sua membrana. Forse il lettore si è reso conto che l’equazione del lattato è una equazione differenziale ordinaria del primo ordine. Si dice “primo ordine” perché nell’equazione compare solo la derivata prima della funzione che dobbiamo determinare (lattato, in questo caso); “ordinario” si riferisce al fatto che il lattato è funzione di una sola variabile (il tempo, in questo caso). Si vede immediatamente che un’equazione come questa può essere scritta come segue:

equation bis.PNG

Supponiamo ora di avere altre due equazioni differenziali di questo tipo, una per il piruvato e una per i protoni (le altre due funzioni del tempo che sono presenti nell’equazione):

equations.PNG

Allora avremmo un sistema di tre equazioni differenziali ordinarie come questo:System.PNG

I valori iniziali delle funzioni che dobbiamo determinare sono raccolti nell’ultima riga: questi sono i valori che le funzioni incognite assumono all’inizio della simulazione (t = 0). In questo caso, questi valori sono le concentrazioni di lattato, piruvato e protoni nel citosol, a riposo. Le tre funzioni del tempo sono chiamate la soluzione del sistema. Questo tipo di sistema di equazioni è un esempio di problema di Cauchy, e sappiamo dalla teoria matematica che non solo ha una soluzione, ma che questa soluzione è unica. Inoltre, mentre questa soluzione  può non essere sempre facilmente trovata con metodi rigorosi, è abbastanza facile risolvere il problema con metodi approssimati, come il  metodo di Runge-Kutta o il metodo di Heun. Detto questo, il sistema di equazioni differenziali ordinarie proposto da Lengert e Drossel per il metabolismo energetico è proprio come quello qui sopra, con l’eccezione che comprende 15 equazioni anziché tre. Quindi, la principale difficoltà in questo tipo di simulazione non è l’aspetto computazionale, ma la determinazione dei parametri (come quelli enzimatici) e dei valori iniziali, che devono essere raccolti dalla letteratura medica o devono essere determinati sperimentalmente, se non sono già disponibili. L’altro problema è come progettare le equazioni: esistono spesso diversi modi per costruire un modello matematico di una reazione chimica o di qualsiasi altro processo biologico.

Il modello matematico della ME/CFS

Come adattiamo ai pazienti ME/CFS un modello del metabolismo energetico che è stato impostato con parametri presi da esperimenti condotti su soggetti sani? Questa è un’ottima domanda, e abbiamo visto che Robert Phair ha dovuto usare i dati genetici dei pazienti ME/CFS relativi agli enzimi chiave del metabolismo energetico, al fine di impostare il suo modello. Ma questi dati non erano disponibili quando Lengert e Drossel hanno progettato le loro equazioni. E allora? I due fisici hanno cercato studi sulla fosforilazione ossidativa nei pazienti ME/CFS e hanno scoperto che qusto processo cellulare era stato misurato con diverse impostazioni sperimentali e da diversi gruppi e che il denominatore comune di tuti gli studi era una riduzione di funzione che andava da circa il 35% (Myhill S et al. 2009), (Booth, N et al 2012), (Argov Z. et al. 1997), (Lane RJ. et al. 1998) a circa il 20% (McCully KK. et al. 1996), (McCully KK. et al. 1999). Quindi l’idea degli autori è stata di moltiplicare i parametro enzimatici di ciascuna reazione appartenente alla fosforilazione ossidativa per un numero compreso tra 0,6 (grave ME / CFS) a 1,0 (persona sana). In particolare, i due fisici hanno scelto un valore di 0,7 per la ME/CFS, nei loro esperimenti in silico (cioè esperimenti virtuali condotti nel processore di un computer).

Previsioni del modello matematico

Il modello matematico è stato utilizzato per eseguire prove di esercizio in silico con varie lunghezze e intensità. Quello che Lengert e Drossel hanno trovato è stato che il tempo di recupero nel paziente ME/CFS medio era sempre maggiore se confrontato con quelli di una persona sana. Il tempo di recupero è definito come il tempo necessario affinché una cellula ripristini il suo contenuto di ATP (97% del livello in stato di riposo) dopo lo sforzo. Nella figura 3 si vedono i risultati della simulazione per un esercizio molto breve (30 secondi) e molto intenso. Come potete vedere, nel caso di una cellula sana (a sinistra) il tempo di recupero è di circa 600 minuti (10 ore) mentre una cellula di una persona con ME/CFS (a destra) richiede più di 1500 minuti ( 25 ore) per recuperare.

half minute 1.png
Figura 3. Risultati della simulazione per un esercizio con una durata di 30 secondi e un’intensità elevata (consumo iniziale di ATP 300 volte il valore di riposo). A sinistra, il caso di una cellula muscolare scheletrica sana, a destra il caso di una cellula di una persona con ME/CFS le cui reazioni mitocondriali hanno una velocità ridotta al 70% della velocità del controllo sano. I grafici li ho ottenuti utilizzando la versione online del software, disponibile qui.

Un altro risultato interessante della simulazione è un aumento di AMP nei pazienti rispetto al controllo (figura 3, linea arancione). Ciò è dovuto all’uso compensativo delle due vie metaboliche in figura 4: la reazione catalizzata dall’adenilato chinasi, in cui due molecole di ADP sono utilizzate per produrre una molecola di ATP e una molecola di AMP; e la reazione catalizzata dalla deaminasi AMP, che degrada AMP in IMP (che viene quindi convertito in inosina e ipoxantina). Queste due reazioni sono utilizzate dai pazienti ME/CFS più che dal controllo sano, al fine di aumentare la produzione di ATP al di fuori dei mitocondri.

adenylate-kinase-and-amp-deaminase1.png
Figura 4. La via metabolica a sinistra è utilizzata dai pazienti ME/CFS più che nel controllo per aumentare la produzione di ATP al di fuori dei mitocondri, secondo questo modello matematico. Il percorso sulla destra degrada l’AMP in IMP.

Se diamo un’occhiata più da vicino alle concentrazioni di AMP e IMP nelle 4 ore successive allo sforzo (figura 5), vediamo effettivamente una maggiore produzione di IMP (linea verde) e AMP (linea arancione) nei muscoli scheletrici dei pazienti (destra) rispetto ai controlli (sinistra).

half minute 3.png
Figura 5. Lo stesso della figura 3, ma ingrandito per dare uno sguardo più da vicino alle concentrazioni durante le 4 ore successive allo sforzo. La cellula sana è a sinistra, mentre le cellule di una persona con ME/CFS sono sulla destra.

Un’ulteriore via di compensazione utilizzata dai pazienti (secondo questo modello) è la produzione di ATP da ADP da parte dell’enzima creatina chinasi (figura 6). Questo è un altro modo che abbiamo per produrre ATP nel citosol senza l’aiuto dei mitocondri. In questo modello di ME/CFS, vi è un aumento nell’uso di questo percorso, che porta a una diminuzione della concentrazione cellulare di fosfocreatina e un aumento della concentrazione cellulare di creatina (figura 7).

creatine kinase
Figura 6. La reazione catalizzata dalla creatina chinasi: una molecola di ADP viene convertita in ATP grazie al gruppo fosfato trasportato dalla fosfocreatina.
half minute 4.png
Figura 7. La concentrazione di fosfocreatina nel citosol delle cellule muscolari scheletriche è inferiore nella ME/CFS (a destra) rispetto al controllo (a sinistra) durante e dopo l’esercizio. Ciò è dovuto al maggiore uso di questa molecola per produrre ATP in modo anaerobico nel metabolismo ME/CFS rispetto al controllo. I parametri per questa simulazione sono gli stessi descritti nella figura 3.

Confronto con i dati metabolici disponibili

Sono curioso di vedere se i dati dei vari studi metabolomici condotti dopo la pubblicazione del modello di Lengert e Drossel sono coerenti con le previsioni del modello stesso. Discuterò questo argomento in un altro articolo perché devo ancora studiare questo aspetto. Vorrei solo sottolineare che se ritenessimo vero l’alto tasso di degradazione dell’IMP proposto in questo modello, probabilmente troveremmo un alto livello di ipoxantina nel sangue dei pazienti, rispetto ai controlli, mentre questo metabolita è diminuito nei pazienti, secondo uno studio (Armstrong CW et al. 2015).

Annunci

Multivariate normal distribution, a proof of existence

Multivariate normal distribution, a proof of existence

For my brother who was

decades ahead of me

in terms of knowledge

albeit being only two years older.

 

Introduction

In a previous blog post, I have discussed the joint density of a random vector of two elements that are normally distributed. I was able to prove the expression for the joint probability, not without fighting against some nasty integrals. In the end, I introduced the expression of the joint probability for a random vector of m normally distributed elements and I left my four readers saying “I have no idea about how it could be proved“. We were in June, I was in the North of Italy back then, hosted by friends but mainly alone with my books in a pleasant room with pink walls and some dolls listening to my speculations without a blink; a student of engineering was sharing with me, via chat, her difficulties with the surprisingly interesting task of analyzing the data from some accelerometers put in the mouth of fat people while they stand on an oscillating platform; my birthday was approaching and I was going to come back in Rome after a short stop in Florence, where I was for the first time fully aware of how astonishingly beautiful a woman in a dress can be (and where I saw the statues that Michelangelo crafted for the Tomb of the Medicis, the monument to Lorenzo in particular, which is sculpted in our culture, more profoundly than we usually realize).

But on a day in late December, while I was planning my own cryopreservation (a thought I often indulge in when my health declines even further), I realized that the covariance matrix is a symmetrical one so it can be diagonalized, and this is the main clue in order to prove the expression of this density. As obvious as it is, I couldn’t think of that when I first encountered the multivariate normal distribution, and the reason for this fault is my continuous setbacks, the fact that for most of the last 20 years I have not only been unable to study but even to think and read. And this is also the reason why I write down these proofs in my blog: I fear that I will leave only silence after my existence, because I have not existed at all, due to my encephalopathy. I can’t do long term plans, so as soon as I finish a small project, such as this proof, I need to share it because it might be the last product of my intellect for a long time. So, what follows is mainly a proof of my own existence, more than it is a demonstration of the multivariate normal distribution.

Before introducing the math, two words about the importance of the multivariate normal distribution. Many biological parameters have a normal distribution, so the normal density is the most important continuous distribution in biology (and in medicine). But what happens when we are considering more than one parameter at the time? Suppose to have ten metabolites that follow a normal distribution each, and that you want to calculate the probability that they are all below ten respective maximal values. Well, you have to know about the multivariate normal distribution! This is the reason why I believe that anyone who is interested in biology or medicine should, at least once in her lifetime, go through the following mathematical passages.

Can’t solve a problem? Change the variables!

In this paragraph, I present a bunch of properties that we need in order to carry out our demonstration. The first one derives directly from the theorem of change of variables in multiple integrals. The second and the third ones are a set of properties of symmetrical matrices in general, and of the covariance matrix in particular. Then, I collect a set of integrals that have been introduced or calculated in the already cited blog post about the bivariate normal distribution. The last proposition is not so obvious, but I won’t demonstrate it here, and those who are interested in its proof, can contact me.

Prop 1_2.PNG
Figure1. The domains of the bijective function Y = Φ(X).

PROPOSITION 1 (change of variables). Given the continuous random vector X = (X_1, X_2, …, X_m) and the bijective function Y = Φ(X) (figure 1), where Y is a vector with m dimensions, then the joint density of Y can be expressed through the joint density of X:

Prop 1_1.PNG

 

 

1)

 

PROPOSITION 2 (symmetrical matrix). Given the symmetrical matrix C, we can always write:

Prop 1_3.PNG

 

2)

 

where λ_1, λ_2, …, λ_m are the eigenvalues of matrix C and the columns of P are the respective eigenvectors. It is also easy to see that for the inverse matrix of C we have:

Prop 1_4.PNG

 

3)

 

Moreover, the quadratic form associated with the inverse matrix is

Prop 2_3.PNG

 

4)

where

Prop 2_4.PNG

 

5)

 

PROPOSITION 3 (covariance matrix). If C is the covariance matrix of the random vector X = (X_1, X_2, …, X_m), which means that

Prop 3_1.PNG

 

6)

 

then, with the positions made in Prop. 2, we have

Prop 3_2.PNG

 

7)

 

where σ_j is the standard deviation of X_j and ρ_i,j is the correlation coefficient between X_i and X_j.

PROPOSITION 4 (some integrals). It is possible to calculate the integrals in the following table. Those who are interested in how to calculate the table can contact me.

exponential integral 4

PROPOSITION 5 (other integrals). It is possible to calculate the two following integrals from the table. Those who are interested in how to calculate them can contact me.

exponential integral 3

 

8)

PROPOSITION 6 (sum of normal random variables). Given the random vector X = (X_1, X_2, …, X_m) whose components are normally distributed, then the density of the random variable Y = X_1 + X_2 + … + X_m is a normal law whose average and standard deviations are respectively given by:

extra.PNG

Multivariate normal distribution

PROPOSITION 7. The joint probability density in the case of a random vector whose m components follow a normal distribution is:Prop 6_1.PNG

 

 

9)

 

Demonstration, first part. The aim of this proof will be to demonstrate that if we calculate the marginal distribution of X_i from the given joint distribution, we obtain a normal distribution with an average given by μ_i. Moreover, we will prove that if we use this joint distribution to calculate the covariance between X_i and X_j, we obtain σ_iσ_jρ_i,j (I have to apologize with the reader for this weird way of writing subscripts, but WordPress doesn’t provide an equation editor). We start operating the following change of variables:

10)  Prop 6_2.PNG

 

whose Jacobian is the identity matrix. So we obtain for the joint density in Eq. 9 the expression:

Prop 6_5.PNG

 

 

11)

 

 

We then consider the substitution Prop 6_3.PNG

 

12)

 

whose Jacobian is the determinant of P which is again the identity matrix, since P is an orthogonal matrix (P is the matrix introduced in Prop. 2, whose columns are eigenvectors of the covariance matrix). Then we have

Prop 6_6.PNG

And, according to Prop. 1, we obtain for the joint distribution in Eq. 9 the expression:

prop 6_4

 

13)

 

So, the marginal distribution of the first random variable is

prop 6_7

We recognize the integrals in Prop. 4, for n = 0. So we have for the marginal distribution:

Prop 6_9.PNG

 

14)

while the joint distribution becomes

prop 6_8

 

15)

 

Let’s now consider another change of variable, the following one:Prop 6_10.PNG

 

16)

 

whose Jacobian is given by:

Prop 6_11.PNG

Then, according to Prop. 1, we have

Prop 6_12.PNG

This proves that the variables X_1”, X_2”, … ,X_m” are independent. But they are also normally distributed random variables whose average is zero and whose standard deviation is

Prop 6_13.PNG

for i that goes from 1 to m. Since we have

Prop 6_14.PNG

we can calculate the marginal distribution of ξ_j according to Prop. 6:

Prop 6_15.PNG

 

17)

 

Remembering the very first substitution (Eq. 10) we then draw the following conclusion:Prop 6_16.PNG

 

18)

 

Now, if you remember Prop. 3, you can easily conclude that the marginal density of X_j is, in fact, a normal distribution with average given by μ_j and standard deviation given by σ_j. This concludes the first part of the demonstration. It is worth noting that we have calculated, in the previous lines, a very complex integral (the first collected in the following paragraph), and we can be proud of ourselves.

Demonstration, second part. We have now to prove that the covariance coefficient of between X_i and X_j, is given by ρ_i,j. In order to do that, we’ll calculate the covariance between X_i and X_j with the formula

19)                     Cov[X_i, X_j] = E[X_i×X_j]E[X_i]×E[X_j] = E[X_i×X_j]μ_i×μ_j

For E[X_i×X_j] we have

Prop 6_18.PNG

Considering the substitution in Eq. 10 we have

prop 6_19

prop 6_20

To simplify the writing, let’s assume i=1 and j=2. For I_1 we have:

prop-6_21.png

Now, considering again Prop. 4, we easily recognize that:

Prop 6_22.PNG

 

 

20)

 

 

So, the integral I_1 becomes:

Prop 6_23.PNG

 

21)

 

For I_2 we have:

Prop 6_24.PNG

So, I_2 is zero and the same applies to I_3, as the reader can easily discover by herself, using Eq. 20. Hence, we have found:

Prop 6_25.PNG

Now, just consider Eq. 7 (the second one) in Prop. 3, and you will recognize that we have found

Prop 6_26.PNG

which is exactly what we were trying to demonstrate. The reader has likely realized that we have just calculated another complex integral, the second one in the following paragraph. It can be also verified that the joint density is, in fact, a density: in order for that to be true it must be

prop 6_29

Now, if we use the substitutions in Eq. 10 and in Eq. 12 we obtain:

Prop 6_30.PNG

And our proof is now complete.

Integrals

Prop. 7 is nothing more than the calculation of three very complex integrals. I have collected these results in what follows. Consider that you can substitute the covariance matrix with any symmetrical one, and these formulae still hold.

prop 6_17

Prop 6_27.PNGProp 6_28.PNG