Maximum of a normal random vector

In evidenzaMaximum of a normal random vector

Ettore Majorana hadn’t got MATLAB

When Ettore Majorana first met Enrico Fermi, between the end of 1927 and the beginning of 1928, Fermi – who was already an acclaimed scientist in the field of nuclear physics – had just solved an ordinary differential equation of the second-order (whose solution is now commonly named the Thomas-Fermi function) – by numerical integration. It required a week of assiduous work for him to accomplish this task, with the aid of a hand calculator. Fermi showed the results (a table with several numbers) to Majorana, who was a 21 years old student of electrical engineering who had some vague idea of switching from the field of boring ways of providing electrical energy for boring human activities, to the quest for the intimate structure of matter, under the guide of Fermi, the brightest Italian scientific star of that period.

Majorana looked at the numerical table, as I wrote, and said nothing. After two days he came back to Fermi’s lab and compared his own results with the table made by Fermi: he concluded that Fermi didn’t make any mistake, and he decided that it could be worth working with him, so he switched from engineering to physics (Segrè E. 1995, page 69-70).

Only recently it has been possible to clarify what kind of approach to the equation Majorana had in those hours. It is worth mentioning that he not only solved the equation numerically, I guess in the same way Fermi did but without a hand calculator and in less than half the time; he also solved the equation in a semianalytic way, with a method that has the potential to be generalized to a whole family of differential equations and that has been published only 75 years later (Esposito S. 2002). This mathematical discovery has been possible only because the notes that Majorana wrote in those two days have been found and studied by Salvatore Esposito, with the help of other physicists.

I won’t mention here the merits that Majorana has in theoretical physics, mainly because I am very very far from understanding even a bit of his work. But as Erasmo Recami wrote in his biography of Majorana (R), a paper published by Majorana in 1932 about the relativistic theory of particles with arbitrary spin (Majorana E. 1932) contained a mathematical discovery that has been made independently in a series of papers by Russian mathematicians only in the years between 1948 and 1958, while the application to physics of that method – described by Majorana in 1932 – has been recognized only years later. The fame of Majorana has been constantly growing for the last decades.

The notes that Majorana took between 1927 and 1932 (in his early twenties) have been studied and published only in 2002 (Esposito S. et al. 2003). These are the notes in which the solution of the above-mentioned differential equation has been discovered, by the way. In these 500 pages, there are several brilliant calculations that span from electrical engineering to statistics, from advanced mathematical methods for physics to, of course, theoretical physics. In what follows I will go through what is probably the least difficult and important page among them, the one where Majorana presents an approximated expression for the maximum value of the largest of the components of a normal random vector. I have already written in this blog some notes about the multivariate normal distribution (R). But how can we find the maximum component of such a vector and how does it behave? Let’s assume that each component has a mean of zero and a standard deviation of one. Then we easily find  that the analytical expressions of the cumulative distribution function and of the density of the largest component (let’s say Z) of an m-dimensional random vector are

Eq 1.JPG

We can’t have an analytical expression for the integral, but it is relatively easy to use Simpson’s method (see the code at the end of this paragraph) to integrate these expressions and to plot their surfaces (figure 1).

surfaces.JPG
Figure 1. Density (left) and cumulative distribution function (right) for the maximum of the m components of a normal random vector. Numerical integration with MATLAB (by Paolo Maccallini).

Now, what about the maximum reached by the density of the largest among the m components? It is easy, again, using our code, to plot both the maximum and the time in which the maximum is reached, in function of m (figure 2, dotted lines). I have spent probably half an hour in writing the code that gives these results, but we usually forget how fortunate we are in having powerful computers on our desks. We forget that there was a time in which having an analytical solution was almost the only way to get mathematical work done. Now we will see how Majorana obtained the two functions in figure 2 (continuous line), in just a few passages (a few in his notes, much more in mine).

maximum.JPG
Figure 2.  The time in which Z reaches its largest value (left) and the largest value (right) in the function of the dimension of the normal random vector. The two curves in dotted lines are obtained through numerical integration (Simpson’s method) while the continuous lines are the function obtained by Majorana in an analytical way.
% file name = massimo_vettore_normale
% date of creation = 22/05/2019

clear all
delta = 0.01;
n(1) = 0.

for i=2:1:301;
  n(i) = delta + n(i-1);
  n_2(i) = - n(i);
end

for i=1:1:301
  f(i) = 0.39894228*( e^(  (-0.5)*( n(i)^2 )  ) );
end

for i=1:1:3
  sigma(1) = 0.;
  sigma(3) = sigma(1) + delta*( f(1) + ( 4*f(2) ) + f(3) )/3;
  sigma(2) = sigma(3)*0.5;
  for j=2:1:299
    sigma(j+2) = sigma(j) + delta*( f(j) + ( 4*f(j+1) ) + f(j+2) )/3;
  end  
end

for i=1:1:301
  F(i) =  0.5 + sigma(i);
  F_2(i) = 1-F(i);
end

for i=1:1:100;
  m(i) = i;
