On the module of random vectors

On the module of random vectors

“I have hopes of being able to achieve

something of value through my current studies or

with any new ideas that come in the future.”

J. F. Nash

The bridge over the Arno

In 1999 I was wandering in Pisa with a booklet in the pocket of a worn coat, too short for my frame. That coat was dark blue on the outside, green and red inside, with one mended sleeve and an austere cowl: I was much like a young monk, with his holy book (figure 1). I can remember neither its title nor the author, though. It was an introduction to statistical thermodynamics, with beautiful figures, a coloured cover, and less than 100 pages. It contained the work by Maxwell on the kinetic theory of ideal gasses, along with other material. I borrowed it from the University Library because I was fascinated by the way in which Maxwell was able to describe the properties of gasses with just a few hypotheses and some relatively easy mathematical passages. I felt that there was an enormous attraction in these methods, I realized with pleasure that math could give the power to completely understand and hold in hand physical systems and even, I started speculating, biological ones.

Analisi, coraggio, umiltà, ottimismo.jpg
Figure 1. A small self-portrait I drew in 1999 in an empty page of my pocket-sized French dictionary.

My second favourite composer back then was Gustav Mahler (the favourite one being Basil Poledouris): he represented my own way to classical music and I chose him because he wasn’t among the musicians my father and my brother shared a love for. I felt, during my teens, that I had to find my private space, and I met it one day on a used book stand: a cassette tape of Das Lied von The Erde, with a few sheets containing the translation to Italian of the songs. Mahler was born in 1860, a few weeks after Maxwell published his pivotal work about ideal gasses in The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science (R) (figure 2). But in 1999 I was repeatedly listening to a CD with a collection of songs sung by Edith Piaf and Charles Trenet, because I was studying French, and I was having a hard time with pronunciation. So, imagine a secular monk in his prime who listens to old French songs while keeping one hand on a book of statistical thermodynamics hidden in his pocket, wherever he goes, wandering in the streets of Pisa, a city which gave birth to Galileo Galilei. This seems a beautiful story, much like a dream, right? Wrong.

Figure 2. Original cover of the journal in which Maxwell published his work on the kinetic theory of ideal gasses.

I had already started my struggle against the mysterious disease that would have completely erased my life in the following years. In the beginning, it had a relapsing-remitting course, so I could sometimes host the hope that I was recovering, only to find myself caught by the evil curse again. At the end of the year 1999, I was losing my mind, I knew that and I was also aware that my holy book couldn’t save me. I clearly remember one evening, I was walking on Ponte di Mezzo, a simple and elegant bridge above the Arno, and I felt that I couldn’t feel sorrow for the loss of my mind: I realized that not only the functions of my brain assigned to rational thinking were gone, but my feelings couldn’t function properly either. In fact, I noted without a true manifestation of desperation that I had lost my emotions. One day, after spending in vain about eleven hours on a single page of a textbook of invertebrate palaeontology, I accepted that I couldn’t read anymore, at least for the moment.

Had I known for sure that I wouldn’t have recovered in the following twenty years, I would have quite certainly taken my own life, jumping from a building; a fate that I have been thinking about almost every day ever since. I considered this possibility during the endless sequence of days in which there has been nothing other than the absence of my thoughts.

The distribution of velocities of an ideal gas and the one hundred years gap

In the already mentioned paper by Maxwell, he derived the probability density of the speed of the molecules of a gas, granted that the three components of the vector of speed are independent random variables (hyp. 1) and that they share the same density (hyp. 2), let’s say f. Moreover, the density of the speed has to be a function only of its module (hyp. 3). These three hypotheses together say that there is a function Φ such that


This is a functional equation (i.e. an equation in which the unknown is a function) whose solution is not detailed in Maxwell’s work. But it can be easily solved moving to polar coordinates (see figure 3) and deriving with respect to θ both members (the second one gives naught since it depends only to the distance from the origin).

space of velocities.png
Figure 3. The space of velocities and its polar coordinates.

Another way to solve the functional equation is to use the method of Lagrange’s multipliers, searching for the extremes of the density of the velocity, when its module is fixed. In either case, we obtain the differential equation:

differential equation.png

which leads to the density for each component of the speed:


where σ can’t be determined using only the three hypotheses mentioned above. Considering then the well-known law of ideal gasses (pV=nRT) and an expression for p derived from the hypothesis that the collisions between the molecules of gas and the container are completely elastic, Maxwell was able to conclude that:


where m is the mass of the molecule of gas, T is the absolute temperature and K_B is the Boltzmann’s constant. It was 1860, Mahler’s mother was going to deliver in Kaliště, Charles Darwin had just released his masterpiece “On the origin of species”, forced to publish much earlier than what he had planned because of the letter he had received from Wallace, in which he described about the same theory Darwin had been working on for the previous 20 years. In the same point in time, Italy was completing its bloody process of unification, with the Mille expedition, led by Giuseppe Garibaldi.

But the functional equation I have mentioned at the beginning of this paragraph has brought with it a mystery for years, until 1976, when an employee at General Motors Corporation published a short note in the American Journal of Physics [R] in which he showed how Maxwell’s functional equation is, in fact, an example of the well known Cauchy’s functional equation:


In order to prove that, you just have to consider the following definition:


given that


The name of the mathematician who made this observation is David H. Nash, and he has the merit of finding something new in one of the most known equation of physics, an equation that is mentioned in every book of thermodynamics, an equation that has been considered by millions of students in more than a century. It was 1976, my mother was pregnant with my brother; Alma, Gustav Mahler’s wife, had died about ten years before.

Figure 4. A sphere with a radius of z. The integrals of the density of within this sphere has the obvious meaning of the repartition function of Z.

Module of random vectors

Once Maxwell found the density of probability for each component of the speed of the molecules of an ideal gas, he searched for the density of the module of the speed. There is a relatively simple way of doing that. With the following notation


we have that the repartition function of Z is given by the integrals of the density of X within the sphere in figure 4. We have:

repartition function.PNG

The second expression is the same as above but in polar coordinates. Then we can obtain the density of Z by derivation of the repartition function. And this method can be extended to an m-dimensional space. This was the method used by Maxwell in his paper. And yet, there is another way to obtain the expression of the module of a random vector: I have explored it in the last months, during the rare hours in which I could function. By the way, only in the Summer of 2007 I was able to study the kinetic theory by Maxwell, eight years after I borrowed the holy book. Such a waste.

The hard way towards the density of the module of a random vector

When a student starts studying statistics, she encounters a list of densities: the normal distribution, the gamma distribution, the exponential distribution etc. Then there are several derived distributions that arise when you operate sums, roots extractions etc. on random variables. In particular, if f_X is the densty of X and Y = X², then we have


On the other hand, if Y = √X we have


Another important result that we have to consider is that given

linera combination.PNG


linear combination 2.PNG

By using these results I have been able to find that the expression of the density for  the module of an m-dimensional random vector is:


In particular, for m = 3 we have

module 2.PNG

Bessel surface.PNG
Figure 5. A plot of the modified Bessel function.

The case of normal random vectors: the modified Bessel function

In particular, if the random vector has dimension 3 and its components are normal random variables with the same expected value and variance, we have that the density of its module is given by

modulo normale.PNG

where I_0 is the modified Bessel function, which is one solution of the differential equation:

Bessel equation.PNG

whose name is modified Bessel equation. The integral expression of the modified Bessel function is:

Bessel function.PNG

I have coded a script in Matlab which integrates numerically this function (available here for download) which plots the surface in figure 5 and also gives the following table of values for this function.

Bessel table.PNG

The following is the flowchart of the script I coded.

flux chart.png

The case of normal random vectors with a naught expected value: the upper incomplete gamma function

If we consider random variables with an average that is zero (this is the case with the components of speed in ideal gasses), then the density is given by

Chi 3 scaled.PNG

which is a Chi distribution with 3 degrees of freedom, scaled with a scale parameter given by s = 1/σ. In the expression of the repartition function, ϒ is the lower incomplete gamma function, which is defined as follows:

incomplete gamma.PNG

I have written a code for its numerical integration (available here for download), the output of which is in figure 6.

lower incomlete.PNG
Figure 6. Lower incomplete gamma function.
maxwell surface.PNG
Figure 7. Maxwell-Boltzmann distribution for N2 in function of temperature and speed.
maxwell surface 2.PNG
Figure 8. Maxwell-Boltzmann distribution for Helium, Neon, Argon, at room temperature.

Conclusion, twenty years later

The density of the module of the velocity of the molecules of an ideal gas is, in fact, a scaled Chi distribution with 3 degrees of freedom, and it is given by


It can be numerically integrated with the following script I made for Octave/Matlab, which gives the plot in figure 7. Another similar script gives the plot in figure 8. These plots represent the Maxwell-Boltzmann distribution, the centre of the holy book that an unfortunate boy was carrying in his pocket, all alone, some twenty years ago. He could have easily died by his own hand in one of the several thousand days of mental and physical disability that he had to face alone. Instead, he has survived. Had it been only for finding the Maxwell-Boltzmann distribution following another path, it would have been worth it. But he has found much more, including a bright girl, the magnificent next stage of evolution of human beings.

% file name = legge_Maxwell-Boltzmann_2

% date of creation = 22/02/2019

% it plots the density and the distribution function for the

% Maxwell-Boltzmann distribution considered as a function of temperature

% and speed

clear all

% we define some parameters

K_B = 1.381*10^(-23)                                                 % Boltzmann's constant

m = 4.65*10^(-26)                                                     % mass of the molecule N_2

% we define the array of temperature from 0° C to 1250° C

T (1) = 273.15

for i = 2:1:1250

T (i) = T (i-1) + 1.;


% it defines f_gamma in 3/2

f_gamma = sqrt(pi)/2.

% delta of integration

delta = 1.0

% it defines the array for the abscissa

z (1) = 0.;

for i = 2:1:2500

z (i) = z(i-1)+delta;


% it defines the density

for j = 1:1:1250

% it defines a constant

c = ( m/(K_B*T(j)) );

for i = 1:1:2500

f (j,i) = ( c^1.5 )*sqrt(2./pi)*( z(i)^2. )*( e^( -0.5*c*z(i)^2. ) );


% it calculates the ripartition function for He

F (j,1) = 0.;

F (j,3) = F (j,1) + delta*( f(j,1) + ( 4*f(j,2) ) + f(j,3) )/3;

F (j,2) = F (j,3)*0.5;

for k=2:1:2500-2

F (j,k+2) = F (j,k)+delta*( f(j,k)+( 4*f(j,k+1) )+f(j,k+2) )/3;



% It plots f and F

figure (1)

mesh(z(1:100:2500), T(1:100:1250), f(1:100:1250,1:100:2500));


xlabel('speed (m/s)');

ylabel('temperature (Kelvin)');

grid on

figure (2)

mesh(z(1:100:2500), T(1:100:1250), F(1:100:1250,1:100:2500));


xlabel('speed (m/s)');

ylabel('temperature (Kelvin)');

grid on


Un modello matematico per la ME/CFS

Un modello matematico per la ME/CFS

La versione in inglese di questo articolo è disponibile qui.


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.

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:


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


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):


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.

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.

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.



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





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

Prop 1_3.PNG




where \lambda_1, \lambda_2, ..., \lambda_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




Moreover, the quadratic form associated with the inverse matrix is

Prop 2_3.PNG




Prop 2_4.PNG




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




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

Prop 3_2.PNG




where \sigma_j is the standard deviation of X_j and \rho_{ij} 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



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:


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





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 \mu_i. Moreover, we will prove that if we use this joint distribution to calculate the covariance between X_i and X_j, we obtain \sigma_i\sigma_j\rho_{ij}. 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






We then consider the substitution Prop 6_3.PNG




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




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



while the joint distribution becomes

prop 6_8




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




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




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




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 \mu_j and standard deviation given by \sigma_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 \rho_{ij}. 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_iX_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:


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

Prop 6_22.PNG






So, the integral I_1 becomes:

Prop 6_23.PNG




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.


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

Sum of independent exponential random variables with the same parameter

Sum of independent exponential random variables with the same parameter


Let X_1, X_2, …, X_m be independent random variables with an exponential distribution. We have already found the expression of the distribution of the random variable Y = X_1 + X_2 + …+ X_m when X_1>0, X_2>0, …, X_m>0 have pairwise distinct parameters (here). The reader will easily recognize that the formula we found in that case has no meaning when the parameters are all equal to λ. I will derive here the law of Y in this circumstance. This can be done with a demonstration by induction, with no particular effort, but I will follow a longer proof. Why? Perhaps because I like pain? Yes, this might be true; but the main reason is that, to me,  this longer demonstration is quite interesting and it gives us the chance to introduce the following proposition. In this blog post, we will use some of the results from the previous one on the same topic and we will follow the same enumeration for propositions.

PROPOSITION 8 (sum of m independent random variables). Let X_1, X_2, …, X_m be independent random variables. The law of Y = X_1 + X_2 + …+ X_m is given by:

Proof. We know that the thesis is true for m=2 (PROP. 1). Consider now that:

 But we know that X_1, X_2, …, X_m are independent. So we have:

We now operate the substitution

And we obtain:

But this has to be true for any possible interval [a, b], which means that:

This is one of the equations in the thesis. The other ones can be derived in the same way ♦

Guessing the solution

I will solve the problem for m = 2, 3, 4 in order to have an idea of what the general formula might look like. 

PROPOSITION 9 (m=2). Let X_1, X_2 be independent exponential random variables with the same parameter λ. The law of Y_1 = X_1 + X_2 is given by:

for y > 0. It is zero otherwise.

Proof. We just have to substitute f_{X_1}, f_{X_2} in Prop. 1. We obtain:

And we find the thesis by solving the integral ♦

PROPOSITION 10 (m=3). Let X_1, X_2, X_3 be independent exponential random variables with the same parameter λ. The law of Y = X_1 + X_2 + X_3 is given by:

for y>0. It is zero otherwise.

Proof. If we define Y_1 = X_1 + X_2 and Y_2= X_3, then we can say – thanks to Prop. 2 – that Y_1 and Y_2 are independent. This means that – according to Prop. 1 – we have

And the proof is concluded ♦

PROPOSITION 11 (m=4). Let X_1, X_2, X_3, X_4 be independent exponential random variables with the same parameter λ. The law of Y = X_1 + X_2 + X_3 + X_4 is given by:

Proof. We can demonstrate the thesis in the same way we did in Prop. 10. We can draw the same conclusion by directly applying Prop. 8: 

Figure 1. This tetrahedron is defined by the conditions y - x_2 - x_3 - x_4 >0 and x_i >0 for i = 2, 3, 4. It is the domain of integration for the integral in the proof of proposition 11.

The domain of integration is the tetrahedron in Figure 1. The generic point P within the tetrahedron has to belong to the segment AC that, in turns, belongs to the triangle shaded grey. This means that the domain of integration can be written [0, y]×[0, yx_2]×[0, yx_2x_3]. Thus, we calculate:

and the proof is concluded ♦


The reader has now likely guessed what the density of Y looks like when m is the generic integer number. But before we can rigorously demonstrate that formula, we need to calculate an integral.

PROPOSITION 12 (lemma). It is true that:

Proof. We have already found that the thesis is true when m = 4. Let us now assume that it is true for m – 1; if we write the thesis for the m – 2 variables x_3, x_4, … , x_m we have:

Now we put η = yx_2:

By integrating both members in [0, y] we obtain:

and the proof is concluded ♦

PROPOSITION 13. Let X_1, X_2, …, X_m be independent random variables. The law of Y = X_1 + X_2 + …+ X_m is given by:

Proof. By directly applying Prop. 8 we have:


Thus, we can write:

The domain of integration can be written in a more proper way (as we did in Prop. 11) as follows:

But this is the integral calculated in Prop. 12, and the proof is concluded ♦

A numerical application

The reader might have recognized that the density of Y in Prop. 13 is, in fact, a Gamma distribution: it is the distribution Γ(m, λ), also known as Erlang’s distribution. There might even be a reader who perhaps remembers that I have discussed that distribution in a post of mine (here) and that I have coded a program that plots both the density and the distribution of such a law. In Figure 2 we see the density of Y (left) and the distribution function (right) for λ = 1.5 and for m that goes from 1 to 3 (note that m is indicated with α in that picture).

Figure 2. Density (left) and distribution function (right) of Y for λ = 1.5 and for m that goes from 1 to 3 (note that m is indicated with α in this picture).

Sum of independent exponential random variables

Sum of independent exponential random variables

For Chiara,

who once encouraged me

to boldly keep trying.


For the last four months, I have experienced the worst level of my illness: I have been completely unable to think for most of the time. So I could do nothing but hanging in there, waiting for a miracle, passing from one medication to the other, well aware that this state could have lasted for years, with no reasonable hope of receiving help from anyone. This has been the quality of my life for most of the last two decades. Then, some days ago, the miracle happened again and I found myself thinking about a theorem I was working on in July. And once more, with a great effort, my mind, which is not so young anymore, started her slow process of recovery. I concluded this proof last night. This is only a poor thing but since it is not present in my books of statistics, I have decided to write it down in my blog, for those who might be interested.

I can now come back to my awkward studies, which span from statistics to computational immunology, from analysis of genetic data to mathematical modelling of bacterial growth. Desperately searching for a cure.

The problem

Let X_1, X_2, ..., X_m be independent random variables with an exponential distribution with pairwise distinct parameters \lambda_1, \lambda_2, ..., \lambda_m , respectively. Our problem is: what is the expression of the distribution of the random variable Y = X_1 + X_2 + ... + X_m? I faced the problem for m = 2, 3, 4. Then, when I was quite sure of the expression of the general formula of f_Y (the distribution of Y) I made my attempt to prove it inductively. But before starting, we need to mention two preliminary results that I won’t demonstrate since you can find these proofs in any book of statistics.

PROPOSITION 1. Let X_1, X_2 be independent random variables. The distribution of  Y = X_1 + X_2 is given by:

where f_X is the distribution of the random vector [X_1, X_2 ].

PROPOSITION 2. Let X_1, X_2, X_3 be independent random variables. The two random variables Y _1 = X_1 + X_2 + ...+ X_n and Y _2 = X_1 + X_2 + ...+ X_m (with n<m) are independent.

DEFINITION 1. For those who might be wondering how the exponential distribution of a random variable X_i with a parameter \lambda_i looks like, I remind that it is given by:

Guessing the solution

As mentioned, I solved the problem for m = 2, 3, 4 in order to understand what the general formula for f_Y might have looked like.

PROPOSITION 3 (m = 2). Let X_1, X_2 be independent exponential random variables with distinct parameters \lambda_1, \lambda_2, respectively. The law of Y _1 = X_1 + X_2 is given by:

Proof. We just have to substitute f_{X_1}, f_{X_2} in Prop. 1. We obtain:

And the demonstration is complete ♦

PROPOSITION 4 (m = 3). Let X_1, X_2, X_3 be independent exponential random variables with pairwise distinct parameters \lambda_1, \lambda_2, \lambda_3, respectively. The law of Y _1 = X_1 + X_2 + X_3 is given by:

Proof. If we define Y _1 = X_1 + X_2 and Y _2 = X_3, then we can say – thanks to Prop. 2 – that Y _1 and Y _2 are independent. This means that – according to Prop. 1 – we have

The reader will now recognize that we know the expression of  f_{Y_1} because of Prop. 3. So, we have: 

For the first integral we find:

For the second one we have:

Hence, we find:

And the thesis is proved

PROPOSITION 5 (m = 4). Let X_1, X_2, X_3, X_4 be independent exponential random variables with pairwise distinct parameters \lambda_1, \lambda_2, \lambda_3, \lambda_4, respectively. The law of Y = X_1 + X_2 + X_3 + X_4 is given by:

for y>0, while it is zero otherwise.

Proof. Let’s consider the two random variables Y_1 = X_1 +  X_2, Y_2 = X_3 + X_4. Prop. 2 tells us that Y_1, Y_2 are independent. This means that – according to Prop. 1 – we can write:

The reader has likely already realized that we have the expressions of f_{Y_2} and f_{Y_2}, thanks to Prop. 3. So we have:

For the four integrals we can easily calculate what follows:

Adding these four integrals together we obtain:

And this proves the thesis

We are now quite confident in saying that the expression of f_Y for the generic value of m is given by:

for y>0, while being zero otherwise. But we aim at a rigorous proof of this expression.


In order to carry out our final demonstration, we need to prove a property that is linked to the matrix named after Vandermonde, that the reader who has followed me till this point will likely remember from his studies of linear algebra. The determinant of the Vandermonde matrix is given by:

PROPOSITION 6 (lemma). The following relationship is true:

Proof. In the following lines, we calculate the determinant of the matrix below, with respect to the second line. In the end, we will use the expression of the determinant of the Vandermonde matrix, mentioned above:

But this determinant has to be zero since the matrix has two identical lines, which proves the thesis

PROPOSITION 7. Let X_1, X_2, ... , X_m be independent exponential random variables with pairwise distinct parameters \lambda_1, \lambda_2, ... , \lambda_m, respectively. The law of Y = X_1 + X_2 + ... + X_m is given by:

for y > 0, while being zero otherwise.

Proof. We already know that the thesis is true for m = 2, 3, 4. We now admit that it is true for m-1 and we demonstrate that this implies that the thesis is true for (proof by induction). Let’s define the random variables Y_1 = X_1 + X_2 + ... + {X_{m-1}} and Y_2 = X_m. These two random variables are independent (Prop. 2) so – according to Prop. 1 – we have:

Now, f_{Y_1} is the thesis for m-1 while f_{Y_2} is the exponential distribution with parameter \lambda_m. So we have:

For the sum we have:

The sum within brackets can be written as follows:

So far, we have found the following relationship:

or, equivalently:

prop 7

In order for the thesis to be true, we just need to prove that

prop 7_1

which is equivalent to the following:

prop 7_3

Searching for a common denominator allows us to rewrite the sum above as follows:

prop 7_4

So, we just need to prove that:

It is possible to show, with some manipulations, that this is equivalent to the thesis of Prop. 6 and this proof is concluded ♦

References. A paper on this same topic has been written by Markus Bibinger and it is available here.


Mark Davis e il test immunitario universale

Mark Davis e il test immunitario universale

1. Introduzione

Queste sono solo alcune note raccolte dal discorso che Mark Davis ha pronunciato in occasione del Community Symposium tenutosi nell’agosto scorso (2017) a Stanford (video). Nei paragrafi 2, 3, 4 e 5 introdurrò alcune nozioni di base sui recettori delle cellule T (T cell receptors: TCR); il paragrafo 6, attraverso riferimenti  al video già menzionato e a tre articoli pubblicati da Davis et al. nel corso degli ultimi quattro anni, descrive una  nuova tecnologica sviluppata da Mark Davis e colleghi. Questi cenni preliminari dovrebbero auspicabilmente fornire i mezzi per comprendere a pieno la portata dei dati pilota presentati da Mark Davis a proposito dell’attività delle cellule T nella ME/CFS (paragrafo 7) e nella malattia di Lyme cronica (paragrafo 8), mostrando perché tale tecnologia prometta di divenire una sorta di test universale per qualsiasi tipo di infezione o patologia autoimmune, nota o sconosciuta.

2. Cellule T

I linfociti T sono una tipologia di leucociti (o globuli bianchi), vale a dire la componente cellulare del nostro sistema immunitario. La gran parte dei linfociti T in circolo è rappresentata da linfociti T helper (T helper cells: Th cells)  e da linfociti T citotossici (cytotoxic T lymphocytes: CTL). Mentre la funzione dei linfociti T helper è quella di regolare l’attività degli altri leucociti attraverso la produzione di un’ampia gamma di trasmettitori chimici (le citochine, cytokines), le CTL sono coinvolte direttamente nella soppressione delle cellule ospiti infette. I linfociti T appartengono al ramo cosiddetto adattivo del sistema immunitario, assieme alle cellule B (le fabbriche di anticorpi), e, in quanto tali, il loro compito è quello di garantire una difesa specifica, su misura, contro gli agenti patogeni: per contrastare uno specifico agente patogeno, il nostro sistema immunitario può schierare in campo non solo anticorpi specifici ma anche specifici linfociti T (Th cells e CTL). Il ramo innato del sistema immunitario, invece, (nel quale rientrano le cellule natural killer, i macrofagi, le cellule dendritiche, i mastociti…) è in grado di fornire soltanto una difesa aspecifica, una prima linea di risposta immunitaria.

3. Recettori dei linfociti T

I linfociti T sono in grado di andare alla ricerca di specifici patogeni grazie a una molecola espressa sopra la propria superficie, chiamata recettore del linfocita T (TCR). Nella figura 1 si può vedere una schematica rappresentazione del TCR e del meccanismo in virtù del quale il linfocita T riconosce il proprio target. Gli antigenti (proteine) degli agenti patogeni vengono indicati ai linfociti T da altre cellule del nostro corpo: vengono esposte sopra molecole chiamate Complesso Maggiore di Istocompatibilità (MHC), che si trova espresso sulla membrana esterna. Se un dato antigene mostra compatibilità con il TCR di uno specifico linfocita T, tale linfocita T si attiva e comincia a proliferare (espansione clonale, clonal expansion). Le due catene principali (α and β) sono assemblate combinando la trascrizione di segmenti di gene, ognuno dei quali ha copie multiple, leggermente diverse fra loro: in altre parole, i TCR vengono assemblati a partire da peptidi scelti a caso da un insieme di diverse alternative possibili. Questo comporta un repertorio di 10^15 diversi possibili TCR  (Mason DA 1998). Ogni linfocita T mostra un solo tipo di TCR.


4. Cellule T helper 

Le cellule T helper sono programmate per riconoscere esclusivamente antigeni esposti dalle molecole MHC di seconda classe (II): questa classe di MHC viene espressa sulla membrana esterna di alcuni leucociti, principalmente le cellule dendritiche, le cellule B e i macrofagi (tutte assieme dette “cellule che presentano l’antigene”, antigen presenting cells: APC). Le molecole MHC II legano il TCR delle cellule T helper grazie al peptide CD4 (espresso unicamente dalle cellule T helper). L’antigene presentato dalle molecole MHC è un peptide lungo 13-17 amminoacidi (Rudensky, et al., 1991) (figura 2).


5. Linfociti T citotossici 

I TCR espressi dai linfociti T citotossici (CTL) possono legare solo antigeni esposti dalle molecole MHC di prima classe (I), che si trovano nella membrana esterna di qualunque cellula del nostro corpo. La glicoproteina CD8 è la molecola che rende i TCR espressi dalle CTL specifici per il MHC I. Mentre gli antigeni esposti dalle APC appartengono a patogeni raccolti sul campo di battaglia di passate infezioni, i peptidi esposti dal MHC I di una specifica cellula appartengono a patogeni che hanno fatto ingresso nella cellula stessa, e pertanto costituiscono la prova di un’infezione intracellulare ancora in atto (figura 3). Nel momento in cui un CTL riconosce un antigene che combacia con il proprio TCR, il CTL iduce l’apoptosi (morte programmata) della cellula che mostra l’antigene. Gli antigeni esposti dal MHC I sono peptidi che vanno dagli 8 ai 10 amminoacidi (Stern, et al., 1994).MHC I.JPG

Figure 3. Una cellula infetta espone un antigene virale sul proprio MHC I. Il TCR di un CTL si lega a questo peptide ed invia un segnale interno diretto al suo proprio nucleo, il quale risponde attivando l’apoptosi (attraverso il rilascio di granzimi, ad esempio) della cellula infetta (disegno di Paolo Maccallini).6. Il test immunologico universale

Nel corso del suo discorso, Mark Davis illustra alcuni concetti base sul sistema immunitario, prima di passare a introdurre i nuovi, entusiasmanti dati riguardo alla ME/CFS e alla Lyme cronica (o post-treatment Lyme disease syndrome: PTLDS). Contestualmente, però, dedica alcuni minuti alla descrizione di un complesso nuovo test che teoricamente renderebbe possibile estrapolare tutte le informazioni contenute nel repertorio di TCR presenti – in un dato momento – nel sangue di un essere umano. Un test del genere – che chiamerei “test immunologico universale” – avrebbe la capacità di determinare se un paziente presenta un’infezione in corso (e, nel caso, indicare il patogeno coinvolto) oppure una malattia autoimmune (anche qui specificando la natura dell’autoantigene, ossia il tessuto attaccato dal sistema immunitario). A quanto mi è dato di comprendere, il test richiede tre passaggi, che elenco nelle sezioni seguenti.

6.1. Primo step: sequenziamento del TCR