end

for i=1:1:301
  for j=1:1:100
    F_Z (i,j) = F(i)^j;
    F_Z_2 (i,j) = F_2(i)^j;
    f_Z (i,j) = 0.39894228*j*( F(i)^(j-1) )*( e^(  (-0.5)*( n(i)^2 )  ) );
    f_Z_2 (i,j) = 0.39894228*j*( F_2(i)^(j-1) )*( e^(  (-0.5)*( n(i)^2 )  ) );
  endfor
endfor

figure (1)
mesh(m(1:2:100),n(1:10:301),F_Z(1:10:301,1:2:100));
grid on
hold on
mesh(m(1:2:100),n_2(2:10:301),F_Z_2(2:10:301,1:2:100));
xlabel('m');
ylabel('t'); 
legend('F',"location","NORTHEAST");

figure (2)
mesh(m(1:2:100),n(1:10:301),f_Z(1:10:301,1:2:100));
grid on
hold on
mesh(m(1:2:100),n_2(2:10:301),f_Z_2(2:10:301,1:2:100));
xlabel('m');
ylabel('t'); 
legend('f',"location","NORTHEAST");

Asymptotic series

I have always been fascinated by integrals since I encountered them a lifetime ago. I can still remember the first time I learned the rule of integration by parts. I was caring for my mother who was dying. That night I was in the hospital with her, but she couldn’t feel my presence, she had a tumour in her brain and she was deteriorating. And yet I was not alone, because I had my book of mathematics and several problems to solve. But when my mind was hit by the disease for the first time, about a year later, and I lost the ability to solve problems, then real loneliness knocked at my door.

Now, why am I talking about the integration by parts? Well, I have discovered a few days ago, while studying Majorana’s notes, that integration by parts – well known by students to be a path towards recursive integrations that usually leads to nowhere – is, in fact, a method that can be useful for developing series that approximate a function for large values of x (remember that Taylor’s polynomials can approximate a function only for values of x that are close to a finite value x_0, so we can’t use them when x goes to ∞). Majorana used one such a series for the error function. He developed a general method, which I tried to understand for some time, without being able to actually get what he was talking about. His reasoning remained in the back of my mind for days, while I moved from Rome to Turin, where I delivered a speech about a paper on the measure of electric impedance in the blood of ME/CFS patients; and when I cried, some minutes later, looking at my drawings put on the screen of a cinema, Majorana was with me, with his silence trapped behind dark eyes. A couple of days later, I moved to a conference in London, searching for a cure that could perhaps allow my brain to be normal again and I talked with a famous scientist who worked on the human genome project. Majorana was there too, in that beautiful room (just a few metres from Parliament Square), sitting next to me. I could feel his disappointment, I knew that he would have found a cure, had he had the chance to examine that problem. Because as Fermi once said to Bruno Pontecorvo, “If a problem has been proposed, no one in the world can resolve it better than Majorana” (Esposito S. et al. 2003). Back in Rome, I gave up with the general method by Majorana and I found the way to calculate the series from another book. The first tip is to write the error function as follows:

Eq 2.JPG

Now by integrating by parts, one gets

Eq 3.JPG

But we can integrate by parts one other time, and we get

Eq 4.JPG

And we can go on and on with integration by parts. This algorithm leads to the series

Eq 5.JPG

whose main property is that the last addend is always smaller (in absolute value) than the previous one. And even though this series does not converge (it can be easily seen considering that the absolute value of its generic addend does not go to zero for k that goes to ∞, so the Cauchy’s criteria for convergence is not satisfied) it gives a good approximation for the error function. From this series, it is easy to calculate a series for the Gaussian function (which is what we are interested in):

Eq 6.JPG

A clever way to solve a transcendental equation if you don’t want to disturb Newton

Taking only the first two terms of the series, we have for the cumulative distribution function of Z the expression:

Eq 7.JPG

The further approximation on the right is interesting, I think that it comes from a well-known limit:

Eq 8.JPG

Now we can easily calculate the density of Z by deriving the cumulative distribution function:

Eq 9.JPG

With a further obvious approximation, we get:

Eq 10.JPG

In order to find the value of x in which this density reaches its largest value, we have to search for the value of x in which its derivative is zero. So we have to solve the following equation:

Eq 11.JPG

Which means that we have to solve the transcendental equation:

Eq 12.JPG

Majorana truncated the second member of the equation on the right and proposed as a solution the following one:

Eq 13.JPG

Then he substituted again this solution in the equation, in order to find ε:

Eq 14.JPG

With some further approximations, we have

Eq 15.JPG

So Majorana’s expression for the value of x in which the density of Z reaches its maximum value is

Eq 16.JPG

I have tried to solve the transcendental equation with Newton’s method (see the code below) and I found that Majorana’s solution is a very good one (as you can see from figure 3). Now, If we compare the approximation by Majorana with what I obtained using numerical integration at the beginning (figure 2) we see that Majorana found a very good solution, particularly for the value of x_M. Note: the transcendental equation that has been solved here seems the one whose solution is the Lambert W function, but it is not the same!

trascendental.JPG
Figure 3. The value of x in which the density of Z reaches its largest value. Majorana’s solution (continuous line) and a numerical solution obtained by means of Newton’s method (dotted line).
Vol II.JPG
Figure 4. From the original manuscript by Majorana, which can be found here (page 70). These are some of the passages that I have discussed here.
% file name = tangenti
% date of creation = 08/06/2019

clear all

x(1) = 1;            %the initial guess
for i=1:1:100
  m(i) = i;
end

for i=1:1:100
  for j = 2:1:1000
   f(j-1) = exp( 0.5*( x(j-1)^2 ) ) - ( m(i)/( x(j-1)*sqrt(2*pi) ) );
   f_p(j-1) = x(j-1)*exp( 0.5*( x(j-1)^2 ) ) + ( m(i)/( (x(j-1)^2)*sqrt(2*pi) ) );
   x(j) = x(j-1) - ( f(j-1)/f_p(j-1) );
   if ( abs(x(j)) < 0.001 )
    break;
   endif
   max_t (i) = x(j);
  endfor
endfor

% the aproximations by Majorana

for j=1:1:100
  max_t_M (j) = sqrt(log(j^2)) - ( log(sqrt(2*pi*log(j^2)))/sqrt(log(j^2)) );
endfor

% it plots the diagrams

plot(m(1:1:100),max_t (1:1:100),'.k','Linewidth', 1)
xlabel('m')
ylabel('time for maximum value')
grid on
hold on
plot(m(1:1:100),max_t_M (1:1:100),'-k','Linewidth', 1)

legend('numerical integration',"Majorana's approximation", "location", 'southeast')

Epilogue

From 1934 to 1938 Majorana continued his studies in a variety of different fields (from game theory to biology, from economy to quantistic electrodynamics), but he never published again (R), with the exception for a work on the symmetric theory of electrons and anti-electrons (Majorana E. 1937). But it has been concluded by biographers that the discoveries behind that work were made by Majorana about five years earlier and yet never shared with the scientific community until 1937 (Esposito S. et al. 2003). And in a spring day of the year 1938, while Mussolini was trying his best to impress the world with his facial expressions, Ettore became a subatomic particle: his coordinates in space and their derivatives with respect to time became indeterminate. Whether he had lived in a monastery in the south of Italy or he had helped the state of Uruguay in building its first nuclear reactor; whether he had seen the boundless landscapes of Argentina or the frozen depth of the abyss, I hope that he found, at last, what he was so desperately searching for.

He had given his contribution to humanity, so whatever his choice has been, his soul was already safe. And as I try to save my own soul, going back and forth from mathematics to biology, in order to find a cure, I can feel his presence. The eloquence of his silence trapped behind dark eyes can always be clearly heard if we put aside the noise of the outside world. And it tells us that Nature has a beautiful but elusive mathematical structure which can nevertheless be understood if we try very hard.

In the meanwhile, I write these short stories, like a mathematical proof of my own existence, in case I don’t have further chances to use my brain.

Until time catches me.

 

 

 

On the module of random vectors

In evidenzaOn 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.

cover
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

Cattura.PNG

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:

density.PNG

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:

varianza.png

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:

Cauchy.png

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

Nash.png

given that

Nash2.png

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.

sphere.PNG
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

notazione.png

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

Quadrato.PNG

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

Radice.PNG

Another important result that we have to consider is that given

linera combination.PNG

then

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:

module.PNG

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

maxwell.PNG

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.;

endfor

% 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;

endfor

% 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. ) );

endfor

% 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;

end

endfor

% It plots f and F

figure (1)

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

legend('Density',"location","NORTHEAST");

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

legend('Probability',"location","NORTHEAST");

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

ylabel('temperature (Kelvin)');

grid on

 

Is it that bad?

Is it that bad?

When I got sick, about 20 years ago, for the first time I started thinking about diseases and loss of health. And I remember coming to the conclusion of how fortunate I was, from a physical standpoint. Not of how fortunate I had been in the past (that was obvious), when I could conquer mountains, running on rocks with my 15 Kg bike on the shoulders; no, of how fortunate I was in that very moment, while confined at home, mostly lying horizontally. I was fortunate, I could still move my hands, I had still all my body, even if I couldn’t use it in the outside world, even if most of the usual activities of a 20 years old man were far beyond reaching, I realized how fortunate I was. And after 20 years I still say that, as for the physical functioning, I am a very lucky man. There are even some years, during the core of the summer, in which I can run for some minutes in the sun. I am blessed.

I am aware that this might offend some patients, but I think that those who have ME/CFS in the vast majority of cases are fortunate too, from a physical standpoint. Yes, it is annoying to need help from others for so many things, but with assistance and some arrangements, you can go on with a productive life… unless you have cognitive impairment.