Come già spiegato nel paragrafo 3, quando un linfocita T incontra un peptide a cui è compatibile, comincia a proliferare; pertanto, nel sangue di un paziente con infezione in corso (o con reazione contro il proprio organismo, cioè con reazione autoimmune) è possibile trovare molteplici copie di linfociti T che esprimono il medesimo TCR: a differenza dei controlli sani, nei quali circa il 10% delle CD8 totali è rappresentato da copie di pochi diverse linfociti T (figura 4, prima linea), nei pazienti affetti da Lyme incipiente –  un esempio di infezione attiva – o sclerosi multipla (MS) –  un esempio di malattia autoimmune – abbiamo una massiccia clonazione di alcune linee di CTL (figura 5, seconda e terza riga, rispettivamente). Il primo step del test immunologico universale starà allora nell’identificazione dell’esatta sequenza di TCR espressa dai linfociti T presenti nel sangue, come si legge in Han A et al. 2014, dove troviamo descritto il sistema per sequenziare i geni delle catene α e β di un dato linfocita T. Tale approccio permette di costruire grafici come quello in figura 4 e quindi permette di determinare se il paziente presenti in atto un’attività anomala dei linfociti T oppure no. Qualora si riscontri un fenomeno di espansione clonale, è legittimo ipotizzare che stia avendo luogo o un’infezione o una condizione autoimmune di qualche sorta.

Clonal expansion CD8
Figure 4. Ogni cerchio rappresenta un paziente. Nella prima riga vediamo quattro controlli sani, che non presentano affatto espansione clonale delle cellule CD8 (come nel primo paziente da sinistra) oppure la presentano in maniera assai moderata (come indicato dalle porzioni in blu, bianco e grigio). Nella seconda riga troviamo invece quattro pazienti con malattia di Lyme attiva (fase incipiente) e possiamo ben notare come ciascuno di loro, nessuno escluso, presenti espansione clonale di solo tre diverse T cells (porzioni in rosso, blu e arancione). Nella terza riga, infine, abbiamo quattro pazienti affetti da MS, le cui cellule CD8 sono per maggior parte rappresentate da cloni di una selezione ristretta di T cells.
Fonte: slide proposte da Mark Davis durante il Community Symposium.

6.2. Secondo step: raggruppamento dei TCR 

Mark Davis e colleghi hanno realizzato un software capace di identificare i TCR che condividono il medesimo antigene, sia in un singolo individuo che trasversalmente a un gruppo. L’algoritmo è stato denominato GLIPH (grouping of lymphocyte interaction by paratope hotspots) ed ha dato prova di poter raggruppare – per fare un esempio – i recettori  dei linfociti T CD4 di 22 soggetti con infezione da M. tuberculosis latente in 16 gruppi distinti, ognuno dei quali comprende TCR provenienti da almeno tre individui (Glanville J et al. 2017). Cinque di questi gruppi sono riportati nella figura 5. L’idea sottostante è che TCR che appartengono allo stesso raggruppamento reagiscano allo stesso complesso peptide-MHC (pMHC).

Figure 5.  Cinque gruppi di TCR provenienti da 22 diversi pazienti affetti da turbercolosi latente, raggruppati grazie al GLIPH. La prima colonna da sinistra riporta l’identificativo dei TCR; la seconda l’identificativo dei pazienti. Le CDR per le catene β e α si trovano, rispettivamente, sulla terza e sulla quinta colonna. Tratto da Glanville J et al. 2017.

6.3. Terzo step: ricerca degli epitopi

Come abbiamo visto, questa nuova tecnologia consente di rilevare se sia in atto un’espansione clonale di linfociti T sequenziando i TCR dal sangue periferico. Consente inoltre di raggruppare i TCR presenti in un singolo paziente o condivisi da più pazienti. Il passaggio successivo è quello di identificare a quale/i tipo/i di antigene ognuno di questi raggruppamenti reagisca. Infatti, se potessimo identificare degli antigeni comuni in un gruppo di pazienti dai sintomi comparabili nei quali si riscontri un’espansione clonale in atto e simili TCR, saremmo messi in grado di comprendere se il loro sistema immunitario stia attaccando un patogeno (e di identificare il patogeno) o se stia piuttosto attaccando un tessuto ospite e, qualora fosse questo il caso, di identificare il tessuto. Come già detto, il numero di possibili combinazioni per il materiale genetico dei TCR è calcolato attorno ai 10^15, ma il numero dei possibili epitopi di cellule Th è circa 20^15, che corrisponde a più di 10^19. Ciò implica che i TCR debbano essere in una qualche misura cross-reattivi se vogliono essere in grado di riconoscere tutti i possibili peptidi esposti dai MHC (Mason DA 1998). Il grado di tale cross-reattività e il meccanismo attraverso il quale viene ottenuta sono stati spiegati con esattezza da Mark Davis e colleghi in un recente articolo (Birnbaum ME et al. 2014), che ci fornisce il terzo step del test immunologico universale. Lo scopo di questa fase consiste nel prendere un dato TCR e trovare il repertorio dei suoi specifici antigeni (giova ripetere che, appunto, ogni TCR reagisce a più antigeni). Per comprendere come ciò sia possibile, guardiamo a uno degli esperimenti descritti nell’articolo più sopra citato. I ricercatori si sono concentrati su due TCR ben precisi (chiamati Ob.1A12 e Ob.2F3), clonati da un paziente con MS e noti per essere capaci di riconoscere i pepetidi 85-99 (figura 6) della proteina basica della mielina (MBP) esposti dall’ HLA-DR15. Hanno poi preparato un insieme di cellule di lievito che esprimono molecole HLA-DR15, ognuna caratterizzata da un diverso peptide formato da 14 amminoacidi, con amminoacidi fissi esclusivamente alle posizioni 1 e 4, dove il peptide è ancorato al MHC (figura 6, sinistra). Quando alla coltura di cellule di lievito  che esprimono complessi pMHC vengono aggiunte copie di Ob.1A12, queste legano solo con alcune di quelle e, come è possibile vedere dalla parte destra della figura 6, per ciascuna posizione degli epitopi legati dal Ob.1A12 esiste un amminoacido con maggior tasso di probabilità: ad esempio, l’epitopo Ob.1A12 tipico presenta preferibilmente alanina (A) in posizione -4, istidina (H) in posizione -3, arginina (R) in posizione -2, e così via. Da notare che istidina (H) in posizione 2 e fenilanina (F) in posizione 3 sono amminoacidi obbligatori per un epitopo di  Ob.1A12.

Ob. 1A12
Figure 6. Sulla sinistra: il peptide 85-99 della proteina basica della mielina (myelin basic protein, MPB) è risaputo essere un epitopo per il TCR Ob.1A12. In posizione 1 e 4 possiede due residui che gli consentono di legare con la molecola MHC. In posizione -2, -1, 2, 3 3 5 troviamo invece i residui che legano con il TCR. La seconda riga rappresenta l’epitopo generico della libreria peptidica utilizzata per identificare il grado di cross-reattività tra tutti i possibili epitopi di Ob.1A12. A destra: la probabilità di ciascun amminoacido per ciascuna posizione è rappresentata da sfumature di viola. Come potete vedere, l’istidina (H) in posizione 2 e la fenilalanina (F) in posizione 3 sono amminoacidi obbligatori affinché un epitopo sia reattivo con Ob.1A12. Da (Birnbaum ME et al 2014).

La tabella sulla destra della figura 6 rappresenta, infatti, una tabella di sostituzione (substitution matrix) di dimensioni 14×20, uno strumento impiegato per scansionare il database dei peptidi in modo da trovare, tra tutti i peptidi conosciuti espressi da creature viventi, tutti i possibili epitopi specifici per Ob.1A12. Le matrici di sostituzione vengono solitamente utilizzate nel cosiddetto allineamento di peptidi (peptide alignment), una tecnica che punta all’identificazione di similitudini tra peptidi. Tali matrici sono basate su considerazioni di tipo evoluzionistico (Dayhoff, et al., 1978) o sullo studio delle regioni conservate delle proteine (Henikoff, et al., 1992). Ma la ricerca degli epitopi specifici di un dato TCR richiede (come abbiamo visto per Ob.1A12) una matrice di sostituzione costruita ad hoc per quel TCR: ogni TCR richiede la propria matrice di sostituzione, ottenuta incubando cellule T esprimenti quel TCR con una coltura di lieviti che espongono sui propri MHC una grande varietà di peptidi casuali, e analizzando poi i dati ricavati dall’esperimento. Quindi, un processo piuttosto complesso! Nel caso di Ob.1A12, questo processo ha portato a 2330 peptidi (incluso MBP), mentre la matrice di sostituzione specifica per Ob.2F3 ha trovato 4824 epitopi all’interno dell’intero database di peptidi. Questi peptidi includevano sia proteine non umane (batteriche, virali…) che peptidi umani. Per 33 di loro (26 non umani e 7 proteine umane), questo gruppo di ricercatori ha eseguito un test per verificare direttamente la previsione: 25/26 dei peptidi ambientali e 6/7 dei peptidi umani hanno indotto la proliferazione di cellule T che esprimono il TCR Ob.1A12 e/o il Ob.2F3, e questa è una prova della validità di questa analisi! Questi 33 peptidi sono riportati nella figura 7. Questo è l’ultimo passaggio del test immunitario universale, quello che dal TCR conduce agli epitopi. Come avete visto, un enorme insieme di diversi peptidi da diverse fonti reagisce con un singolo tipo di TCR; in altre parole, la cross-reattività è una proprietà intrinseca del TCR. Ciò significa anche che i test di trasformazione linfocitaria (LTT), ampiamente utilizzati in Europa per l’individuazione di infezioni da Borrelia burgdorferi e altri patogeni, comportano un rischio elevato di risultati falsi positivi e richiedono un processo di validazione sperimentale e teorica, attualmente mancante.

Crossreactive epitopes
Figura 7. Una serie di 33 peptidi che si suppongono essere epitopi specifici sia per Ob.1A12 che per Ob.2F3. Tratto da Birnbaum ME et al. 2014.

Siamo ora pronti ad apprezzare appieno i dati pilota che Mark Davis ha presentato al Symposium sull’espansione clonale delle cellule T CD8 nella ME/CFS e nella Lyme cronica.

7. “We have a hit!”