Actually, I never felt fortunate concerning my cognitive functioning. That has been a true tragedy, I have lost my entire life because of the cognitive damage, by any means because of the limitations of my body. I would have had a meaningful life (according to any standard) even with physical limitations way worse than mine if I had had a functioning brain. But I lost it when I was 20 before I could get the best from it and that is the only real tragedy.

As an example of what I am saying, consider a person like Stephen Hawking: he has been a prolific scientist, he had a stellar career (literally), he shaped our culture with his books for the general public, he was a father of three, and so forth. And yet for most of his life, while he was doing all these things, his physical functioning was way worse than the one of the average ME/CFS patient (worse even of a very severe patient, probably). For many years he could use only the muscles of the head. That’s it. Think about that. Was he missing or invisible? Definitely not!

Was he a disabled man? Yes, technically speaking he was, but in fact, he has never really been disabled, if we judge from his accomplishments. He once said that the disease somehow even helped him in his work, because he could concentrate better on his quantum-relativistic equations¹.

So what I am saying is that if we consider the physical functioning in ME/CFS, it is a negligible problem in most of the cases. The only thing that matters is the impairment in cognition (in the cases in which it is present), especially if it starts at a young age (consider that after you reach your thirties there is very little chance that your brain will give a significant contribution to humanity, even if you are perfectly healthy, so it is not a big loss if you get sick after that age). That is disastrous and there is no wheelchair you can use for it. For all the other symptoms you can find a way to adjust, just as Hawking did in a far worse situation.

This might be one of the reasons for the bad reputation of ME/CFS: there isn’t awareness about cognitive issues, no one talk about them (and in some cases, they are in fact not present at all). But I am pretty sure that the real source of disability in these patients, the lack of productivity, is due to their cognitive problems (when present). Also because in a world like ours, you can work even without using most of your muscles. It is not impossible, I would say it is the rule for a big chunk of the population.

The following one is an interview with Norwegian neurologist Kristian Sommerfelt, in which he points out some analogous considerations. He has done some research on ME/CFS with the group of Fluge and Mella, including the well-known study on pyruvate dehydrogenase. From the subtitles of the video (minute 3:38):

“This [the cognitive problem] is a very typical ME symptom and some of what I believe causes the main limitation. I don’t think the main limitation is that they’re becoming fatigued and exhausted by moving around, walking, running, or having to sit still. If it were just that, I think many ME-patients could have had a much better life. But the problem is that just actively using the mind leads to problems with exactly that, using the mind. It comes to a stop, or slow down, depending on how ill they are.”


¹ It seems that Stephen Hawking had a very unusual presentation of Amyotrophic Lateral Sclerosis (ALS): one half of those with this disease die within 30 months from the first symptom; moreover one out of two ALS patients has a form of cognitive impairment which in some cases can be diagnosed as frontotemporal dementia [R]. So, Stephen Hawking was somehow lucky, in his tragedy, and he doesn’t represent the average ALS patient. I mentioned his case as an example of a person with very severe physical impairment and no apparent cognitive decline, not as an example of the average ALS patient. 

About me

About me

The following video is meant to be a presentation of both the blog and of myself. As I started improving again, some days ago, I decided to record this monologue, so that there could be a video memory of my struggle.

This winter I have spent almost three months in South America, to see if I would have improved during the austral summer, as I usually improve during the Italian summer; and in fact, I did improve. When I came back to Italy (in March) I had a relapse, though. For the last three months, I have been mostly horizontal, without reading or thinking for most of the time.

Now I am climbing the mountain again: I started my rehabilitation reading novels some days ago, then I switched to simple calculations and now I have written my first small code since March. And when I will reach the cognitive level I had about 20 years ago just before I got sick, I will lose everything for months (or years) and I will have to wait without thinking much (despite my best efforts) until I can start all over again…

The myth of Sisyphus has been shaped after me.

 

 

My saviour

One of my short-lived summer improvements (2013). During all these years, as soon as I started feeling better, I opened my books, even before taking a shower and having my hair cut. Happy as a child for most of the time, but also profoundly saddened for the time lost, especially at the beginning of the improvement, when I could realize how much time had passed from the previous positive phase.

97284908_10220774536729711_4708150532923981824_n

Each time I had to start exactly from where I had left many months or even years before (the longest gap has been 5 years without studying). I had to do a cognitive rehabilitation each time, learn again how to read properly, how to do math, how to discipline my thoughts, how to code. It is a hard process each time. And then, a few weeks after, when I recovered enough to function mentally, I relapsed again.

I am pretty sure that only this complete, obsessive devotion to studying has saved me from very bad cognitive disability.

Before getting sick, coding and math had taught me how to think. Then, when I became ill, each equation I wrote, each drawing and code, all those efforts made to bring my soul back from wherever it was, they kept me alive for all these years.

My first time with LaTeX

Each mathematical formula I have used in my blog till now was an image: I wrote them in a Word document using the equation editor by Microsoft (that I like very much), then I captured the screen and saved it as an image and then I uploaded them in my blog.

In the past, I tried to find a way to write them as text, without succeeding. But in fact, it is possible, using the LaTeX syntax and a few HTML commands.

These are the first “true” mathematical formulae present in this blog. I will use this page to test the LaTeX syntax for mathematical expressions (a handbook is available here).

\int_{A_1}\int_{A_2}\frac{\partial \Psi(x,y)}{\partial x}dxdy

i\hbar\frac{\partial}{\partial t}\left|\Psi(t)\right>=H\left|\Psi(t)\right>

A symbol within the paragraph: \pi. A inequality within the paragraph: x_i>0. This is an example of an expression with a pedix that has a pedix: f_{X_1}.

The time machine

I am aware that these are just messages floating in the silence, stored somewhere in the planet as binary numbers. I am writing to myself, mainly, from my remote hiding place.

I have travelled through ages, without really being part of them. All alone with my problem. As a patient with a rare disease that doesn’t even have a proper description, I do not belong to humankind.

But humans have paradoxical behaviours, they care more about a man who lived five thousand years ago in the north of Italy, trapped in the ice of our highest mountains, than of clochards that live right now in pain and loneliness in their community. So it might be that generations from now, someone will find these notes, an archaeologist who will try to build my story, from fragments of what I left behind: drawings and calculations. Mathematics is a universal language, after all, and to some extent, even art is universal; not always but often, good art is forever.

If I fail my mission, history will never record my existence. But it might be that at some point in the future someone will find these notes frozen in the ice of a planet long forgotten by humans themselves, as we now have forgotten Africa, the place we all come from.

 

Why I study my own disease

Why I study my own disease

A lot of patients have asked me why I use my little energies to study my disease, instead of just waiting for science to conquer it.

There are many reasons, the first one being that I am desperate because of the cognitive disability that is worse than death. I am not concerned about the physical limitations, at all, even though I have been mostly housebound for the last 20 years. Another reason is that I like computational biology, and I started studying engineering before getting sick with the idea of switching to bioengineering after graduating. So, this is my job.

The other point is that even though I fit the criteria for ME/CFS, I have a rare disease, granted that ME/CFS is not a rare disease: it has a prevalence of 0.4% according to some studies, so it is relatively common. Then why am I a rare patient? Let’s do some calculations: since the median age at onset is 36.6 years with a standard deviation of 12.3 years, the proportion of males is 19% and the proportion of those who are housebound for most of their disease is 25% (R), the probability for a ME/CFS patient of having my same characteristics is (assuming that these are all independent random variables) given by p=0.0075. So, less than 1 out of 100 ME/CFS patients has my type of illness. If we consider also that pain or aching in muscles is present in 59% of patients and it is mostly absent in my case, the above-mentioned probability is even lower: p=0.0031. Which means that only 3 ME/CFS patients out of 1’000 have my illness.

Taken all these data together, the prevalence of my disease in the general population becomes 1/100’000, which means that I have a rare disease, according to the European definition (where a disease is defined rare if it has a prevalence ≤1/2’000) [R]. In the US a disease is defined rare only if it has a prevalence ≤1/200’000, though [R].

And besides that, it is pretty obvious that I am a unique case. I have never found a patient like me, so far. It can also be noted that in Italy, given a population of 60 million inhabitants, those who have my condition are only 600. This might be the reason why I haven’t met them, yet.

This means that I am probably the only one who is studying my disease, on the planet! This is why I’m doing what I’m doing.

91569815_10220268929849855_9111778109215997952_o

Above, my personal book of immunology, built page by page, paper by paper. I have three other books about this discipline. One of them was a gift from a neurologist that perhaps thought that gift was the only thing she could do to save me. Another one is a very sophisticated text on cutting edge immunology. But this one is the best one because I have selected and read each one of its pages. It is by no means a complete book, it is mainly focused on B cells and B cell autoimmunity, but it has been very useful.

I have built several other books like this, on computational methods in immunology, on metabolism, on neurosciences, on microbiology, and on some diseases: Lyme, ME/CFS, mast cell activation, POTS…

I have learned a great deal, even though outside academia. But I had no choice, I have been too sick and too slow to study at university: I had only a few weeks in which I could study, and then months or years in which I had to wait. This has been my routine. Moreover, given the lack of energy and time, I had to study only what was truly important for my health. Because my goal was to cure myself and save me from a lifetime of cognitive disability.

Testing hypotheses

Testing hypotheses

Introduction

My ME/CFS improves during summer, in the period of the year that goes from May/June to the end of September. I don’t know why. I have several hypotheses. One possible reason for the improvement in summer is an interaction between the light from the Sun and some parts of my physiology, the immune system for instance. We know that ME/CFS tends to have an oscillating course in most of the patients (Chu L. et al. 2019), but the presence of a seasonal pattern in this patient population has not been investigated so far, to my knowledge. And yet, if you ask directly to patients, many of them say that they feel better in summer. Unfortunately, we don’t have scientific data on that, this is an area worth investigating with some carefully done survey.

Seasonal variation of the immune system

The immune system has a high degree of variation for several reasons (Brodin P et Davis MM 2017). In particular, there are studies about the seasonal fluctuations in the expression of some crucial genes of the immune response (Dopico XC et al. 2014).