Mark Davis, insieme a Jacob Glanville e José Montoya, hanno sequenziato TCR dal sangue periferico di 50 pazienti ME/CFS e 49 controlli (primo passo del test immunitario universale, ricordate?), quindi li hanno raggruppati usando l’algoritmo GLIPH (secondo passo). Hanno trovato 28 cluster, ciascuno costituito da più di 2500 sequenze simili, e ogni cluster raccoglie sequenze multiple dallo stesso individuo e sequenze (che sono forse più importanti) da pazienti diversi (figura 8). Il cluster che ho cerchiato in rosso, ad esempio, è una raccolta di oltre 3000 TCR simili. La presenza di questi ampi cluster nei pazienti ME/CFS, rispetto ai controlli sani, rappresenta una prova indiretta di una risposta specifica delle cellule T a un trigger comune in questo gruppo di pazienti, che potrebbe essere un agente patogeno o un tessuto del corpo (o tutti e due).

Clustered TCR
Figura 8. Nella ME/CFS le sequenze di TCR ricavati da 50 pazienti formano 28 raggruppamenti che presentano più di 2500 sequenze simili – cosa che assolutamente non avviene nei controlli sani. Questo fa pensare ad una qualche risposta immunitaria ad un patogeno o ad un tessuto umano (o entrambi). Fonte: slide proposta da Mark Davis durante il Community Symposium.

Tra questi 50 pazienti ME/CFS, Davis e colleghi hanno selezionato 6 pazienti con geni HLA simili (figura 9, sinistra), 5 femmine tra loro, e hanno studiato i loro TCR più in profondità. Nella metà destra della figura 9, è possibile vedere per ciascun paziente il grado di espansione clonale delle CTL. Ricordate che nei controlli sani solo circa il 10% dei CTL è composto da cloni di alcune cellule (figura 4, prima riga), mentre qui vediamo che circa il 50% di tutti i CTL è composto da cloni. Quindi, una “marcata espansione clonale” delle cellule T CD8, come ha detto Davis.

ME subjects CD8
Figura 9. A sinistra: sono stati selezionati 6 pazienti ME/CFS con HLA simili. Sulla prima colonna da sinistra sono stati riportati gli identificativi dei pazienti; la seconda colonna ci informa sull’età di ciascuno; la terza sul genere; la quarta avvisa di eventuali esposizione a citomegalovirus; la quinta riguarda i geni del MHC I. A destra: l’analisi dell’espansione clonale delle cellule T CD8 per ognuno dei pazienti. L’espansione clonale è marcata (circa al 50%), se comparata a quella dei controlli sani (circa al 10%).

Le sequenze delle catene α e β dei TCR di tre dei sei pazienti (pazienti L4-02, L4-10 e L3-20) sono riportate in figura 10 (ho verificato che effettivamente si tratta di catene α e β di TCR umani, inserendole in BLAST).

TCR epitope
Figura 10. Catene β (prima colonna) e rispettive catene α (quinta colonna) provenienti da tre pazienti ME/CFSchains  (L4-02, L4-10, and L3-20, ultima colonna).

Quindi, abbiamo visto finora i primi due passaggi del test immunitario universale. E il terzo passo? Nel suo discorso, Mark Davis non ha presentato alcun particolare epitopo, ha solo mostrato una diapositiva con quella che probabilmente è la selezione degli epitopi dalla libreria discussa nel paragrafo 6.3 da parte di uno dei TCR riportati in figura 10. Questa selezione è riportato in figura 11, ma da quella foto non è possibile raccogliere alcuna informazione sull’identità di questi epitopi. Come probabilmente ricorderete dal paragrafo 6.3, l’analisi dei peptidi selezionati da un TCR nella libreria di peptidi  consente l’identificazione di una matrice di sostituzione che può essere utilizzata per selezionare tutti i possibili epitopi di quel TCR specifico, dal database delle proteine. Quest’ultima fase cruciale deve essere ancora eseguita, o è già stata eseguita, ma Davis non ha comunicato i risultati preliminari durante il suo discorso. Recentemente sono state messe a disposizione nuove risorse dalla Open Medicine Foundation, affinché questa ricerca promettente possa essere ulteriormente perseguita (R). Lo scopo qui, come già detto, è di trovare l’antigene che innesca questa risposta delle cellule T. Come ha detto Mark Davis, potrebbe essere un antigene di un agente patogeno specifico (forse un patogeno comune che va e viene) che suscita una risposta immunitaria anomala che finisce per colpire alcuni tessuti ospiti (microglia, per esempio), portando così attivazione immunitaria che è stata recentemente segnalata da Mark Davis stesso e altri in ME/CFS (Montoya JG et al. 2017). L’idea di un patogeno comune che innesca una risposta immunitaria patologica non è nuova in medicina, e la febbre reumatica (RF) è un esempio di una tale malattia: la RF è una malattia autoimmune che attacca il cuore, il cervello e le articolazioni ed è generalmente innescata da uno streptococco che infetta la gola (Marijon E et al. 2012). L’altra possibilità è, naturalmente, quella di un’infezione in corso di qualche tipo, che deve ancora essere rilevata. Come detto (vedi par. 6.1), l’espansione clonale delle cellule T CD8 è presente sia nelle infezioni acute (come la malattia di Lyme) che nelle malattie autoimmuni (come la SM) (figura 4), quindi dobbiamo aspettare l’identificazione dell’antigene se vogliamo capire se l’attività del CTL è contro un agente patogeno e/o contro un tessuto del nostro corpo.

peptide library
Figura 11. Nella figura possiamo osservare la selezione, che avviene in più momenti, di una serie di peptidi da parte di un particolare TCR proveniente da un paziente ME/CFS. La selezione ha luogo tra una enorme quantità di peptidi esposti dall’ HLA-A2 (MHC I) espresso da cellule di lievito. Ad ogni passaggio il numero di possibili peptidi si riduce.

8. La Lyme cronica esiste

È stato probabilmente trascurato il fatto che nel suo discorso, Mark Davis ha riportato anche dati molto interessanti sulla sindrome della malattia di Lyme post-trattamento (PTLDS, nota anche come malattia di Lyme cronica). In particolare, ha trovato un’espansione clonale marcata nelle cellule T CD8 di 4 pazienti PTLDS (circa il 40% dei CTL totali) come riportato nella figura 12: si consideri che in questo caso le fette blu rappresentano cellule T uniche, mentre tutte le altre fette rappresentano cloni! Tutto ciò che è stato detto sull’espansione clonale CD8 nella ME/CFS si applica anche in questo caso: potrebbe essere la prova di un’infezione in corso – forse la stessa B. burgdorferi, come suggerito da diversi modelli animali (Embers ME et al. 2017), (Embers ME et al. 2012), (Hodzic E et al. 2008), (Yrjänäinen H et al. 2010) –  o una coinfezione (un virus?). Oppure potrebbe essere l’espressione di una reazione autoimmune innescata dalla infezione iniziale. Questo deve ancora essere scoperto, eseguendo il test immunitario universale completo, ma ciò che è già chiaro dalla figura 12 è che la PTLDS è una condizione reale, con qualcosa di veramente anomalo nella risposta immunitaria: la Lyme cronica esiste.

Figura 12. Espansione clonale di cellule T CD8 in quattro pazienti affetti da Lyme cronica. L’espansione clonale, che indica l’attività delle cellule T contro un patogeno o un tessuto ospite, è assai marcata.

9. Conclusioni

Mark Davis e altri ricercatori hanno sviluppato un test complesso che è in grado di sequenziare i TCR dai pazienti, raggrupparli in gruppi di TCR che reagiscono agli stessi antigeni e scoprire gli antigeni che hanno attivato quella particolare risposta delle cellule T. Questo test è una sorta di test immunitario universale che è teoricamente in grado di riconoscere se una persona (o un gruppo di pazienti) presenta una risposta immunitaria contro un agente patogeno o contro uno dei loro stessi tessuti (o entrambe le cose). Questo approccio ha già fornito dati pilota su una attivazione anomala delle cellule T CD8 nei pazienti ME/CFS e nei pazienti PTLDS e, si spera, identificherà il trigger di questa risposta immunitaria nel prossimo futuro. Se la ME/CFS è causata da un’infezione attiva, da una malattia autoimmune o da entrambe le cose, il test immunologico universale potrebbe essere in grado di dircelo. Questa nuova tecnologia è per l’immunologia, ciò che il sequenziamento dell’intero genoma è per la genetica, o la metabolomica è per le malattie molecolari: non cerca un particolare agente patogeno o una particolare malattia autoimmune. No, cerca tutte le possibili infezioni e disturbi immunitari, anche quelli che devono ancora essere scoperti.

Il teorema di Eulero

Il teorema di Eulero

Quello che segue è un capitolo introduttivo del manuale di Meccanica Razionale che scrissi molti anni fa, in un periodo di miglioramento, allo scopo di recuperare nozioni apprese prima della malattia e – allo stesso tempo – di riabilitare la mente annichilita dalla misteriosa infermità che è oggetto di questo blog.

Molti anni dopo, i miei tentativi disperati di analisi della patologia che mi ha tolto dal mondo, mi portarono a studiare – tra le altre cose – la meravigliosa struttura delle proteine e le loro interazioni con gli anticorpi. E ho ritrovato la meccanica del corpo rigido dove non osavo cercarla: le proteine sono spesso approssimabili come corpi rigidi e lo studio degli autoepitopi si avvale di un oggetto matematico centrale nella meccanica razionale, l’ellissoide di inerzia. Questa è comunque un’altra storia, su cui tornerò; una storia d’amore perduto e ritrovato.

Quello che segue è il capitolo in cui definisco il corpo rigido e lo spostamento rigido. Segnalo in particolare la dimostrazione del teorema di Eulero, che sancisce la natura rotatoria di ogni spostamento sferico. Il mio volume completo è disponibile qui.