How does this regulation happen? Different mechanisms are possible, some of them might be related to changes in the light we receive from the Sun as the Earth rotates around it. We know that the length of the day has an effect on innate immunity: the more the hours of light, the lower the power of the innate immune system (Pierre K. et al. 2016). We also know that ultraviolet radiation, particularly UVB, is an agonist for the aryl hydrocarbon receptor (AhR) (Navid F. et al. 2013). This receptor seems to reduce the expression of the major histocompatibility complex II (MHC II) in dendritic cells (DCs), thus reducing their antigen-presenting activity (Rothhammer V. et Quintana F.J. 2019). UVB might be able to reach dendritic cells when they are circulating near the skin, during summer, thus inhibiting their antigen-presenting activity. Infrared radiation, on the other hand, seems to have an effect on energy metabolism: in Fall we lose a significant amount of infrared radiation in a wavelength range (0.7-1.0 nm) that is known to have an effect on mitochondrial activity (Nguyen L.M. et al. 2013) and it might perhaps have an indirect effect on immunity too.

As further proof of seasonal fluctuation in immunity, some immunological diseases have this kind of seasonality: Rheumatoid arthritis (Nagamine R. et al. 2014) and Rheumatic fever (Coelho Mota C.C. et al. 2010) are two examples. Moreover, the prevalence of Multiple Sclerosis is directly proportional to the latitude (Simpson S. et al. 2011). We also know that there is seasonal fluctuation in serum autoantibodies (Luong T.H. et al. 2001).

Of course, sunlight might be just one of the variables into play. The other aspect I am considering is the seasonal distribution of some common pathogens. Streptococcus, Enteroviruses and Fungi of the genus Penicillium are known to have a seasonal distribution with a peak in Fall and/or Winter (Ana S.G. et al. 2006), (Salort-Pons M et al. 2018), (Coelho Mota C.C. et al. 2010). Common influenza has this pattern too. Rheumatic fever, a disease due to an abnormal immune response to Streptococcus, has its flares in Fall because Streptococcus is more common in that period of the year (Coelho Mota C.C. et al. 2010). Even the composition of the gut microbiota has a seasonal pattern (Koliada A. et al. 2020). I am currently investigating my immunosignature, measured with an array of 150.000 random peptides, to see if I can find some relevant pathogen in my case. You can find this study here.

An experiment

I moved from Rome (Italy) to Rosario (Argentina) at the beginning of January. I was very sick and I steadily improved after about 40 days. I became a less severe ME/CFS patients and I could work several hours a day and care for myself, granted that I did not exceed with aerobic exercise. At the end of March, I started deteriorating as it usually happens at the end of September, when I am in Rome. In order to study this phenomenon, I have built a complete model of solar radiation at sea level, which considers the inclination of sunrays in function of the latitude and of the day of the year. It takes into account the effect of the atmosphere (both diffusion and absorption) and the eccentricity of the orbit (Maccallini P. 2019). If you look at the figure below (a byproduct of my mathematical model) you can see that when I started deteriorating in Rosario, the power of sunrays at noon in that city was still as high as it is in Rome during the summer solstice (this is due to the fact that the Earth is closer to the Sun in this period and to the fact that Rosario is closer to the Equator than Rome is).

Rome_vs_Rosario

So I have to discard the original idea that the power within the infrared range, or the ultraviolet radiation, or the visible one is responsible for my improvement in summer. If I still have to consider that sunlight has something to do with my improvement, I must conclude that it is the length of the day the relevant parameter: I may need more than 12 hours of light to feel better. Why? Because the longer the day, the lower the strength of the innate immunity. This is now my working hypothesis and I will start from the following mathematical model to pursue this research: (Pierre K. et al. 2016).

I will also use full-spectrum lamps early in the morning and in the evening to reproduce a 15 hours day, so to dampen down my innate immune system in a safe, drug-free way. I have to reproduce a day of 15 hours and see what happens. In the figure below the hours of the day at dawn and at dusk and the length of the day for Rome, for each day of the year (this is also a plot from my model).

ore di luce

What follows is the script I have coded in order to plot the first figure of this post. More details on this model of solar radiation are here: (Maccallini P. 2019). As a further note, I would like to acknowledge that I started pursuing this avenue in the summer of 2009: I was building the mathematical model of solar radiation (see figure below, made in 2009) but as the summer finished, I turned into a statue and I had to stop working on it. When I improved, about a year later I started working on the systematic analysis of the mechanical equilibrium of planar structures (it is a chapter of this book). I am proud of that analysis, but it has not been very useful for my health…

trigonometria sferica 2

% file name = sun emissive power sea level Rosario vs Roma
% date of creation = 07/11/2019
% sun emissive power per unit area, per unit wavelength at sea level
clear all
% three parameters of the orbit
A = 6.69*( 10^(-9) ); % 1/km
B = 1.12*( 10^(-10) ); % 1/km
delta = pi*313/730;
% the two parameters of Plunk's law
C_1 = 3.7415*( 10^(-16) ); % W*m^2
C_2 = 1.4388*( 10^(-2) ); % mK
% Stefan-Boltzmann parameter ( W/( (m^2)*(K^4) ) )
SB = 5.670*( 10^(-8) );
% radius of the photosphere (m)
R_S = 696*(10^6); % m
% temperature of the photosphere (K)
T_S = 5875;
% conversion of units of measurments
N = 20; % dots for the equator
R = 3.8; % radius of the orbit
ro_E = 1.3; % radius of the earth
lambda_Rosario = -32*pi/180; % latitude of Rosario (radiants)
lambda_Roma = 41*pi/180; % latitude of Rome (radiants)
delta = 23.45*pi/180; % tilt angle
% the array of theta
theta(1) = 0; % winter solstice (21/22 December)
i_ws = 1;
day = 2*pi/365;
days = [1:1:366];
for i = 2:366
theta(i) = theta(i-1) + day;
if ( abs( theta(i) - (pi/2) ) <= day )
i_se = i; % spring equinox (20 March)
endif
if ( abs( theta(i) - pi ) <= day )
i_ss = i; % summer solstice (20/21 June)
endif
if ( abs( theta(i) - (3*pi/2) ) <= day )
i_ae = i; % autumn equinox (22/23 September)
endif
endfor
% the array of the radius (m)
for i=1:1:366
o_omega (i) = (10^3)/[ A + ( B*sin(theta(i) + delta ) ) ]; % m
endfor
% the array of the wavelength in micron
N = 471;
L(1) = 0.3;
L(N) = 5.0;
delta_L = ( L(N) - L(1) )/(N-1);
for j = 2:N-1
L (j) = L(j-1) + delta_L;
endfor
% the array of beta*L
load beta_int.mat S
% the array of L in metres
L_m = L*( 10^(-6) );
% angle psi
psi(1) = 0;
minute = pi/(12*60);
for i = 2:(24*60)+1
psi(i) = psi(i-1) + minute;
endfor
% -----------------------------------------------------------------------------
% Rosario
lambda = lambda_Rosario
% angle between n and r at noon in Rosario
for i= [i_ws, i_ae, i_ss, i_se]
for j=1:(24*60) + 1
% scalar product between n and r
scalar_p(j) = [cos(lambda)*sin(psi(j))*cos(delta) + sin(lambda)*sin(delta)]*( -cos(theta(i)) )+ [(-1)*cos(lambda)*cos(psi(j))]*( -sin(theta(i)) );
endfor
% value of psi at noon
for j=1:(24*60) + 1
if ( ( scalar_p(j) ) == ( max( scalar_p ) ) )
j_noon = j;
psi_noon (i) = psi(j);
endif
endfor
% angle between n and r at noon
cos_gamma (i) = scalar_p(j_noon);
endfor
% the array of the emissive power (W/(m^2)*micron) in Rosario
for i = i_se:i_se
for j=1:N
num = C_1*( (R_S)^2 );
den = ( (L_m(j)^5)*( (e^(C_2/( L_m(j)*T_S ))) - 1)*( (o_omega(i))^2 ) )*10^6;
power(j,i) = ( num/den )*( e^(-S(j)/cos_gamma (i)) );
endfor
% plotting
plot (L (1:N), power(1:N,i), '-r', "linewidth", 2)
xlabel('wavelenght ({\mu})');
ylabel('W/m^{2}{\mu}');
axis ([0.3,5,0,1500])
grid on
endfor
hold on
% -----------------------------------------------------------------------------
% Rome
lambda = lambda_Roma
% angle between n and r at noon in Rosario
for i= [i_ws, i_ae, i_ss, i_se]
for j=1:(24*60) + 1
% scalar product between n and r
scalar_p(j) = [cos(lambda)*sin(psi(j))*cos(delta) + sin(lambda)*sin(delta)]*( -cos(theta(i)) )+ [(-1)*cos(lambda)*cos(psi(j))]*( -sin(theta(i)) );
endfor
% value of psi at noon
for j=1:(24*60) + 1
if ( ( scalar_p(j) ) == ( max( scalar_p ) ) )
j_noon = j;
psi_noon (i) = psi(j);
endif
endfor
% angle between n and r at noon
cos_gamma (i) = scalar_p(j_noon);
endfor
% the array of the emissive power (W/(m^2)*micron) in Rosario
for i = [i_ae, i_ss]
for j=1:N
num = C_1*( (R_S)^2 );
den = ( (L_m(j)^5)*( (e^(C_2/( L_m(j)*T_S ))) - 1)*( (o_omega(i))^2 ) )*10^6;
power(j,i) = ( num/den )*( e^(-S(j)/cos_gamma (i)) );
endfor
endfor
hold on
plot (L (1:N), power(1:N,i_ae), '-k', "linewidth", 2)
plot (L (1:N), power(1:N,i_ss), '--k', "linewidth", 2)
legend ('spring equinox in Rosario', 'autumn equinox in Rome', 'summer solstice in Rome', "location",'NORTHEAST')
hold on
plot ([0.4,0.4], [0,1500], '--k', "linewidth", 1)
plot ([0.7,0.7], [0,1500], '--k', "linewidth", 1)