Capitolo 2. Spostamenti-001

Capitolo 2. Spostamenti-002

Capitolo 2. Spostamenti-003

Capitolo 2. Spostamenti-004

Capitolo 2. Spostamenti-005

Capitolo 2. Spostamenti-006

Capitolo 2. Spostamenti-007

Capitolo 2. Spostamenti-008

Capitolo 2. Spostamenti-009

Capitolo 2. Spostamenti-010

Capitolo 2. Spostamenti-011

Capitolo 2. Spostamenti-012

Capitolo 2. Spostamenti-013

Capitolo 2. Spostamenti-014

Capitolo 2. Spostamenti-015

Capitolo 2. Spostamenti-016

Capitolo 2. Spostamenti-017

Capitolo 2. Spostamenti-018

Why we can’t use LTTs, yet

A line of T cells (called Ob.2F3) expressing the same T cell receptor (TCR) from an MS patient was studied in 2014 and it was found to proliferate when incubated with 4824 different peptides. Thirty-three of them were further studied (see figure) and found to belong to both Homo sapiens and several different, unrelated microbes (Birnbaum ME et al. 2014). The taking home message here is that T cells are not specific to a single pathogen, they are highly cross-reactive, as it was already pointed out in this pivotal study: (Mason DA 1998). And this means that we can’t use lymphocyte transformation tests (LTTs) the way we do now. 

I feel really frustrated when patients send me their LTTs and ask me to comment the results. I have to say that they have vasted their money and that these results are useless. I do hope that my blog can make a difference and stop this unfair commerce at the expenses of desperate folks.

Crossreactive epitopes
Figure. A set of 33 peptides (both human and environmental) predicted to be specific epitopes for both Ob.1A12 and Ob.2F3. From (Birnbaum ME et al. 2014).






Integrals for Statistics

In the following four tables, I present a collection of integrals and integral functions that are of common use in Statistics. The values for the Beta function have been calculated by a code of mine that uses – as a recursive function – another code that I wrote for the calculation of the Gamma Function (see the related blog post). Below these tables, you find a representation of the Beta function, plotted by my code. I have demonstrated each one of these results and these calculations will be included in an upcoming book.

Immagine.pngImmagine 2Immagine 3

ME/CFS, a mathematical model

ME/CFS, a mathematical model

Robert Phair and the trap

Many of my readers are probably aware of the attempts that are currently being made to mathematically simulate the energy metabolism of ME/CFS patients, integrating metabolic data with genetic data. In particular, dr. Robert Phair has developed a mathematical model of the main metabolic pathways involved in energy conversion, from energy stored in the chemical bonds of big molecules like glucose, fatty acids, and amino acids, to energy stored in adenosine triphosphate (ATP), ready to use. Phair, who is an engineer, determined the differential equations that rule this huge amount of chemical reactions and adapted them to the genetic profile found in ME/CFS patients. Genetic data are the result of the analysis of the exomes of more than 40 patients (I am among them). He found, in particular, that the enzyme ATP synthase (one step of the mitochondrial electron transport chain, called complex V) presents a damaging variant that is found in less than 30% of healthy controls while being detected in more than 60% of ME/CFS patients. He found other enzymes (I don’t know their exact number) of energy metabolism meeting these same criteria. Setting a reduced activity for these enzymes in his mathematical model, he found that, after particularly stressing events (such as infections), ME/CFS patients metabolism can fall into a low-energy state without being able to escape from it (this theory has been called the “metabolic trap hypothesis”) (R, R, R). Phair’s hypothesis is currently being experimentally tested and has not been published yet, but some years ago two physicists published an interesting mathematical model of energy metabolism during and after exercise, in ME/CFS patients compared to healthy controls (Lengert N. et Drossel B. 2015). In what follows I will describe this model and its predictions and we will also see how these differential equations look like. 

Metabolic pathways that have been analyzed

The model by Lengert and Drossel extends two previously published systems of differential equations that describe the behaviour of glycolysis, Krebs cycle (hugely simplified as a single reaction!), mitochondrial electron transport chain (described in detail), creatine kinase system, and of conversion of adenosine diphosphate (ADP) in ATP, in skeletal muscles (Korzeniewski B. et Zoladz JA. 2001), (Korzeniewski B. et Liguzinski P. 2004). They added equations for lactate accumulation and efflux out of the cell, for de novo synthesis of inosine monophosphate (IMP) during recovery, for degradation of adenosine monophosphate (AMP) into IMP, for degradation of IMP to inosine and hypoxanthine. All the pathways involved are collected in figure 1. These reactions are described by 15 differential equations and the solution is a set of 15 functions of time that represent the concentration of the main metabolites involved (such as lactate, pyruvate, ATP etc). Let’s now take a closer look at one of these equations and at the general structure of the whole system of equations.

Figure 1. This is a schematic representation of the metabolic pathways described by the mathematical model developed by Lengert and Drossel. In detail: cytosolic synthesis and degradation of ADP, AMP and IMP (left), protein kinase pathway and glycolysis (centre), electron transport chain and TCA cycle (right). From Lengert N. et Drossel B. 2015.

lactate dehydrogenase.PNG
Figure 2. Lactate dehydrogenase is the enzyme involved in the catalysis of the conversion of lactate to pyruvate. This reaction goes in both directions.

Differential equations for chemical reactions

Let’s consider the equation used by the Authors for the reaction catalyzed by lactate dehydrogenase (the transformation of pyruvate into lactate, figure 2) where they also took into account the efflux of lactate from the cytosol. The differential equation is the following one:


where the three parameters are experimentally determined and their values are


The first one describes the activity of the enzyme lactate dehydrogenase: the higher this parameter is, the more active the enzyme is. The second one describes the inverse reaction (from lactate to pyruvate). The third one is a measure of how much lactate the cell is able to transport outside its membrane. You have probably recognized that the equation of lactate is a first-order ordinary differential one. We say “first-order” because in the equation there is only the first-order derivative of the function that we have to determine (lactate, in this case); “ordinary” refers to the fact that lactate is a function only of one variable (time, in this case). You can easily realize that an equation like that can be written as follows:

equation bis.PNG

Suppose now that we had other two differential equations of this type, one for pyruvate  and one for protons (the other two functions of time that are present in the equation):


We would have a system of three ordinary differential equations like this oneSystem.PNG

The initial values of the functions that we have to determine are collected in the last row: these are the values that they have at the beginning of the simulation (t=0). In this case, these values are the concentrations of lactate, pyruvate and protons in the cytosol, at rest. The three functions of time are called the solution to the system. This kind of system of equations is an example of a Cauchy’s problem, and we know from mathematical theory that not only it has a solution, but that this solution is unique. Moreover, whereas this solution can’t be always easily found with rigorous methods, it is quite easy to solve the problem with computational methods, like the Runge-Kutta method or the Heun’s method. All that said, the system of ordinary differential equations proposed by Lengert and Drossel for energy metabolism is just like this one, with the exception that it comprises 15 equations instead of three. So, the main difficulty in this kind of simulation is not the computational aspect but the determination of the parameters (like the enzymatic ones) and of the initial values, that have to be collected from the medical literature or have to be determined experimentally, if not already available. The other problem is how to design the equations: there are several ways to build a model for a chemical reaction or for any other biological process.

The mathematical model of ME/CFS

How do we adapt to ME/CFS patients a model of energy metabolism that has been set with parameters taken from experiments performed on healthy subjects? This is a very good question, and we have seen that Robert Phair had to use genetic data from ME/CFS patients on key enzymes of energy metabolism, in order to set his model. But this data was not available when Lengert and Drossel designed their equations. So what? They looked for studies about the capacity of oxidative phosphorylation in ME/CFS patients in comparison with healthy subjects, and they found that it had been measured with different experimental settings by various groups and that the common denominator was a reduction ranging from about 35% (Myhill S et al. 2009), (Booth, N et al 2012), (Argov Z. et al. 1997), (Lane RJ. et al. 1998) to about 20% (McCully KK. et al. 1996), (McCully KK. et al. 1999). So the idea of the Authors was to multiply the enzymatic parameter of each reaction belonging to the oxidative phosphorylation by a number ranging from 0.6 (severe ME/CFS) to 1.0 (healthy person). In particular, they used a value of 0.7 for ME/CFS, in their in silico experiments.

Predictions of the mathematical model

The mathematical model was used to perform in silico exercise tests with various length and intensities. What they found was that the time of recovery in the average ME/CFS patient was always greater if compared to a healthy person. The time of recovery is defined as the time that a cell needs to replenish its content of ATP (97% of the level in resting state) after exertion. In Figure 3 you see the results of the simulation for a very short (30 seconds) and very intense exercise. As you can see, in the case of a healthy cell (on the left) the recovery time is of about 600 minutes (10 hours) whereas a cell from a person with ME/CFS (on the right) requires more than 1500 minutes (25 hours) to recover.

half minute 1.png
Figure 3. Results of the simulation for an exercise with a duration of 30 seconds and a high intensity (initial ATP consumption 300 times the resting value). On the left, the case of a healthy skeletal muscle cell, on the right the case of a cell from a person with ME/CFS whose mitochondrial reactions have a velocity reduced to 70% of the velocity of healthy control. The plot has been obtained by using the online version of the software, available here.

Another interesting result of the simulation is an increase in AMP in patients vs control (figure 3, orange line). This is due to the compensatory use of the two metabolic pathways in figure 4: the reaction catalyzed by adenylate kinase, in which two molecules of ADP are used to produce a molecule of ATP and a molecule of AMP; and the reaction catalyzed by AMP deaminase, that degrades AMP to IMP (that is then converted to inosine and hypoxanthine). These two reactions are used by ME/CFS patients more than in healthy control, in order to increase the production of ATP outside mitochondria.

adenylate kinase and amp deaminase.PNG
Figure 4. The metabolic pathway on the left is used by ME/CFS patients more than in control to increase the production of ATP outside mitochondria, according to the present model. The pathway on the right then degrades AMP to IMP.

If we give a closer look at the concentrations of AMP and IMP in the 4 hours following the exertion (figure 5), we actually see an increased production of IMP (green line) and AMP (orange line) in skeletal muscles of ME/CFS patients (on the right) vs controls (left).

half minute 3.png
Figure 5. The same as figure 3, but magnified in order to have a closer look at concentrations during the 4 hours following the exertion. The healthy cell is on the left, whereas the cells from a person with ME/CFS is on the right.

A further compensatory pathway used by patients (according to this model) is the production of ATP from ADP by the enzyme creatine kinase (figure 6). This is another way that we have on our cells to produce ATP in the cytosol without the help of mitochondria. In this model of ME/CFS, there is an increase in the use of this pathway, which leads to a decrease in the cellular concentration of phosphocreatine and an increase in the cellular concentration of creatine (figure 7).

creatine kinase
Figure 6. The reaction catalyzed by creatine kinase: a molecule of ADP has converted to ATP thanks to the phosphate group carried by phosphocreatine.

half minute 4.png
Figure 7. The concentration of phosphocreatine in the cytosol of skeletal muscle cells is lower in ME/CFS (right) versus control (left) during and after exercise. This is due to the greater use of this molecule to produce ATP anaerobically in ME/CFS metabolism vs control. Parameters for this simulation are the same as described in figure 3.

Comparison with available metabolic data

I am curious to see if data from the various metabolomic studies done after the publication of the model by Lengert and Drossel are in agreement with it. I will discuss this topic in another article because I still have to study this aspect. I would just point out that if we assumed true the high rate of IMP degradation proposed in this model, we would probably find a high level of hypoxanthine in the blood of patients, compared to controls, whereas this metabolite is decreased in patients, according to one study (Armstrong CW et al. 2015).

Comparison with the model by Phair

The metabolic model by Robert Phair will probably give a more accurate simulation of the energy metabolism of patients if compared with the system of ordinary differential equations that we have discussed in this article. And there are two reasons for that. The first one is that Phair has included equations also for fatty acid beta-oxidation, pentose phosphate pathway, and NAD synthesis from vitamin niacin. The other one is that, whereas the two German physicists reduced the velocities of all the enzymatic reactions that happen in mitochondria, Phair has genetic data for every enzyme involved in these reactions for a group of ME/CFS patients and thus he can determine and set the actual degree of activity for each enzyme. But there is a further level required in order to bring the mathematical simulation closer to the reality: gene expression. We know, for instance, that in ME/CFS patients there is a higher than normal expression of aconitase (an enzyme belonging to the TCA cycle) and of ATP synthase (Ciregia F et al 2016) and this should be taken into account in a simulation of ME/CFS patients energy metabolism. Note that ATP synthase is exactly the enzyme that Phair has found to be genetically damaged in patients, and this makes perfect sense: if an enzyme has reduced activity, its reaction can be speeded up by expressing more copies of the enzyme itself.

One could expect that, in a near future, genetic data and gene expression data from each of us will be used to set mathematical models for metabolic pathways, in order to build a personalized model of metabolism that might be used to define, study and correct human diseases in a personalized fashion. But we would need a writer of sci-fiction in order to tell this chapter of the future of medicine.


Bivariate normal distribution

Bivariate normal distribution

In a previous blog post, I have discussed the case of normally distributed random variables. But what happens if we have a vector whose components are random variables with this kind of distribution? In particular, what is the expression of the joint probability density? I don’t know how to derive this formula from the marginal distributions, but I have been able to demonstrate, after several hours of work, the following proposition.

Proposition 1. The joint probability density in the case of a random vector whose two components follow a normal distribution is:

probability density.png

where σ_1 and σ_2 are the standard deviations, μ_1 and μ_2 are the expectations of the two random variables and ρ is the correlation coefficient between the two variables.

Proof. The marginal distribution of x_1 is given by

marginal density.png

By operating the substitutions

substitutions.pngwe obtain

marginal density 2.png

This demonstrates that x_1 is in fact a normally distributed random variable and that its expectation and standard deviation are σ_1 and μ_1, respectively. The same applies to x_2. In order to demonstrate that ρ is the correlation coefficient between x_1 and x_2 we have to calculate the covariance of f_X, which is given by E[X_1,X_2]-E[X_1]E[X_2]. The first addend is given by:

covariance 1.png

Let’s consider the inner integral and the following substitutions

substitutions 2.png

We have:covariance 2.png

The two integrals inside the root square are not the easiest ones, and I will not discuss them here. Just consider that they can be solved using the well known results:

exponential integral 4

We have exponential integral 3.png

 If we put these two integrals inside I, we obtain

covariance 3.png

The two inner integrals are two difficult ones, again. The first one can be solved using the same exponential integral used above:

covariance 4

For the other inner integral we have

covariance 5.png

The first and the second integrals can be solved using the table previously reported, while the third one is zero. At the end of all this work, we have:

covariance 6.png

And this proves that ρ is the correlation coefficient between x_1 and x_2. This is a complete demonstration of the expression of the joint probability density of a bivariate normal distribution ■

This formula can be extended to the case of m components.

Proposition 2. The joint probability density in the case of a random vector whose m components follow a normal distribution is:

probability density 2.png

where C is the covariance matrix, given by

covariance matrix.png

Proof. This formula is easy to prove in the case of m=2, using proposition 1. I have no idea about how it could be proved when m>2 ■

I have written a code in Fortran (download) that asks you for expectations, standard deviations and correlation coefficient and plots the joint density of the bivariate normal distribution. I have written a code in Octave that does the same job (download). In the picture below you find the joint density for two standard normally distributed random variables with ρ=0 (left), ρ=0.7 (right, upper half) and ρ=-0.7 (right, lower half). These pictures have been plotted by the code in Fortran. The other image is a comparison between the two programs, for another bivariate normal distribution.






Three new possible autoepitopes in ME/CFS

Paolo Maccallini


I have performed a set of analysis of experimental data already published about autoimmunity to muscarinic receptors in ME/CFS. My predictions are that extracellular loop 2 and 3, and also transmembrane helix 5 of both muscarinic cholinergic receptors 4 and 3, are main autoantigens in a subset of ME/CFS patients. Moreover, I have found that autoimmunity to M4 and M3 ChR is independent from autoimmunity to beta 2 adrenergic receptor, also reported in ME/CFS patients.  


Myalgic encephalomyelitis/chronic fatigue syndrome (ME/CFS) is a debilitating disease characterized by cognitive deficits, fatigue, orthostatic intolerance with symptoms exacerbated after exertion (post-exertional malaise, PEM) (IOM, 2015). This disease has no known cause but several abnormalities have been observed in energy metabolism (Tomas C. and Newton J. 2018), immune system, gut flora (Blomberg J. et al. 2018), brain (Zeineh MM. et al. 2014). A possible role for autoantibodies in the pathogenesis of the disease has been suggested by the finding of reactivity of patient sera to two nuclear antigens (Nishikai, et al., 1997), (Nishikai, et al., 2001), to cardiolipin (Hokama, et al., 2009), to HSP60 (Elfaitouri A. et al. 2013), and to muscarinic cholinergic (M ChR) and beta adrenergic receptors (ß AdR) (Tanaka S et al. 2003), (Loebel M et al. 2016); reactivity that was significantly elevated when compared to healthy contols. Reactivity to adrenergic and muscarinic Ch receptors has been confirmed by two independent groups, but these results have not been published yet (R). A role for autoantibodies in at least a subgroup of patients has also been suggested by a response to rituximab, a CD20 B cells depleting agent (Fluge Ø. et al. 2011), (Fluge Ø. et al. 20115), and to immunoadsorption (Scheibenbogen C. et al 2018). Sera response to muscarinic cholinergic receptors is confirmed in two studies but both of them used an immune assay with proteins coated on a plate. This kind of test does not allow to identify the exact autoepitope on the receptor and – even more importantly – it is subjected to false positive results because it exposes to sera surfaces of receptors that are hidden when they are in their physiological position (Ramanathan S et al 2016). Nevertheless, the amount of data provided in the study by Loebel et al. where reactivity of sera to 5 subtypes of muscarinic cholinergic receptors have been measured simultaneously, has – in our opinion – the potential to unveil the exact autoepitope(s). Thus, I performed a bioinformatical analysis on experimental data from this study in order to extract hidden information. I used a software for the in silico study of B cell epitope cross-reactivity (Maccallini P. et al. 2018) and a software for amino acid protrusion index calculation (Ponomarenko J. et al., 2008).  Our prediction is that patients sera mainly react to three epitopes that belong to the second and third extracellular loop of M3 and M4 ChR, but also to a hidden epitope of the same two receptors, leading to possible false positive results of this test. I have also found that the reactivity to beta 2 adrenergic receptor (ß2 AdR) found in the study by Loebel et al. is not due to the same antibody that reacts to muscarinic cholinergic receptors.


Search for cross-reactive epitopes. Cross-reactivity between muscarinic cholinergic receptors M4 and M3, and between M4 and M1 has been studied in silico using EPITOPE, a software already described (Maccallini P. et al. 2018). Briefly, EPITOPE searches for cross-reactive epitopes shared between two proteins (let’s say protein A and protein B) by comparing each possible 7-mer peptide of A with each possible 7-mer peptide of B. The comparison is made using the algorithm by Needleman and Wunsch (Needleman SB. and Wunsch CD. 1970)  with a gap model a + b·x, where a is the opening gap penalty, b is the extending one, and x is the extension of the gap. A penalty for gaps at the end of the alignment was also assumed. The choice for gap penalties and substitution matrix were done according to the theory already developed for peptide alignments (Altschul SF. 1991), (Karlin S. and Altschul SF. 1990). Available experimental data on cross-reactivity between γ enolase and α enolase (McAleese SM. et al. 1988)  have been used for EPITOPE calibration: a score >60 was considered the cut-off for cross-reactivity, a score below 50 indicates non-cross-reactive epitopes; a score between 50 and 60 defines a borderline result. A simpler version of EPITOPE has been used for single local alignments. The main program used for M4-M3 comparison, its subroutine NeWalign and the substitution matrix employed are available for download. Primary structures used in this work have been downloaded from UniProt and are the following ones: M1 ChR (P11229), M3 ChR (P20309), M4 ChR (P08173), B2 AdR (P07550).