Back home

Back home

Continuation of this post.

Forty-four hours of travelling, in total, from Rosario to Rome, by pullman, by plane, and by train. With 40 kilos of books and papers.

I had a flight for Rome that was programmed to take off from Ezeiza, the International airport of Buenos Aires, on April 13th, but I decided to take the one organized by the Italian government for March 23th, a special flight set up to bring back home Italian citizens abroad, before a complete shut down of international flights from Argentina to our country. There were no flights from Rosario, my city, to Buenos Aires, though, but I managed to find a company that organizes transportations by pullman from one city to the other, in Argentina: Tienda Leon.

So, on March 22nd, I moved to Ezeiza where I waited several hours before sitting on my chair, on a brand new Boeing 787 bearing the colours of the Italian company Neos.

While at the airport, I met most of the Italians that were going to get the same flight, all wearing their masks. Some of them with some very fancy models, that made them look like a Star Wars character. I was there, well aware that I was going far beyond the limit set by my disease. I had to lay down continuously and I could see how frail I was in comparison with the other passengers waiting for the flight. No one knew how sick I was, I told nobody. No one knew that I have been living in my bedroom for most of the last 20 years. And that this was the very first long travel abroad for me.

I have just received the notification that my flight for April 13th has been cancelled, so my choice to come back as soon as possible has been a wise one. I took that decision also because of the advice from the diplomatic offices of the Italian Consulate in Rosario.

A friend has crafted the picture above, not knowing how much Indiana Jones has meant for me when I was a teenager. But, even though an appealing adventure, the tragedy behind it is real, it is not a movie. Once in Milan, I could start seeing the effects of the pandemic in the eyes of the staff of the airport of Malpensa: the fear and the concern. Then I moved from Milan to Fiumicino, where I found a train for Rome, my city. A city that I left two and a half months ago full of life and noise, now empty as in a dream.

 

From Argentina to Italy, during a pandemic

From Argentina to Italy, during a pandemic

It has been a great ride, my almost-three-month period here in Rosario, next to the huge slowly flowing river of Paranà. This is a city full of life, embraced by the warmest summer I have ever seen. Populated by wonderful citizens.

WhatsApp Image 2020-03-09 at 09.37.08

I have been living in an apartment where the sun awakes me very early in the morning, through a wide window next to my bed. I could see the roofs of the centre of the city as I opened my eyes, including the top of the Monumento National a la Bandera, a gigantic building that celebrates this great nation. For the second half of the day, I had the light from the opposite window and I could follow its changes, while I was working, as the hours passed by; I saw every day the same magic ritual: as the photons from our star went through thicker layers of the atmosphere, they changed their frequency, turning redder and redder, culminating in a warm explosion, just before the night.

And in the meanwhile, news from Italy was scarier and scarier and the hypothesis that the new coronavirus could reach this continent was more obvious as the weeks passed by. Now we have the virus here, and president Alberto Fernandez has declared the state of quarantine from March 20th.

At that point, the connection between Argentina and other countries (including Italy) has become uncertain; my flight planned for March 28th has been cancelled and I have decided to get one of the special flights organized by the Italian government to bring back its citizens from Argentina, before a complete shut down of international travels. So I had about 48 hours to find a means of transport from Rosario to the airport in Buenos Aires, where the flight will take off tomorrow, at 1:00 AM.

But there was no way I could find a flight from Rosario to Buenos Aires in such a short time, also because of the shut down of Argentina, and no trains were available. Fortunately enough I have found a Pullman, and I am going to leave this apartment in a few hours.

My come back to Italy is becoming more and more adventurous also because I will land in Milan, one of the places most hit by the infection in the whole planet. There I have to reach Rome. I have been able to find a plane from Milan to Rome, so it will be possible to be at home on the evening of March 23rd. I have to avoid to get the virus though, during this travel. I will be exposed to it for sure, so I am taking any possible measure to ensure my safeness.

This travel to Argentina has been a success. My health has improved, even though now I am deteriorating again, as was expected, as the light of the summer of the southern hemisphere becomes weaker. But I have been able to use my new energies to write and submit a paper on the cingulate cortex in ME/CFS, I have gone further with my studies on the analysis of the immunosignature (measured using random peptides) in my own serum (R), I have started the study of a mathematical model for the diffusion of Coronavirus 19 among the Italian population (R). I have learnt a great deal about computational neuroanatomy (R) and neurosciences in general. I have finished a complete model for solar radiation at sea level (R) and I might have found one of the environmental parameters that determine my improvement during summer. And yes, I have also been able to draw a portrait.

WhatsApp Image 2020-03-18 at 10.07.35

The adventure in the realm of science and art has been great, now I have to live the adventure of coming back home going through a world that is facing one of the greatest health challenges of the last century.