Surface exposure. In order to select only those 7-mer peptides that are on the surface of proteins, I have considered their mean protrusion indexes. A protrusion index of at least 0.5 has been considered the cut-off for surface exposure. Protrusion indexes of single amino acids have been calculated with ElliPro. A protrusion index of 0.5 means that the amino acid is outside the ellipsoid of inertia which includes 50% of the centers of mass of all the amino acids of the protein (Ponomarenko J. et al., 2008). For M4 ChR I have used the crystal structure 5DSG (Thal DM. et al. 2016). The 3D structure of human M3 ChR has not been experimentally determined yet, so I have used a theoretical model built using murine M3 ChR (PDB ID: 4DAJA) as a template, provided by ModBase.

M ChR plot
Figure 1. The position of the first amino acid of each possible 7-mer peptide of M4 ChR is reported on the abscissa, the score for the comparison of each of these peptides with M1 ChR (blue line) and M3 ChR (orange line) is reported on the ordinate. N terminus, extracellular loop 1, 2 and 3 are also indicated. Scores above the yellow line indicate cross-reactivity, scores below the blue line indicate a lack of cross-reactivity.

Selection criteria. Our purpose is to predict whih epitopes of M3 and M4 ChRs sera from ME/CFS patients react to. So, we search for M4 ChRs 7-mer peptides that are cross-reactive to M3 ChR, but non-cross-reactive to M1 ChR. Moreover, they must present surface exposure both on M4 and on M3 ChR (otherwise antibodies can’t reach them). So, selection criteria for M4 ChR epitopes are as follows:

  1. they must be cross-reactive to M3 ChR;
  2. they must be non-cross reactive to M1 ChR or borderline;
  3. they must present a mean protrusion index ≥0.5;
  4. M3 ChR peptides they cross-react to have to present a mean protrusion index ≥0.5.

We will refer to strict criteria when we assume only non-cross-reactivity in 2, while weak selection criteria are fulfilled when M4 ChR epitopes have borderline reactivity to M3 Chr peptides.

M4 vs M1, M3
Figure 2. Distribution of the scores from the comparison of M4 ChR with M1 ChR (left) and with M3 ChR (right). M3 ChR presents a slightly higher mean score.


The search for 7-mer peptides of M4 ChR that are cross-reactive to M3 ChR found 108 sequences. We then studied cross-reactivity to M1 ChR for each of these peptides and we found that 11 of them are non-cross-reactive and that other 9  peptides have borderline reactivity. None of these 20 peptides presented a cross-reactivity to B2 AdR (Table 1S, column 1). Scores between peptides of M4 ChR and the other two muscarinic cholinergic receptors are plotted in Figure 1. The distribution of scores from the comparison of M3 ChR with M1 ChR and with M3 ChR are reported in Figure 2. For the M4 ChR 20 epitopes mentioned above, we calculated the mean protrusion indexes and we did the same calculation for their cross-reactive peptides on M3 ChR. We also indicated their position with respect to the plasma membrane. All these data are collected in Table 2S. Once we apply selection criteria on these 20 peptides, we obtain 9 epitopes (Table 1). Of these selected epitopes, one belongs to a transmembrane helix: peptide 186-192 of M4 ChR, which cross-reacts to peptide 231-237 of M3 ChR. Peptide 418-431 of M4 ChR is partially immersed in the plasmatic membrane, even though its cross-reactive peptide of M3 ChR is entirely exposed to the extracellular space, and the same applies to the other two epitopes found (figure 1). Peptide 175-181 of M4 ChR cross-reacts to peptide 211-217 of M3 ChR; peptide 186-192 of M4 Chr cross-reacts to peptide 222-228 of M3 ChR; peptide 418-431 of M4 Chr cross-reacts to peptide 513-522 of M3 ChR. Sequences that fulfill selection criteria and their respective inverted sequences are collected in  Table 2.

Table 1
Table 1. This is the collection of M4 Chr 7-mer peptides that are cross-reactive to M3 ChR; are not cross-reactive or borderline with M1 ChR; have a mean protrusion index higher than 0.5; are cross-reactive with epitopes of M3 ChR with a protrusion index higher than 0.5.


B cells autoimmunity to muscarinic cholinergic receptors in ME/CFS has been reported in two studies (Tanaka S et al. 2003), (Loebel M et al. 2016) and this finding has been recently confirmed by two other independent groups who have not published yet (R). The two studies mentioned used full-length proteins coated on a plate in order to perform the immune assay. With this kind of technique we may have both false positives (due to the fact that sera react with peptides that are not in the extracellular domain) and false negatives (due to protein denaturation, which leads to the formation of epitopes that would not be present if the protein were correctly folded) as it has been reported in the case of anti-MOG antibodies (Ramanathan S et al 2016). A way to solve the possible inaccuracy of these data would thus be to measure sera reactivity with a cell-based assay (CBA) which is a test where receptors are expressed by eukaryotic cells and thus they are held in their physiological position.

Figure 1. Peptides of table 1 that belong to the extracellular domain of M3 and M4 ChR are here highlighted directly on the 3D structures of their respective receptors.

Nevertheless, we can still try to extract hidden information from experimental data and predict the position of the epitope(s) ME/CFS patients sera react to. Knowing that sera from patients react to M4, M3 ChRs and that there is a low correlation between reactivity to M4 ChR and reactivity to M1 ChR (Loebel M et al. 2016), we selected 7-mer peptides of M4 ChR that cross-react (in silico) to M3 ChR but not to M1 ChR (Table 2S). We then selected, among them, only those peptides that have surface exposure on their respective proteins (Table 1). The result is that patient sera react to extracellular loops 2 and 3 of both M3 and M4 ChRs (Figure 1), but also to a hidden antigen, a peptide of transmembrane helix 5 of both M3 and M4 ChR.

Our results are of interest because extracellular loops 2 and 3 of M3 ChR are known autoepitopes in Sjögren’s syndrome (Ss) (Deng C. et al. 2915). Moreover, sera from patients with orthostatic hypotension (OH) react to extracellular loop 2 of M3 ChR, where they show an agonistic effect, thus acting as vasodilators (Li H. et al. 2012). OH, a form of orthostatic intolerance has been reported in ME/CFS patients (Bou-Holaigah et al. 1995) while fatigue similar to post-exertional malaise has been described in Ss (Segal B. et al. 2008). A pathogenic role of these antibodies in fatigue for both ME/CFS and Sjögren syndrome could perhaps be due to their vasodilatory effect.

Our analysis unveiled reactivity to a hidden autoepitope, which belongs to transmembrane helix 5 of M3 and of M4 ChR. This epitope is buried inside the plasma membrane when these two receptors are in their physiological position, so this reactivity can’t contribute to the pathogenesis of ME/CFS.

None of the 7-mer peptides of M4 ChR that cross-react to M3 ChR and at the same time don’t cross-react to M1 ChR presents in silico reactivity to B2 AdR. This means that in those patients whose sera present reactivity to both M4-M3 ChR and B2 AdR, there are two distinct autoantibodies. This prediction of our model is consistent with the low correlation found by Loebel and colleagues between anti-M4 ChR and anti-B2 AdR antibodies (Loebel M et al. 2016).

Most B cells epitopes on non-denaturated proteins (i.e. proteins that conserve their tertiary structure) are believed to be conformational (Morris, 2007), so a significant limitation of this study is due to the fact that our analysis considers only linear epitopes. Nevertheless, the main limitation of this study remains by far my encephalopathy.


This analysis of previously published data suggests a role for the second and the third extracellular loop of M4 and M3 ChR as autoantigens in ME/CFS. It also predicts the presence of a hidden autoantigen and thus a risk of false-positive results with standard ELISA.  The eight peptides found by this analysis and their inverse sequences (Table 2) should be employed as query sequences for the search for possible triggering pathogens and for other autoantigens. These predictions should be tested using both cell-based assays and ELISA tests with these 8 peptides coated on the plate.

Table 2.PNG
Table 2. Peptides belonging to M4 and M4 ChR that fulfill our selection criteria are collected on the left. On the right, their reverse sequences. These 16 peptides can be used in BLAST in order to serach for triggering pathogens and for other possible autoepitopes.

Supplementary material. The following two tables represent the first two steps of the analysis presented in this paper. M4 ChR 7-mer peptides that are cross-reactive to M3 ChR are collected in Table 1S, while those of them that are non-cross-reactive (or borderline) to M1 ChR are collected in Table 2S.

Table 1S. Peptides of M4 ChR that are cross-reactive to M3 ChR are collected in the first column. In the second column are collected the scores of these 7-mer peptides obtained from the comparison with M1 ChR. For those that obtained a score below 60, the score from the comparison with B2 AdR is reported in column 5. Positions of peptides of interest that belong to M3ChR and B2 AdR are collected in columns 4 and 6 respectively.

Table 2S.PNG
Table 2S. These 20 peptides are those M4 ChR peptides that cross-react to M3 ChR and at the same time are non-cross-reactive or borderline when compared to M1 ChR. Reactivity to B2 AdR is also indicated, as well as positions with respect to the plasma membrane and mean protrusion indexes. On the left are indicated those peptides of M4 ChR that pass the selection according to our criteria. Both a strict selection and a selection with more weak criteria are reported.