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.
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:
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 one
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.
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.
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).
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).
Figure 6. The reaction catalyzed by creatine kinase: a molecule of ADP has converted to ATP thanks to the phosphate group carried by phosphocreatine.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.
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:
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
By operating the substitutions
we obtain
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:
Let’s consider the inner integral and the following substitutions
We have:
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:
We have
If we put these two integrals inside I, we obtain
The two inner integrals are two difficult ones, again. The first one can be solved using the same exponential integral used above:
For the other inner integral we have
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:
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:
where C is the covariance matrix, given by
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.
I have performed a set of analyses of experimental data previously 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 of autoimmunity to beta 2 adrenergic receptors, also reported in ME/CFS patients.
Introduction
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.
Methods
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.
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, and scores below the blue line indicate a lack of cross-reactivity.
Selection criteria. Our purpose is to predict which 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, the selection criteria for M4 ChR epitopes are as follows:
they must be cross-reactive to M3 ChR;
they must be non-cross reactive to M1 ChR or borderline;
they must present a mean protrusion index ≥0.5;
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.
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.
Results
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 to 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. 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.
Discussion
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-denatured 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.
Conclusion
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. 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. 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 weaker criteria are reported.
We present an attempt at exome analysis in two ME/CFS patients. Pt. 1 presents a mild form of carboxypeptidase N (CPN1) deficiency (a missense in exon 3) while Pt. 2 revealed two rare intronic variants in the same gene. CPN1 is an enzyme that inactivates kinins and complement proteins split products (such as C4a, a known anaphylatoxin). Therefore, CPN1 deficiency could explain C4a increase after exercise and mast cell abnormalities previously reported in ME/CFS. It could also explain the high prevalence of POTS in ME/CFS since kinins are vasodilators.
Introduction
Myalgic encephalomyelitis/chronic fatigue syndrome (ME/CFS) is a debilitating disease characterized by cognitive deficits, fatigue, orthostatic intolerance with symptoms exacerbated after exertion (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 and gut flora (Blomberg J. et al. 2018), brain (Zeineh MM. et al. 2014). In this population of patients, several abnormalities have been found to be triggered by exercise, such as abnormal aerobic performance (Snell C. et al. 2013), enhanced gene expression of specific receptors (White AT. et al. 2012), abnormal gut flora translocation (Shukla SK et al. 2015) and failure in blood clearance of complement protein 4 split product A (C4a) (Sorensen B et al. 2003). An increase in C4a is part of the human physiologic response to physical exercise, but these levels return to baseline within 30 minutes to 2 hours (Dufaux B et al. 1991) while in ME/CFS there is a peak in serum C4a six hours after exertion. A possible explanation for slow C4a inactivation could be a problem in carboxypeptidase N (CPN1), an enzyme involved in the inactivation of C3a, C4a, C5a. CPN1 is required for kinins inactivation too, such as bradykinin, kalladin (Hugli T. 1978), (Plummer TH et Hurwitz MY 1978), that are vasodilators. We report on the case of a ME/CFS patient (Pt. 1) with a missense variant in CPN1 gene that is linked to reduced function of the enzyme and of another ME/CFS patient (Pt. 2) with rare variants in introns 1 and 6 of the same gene with uncertain significance (table 1, figure 1).
Table 1. This is a collection of the variants within gene CPN1 found in Pt. 1 and 2. The verdict is the one given by VarSome. A damaging SNP present in exon 3 of CPN1 from Pt. 1 and two rare SNPs found in intron 1 and intron 6 of CPN1 from Pt. 2 are highlighted in orange. Figure 1. SNPs and InDels along gene CPN1 found in Pt. 1 (first row) and Pt. 2 (second row) as reported by IGV.
Materials and Methods
Whole exome sequencing (WES) has been performed on cells from the saliva of two ME/CFS patients, with an average 100X coverage (Dante Labs). The first search for pathogenic variants and insertions/deletions was performed with the software EVE, provided by Sequencing.com. A further refinement of the search was conducted by manual insertion of these SNPs in VarSome. The search for possible unknown pathogenic variants within the gene for CPN1 has been performed using Integrative Genome Viewer (IGV), an opensource tool for genetic data analysis.
Results
Results from the analysis of the two exomes performed with EVE and refined with VarSome are collected in table 2 (Pt. 1) and table 3 (Pt. 2).
Pt. 2 is carrier of a mitochondrial disease (table 3, first line): a missense in gene for medium-chain acyl-CoA dehydrogenase (MCAD) which leads to mild functional impairment of the enzyme involved in the oxidation of fatty acids (44% residual activity) (Koster KL. et al. 2014).
Pt. 2 is also homozygous for a variation in gene arylsulfatase A (ARSA) that is linked to a residual activity of only 10% of normal (Gomez-Ospina N. 2010). Arylsulfatase A deficiency (also known as metachromatic leukodystrophy or MLD) is a disorder of impaired breakdown of sulfatides (cerebroside sulfate or 3-0-sulfo-galactosylceramide), sulfate-containing lipids that occur throughout the body and are found in greatest abundance in nervous tissue, kidneys, and testes. Sulfatides are critical constituents in the nervous system, where they comprise approximately 5% of the myelin lipids. Sulfatide accumulation in the nervous system eventually leads to myelin breakdown (leukodystrophy) and a progressive neurologic disorder (Von Figura et al 2001). Nevertheless, this genotype does not cause MLD, and this benign condition of reduced ARSA activity is called ARSA pseudodeficiency. There are about 4 homozygotes in 1000 persons among non-Finnish Europeans (VarSome)
Pt. 1 is a carrier of a missense in gene CPN1 (table 2, first line) which leads to a loss of more than 60% of activity, according to a study on a single patient (Mathews KP. et al. 1980), (Cao H. et Hegele RA. 2003). The study of gene CPN1 in both patients (using IGV) has led to the identification of two rare variants (frequency less than 0.002) in intron 1 and 6 of one allele from Pt. 2 (table 1, figure 1). In MCAD no other damaging variations have been identified in these two patients by direct inspection with IGV (data not shown).
Table 2. Possible pathogenic variants found in exome from Pt. 1.Table 3. Possible pathogenic variants found in exome from Pt. 2.
Discussion
Whole exome sequencing (WES) is a technique that aims at the sequencing of the fraction of our genome that encodes for proteins: about 30 million base pairs (1% of the all the human DNA) divided into about 20 thousand genes (Ng SB et al. 2009). It has become increasingly clear that the use of WES can positively improve the rate of diagnosis and decrease the time needed for a definitive diagnosis in patients with rare genetic diseases (Sawyer SL et al. 2016). WES also positively impacts the ability to discover new pathogenic variants in known disease genes (Polychronakos C. et Seng KC. 2011) and the discovery of completely new disease genes (Boycott KM 2013). ME/CFS seems to have a genetic component: a US study found clear evidence of familial clustering and elevated risk for the disease among relatives of ME/CFS cases (Albright F et al. 2011) and several SNPs in various genes have been reported as more prevalent in ME/CFS patients versus healthy controls (Wang T et al. 2017). And yet, no studies that analyzed whole exomes of ME/CFS patients have been published, to my knowledge.
In this study, we searched for known genetic diseases in the exomes of two ME/CFS patients who fit the IOM criteria for SEID (IOM, 2015), with postural orthostatic tachycardia syndrome (POTS) identified by positive tilt table test. We detected a missense variant in CPN1 (rs61751507) in Pt. 1 (heterozygosis) that has been associated to a loss of activity of the enzyme of at least 60% in a previous study (Mathews KP. et al. 1980), (Cao H. et Hegele RA. 2003). We then found that, although Pt. 2 was not a carrier of this SNP, she had two rare SNPs in intron 1 (rs188667294) and 6 (rs113386068) of gene CPN1 (both present in less than 1/500 alleles, table 1, figure 1). These intronic variations have not been studied, to our knowledge, so their pathogenicity can’t be excluded at present. Variations in introns can be damaging just as missense and nonsense mutations in exons; suffice to say that the main known pathogenic SNP of gene CPN1 is a substitution in intron 1 (Cao H. et Hegele RA. 2003).
Carboxypeptidase N (CPN1) is an enzyme involved in the inactivation of C3a, C4a, C5a, and of kinins (bradykinin, kalladin) (Hugli T. 1978), (Plummer TH et Hurwitz MY 1978). In ME/CFS the physiologic increase in blood of C4a (the split product of the complement protein C4) after exercise is significantly more pronounced than in healthy controls as if there was a defect in C4a inactivation (Sorensen B et al. 2003). Such a defect could very well be a loss of function in CPN1, as found in Pt 1. Moreover, CPN1 is involved in inactivation of bradykinin, which is known to induce vasodilatation (Siltari A. et al. 2016), therefore CPN1 deficiency could play a role in POTS and in orthostatic intolerance in general. Both patients have a tilt table test positive for POTS. C4a has been recently considered to play a causal role in the cognitive deficit of schizophrenia, because of its role in synapsis pruning (Sekar, A et al, 2016); therefore a failure in its inactivation could be implicated in the incapacitating cognitive defects lamented by ME/CFS patients.
Only two patients with CPN1 deficiency have been reported so far in medical literature (Mathews KP. et al. 1980), (Willemse Jl et al. 2008), and the enzymatic defect has been associated to angioedema that most often involved the face and tongue, urticaria, and hay fever and asthma precipitated by exercise. This clinical presentation could be due, at least in part, to mast cell activation: in fact, C4a is a known anaphylatoxin that induces mast cells degranulation and release of histamine (Erdei A. et al. 2004). That said, we can observe that even if the clinical presentation of the only two known cases of CPN1 deficiency doesn’t fit the clinical picture of ME/CFS, mast cell activation syndrome (MCAS) has some commonalities with ME/CFS (Theoharides, TC et al. 2005), and mast cell abnormalities have been reported among ME/CFS patients (Nguyen T. et al. 2016). So we can’t exclude that activation of mast cells by a failure in C4a inactivation may lead to ME/CFS symptoms. The role of exercise as a trigger for symptoms in CPN1 deficiency is also highly suggestive because this is a pathognomonic feature of ME/CFS.
Conclusion
CPN1 deficiency is present (even if in a mild form) in Pt. 1, while Pt. 2 presents two rare intronic variants whose pathogenic role can’t be excluded. CPN1 deficiency could explain the abnormal increase of C4a after exercise and might be a contributing factor to post-exertional malaise and cognitive symptoms in ME/CFS. A search for pathogenetic SNPs in gene CPN1 among ME/CFS patients would clarify the role (if any) of this gene.
Acknowledgments. I would like to thank Chiara Scarpellini for her careful collection of annotations for each of the 2 hundred or so variants found by EVE within the exomes of Pt. 1 and Pt. 2 (table 2 and table 3).
In what follows I will present a set of codes that I have written in Octave and in Fortran for the study of random variables that follow a gamma distribution. In another post, we will discuss the beta distribution.
2_The gamma distribution
A random variable is said to have a gamma distribution with parameters (α, λ) if its density is given by
with α, λ>0 and where Γ(α) is the gamma function, which is quite a pain in the lower back, as you will see. This random variable can be used to describe the probability that the duration of life of an entity (whether it is a mechanical component or a bacterium) depends on time. The gamma distribution allows to describe the following situations:
for α>1 the probability of failure of a mechanical component (or death of an organism) increases as the time passes by (increasing hazard rate);
for α=1 the probability of failure is constant (constant hazard rate);
for α<1 the failure is less likely as time passes by (decreasing hazard rate).
The hazard rate is the limit below:
where F_X(t) is the distribution function of the gamma distribution. I have not been able to find an analytic expression for the hazard rate of the gamma distribution, but I have been able to obtain it with the computer-aided calculation (see below). But in order to calculate F_X(t), we need to calculate the gamma function…
3_The gamma function
Γ(α) (the gamma function) is an integral function that cannot be calculated analytically:
It has to be calculated with numerical integration, but a problem is that Γ(α) is an indefinite integral. One way to face this problem is to stop the numerical integration when the addend is below a certain cut off (let’s say 0.0001). It is worth noting that an important feature of this function is that if you know its value in an interval such as [x, x+1], then you can derive Γ(α) for any other possible value of α. This is due to the recursive formula: Γ(α+1) = αΓ(α). That said, I have calculated Γ(α) in [1, 2] and then I have used these values to define Γ(α) in [0, 4] with the .exe file (written in Fortran) available here. Download it and put it in a folder. When you run it, it gives a .txt file with the values of Γ(α) in [0, 4] and two .mat files that can be used to build the plot of Γ(α) (see figure 1). Values of Γ(α) in [0, 2] are tabulated below, while those of you who like the very long list of numbers can find a huge table with values of the gamma function in [0, 4] at the end of this article.
Table 1. Values of the gamma function in [0.1, 1.99], with a step of 0.01.Figure 1. The gamma function between 0 and 4, as calculated by my code in Fortran. Plotted by Octave.
You can find below the script of main_f_gamma.exe and its module.
!Data: 1/03/2018 PROGRAM main_f_gamma
USE DISLIN !libreria grafica USE mod_f_gamma !modulo con procedure
IMPLICIT NONE
REAL(KIND=num), DIMENSION (1:400):: f_gamma !valori della funzione gamma fra 1 e 4 REAL(KIND=num), DIMENSION (1:400):: alpha !valori della ascissa fra 1 e 4 CHARACTER(len=10):: chiusura !serve a chiudere il programma INTEGER:: i,j !indici di ciclo
CALL sub_f_gamma_1_2 (f_gamma,alpha) !calcola la funzione gamma fra 1 e 2
!calcolo della funzione gamma e della ascissa fra 2 e 4
ciclo: DO i=201,400,1 alpha(i) = 1. + alpha(i-100) f_gamma(i) = alpha(i-100)*f_gamma(i-100) END DO ciclo
!calcolo della funzione gamma e della ascissa fra 0 e 1
ciclo: DO i=2,100,1 alpha(i) = alpha(100+i)-1. f_gamma(i) = f_gamma(100+i)/alpha(i) END DO ciclo
!scrivo la funzione gamma tra 0 e 4 su un file di testo
WRITE(*,*) “To close this program press any key, then press RETURN.” READ(*,*) chiusura
END PROGRAM main_f_gamma
!Data: 1/03/2018
MODULE mod_f_gamma
!sezione dichiarativa
IMPLICIT NONE
!costanti
REAL,PARAMETER:: e = 2.71828182 !numero di Nepero REAL,PARAMETER:: d_alpha = 0.01 !incremento delle ascisse della funzione gamma INTEGER,PARAMETER:: num = SELECTED_REAL_KIND (p=13) !aumento il numero di cifre significative
CONTAINS
!————————————————————————————-
SUBROUTINE sub_f_gamma_1_2 (f_gamma, alpha)
!Esegue l’integrale necessario al calcolo della funzione gamma !per valori compresi fra 1 e 2 e registra questi valori in un array. !L’integrazione viene eseguita con il metodo di Cavalieri-Simpson
IMPLICIT NONE
!dichiaro gli argomenti fittizi
REAL(KIND=num),INTENT(OUT),DIMENSION(1:400):: f_gamma !valori della funz. gamma in [1,2] REAL(KIND=num),INTENT(OUT),DIMENSION(1:400):: alpha !contiene l’ascissa della funzione gamma
!dichiaro le variabili locali
REAL(KIND=num),DIMENSION(1:1000000):: funzione !è la funzione integranda REAL(KIND=num),DIMENSION(1:1000000):: primitiva !è l’integrale illimitato REAL(KIND=num),PARAMETER:: d_x = 0.001 !intervallo di integrazione INTEGER:: i,j !indice del ciclo INTEGER:: index = 1 !contiene l’estremo superiore della integrazione numerica
!sezione esecutiva
!assegno i valori di funzione e ascissa, eseguo l’integrazione numerica, assegno i valori di f_gamma_1_2
alpha(101) = 1. !assegna il primo valore di alpha alpha(102) = alpha(101)+d_alpha !assegna il primo valore di alpha f_gamma(101) = 1. !assegna il primo valore di gamma
ciclo_esterno: DO j=102,200,1 !assegna i valori a f_gamma_1_2
!assegna i valori a funzione
ciclo_funzione: DO i=1,1000000,1
funzione (i) = ( (d_x*REAL(i))**(alpha(j)-1.) )*( e**((-1)*d_x*REAL(i)) ) index = i IF (funzione(i)<=0.0001) EXIT ciclo_funzione
f_gamma(j) = primitiva(index) !assegna il valore j-mo a f_gamma_1_2
alpha (j+1) = alpha (j) + d_alpha
END DO ciclo_esterno
WRITE(*,*) alpha WRITE(*,*) f_gamma
RETURN
END SUBROUTINE sub_f_gamma_1_2
!————————————————————————————-
END MODULE mod_f_gamma
4_The distribution function
The distribution function of the gamma density cannot be calculated analytically so, again, I had to make good use of my poor coding skills and I wrote a code in Octave (here, it uses the two .mat files generated by the code in Fortran mentioned above!) and one in Fortran (available here). In figure 2 (code in Octave) you can find the density and the distribution function for different values of α.
Figure 2. Densities (left) and distribution functions (right) for α = 1.5, 2.0, 3.5 and λ = 1.5.
The code in Fortran, on the other hand, asks the user for the values of α and λ and gives a plot for the density and one for the distribution function. Moreover, it gives a .txt file with the values of the distribution function. Note that you can only place values for α that are below 25, otherwise you generate an overflow condition. I am not sure why you have overflow condition with a value of α that is below 34, but this is another story. Nevertheless, for those who are interested, Γ(n) = (n-1)! and the biggest factorial that can be handled by a computer in single precision mode (where 8 bits are used for the exponent) is 33!. So, we have a little mystery here… The reason why I have stuck with single precision is that DISLIN – the graphics library – seems to be unable to handle double precision numbers. If you are interested in the gamma function, the density and the distribution function for α = 20 and λ=2, then the code in Fortran gives for Γ(20) a value of 1.216451E+17 and the two plots in figure 3.
Figure 3. Density and distribution function for α = 20 and λ=2.
In what follows you find the scripts of the Fortran .exe mentioned in this paragraph.
PROGRAM main_legge_gamma USE DISLIN !libreria grafica USE mod_legge_gamma !modulo con procedure
IMPLICIT NONE
REAL(KIND=num):: alph, lamb !valore di alpha e di lambda REAL(KIND=num):: f_g !il valore della funzione gamma REAL(KIND=num):: estremo !estremo superiore di x REAL(KIND=num):: delta !intervallo di integrazione numerica REAL(KIND=num), DIMENSION (1:500000):: f_gamma !valori della funzione gamma fra 0 e 5000 REAL(KIND=num), DIMENSION (1:500000):: alpha !valori della ascissa fra 0 e 5000 REAL(KIND=num), DIMENSION (1:1000):: d_g !contiene la densita’ REAL(KIND=num), DIMENSION (1:1000):: f_r !funzione di ripartizione REAL(KIND=num), DIMENSION (1:1000):: x !ascissa CHARACTER(len=10):: chiusura !serve a chiudere il programma INTEGER:: i,j,k !indici di ciclo INTEGER:: ind !il valore di alpha in cui calcolo la funz. gamma
!calcolo della funzione gamma fra 1 e 2
CALL sub_f_gamma_1_2 (f_gamma,alpha) !calcola la funzione gamma fra 1 e 2
WRITE(*,*) “This program calculates the gamma distribution, for alpha and lambda specified” WRITE(*,*) “by the user. It plots the diagram of the density and distribution function, and” WRITE(*,*) “writes a file .txt with some values of the distribution function.” WRITE(*,*) ” “ WRITE(*,*) “Please, type the value of alpha (a positive number between 0.01 and 25.0).” WRITE(*,*) “Use dot as decimal separator.” READ(*,*) alph WRITE(*,*) ” “ WRITE(*,*) “Please, type the value of lambda (a positive number).” WRITE(*,*) “Use dot as decimal separator.” READ(*,*) lamb
!calcolo della funzione gamma e della ascissa fra 2 e 25
ciclo_2_5000: DO i=201,2500,1 alpha(i) = 1. + alpha(i-100) f_gamma(i) = alpha(i-100)*f_gamma(i-100) END DO ciclo_2_5000
!calcolo della funzione gamma e della ascissa fra 0 e 1
ciclo_0_1: DO i=2,100,1 alpha(i) = alpha(100+i)-1. f_gamma(i) = f_gamma(100+i)/alpha(i) END DO ciclo_0_1
!cerco il valore di alpha specificato dall’utente
ciclo_alph: DO i=1,2500,1 ind = i IF (ABS(alpha(i)-alph)<=0.01) EXIT ciclo_alph END DO ciclo_alph
!assegno il valore alla funzione gamma
f_g = f_gamma(ind)
WRITE(*,*) “The value of the gamma function in”, alph, “is”, f_g
!calcolo la densita’ gamma
IF (alph <= 1) THEN estremo = 2. delta = estremo/1000. ciclo_densita_1: DO k=1,1000,1 x(k)= delta*k d_g(k) = ( ( lamb**alph )/( f_g ) )*( x(k)**( alph – 1. ) )*( e**( (-1)*lamb*x(k) ) ) END DO ciclo_densita_1 ELSE estremo = ( (alph-1.)/lamb )*3. delta = estremo/1000. ciclo_densita_1: DO k=1,1000,1 x(k)= delta*(k-1.) d_g(k) = ( ( lamb**alph )/( f_g ) )*( x(k)**( alph – 1. ) )*( e**( (-1)*lamb*x(k) ) ) END DO ciclo_densita_1 END IF
WRITE(20,*) “The gamma distribution for alpha =”, alph, “and lambda =”, lamb WRITE(20,*) ” “ WRITE(20,*) “The value of the gamma function in”, alph, “is”, f_g WRITE(20,*) ” “ WRITE(20,*) “distribution function” WRITE(20,*) ” “ WRITE(20,*) ” x “,”|”,”distribution function” WRITE(20,*) “_________________________________”
ciclo_write: DO j=1,1000,1 WRITE(20,100) x(j),”|”,f_r(j) END DO ciclo_write 100 FORMAT (1X,1F8.4,A,1F8.4)
!disegno il grafico della densita’
CALL diagramma (x, d_g, ‘Density ‘)
!disegna il grafico della funzione di ripartizione
CALL diagramma (x, f_r, ‘distribution function’)
WRITE(*,*) “To close this program press any key, then press RETURN.” READ(*,*) chiusura
END PROGRAM main_legge_gamma
!Data: 1/03/2018
MODULE mod_legge_gamma
!sezione dichiarativa
IMPLICIT NONE
!costanti
REAL,PARAMETER:: e = 2.71828182 !numero di Nepero REAL,PARAMETER:: d_alpha = 0.01 !incremento delle ascisse della funzione gamma INTEGER,PARAMETER:: num = 1 !singola precisione
CONTAINS
!————————————————————————————-
SUBROUTINE sub_f_gamma_1_2 (f_gamma, alpha)
!Esegue l’integrale necessario al calcolo della funzione gamma !per valori compresi fra 1 e 2 e registra questi valori in un array. !L’integrazione viene eseguita con il metodo di Cavalieri-Simpson
IMPLICIT NONE
!dichiaro gli argomenti fittizi
REAL(KIND=num),INTENT(OUT),DIMENSION(:):: f_gamma !valori della funz. gamma in [1,2] REAL(KIND=num),INTENT(OUT),DIMENSION(:):: alpha !contiene l’ascissa della funzione gamma
!dichiaro le variabili locali
REAL(KIND=num),DIMENSION(1:1000000):: funzione !è la funzione integranda REAL(KIND=num),DIMENSION(1:1000000):: primitiva !è l’integrale illimitato REAL(KIND=num),PARAMETER:: d_x = 0.001 !intervallo di integrazione INTEGER:: i,j !indice del ciclo INTEGER:: index = 1 !contiene l’estremo superiore della integrazione numerica
!sezione esecutiva
!assegno i valori di funzione e ascissa, eseguo l’integrazione numerica, assegno i valori di f_gamma_1_2
alpha(101) = 1. !assegna il primo valore di alpha alpha(102) = alpha(101)+d_alpha !assegna il primo valore di alpha f_gamma(101) = 1. !assegna il primo valore di gamma
ciclo_esterno: DO j=102,200,1 !assegna i valori a f_gamma_1_2
!assegna i valori a funzione
ciclo_funzione: DO i=1,1000000,1
funzione (i) = ( (d_x*REAL(i))**(alpha(j)-1.) )*( e**((-1)*d_x*REAL(i)) ) index = i IF (funzione(i)<=0.0001) EXIT ciclo_funzione
f_gamma(j) = primitiva(index) !assegna il valore j-mo a f_gamma_1_2
alpha (j+1) = alpha (j) + d_alpha
END DO ciclo_esterno
RETURN
END SUBROUTINE sub_f_gamma_1_2
!————————————————————————————-
SUBROUTINE diagramma (x, funzione,label_y)
!disegno il grafico di funzione
!dichiaro gli argomenti fittizi
REAL(KIND=num),INTENT(IN),DIMENSION(1:1000):: x !ascissa REAL(KIND=num),INTENT(IN),DIMENSION(1:1000):: funzione !funzione CHARACTER (len=20), INTENT(IN):: label_y !nome delle ordinate
!dichiaro le variabili locali
REAL(KIND=num),DIMENSION(1:1000):: max_x !estremo ascissa REAL(KIND=num),DIMENSION(1:1000):: max_funz !estremo funzione
CALL METAFL (‘PDF’) !scelgo il formato dell’output CALL SETPAG (‘DA4P’) !formato della pagina CALL SCRMOD (‘revers’) !scritta nera su fondo bianco CALL DISINI !inizializza DISLIN impostando parametri di default
CALL PAGERA !chiedo una linea di contorno per il piano xy CALL DUPLX !chiedo un font a doppio spessore
CALL AXSPOS (300,1400) !coordinate dell’angolo in basso a sinistra CALL AXSLEN (1500,1000) !la lunghezza, in pixel, dei due assi
CALL LABDIG (2,’X’) !numero di cifre decimali negli assi CALL LABDIG (3,’Y’) !numero di cifre decimali negli assi
CALL NAME (‘X’,’x’) !assegno un nome all’asse delle x CALL NAME (label_y,’y’) !assegno un nome all’asse delle y
CALL GRAF (0.,max_x,0.,max_x/5.,0.,max_funz+max_funz/5.,0.,max_funz/5.) !ambiti degli assi CALL TEXMOD (‘ON’) !modalita’ TeX CALL GRID (2,2) !impone una griglia CALL CURVE (x,funzione,1000)
CALL DISFIN !termina DISLIN
END SUBROUTINE diagramma
!————————————————————————————-
END MODULE mod_legge_gamma
5_The hazard rate
As mentioned in paragraph 2, the hazard rate of the gamma distribution cannot be easily studied. I forgot to mention that it is relatively easy to prove that it is constant for α=1. But what happens for α>1 and for α<1? In order to answer, I have written a code (download) that plots it and it seems to be possible to say that the hazard rate increases for α>1 and decreases for α<1. This means that the gamma distribution is suitable for modeling each of the three kinds of patterns for the probability of failure as the time passes by. In figure 4 you find an example of the hazard rate for α<1 (left) and an example for α>1 (right).
Figure 4. Hazard rate for α = 0.8, λ = 1.5 (left) and α = 3.0, λ=1.5 (right).
This is the script of the Fortran exe program mentioned in this paragraph:
!Data: 7/03/2018
PROGRAM main_legge_gamma_avaria
USE DISLIN !libreria grafica USE mod_legge_gamma_avaria !modulo con procedure
IMPLICIT NONE
REAL(KIND=num):: alph, lamb !valore di alpha e di lambda REAL(KIND=num):: f_g !il valore della funzione gamma REAL(KIND=num):: estremo !estremo superiore di x REAL(KIND=num):: delta !intervallo di integrazione numerica REAL(KIND=num), DIMENSION (1:500000):: f_gamma !valori della funzione gamma fra 0 e 5000 REAL(KIND=num), DIMENSION (1:500000):: alpha !valori della ascissa fra 0 e 5000 REAL(KIND=num), DIMENSION (1:1000):: d_g !contiene la densita’ REAL(KIND=num), DIMENSION (1:1000):: f_r !funzione di ripartizione REAL(KIND=num), DIMENSION (1:1000):: x !ascissa REAL(KIND=num), DIMENSION (1:1000):: h !tasso di avaria CHARACTER(len=10):: chiusura !serve a chiudere il programma INTEGER:: i,j,k !indici di ciclo INTEGER:: ind !il valore di alpha in cui calcolo la funz. gamma
!calcolo della funzione gamma fra 1 e 2
CALL sub_f_gamma_1_2 (f_gamma,alpha) !calcola la funzione gamma fra 1 e 2
WRITE(*,*) “This program calculates the gamma distribution, for alpha and lambda specified” WRITE(*,*) “by the user. It plots the diagram of density, distribution function, failure rate” WRITE(*,*) “and writes a file .txt with some values of the distribution function.” WRITE(*,*) ” “ WRITE(*,*) “Please, type the value of alpha (a positive number between 0.01 and 25.0).” WRITE(*,*) “Use dot as decimal separator.” READ(*,*) alph WRITE(*,*) ” “ WRITE(*,*) “Please, type the value of lambda (a positive number).” WRITE(*,*) “Use dot as decimal separator.” READ(*,*) lamb
!calcolo della funzione gamma e della ascissa fra 2 e 25
ciclo_2_5000: DO i=201,2500,1 alpha(i) = 1. + alpha(i-100) f_gamma(i) = alpha(i-100)*f_gamma(i-100) END DO ciclo_2_5000
!calcolo della funzione gamma e della ascissa fra 0 e 1
ciclo_0_1: DO i=2,100,1 alpha(i) = alpha(100+i)-1. f_gamma(i) = f_gamma(100+i)/alpha(i) END DO ciclo_0_1
!cerco il valore di alpha specificato dall’utente
ciclo_alph: DO i=1,2500,1 ind = i IF (ABS(alpha(i)-alph)<=0.01) EXIT ciclo_alph END DO ciclo_alph
!assegno il valore alla funzione gamma
f_g = f_gamma(ind)
WRITE(*,*) “The value of the gamma function in”, alph, “is:”, f_g
!calcolo la densita’ gamma
IF (alph <= 1) THEN estremo = 2. delta = estremo/1000. ciclo_densita_1: DO k=1,1000,1 x(k)= delta*k d_g(k) = ( ( lamb**alph )/( f_g ) )*( x(k)**( alph – 1. ) )*( e**( (-1)*lamb*x(k) ) ) END DO ciclo_densita_1 ELSE estremo = ( (alph-1.)/lamb )*3. delta = estremo/1000. ciclo_densita_1: DO k=1,1000,1 x(k)= delta*(k-1.) d_g(k) = ( ( lamb**alph )/( f_g ) )*( x(k)**( alph – 1. ) )*( e**( (-1)*lamb*x(k) ) ) END DO ciclo_densita_1 END IF
WRITE(20,*) “The gamma distribution for alpha =”, alph, “and lambda =”, lamb WRITE(20,*) ” “ WRITE(20,*) “The value of the gamma function in”, alph, “is”, f_g WRITE(20,*) ” ” WRITE(20,*) “distribution function” WRITE(20,*) ” “ WRITE(20,*) ” x “,”|”,”distribution function” WRITE(20,*) “__________________________________”
ciclo_write: DO j=1,1000,1 WRITE(20,100) x(j),”|”,f_r(j) END DO ciclo_write 100 FORMAT (1X,1F8.4,A,1F8.4)
!calcolo il tasso di avaria
ciclo_avaria: DO k=1,1000,1 h(k) = d_g(k)/( 1. – f_r(k) ) END DO ciclo_avaria
!disegno il grafico della densita’
CALL diagramma (x, d_g, ‘Density ‘)
!disegna il grafico della funzione di ripartizione
CALL diagramma (x, f_r, ‘distribution function’)
!disegno il grafico del tasso di avaria
CALL diagramma (x, h, ‘Failure rate ‘)
WRITE(*,*) “To close this program press any key, then press RETURN.” READ(*,*) chiusura
END PROGRAM main_legge_gamma_avaria
!Data: 1/03/2018
MODULE mod_legge_gamma_avaria
!sezione dichiarativa
IMPLICIT NONE
!costanti
REAL,PARAMETER:: e = 2.71828182 !numero di Nepero REAL,PARAMETER:: d_alpha = 0.01 !incremento delle ascisse della funzione gamma INTEGER,PARAMETER:: num = 1 !aumento il numero di cifre significative
CONTAINS
!————————————————————————————-
SUBROUTINE sub_f_gamma_1_2 (f_gamma, alpha)
!Esegue l’integrale necessario al calcolo della funzione gamma !per valori compresi fra 1 e 2 e registra questi valori in un array. !L’integrazione viene eseguita con il metodo di Cavalieri-Simpson
IMPLICIT NONE
!dichiaro gli argomenti fittizi
REAL(KIND=num),INTENT(OUT),DIMENSION(:):: f_gamma !valori della funz. gamma in [1,2] REAL(KIND=num),INTENT(OUT),DIMENSION(:):: alpha !contiene l’ascissa della funzione gamma
!dichiaro le variabili locali
REAL(KIND=num),DIMENSION(1:1000000):: funzione !è la funzione integranda REAL(KIND=num),DIMENSION(1:1000000):: primitiva !è l’integrale illimitato REAL(KIND=num),PARAMETER:: d_x = 0.001 !intervallo di integrazione INTEGER:: i,j !indice del ciclo INTEGER:: index = 1 !contiene l’estremo superiore della integrazione numerica
!sezione esecutiva
!assegno i valori di funzione e ascissa, eseguo l’integrazione numerica, assegno i valori di f_gamma_1_2
alpha(101) = 1. !assegna il primo valore di alpha alpha(102) = alpha(101)+d_alpha !assegna il primo valore di alpha f_gamma(101) = 1. !assegna il primo valore di gamma
ciclo_esterno: DO j=102,200,1 !assegna i valori a f_gamma_1_2
!assegna i valori a funzione
ciclo_funzione: DO i=1,1000000,1
funzione (i) = ( (d_x*REAL(i))**(alpha(j)-1.) )*( e**((-1)*d_x*REAL(i)) ) index = i IF (funzione(i)<=0.0001) EXIT ciclo_funzione
f_gamma(j) = primitiva(index) !assegna il valore j-mo a f_gamma_1_2
alpha (j+1) = alpha (j) + d_alpha
END DO ciclo_esterno
RETURN
END SUBROUTINE sub_f_gamma_1_2
!————————————————————————————-
SUBROUTINE diagramma (x, funzione,label_y)
!disegno il grafico di funzione
!dichiaro gli argomenti fittizi
REAL(KIND=num),INTENT(IN),DIMENSION(1:1000):: x !ascissa REAL(KIND=num),INTENT(IN),DIMENSION(1:1000):: funzione !funzione CHARACTER (len=20), INTENT(IN):: label_y !nome delle ordinate
!dichiaro le variabili locali
REAL(KIND=num),DIMENSION(1:1000):: max_x !estremo ascissa REAL(KIND=num),DIMENSION(1:1000):: max_funz !estremo funzione
CALL METAFL (‘PDF’) !scelgo il formato dell’output CALL SETPAG (‘DA4P’) !formato della pagina CALL SCRMOD (‘revers’) !scritta nera su fondo bianco
CALL DISINI !inizializza DISLIN impostando parametri di default
CALL PAGERA !chiedo una linea di contorno per il piano xy CALL DUPLX !chiedo un font a doppio spessore
CALL AXSPOS (300,1400) !coordinate dell’angolo in basso a sinistra CALL AXSLEN (1500,1000) !la lunghezza, in pixel, dei due assi
CALL LABDIG (2,’X’) !numero di cifre decimali negli assi CALL LABDIG (3,’Y’) !numero di cifre decimali negli assi
CALL NAME (‘X’,’x’) !assegno un nome all’asse delle x CALL NAME (label_y,’y’) !assegno un nome all’asse delle y
CALL GRAF (0.,max_x,0.,max_x/5.,0.,max_funz+max_funz/5.,0.,max_funz/5.) !ambiti degli assi CALL TEXMOD (‘ON’) !modalita’ TeX per le formule CALL GRID (2,2) !impone una griglia CALL CURVE (x,funzione,1000)
CALL DISFIN !termina DISLIN
END SUBROUTINE diagramma
!————————————————————————————-
END MODULE mod_legge_gamma_avaria
A table for the gamma function in [0.01, 3.99] is reported below. In each column, the value of α is indicated on the left while the corresponding value for Γ(α) is indicated on the right.
I am currently trying to learn as much as I can about basic statistics; and one of the reasons for that is the need to gain access to the tools useful to better understand statistical analyses of medical data that I find in all the research papers on ME/CFS and related topics, which I read almost on a daily basis. Hopefully, I should also be able, at the end of this path, to perform my own statistical analyses of biological data.
When I study a subject where it is required the use of computer-aided calculation, I usually code my own software, even if there are programs ready to use for that purpose. The main reason for that is that in order to write a program that really works, you have to reach a complete understanding of the argument. So, it is a test in some way: if your program works properly, then you have learned something new. The other reason is that I like this activity, even though I can’t perform it for most of the time, because of the cruel encephalopathy due to my mysterious disease.
To give you an example of this attitude that I have, when I decided to study molecular mimicry between B-cell epitopes from a mathematical point of view, I found several programs for peptides alignments; nevertheless I decided to write my own software, and I spent several months of pure joy (but also of profound frustration because I was often too sick to write even a single line of code), learning about substitution matrices, dynamic programming and performing my own in silico experiments of cross-reactivity between enolases (a glycolitic enzyme evolutionary conserved across species) from different organisms. This work was then used in a paper recently published in the scientific journal Medical Hypotheses (Maccallini P et al. 2018). In figure 1 the flow diagram that I drew, while here and here you can download the main program and the subroutine, respectively, that I coded in Octave in order to perform comparisons between proteins. But this is another story… So, let’s come back to statistics.
Figure1. The flow diagram of my own software for the search of cross-reactive B cell epitopes in two proteins.
Many measurements of physical or biological phenomena obey to what is known as normal distribution: indicated with x the value of a specific measure, we have that values near the mean (also known as expectation and indicated with the letter μ) are more frequently assumed by x. In particular, if we indicate f(t) a function (named density) such that f(t)dt represents the probability that x takes a value between t and t+dt, then f(t) is the well known bell-shaped curve in figure 2 (on the left). As you can see, the higher the standard deviation σ is, the lower the bell: this means that standard deviation represents a measure of how much the values taken by x are spread on the x-axis.
Figure 2. Density (left) and distribution function (right) for some values of expectation and standard deviation for a normally distributed random variable.
Figure 2 has been obtained with a code that I have written in Octave, which can be downloaded here. So, what is the probability that x takes a value below a? This probability is given by the integral of f(t) between -∞ and a. This function is named distribution function, Φ(x). In figure 2, on the right, Φ(x) is plotted for different values of expectation and standard deviation. I have calculated these integrals using the numerical method known as Simpson’s rule. For μ=0 and σ=1 we have a particular type of normal distribution, usually called standard normal distribution. Some values of its distribution function are collected in the table in figure 3, that I have derived with this code in Octave. A table like this is available in any book of statistics, but I found useful to calculate my own table, for the reasons we have already mentioned.
Figure 3. The distribution function for a standard normal distribution.
Now, what about the distribution function for a normal non-standard distribution? Of course, you can calculate Φ(x) for any given values of μ and σ using available tools (like, for instance, function normcdf(x, m, s) of Octave). But, as always, I have coded my own software (in Fortran in this case). You can download the executable file here. Put it in a folder, then run it: it will ask you to type the values of μ, σ and it will give as output two diagrams (one for f, the other for Φ) and a .txt file which contains several values of Φ that can be used for applications. For μ=2, σ=0.5 we obtain the diagrams in figure 4, and a table whose first lines are in figure 5.
Figure 4. Density (left) and distribution function (right) for μ=2, σ=0.5 for a random variable with normal distribution.
A traslation of this blog post to Spanish can be downloaded here. I would like to thank Humbert.Cat for the translation.
1. Introduction
These are some notes about the talk that Mark Davis gave during the Community Symposium held in August at Stanford (video). I will introduce some basic notions about T cell receptors (TCR) in paragraphs 2, 3, 4, and 5. Paragraphs 6 is a description of an innovative technology developed by Mark Davis and his colleagues, based on information gathered from the video itself and three research papers published by Davis and others in the last 4 years. This background should be hopefully enough to allow a good understanding of the exciting pilot data presented by Mark Davis on T cell activity in ME/CFS (paragraph 7), and in chronic Lyme (paragraph 8), and to realize why this technology promises to be some sort of universal test for any kind of infectious and autoimmune diseases, known or unknown.
2. T cells
T cells are a type of leukocytes (also known as white blood cells), the cellular component of our immune system. Most of our circulating T cells are represented by T helper cells (Th cells) and cytotoxic T lymphocytes (CTL). While the function of Th cells is to regulate the activity of other leukocytes through the production of a wide range of chemicals (cytokines), CTLs are directly involved in the killing of host cells infected by pathogens. T cells belong to the adaptive arm of the immune system, along with B cells (the factories of antibodies), and as such, they are meant to provide a defence tailored to specific pathogens: our immune system can provide not only antibodies specific for a given pathogen but also specific T cells (both Th cells and CTLs). The innate arm of the immune system (which includes natural killer cells, macrophages, dendritic cells, mast cells…) on the other hand can provide only a one-fits-all type of defense, which represents the first line of immune response, during an infection.
3. T cell receptor
T cells search for their specific pathogens thanks to a molecule expressed on their surface, called T cell receptor (TCR). In figure 1 you can see a schematic representation of the TCR and of the mechanism by which T cells recognize their targets. Antigens (proteins) from pathogens are presented to T cells by other cells of our body: they are displayed on molecules called major histocompatibility complex (MHC), expressed on the outer membrane; if the antigen fits the TCR of a specific T cell, then this T cell is activated and proliferates (clonal expansion). The two chains (α and β) are assembled using the transcription of gene segments with several copies each: in other words, TCRs are assembled with peptides chosen randomly within a set of several possible choices. This leads to a repertoire of 10^15 possible different TCRs (Mason DA 1998). Each T cell displays only one type of TCR.
Figure 1.Upper half. Th cells and CTLs share the same TCR: in both cases this molecule is the assembly of two peptides (chain α and chain β), but while the TCR of Th cells (on the right) is expressed next to the molecule CD4 (which binds to class II MHC), the TCR of CTL is associated with the molecule CD8 (on the left), which is specific for MHC I. Black bars represent four chains (a complex called CD3) that are involved in the signaling of the TCR with the nucleus of the cell (by Paolo Maccallini). Lower half. A beautiful structural representation of the TCR, bound to the peptide-MHC complex (pMHC), from (Gonzàlez PA et al. 2013). In green the peptide, in blue the β chain, in dark green the α chain. CDRs (complementarity determining regions, orange) are composed of those residues of the α and β chains that directly bind the pMHC.
4. T helper cells
Th cells can recognize only antigens presented by class II MHC: this class of MHC is expressed on the outer membrane of some leukocytes, mainly dendritic cells, B cells, and macrophages (referred to as antigen presenting cells, APCs). MHC II engages the TCR of Th cells thanks to peptide CD4 (expressed exclusively by Th cells). The antigen presented by MHC II is a peptide with a length of 13-17 amino acids (Rudensky, et al., 1991) (figure 2).
Figure 2. The TCR expressed by a Th cell binds an epitope presented by a class II MHC expressed on the plasma membrane of an APC. Chains α and β of MHC II are also represented (by Paolo Maccallini).
5. Cytotoxic T lymphocytes
TCRs expressed by CTLs can bind only antigens displayed by class I MHC, which can be found on the outer membrane of any cell of our body. CD8 is the molecule that makes the TCR expressed by CTLs specific for MHC I. While antigens presented by APCs belongs to pathogens that have been collected on the battlefield of the infection, peptides displayed by class I MHC of a specific cell belong to pathogens that have entered the cell itself, therefore they are the proof of an ongoing intracellular infection (figure 3). When a CTL recognizes an antigen that fits its TCR, then the CTL induces apoptosis (programmed death) of the cell that displays it. Antigens presented by MHC I are peptides in the range of 8 to 10 amino acids (Stern, et al., 1994).
Figure 3. An infected cell displays a viral antigen on its MHC I. The TCR of a CTL binds this peptide and send a signal to the nucleus of the CTL itself, that responds with the induction of apoptosis (releasing granzymes, for instance) of the infected cell (by Paolo Maccallini).
6. The universal immune testing
In his talk, Mark Davis presents an overview of some basic concepts about the immune system, before introducing his exciting new data about ME/CFS and post-treatment Lyme disease syndrome (PTLDS, also known as chronic Lyme). But he also describes a few details of a complex new assay that is theoretically able to read all the information packed in the repertoire of TCRs present – in a given moment – in the blood of a human being. As such, this test – that I have named the universal immune testing – seems to have the potential to determine if a given patient has an ongoing infection (and the exact pathogen) or an autoimmune disease (and the exact autoantigen, i.e. the tissue attached by the immune system). To my understanding, this assay requires three steps, described in the following sections.
6.1. First step: TCR sequencing
As said in paragraph 3, when T cells encounter their specific peptide presented by MHC, they proliferate so that in blood of patients with an ongoing infection (or with a reaction against self, i.e. autoimmunity) we can find several copies of T cells expressing the same TCR: while in healthy controls about 10% of total CD8 T cells is represented by clones of a few different T cells (figure 4, first line), in early Lyme disease – an example of active infection – and in multiple sclerosis (MS) – an example of autoimmune disease – we have a massive clonation of a few lines of CTLs (figure 5, second and third line, respectively). The first step of the universal immune testing is represented by the identification of the exact sequence of TCRs expressed by T cells in blood, as reported in (Han A et al. 2014) where it is described how to sequence genes for the α and the β chain of a given T cell. This approach allows to build graphs of the kind in figure 4, and therefore to determine whether the patient has an abnormal ongoing T cell activity or not. If a clonal expansion is found, then we can speculate that either an active infection is present or some sort of autoimmune condition.
Figure 4. Each circle represents a patient. In the first line, we have four healthy controls, with no clonal expansion of CD8 T cells (the first one, left) or with only a low-level clonal expansion (slices in blue, white, and grey). In the second line, we have four patients with active Lyme disease (early Lyme) and all of them present a massive expansion of only three different T cells (slices in red, blue and orange). In the third line, we have four MS patient with most of their CD8 T cells represented by clones of a bunch of T cells. From the talk by Mark Davis.
6.2. Second step: TCR clustering
Mark Davis and his group have been able to code a software that allows to identify TCRs that share the same antigen, either within an individual or across different donors. This algorithm has been termed GLIPH (grouping of lymphocyte interaction by paratope hotspots) and has been found capable – for instance – to cluster T CD4 cell receptors from 22 subjects with latent M. tuberculosis infection into 16 distinct groups, each of which comprises TCRs from at least 3 different donors (Glanville J et al. 2017). Five of these groups are reported in figure 5. The idea here is that TCRs that belong to the same cluster, react to the same peptide-MHC complex (pMHC).
Figure 5. Five group of TCRs from 22 different donors with latent tuberculosis, clustered by GLIPH. The first column on the left has TCRs IDs, the second one reports donors IDs. Complementarity determining regions (CDR) for the β and the α chains are reported in the third and fifth column, respectively. From (Glanville J et al. 2017).
6.3. Third step: quest for the epitope(s)
As we have seen, this new technology allows to recognize if T cell clonal expansion is an issue in a given patient, by sequencing TCRs from his peripheral blood. It also allows to cluster TCRs either within an individual or across different patients. The next step is to identify what kind of antigen(s) each cluster of TCRs reacts to. In fact, if we could recognize these antigens in a group of patients with similar symptoms, with T cell clonal expansion and similar TCRs, we would be able to understand whether their immune system is fighting a pathogen (and to identify the pathogen) or if it is attacking host tissues and, if that was the case, to identify what tissue. As mentioned, the number of possible TCR gene rearrangement is supposed to be about 10^15, but the number of possible Th cell epitope is about 20^15 which is more than 10^19. This implies that TCRs have to be cross-reactive to some extent, in order to recognize all possible peptides presented by MHCs (Mason DA 1998). The exact extent of this cross-reactivity and the mechanism by which it is obtained has been elucidated by Mark Davis and his colleagues in a recent paper (Birnbaum ME et al. 2014) that gives us the third step of the universal immune testing. The aim of this phase is to take a given TCR and to find the repertoire of his specific antigens (as said, one TCR reacts to several antigens). In order to understand how this is possible let’s consider one of the experiments described in the paper mentioned above. The researchers considered two well-defined TCRs (named Ob.1A12 and Ob.2F3), cloned from a patient with MS and known to recognize peptide 85-99 (figure 6) of myelin basic protein (MBP) presented by HLA-DR15. They then prepared a set of yeast cells expressing HLA-DR15 molecules, each presenting a different peptide of 14 amino acids, with fixed residues only at position 1 and 4, where the peptide is anchored to MHC (figure 6, left). When copies of Ob.1A12 are added to this culture of yeast cells expressing pMHC complexes, they bind only some of them, and as you can see in the right half of figure 6, for each position of the epitopes bound by Ob.1A12, there is an amino acid that is more likely: for instance, the typical epitope of Ob.1A12 preferentially has alanine (A) at position -4, histidine (H) at position -3, arginine (R) at position -2, and so forth. As you can see, histidine (H) at position 2 and phenylalanine (F) at position 3 are obligate amino acids for a Ob.1A12 epitope.
Figure 6. On the left: peptide 85-99 of myelin basic protein (first row) is known to be an epitope for Ob.1A12. At position 1 and 4 it has two residues that allow its binding to the MHC molecule. At position -2, -1, 2, 3, and 5 we find those residues that bind the TCR. The second row represents the generic epitope of the peptide library used to identify the degree of crossreactivity between all the possible Ob.1A12 specific epitopes. On the right: the likelihood of amino acids for each position of Ob.1A12 specific epitopes is represented by shades of violet. As you can see, histidine (H) at position 2 and phenylalanine (F) at position 3 are obligate amino acids for a Ob.1A12 epitope. From (Birnbaum ME et al. 2014).
The table on the right side of figure 6 is, in fact, a substitution matrix with dimension 14×20, a tool that can be used to scan the peptide database in order to find, among all the known peptides expressed by living creatures, all the possible Ob.1A12 specific epitopes. Substitution matrices are commonly used for the so-called peptide alignment, a technique that aims at the identification of similarities between peptides. These matrices are based on evolutionary considerations (Dayhoff, et al., 1978) or on the study of conserved regions in proteins (Henikoff, et al., 1992). But the search for specific epitopes of a given TCR requires (as we have seen here for Ob.1A12) a substitution matrix built ad hoc for that TCR: each TCR requires its own substitution matrix that is obtained adding clones of that TCR on a culture of yeast cells presenting a huge peptide library on their MHCs, and analyzing data from this experiment. So, quite a complex process! In the case of Ob.1A12, this process led to 2330 peptides (including MBP), while the Ob.2F3 specific substitution matrix found 4824 epitopes within the whole peptide database. These peptides included both non-human proteins (bacterial, viral…) and human peptides. For 33 of them (26 non human and 7 human proteins), this group of researchers performed a test in order to directly verify the prediction: 25/26 of environmental peptides and 6/7 of the human peptides induced proliferation of T cells expressing Ob.1A12 and/or Ob.2F3, and this is a huge proof of the validity of this analysis! These 33 peptides are reported in figure 7. This is the last step of the universal immune testing, the one that from the TCR leads to the epitopes. As you have seen, a huge set of different peptides from different sources is linked to each single TCR, in other words, crossreactivity is an intrinsic property of TCR. This also means that lymphocyte transformation tests (LTTs), widely used in Europe for the detection of infections like Borrelia burgdorferi and others, bear a high risk of false-positive results and require a process of experimental and theoretical validation, that is currently lacking (see also this post on this issue).
Figure 7. 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).
We are now ready to fully appreciate the pilot data that Mark Davis presented at the Symposium on CD8 T cell clonal expansion in ME/CFS and in chronic Lyme.
7. We have a hit!
Mark Davis, along with Jacob Glanville and José Montoya, has sequenced TCRs from the peripheral blood of 50 ME/CFS patients and 49 controls (first step of the universal immune testing, remember?), then they have clustered them using the GLIPH algorithm (second step). They have found 28 clusters with more than 2500 similar sequences each, where each cluster collects multiple sequences from the same individual as well as (which is perhaps more important) sequences from different patients (figure 8). The cluster that I have circled in red, for instance, is a collection of more than 3000 similar TCRs. The presence of this wide clusters in ME/CFS patients, compared to healthy controls, represents an indirect proof of a specific T cell response to some common trigger in this group of patients, which might be a pathogen or a tissue of the body (or both).
Figure 8. In ME/CFS, TCRs sequences from 50 patients form 28 clusters with more than 2500 similar sequences, and this is clearly not the case in healthy controls. This point to some specific immune response to a pathogen or to a human tissue (or both). This slide is from the talk by Mark Davis.
Among these 50 ME/CFS patients, Davis and colleagues selected 6 patients with similar HLA genes (figure 9, left), 5 females among them, and studied their TCRs deeper. In the right half of figure 9, you can see for each patient the degree of CTL clonal expansion. Remember that in healthy controls only about 10% of CTLs is composed by clones of a few cells (figure 4, first raw), while here we see that about 50% of all CTLs is composed by clones. So, a “marked clonal expansion” of CD8 T cells, as Davis said.
Figure 9. On the left: 6 MECFS patients with similar HLA genes have been selected. Patient ID is reported in the first column on the left, the second column indicates the age of each patient, the third indicates the gender, the fourth column is about exposure to cytomegalovirus, the third one is on MHC I genes. On the right: analysis of clonal expansion of CD8 T cells for each of the six patients. There is a marked clonal expansion (about 50%) compared to healthy controls (about 10%).
Sequences of α and β chains of TCRs from three of the six patients (patients L4-02, L4-10, and L3-20) are reported in figure 10 (I have verified that in fact these are sequences of α and β chains of human TCRs using them as query sequences in standard protein BLAST).
Figure 10. Beta chains (first column) and respective α chains (fifth column) from 3 ME/CFS patients (L4-02, L4-10, and L3-20, last column).
So, we have seen so far the first two steps of the universal immune testing in ME. What about the third step? In his talk, Mark Davis didn’t present any particular epitope, he just showed a slide with what likely is the selection of the epitopes from the peptide library (see paragraph 6.3) by one of the TCRs reported in figure 10. This selection is reported in figure 11, but from that picture, it is not possible to gather any information about the identity of these epitopes. As you probably remember from paragraph 6.3, the analysis of the peptides selected by a TCR among the peptide library allows the identification of a substitution matrix that can be used to select all the possible epitopes of that specific TCR, from the peptide database. This last crucial step has to be performed yet, or it has been already performed, but Davis has not communicated the preliminary results during his talk. Recently new resources have been made available by Open Medicine Foundation, for this promising research to be further pursued, among other projects (R). The aim here, as already said, is to find the antigen that triggers this T cell response. As Mark Davis said, it might be an antigen from a specific pathogen (perhaps a common pathogen that comes and goes) that elicits an abnormal immune response which ends targeting some host tissue (microglia, for instance), thus leading to the kind of immune activation that has been recently reported by Mark Davis himself and others in ME/CFS (Montoya JG et al. 2017). The idea of a common pathogen triggering a pathologic immune response is not new in medicine, and rheumatic fever (RF) is an example of such a disease: RF is an autoimmune disease that attacks heart, brain and joints and is generally triggered by a streptococcal throat infection (Marijon E et al. 2012). The other possible avenue is, of course, that of an ongoing infection of some kind, that has yet to be detected. As said (see par. 6.1), CD8 T cell clonal expansion is present in both acute infections (like early Lyme disease) and autoimmune diseases (like MS) (figure 4), so we have to wait for the antigen identification if we want to understand if the CTLs activity is against a pathogen and/or against a host tissue.
Figure 11. In this picture, we can see the selection, through several rounds, of a bunch of peptides by a particular TCR from a ME patient. The selection takes place among a huge collection of peptides presented by HLA-A2 (MHC I) expressed by yeast cells. At each round the number of possible peptides is smaller.
8. Chronic Lyme does exist
It has probably been overlooked that in his talk, Mark Davis reported also very interesting data on post-treatment Lyme disease syndrome (PTLDS, also known as chronic Lyme disease). In particular, he found a marked clonal expansion in CD8 T cells of 4 PTLDS patients (about 40% of total CTLs) as reported in figure 12: consider that in this case, blue slices represent unique T cells, while all the other slices represent clones! All that has been said about CD8 clonal expansion in ME/CFS does apply in this case too: it might be the proof of an ongoing infection – perhaps the same B. burgdorferi, as suggested by several animal models (Embers ME et al. 2017), (Embers ME et al. 2012), (Hodzic E et al. 2008), (Yrjänäinen H et al. 2010) – or a coinfection (a virus?) or it could be the expression of an autoimmune reaction triggered by the initial infection. This has still to be discovered, running the complete universal immune testing, but what is already clear from figure 12 is that PTLDS is a real condition, with something really wrong going on within the immune response: chronic Lyme does exist.
Figure 12. CD8 T cells clonal expansion in four chronic Lyme patients: there is a marked clonal expansion that stands for T cell activity against a pathogen or against host tissue.
9. Conclusions
Mark Davis and other researchers have developed a complex assay that is able to sequence TCRs from patients, cluster them into groups of TCRs that react to the same antigens, and discover the antigens that triggered that particular T cell response. This assay is a kind of universal immune testing that is theoretically able to recognize if a person (or a group of patients) presents an immune response against a pathogen or against one of his own tissues (or both). This approach has already given pilot data on an ongoing CD8 T cell activity in ME/CFS patients and in chronic Lyme patients and will hopefully identify the trigger of this immune response in the near future. Whether ME/CFS is an ongoing infection, an autoimmune disease or both, the universal immune testing might be able to tell us. This new technology is for immunology, what whole genome sequencing is for genetics, or metabolomics is for molecular diseases: it doesn’t search for a particular pathogen or a particular autoimmune disease. No, it searches for all possible infections and immune disorders, even those that have yet to be discovered.
Le misure metaboliche in vitro dello studio norvegese in cui fu evidenziato un possibile blocco al livello del piruvato deidrogenasi nei pazienti ME/CFS, sono state effettuate con il dispositivo Seahorse XFe96 della Agilent. Questo apparecchio (grande come una stampante da tavolo) permette di misurare in tempo reale – in vitro – il metabolismo energetico di cellule prelevate da pazienti (ad esempio linfociti). Come è spiegato nel video che segue, il tutto si riduce a due misure:
una misura del consumo di ossigeno, che fornisce una stima del funzionamento mitocondriale;
una misura della concentrazione di protoni, che fornisce una stima del funzionamento della glicolisi.
Mi risulta che l’Università degli Studi di Firenze sia in possesso di questo apparecchio (vedi qui).
Il Seahorse è attualmente impiegato nello studio di Avindra Nath (NIH) su 40 pazienti con ME/CFS post-infettiva (PI-ME/CFS). Come gruppo di controllo per questa ricerca, oltre a 20 persone sane, sono state selezionate anche 20 persone che hanno avuto la Lyme e sono completamente guarite.
In tre miei articoli precedenti (qui, qui e qui) – a commento del recente lavoro pubblicato da Fluge e colleghi (Fluge O et al. 2016) – ho descritto uno scenario in cui una funzione ridotta dell’enzima piruvato deidrogenasi (PDH) nei pazienti ME/CFS porta a una inefficiente sintesi di energia in seno al ciclo del TCA. La ridotta funzione del PDH è stata dedotta dalla presenza di un fenomeno di catabolismo degli amminoacidi, e dalla sovra espressione degli enzimi piruvato deidrogenasi kinasi (PDK), in particolare le isoforme 1, 2, 4. Ora la domanda è: cosa causa questa alterazione metabolica? Gli Autori, sulla scorta del loro successo terapeutico con il Rituximab, ipotizzano che un autoanticorpo possa – in alcuni pazienti – attivare o disattivare dei circuiti legati alla regolazione del metabolismo energetico. In quanto segue propongo uno scenario alternativo, basato su uno studio su topi con l’influenza.
Primo atto: influenza A e piruvato deidrogenasi
Nel 2014 un gruppo giapponese (Yamane K et al. 2014) ha inoculato il virus della influenza A (IAV) in topi da laboratorio, e nei 7 giorni successivi ha condotto uno studio sulle cavie malcapitate, simile a quello eseguito da Fluge e Mella sui pazienti ME/CFS, se non per il fatto che i topi sono stati sacrificati in modo da poter effettuare le misure direttamente nei tessuti. Come si vede in figura 1.A, l’attività del piruvato deidrogenasi si riduce col passare dei giorni nei vari tessuti esaminati, con la sola eccezione del cervello. Contestualmente (figura 1.B) anche il livello di ATP scende (ovunque, tranne che nel cervello). Questa prima parte dell’esperimento si può considerare equivalente alla prima parte dello studio di Fluge e Mella, quella che ho discusso qui. Cambia il tipo di misure effettuate, ma il risultato è lo stesso: il metabolismo energetico è depresso e si ha una perdita di attività del piruvato deidrogenasi.
Figura 1. Attività dell’enzima piruvato deidrogenasi in vari tessuti (A) e concentrazione di ATP nei medesimi tessuti (B).
Secondo atto: piruvato deidrogenasi kinasi, il solito sospetto
Esattamente come Fluge e Mella, anche i ricercatori giapponesi si sono chiesti se una espressione genica insolitamente alta degli enzimi piruvato deidrogenasi kinasi (ce ne sono quattro, indicati PDK1, PDK2…) potesse essere responsabile della ridotta attività del piruvato deidrogenasi. Infatti questi quattro enzimi hanno proprio la funzione di inibire il pirvato deidrogenasi. Come si vede in figura 2, il PDK4 aumenta rapidamente col passare dei giorni sia nel cuore, che nei polmoni, così come nel fegato e nei muscoli scheletrici.
Figure 2. Espressione del PDK4 in vari tessuti, in funzione dei giorni contati a partire dal momento in cui è stato inoculato il virus dell’influenza.
Questo secondo esperimento è simile agli esperimenti sulla espressione genica nelle cellule mononucleari del sangue periferico, effettuati dal gruppo di Fluge e Mella sui pazienti ME/CFS (qui). Anche in quel caso si è trovata una sovra espressione del PDK4, ma anche del PDK1 e 2, che invece nei topi sono normali. Comunque esiste una sovrapponibilità fra i due risultati.
Uomini e topi
In base a quanto visto, durante i primi 7 giorni dalla inoculazione del virus della influenza A, i topi cominciano a sviluppare una disfunzione metabolica simile a quella descritta da Fluge e Mella nei pazienti ME/CFS: un aumento del piruvato deidrogenasi kinasi si associa a una perdita di funzione del piruvato deidrogenasi e a un complessivo decadimento del metabolismo energetico. Questo cosa significa? E’ difficile trarre conclusioni, ma potremmo forse azzardare l’ipotesi che:
la alterazione metabolica descritta nella ME/CFS altro non è che quella che si verifica durante una infezione, a partire già dal primo giorno.
Poichè i primi anticorpi (classe IgM) si formano solo una o due settimane dopo l’inizio della infezione, possiamo escludere che le alterazioni osservate nei topi siano da imputare agli anticorpi. Gli stessi autori le attribuiscono a varie citochine (vei figura 3).
Figura 3. Il virus della influenza induce la sintesi di citochine che, a loro volta, attivano la sovra espressione di PDK che inibisce il piruvato deidrogenasi.
Questo significa che una possibile ipotesi per il difetto del piruvato deidrogenasi nella ME/CFS può essere semplicemente la presenza di una infezione cronica, in accordo con quanto suggerito da Antony Komaroff, fra gli altri (vedi qui). Ovviamente questa è solo una fra le tante ipotesi possibili.
E il rituximab?
Se gli anticorpi non c’entrano, allora perché il farmaco rituximab – che uccide le cellule B che esprimono il CD20 – ha un effetto terapeutico in più di metà dei pazienti ME/CFS? Questa è un’ottima domanda, se si potesse rispondere alla quale saremmo più vicini alla soluzione. Tuttavia si consideri che le cellule B non sono solo fabbriche di anticorpi, ma sono anche prensentatori di antigeni, produttrici di citochine (Frances E. Lund 2009) e rilasciano DNA mitocondriale (dati non pubblicati, anticipati da Anders Rosén in questo video, minuto 14:16) esattamente come i mastociti (Zhang et al. 2012). Si ritiene che il DNA mtocondriale sia fortemente infiammatorio (infatti assomiglia a quello di un batterio) e quindi potrebbe essere la causa di diversi disturbi (Zhang et al. 2012), tra cui anche magari la disregolazione del piruvato deidrogenasi. Quindi l’effetto del rituximab nella ME/CFS non è necessariamente legato alla presenza di autoanticorpi.
Conclusione
Abbiamo visto che l’alterazione metabolica recentemente ipotizzata nella ME/CFS (Fluge O et al. 2016) è presente anche durante i primi 7 giorni di una infezione virale. Quindi il detto comune secondo il quale “la CFS è come una influenza che non passa mai” sembra corretto anche da un punto di vista metabolico, e potrebbe avere un potere descrittivo ben più profondo di quanto si sia potuto immaginare fino ad oggi.
Donazione
Considera una donazione per sostenere questo blog.
Per una lettura veloce andare direttamente ai paragrafi 7 e 8.
1.Introduzione
Nei precedenti due articoli (qui e qui) sul lavoro pubblicato da Fluge, Mella e collaboratori (Fluge et al. 2016) abbiamo visto che:
qualcosa induce l’espressione di una serie di meccanismi che riducono la funzione del piruvato deidrogenasi, costringendo i pazienti ME/CFS a bruciare amminoacidi al posto dello zucchero. Ma dei sistemi di compenso intervengono per cercare di riportare il metabolismo energetico alla normalità, senza riuscirci. I sistemi di compenso sono diversi fra maschi e femmine, ma il difetto metabolico a monte è lo stesso nei due sessi.
Allo scopo di individuare la molecola (o le molecole) a cui ricondurre l’origine di questa disfunzione, il gruppo di scienziati ha preparato delle colture di cellule umane provenienti da muscoli scheletrici, e le ha esposte al siero di 12 pazienti ME/CFS e di altrettanti controlli sani. Le colture cellulari sono state poi sotToposte a due tipi di misure:
il livello di consumo di ossigeno (OCR, oxygen consumption rate), che indica il livello di attività dei mitocondri;
il livello di acidificazione dello spazio extracellulare (ECAR, extracellular acidification rate) che è una misura surrogata del livello di acido lattico prodotto.
Entrambi i parametri sono stati misurati sia in presenza di glucosio, che in presenza di amminoacidi. Nel complesso sono state fatte misure in quattro condizioni sperimentali.
2.Primi due esperimenti: amminoacidi e glucosio
Le cellule muscolari a riposo – esposte al siero dei pazienti – e coltivate in presenza di amminoacidi, presentano un consumo di ossigeno maggiore delle cellule esposte al siero di controlli sani (figura1, B.I). Aggiungendo il glucosio (figura 1, B.II), la situazione non cambia: ancora la respirazione delle cellule muscolari esposte al siero dei pazienti è maggiore di quella delle cellule muscolari esposte al siero di controlli sani. La sintesi di acido lattico non differisce fra le cellule esposte a siero di pazienti e cellule esposte a siero di controlli sani (figura 1, D.I e D.II).
Figura 1. Consumo di ossigeno e nei quattro esperimenti (B) e sintesi di acido lattico (D).
3.Terzo esperimento: blocco dell’enzima ATP sintasi
Il blocco dell’enzima ATP sintasi, il quale si occupa di convertire ADP in ATP alla fine della catena respiratoria, riduce drammaticamente il consumo di ossigeno in entrambe le colture, ma quella con siero di pazienti è meno colpita (figura 1, B.III). Vale la pena ricordare che la subunità beta dell’enzima ATP sintasi è sovra espressa nei pazienti ME/CFS (vedi qui), e questo potrebbe essere legato all’effetto benefico del siero dei pazienti in questo esperimento. La produzione di lattato in questo esperimento è maggiore nelle cellule esposte a siero di pazienti (figura 1, D.III).
4.Quarto esperimento: blocco della catena respiratoria
I ricercatori hanno esposto le colture cellulari a una molecola (la CCCP) che inibisce la catena respiratoria. Come forma di compenso le colture cellulari aumentano drasticamente il consumo di ossigeno ma la coltura esposta al siero ME/CFS consuma più ossigeno (figura 1, B.IV) e produce più lattato (figura 1, D.IV).
Figura 2. Sottrazioni fra misure effettuate negli esperimenti su coltura cellulare.
5.Sottrazioni
I ricercatori hanno poi effettuato i seguenti calcloli, i cui risltati sono riportati in figura 2:
consumo di ossigeno dell’esperimento II meno quello dell’esperimento III (figura 2, E);
consumo di ossigeno dell’esperimento IV meno quello dell’esperimento III (figura 2, E);
produzione di acido lattico dell’esperimento II meno quello dell’esperimento I (figura 2, F);
produzione di acido lattico dell’esperimento III meno quello dell’esperimento II (figura 2, F);
produzione di acido lattico dell’esperimento IV meno quello dell’esperimento II (figura 2, F).
6.Cosa ci dicono questi esperimenti?
Non è immediato dedurre il significato di questi esperimenti, almeno per me. Tuttavia possiamo dire quanto segue.
Il consumo di ossigeno delle cellule muscolari esposte al siero dei pazienti è maggiore di quello del gruppo di controllo (esperimenti I-IV), e questo è apparentemente in disaccordo con quanto risulta dai test ergosirometrici nella ME/CFS, in cui il consumo di ossigeno sistemico – per watt erogato – è minore nei pazienti (Vanness, 2007), (Snell, 2013).
Stressando chimicamente la catena respiratoria (esperimenti III e IV) la produzione di acido lattico aumenta nella coltura esposta a siero di pazienti più di quanto non aumenti nel contollo. Forse questo è il dato più interessante, che potrebbe rispecchiare un blocco nel piruvato deidrogenasi indotto dal siero dei pazienti nella coltura cellulare.
7.Conclusioni
Gli esperimenti in vitro di Fluge e colleghi dimostrano che nel siero dei pazienti ME/CFS è presente un fattore X (non noto) che:
aumenta la capacità delle cellule di consumare ossigeno;
aumenta la produzione di lattato in condizioni di aumentato fabbisogno energetico (simulate in provetta con inibitori della catena respiratoria).
La prima osservazione può indicare – secondo gli autori – la presenza di meccanismi di compenso attivati da messaggeri chimici contenuti nel sangue dei pazienti, che potenziano l’attività mitocondriale. La seconda osservazione sperimentale è in accordo con lo studio di Armstrong, che non rileva un aumento del lattato basale (Armstrong CW et al. 2015), e con osservazioni del gruppo norvegese (non ancora pubblicate) che indicano un aumento significatico di lattato dopo esercizio, rispetto ai controlli sani. L’aumento di lattato a seguito di esercizio è ciò che ci si aspetterebbe in presenza di un blocco del piruvato deidrogenasi. Quindi il fattore X (o i fattori?) contenuto nel siero dei pazienti è responsabile di due azioni, apparentemente opposte: da un lato potenzia i mitocondri, dall’altro fa aumentare la produzione di lattato. Questo potrebbe voler dire che:
nel siero dei pazienti è presente sia la causa dell’inibizione del piruvato deidrogenasi, che un fattore di compenso, che cerca di porre rimedio al difetto metabolico.
Questa è solo una delle interpretazioni possibili, ovviamente. Aggiungo che un aumento significativo del lattato dopo attività fisiche blande era stato segnalato in un precedente studio su un solo paziente, realizzato dal paziente stesso (Mark Vink 2015).
8.Un possibile test
Come abbiamo visto, nella ME/CFS il lattato basale sembra normale (Armstrong CW et al. 2015) e lo studio di Fluge e Mella – qui discusso – conferma questa dinamica anche con esperimenti in vitro. Tuttavia, a seguito di stress energetici, la produzione di lattato aumenta più di quanto ci si aspetterebbe normalmente, come già dimostrato da un paziente/ricercatore (Mark Vink 2015) e come evidenziato anche nello studio di Fluge e Mella, in vitro. Allora ho pensato che un possibile test in grado di rilevare questa dinamica metabolica potrebbe essere la misura di lattato e ammonio nel test da sforzo ischemico dell’avambraccio. In questo test viene fatta contrarre ripetutamente una mano del paziente (con una pallina morbida), essendo il flusso sanguigno bloccato con laccio emostatico. Dopo l’esercizio (che dura un minuto), diversi prelievi venosi vengono fatti nei successivi 10 minuti, con il laccio ancora stretto. Questi prelievi forniscono la produzione locale di lattato da parte dei muscoli scheletrici (Livingstone C et al. 2001). In figura 3 potete vedere l’esame di un paziente: il livello basale di lattato è perfettamente nella media, ma dopo 3 minuti la curva (in rosso) si discosta dal valore medio e sale oltre il valore massimo (anche se di poco). Questo tipo di andamento è coerente con quanto evidenziato in vitro da Fluge e Mella, e con la loro ipotesi sul piruvato deidrogenasi.
Figura 3. Curva lattato, dopo sforzo ischemico dell’avambraccio. In rosso il livello di acido lattico del paziente, che da un valore perfettamente normale, sale superando il valore massimo (anche se di poco).
In questo test si misura anche l’ammonio, come verifica del fatto che il paziente si sia ‘impegnato’ nell’esecuzione dello sforzo. Infatti l’ammonio è un prodotto del metabolismo energetico (specialmente anaerobico), secondo le seguenti reazioni:
2ADP –> ATP + AMP
AMP –> IMP + NH3
9.Il cavaluccio marino
Le misure di cui ho parlato in questo articolo sono state effettuate con il dispositivo Seahorse XFe96 della Agilent. Questo apparecchio permette di misurare in tempo reale il metabolismo energetico cellulare (ad esempio di linfociti) attraverso una misura del consumo di ossigeno (che fornisce una stima semplice del funzionamento mitocondriale) e la produzione di protoni (che si può ritenere una misura della glicolisi), in varie condizioni sprimentali. Video esplicativo.
There is a European study that has checked for the expression of 1007 mitochondrial proteins in platelets from 2 twins, one with ME/CFS and the other one healthy (Ciregia F et al 2016). Of these proteins, 194 were significantly modified in the sick twin, in comparison with the healthy one. I have checked for differences in pyruvate dehydrogenase complex, ADP/ATP translocase subunits and pyruvate dehydrogenase kinases. This is what I have found in these two twins:
1) Pyruvate dehydrogenase E1 subunits alpha (PDHA) and beta are both increased in the sick twin, which is in partial agreement with the increase in PDHA found by Fluge and Mella (Fluge et al. 2016);
2) ADP/ATP translocase, subunit 2 and 3, are low in the sick twin, compared with the healthy one, which could be in accordance with the study by Myhill and colleagues (Myhill S et al. 2009), (Booth, N et al 2012), if only we assume that the problem with this enzyme found in group defined “HiBlk” is not due to blockage from a molecule or an autoantibody, but is instead due to underexpression of the enzyme itself;
3) Pyruvate dehydrogenase kinases 1 and 3 are overexpressed, which again is in partial agreement with what Fluge and Mella have found in their recent paper, where PDK 1, 2, 4 are overexpressed in ME/CFS patients.
In conclusion, the sick twin does not seem to have any blockage of ADP/ATP translocase, because if that was the case he would present an overexpression of the enzyme, while the enzyme is underexpressed. On the other hand, he does seem to have a problem with his pyruvate dehydrogenase, in fact there is inhibition by overexpressed PDK 1 and 3 and – at the same time – he is expressing more PDHA than his healthy twin.
The same European study (Ciregia F et al 2016) then selected three enzymes from the 194 significantly modified in the sick twin: ACON, ATPB e MDHM. They then evaluated the expression of these enzymes in a cohort of 45 Italian patients with ME/CFS and in 45 matched controls. In this case, they considered mitochondria from saliva. They found that both ACON and ATPB are overexpressed in patients.
ACON stands for aconitase, which is an enzyme of the TCA cycle, that catalyzes the step from citrate to cis-aconitate (see figure). Thus, its over-expression in this cohort of patients is in agreement with the depletion of these two metabolites, found in a recent Japanese study (Yamano E, et al. 2016) and with the study by Fluge and Mella: in fact, if we assume that the TCA cycle is poorly supplied by glycolysis, it would overexpress one or more enzymes in order to increase the energy production from the substrate available. Thus, this finding seems in agreement with both the Norwegian and the Japanese study and seems to complete the picture.
In this picture, I have put together the levels of metabolites of the TCA from the Japanese study and the picture used by Fluge and Mella in their metabolic study. I have indicated the function of aconitase in TCA cycle.
ATPB is subunit beta of ATP synthase, and it is involved in the last step of mitochondrial metabolism, the conversion of ADP into ATP. Again, an overexpression of this enzyme seems to be in agreement with poor energy supply.
A discussion on these observations can be found in this thread in Phoenix Rising.
Nel mio precedente articolo sullo studio pubblicato dal gruppo di Fluge e Mella, mi sono limitato a riportare i loro risultati sulla analisi dei 20 amminoacidi standard nel sangue dei pazienti ME/CFS. Riassumendo:
nelle donne con ME/CFS gli amminoacidi che si trasformano in acetil-CoA e quelli che alimentano direttamente il ciclo dell’acido citrico (ciclo TCA) sono ridotti;
nei maschi non si ha questa riduzione degli amminoacidi nel sangue, ma si ha un aumento della 3-metilistidina (3-MHis), che suggerisce che i muscoli vengono smembrati (catabolizzati) per alimentare il ciclo del TCA. Nelle donne questo fenomeno è assente.
Gli autori hanno formulato allora la seguente ipotesi:
i pazienti ME/CFS (sia maschi che femmine) utilizzano gli amminoacidi per alimentare i mitocondri (il ciclo del TCA avviene dentro i mitocondri) più di quanto facciano i controlli sani.
I mitocondri – in condizioni normali – usano come fonte principale di alimentazione lo zucchero. Lo zucchero – trasformato in piruvato fuori dai mitocondri (nella glicolisi) – viene ulteriormente trasformato in acetil-CoA dal PDH, e quindi alimenta il ciclo del TCA (vedi figura 1).
Figura 1. Rappresentazione semplificata del metabolismo energetico. Lo zucchero è trasformato in piruvato nella glicolisi, il piruvato è trasformato in acetil-CoA dal piruvato deidrogenasi (PDH). L’acetil-CoA entra nei mitocondri dove alimenta il ciclo del TCA. Si vedono anche i vari punti di ingresso degli amminoacidi nei percorsi metabolici.
Ciò posto, una possibile spiegazione per l’aumentato consumo degli amminoacidi da parte dei pazienti ME/CFS è che:
l’attività dell’enzima piruvato deidrogenasi (PDH) sia ridotta – per qualche motivo – e quindi l’approvigionamento di acetil-CoA è deficitario. In questo contesto, i mitocondri bruciano gli amminoacidi per produre energia, al posto dell zucchero.
Espressione genica
Per verificare l’ipotesi di un ridotto funzionamento del piruvato deidrogenasi (PDH), gli autori sono passati a un secondo gruppo di analisi, che sono l’oggetto del presente articolo. In particolare, il gruppo norvegese ha misurato l’RNA messaggero (l’espressione genica) di diversi geni che sono coinvolti nella regolazione della attività del PDK. L’espressione genica è stata misurata in un gruppo di cellule facilmente accessibili, ovvero le cellule mononucleari del sangue periferico (PBMC). In pratica le PBMC sono: linfociti (cellule B, T, NK) e monociti (macrofagi). Ebbene, diversi geni coinvolti nella inibizione del piruvato deidrogenasi risultano sovra espressi, tanto nei maschi che nelle femmine con ME/CFS. Vediamo quali.
Figura 2. Gli enzimi PDK1, PDK2 e PDK 4 (riquadri rossi) sono sovraespressi nella ME/CFS.
Piruvato deidrogenasi chinasi
L’enzima piruvato deidrogenasi chinasi (PDK) ha come principale funzione quella di inibire il PDH. Dunque una buona domanda da farsi è: quanto PDK è espresso nelle cellule delle persone con ME/CFS? Se l’attività del PDH è ridotta, ci si potrebbe aspettare che l’espressione del PDK sia aumentata. In effetti questo è proprio ciò che è stato trovato. Premesso che l’enzima PDK esiste in 4 varianti negli esseri umani (PDK1, PDK2 etc), i norvegesi hanno trovato che PDK1, PDK2 e PDK4 sono sovraespressi nei pazienti ME/CFS (sia maschi che femmine), come indicato in figura 2. La maggiore differenza fra pazienti e controllo sano riguarda il PDK1 (p = 0.0018), quindi gli autori hanno indagato ulteriormente l’espressione di questo enzima, scoprendo che:
il PDK1 è espresso maggiormente nelle donne che negli uomini (figura 3, M);
il PDK1 è tanto più espresso quanto più severa è la patologia (figura 3, N);
il PDK1 è espresso in modo maggiore in chi è malato da più di 10 anni (figura 3, O);
il PDK1 è espresso maggiormente in chi compie meno di 2200 passi al giorno (figura 3, P).
In condizioni di salute, l’aumento della espressione genica dei PDK (con conseguente inibizione del PDH) si verifica durante periodi di mancanza di cibo. Quindi – in un certo senso – i pazienti ME/CFS soffrono le conseguenze della mancanza di cibo, pur vivendo nella disponibilità di risorse alimentari.
Figura 3. Livello del PDK1 in funzione del sesso dei pazienti (M), della severità della patologia (N), della durata della patolgia (O) e del numero di passi giornalieri (P).
PPARD
Fluge e Mella hanno esaminato – sempre nelle PBMC – l’espressione di una serie di altri geni, tra cui i geni PPARD (Peroxisome proliferator-activated receptor) che svolgono moltissimi ruoli nella regolazione del metabolismo cellulare, tra cui la regolazione del metabolismo degli zuccheri, dei carboidrati e delle proteine (Dunning, KR et al. 2014). I PPAR sono tre, ma tra questi l’unica differenza significativa tra pazienti e controllo sano è nel PPARD, che è sovra espresso nei pazienti (vigura 4, F). E’ interessante notare che PPARD controlla l’espressione genica del PDK4: infatti agonisti del recettore PPARD aumentano l’espressione genica del PDK4, nelle cellule muscolari umane (Abbot EL et al. 2005). Questo significa che una possibile catena di eventi è questa:
l’aumentata espressione di PPARD aumenta a sua volta l’espressione di PDK4, che inibisce il piruvato deidrogenasi.
Tuttavia questo è solo uno dei possibili percorsi regolatori coinvolti in questo complessi fenomeni.
Figura 4. Il gene PPARD è sovraespresso.
Espressione di altri geni
Gli Autori norvegesi hanno indagato l’espressione di numerosi altri geni, trovando altre differenze fra controlli sani e pazienti. In particolare, il gene SIRT4 è sovra espresso (figura 5, K), così come il gene MPC1 e il gene per la subunità PDHA del PDH (figura 5, I e G). E’ interessante che questi tre geni siano strettamente legati fra loro, infatti la proteina codificata da SIRT4 inibisce la subunità PDHA del PDH, mentre MPC1 aumenta l’espressione genica di PDHA. Insomma, da un lato si inibisce questo pezzo del PDH e dall’altro si cerca di aumentarne il numero, per compensare!
Figura 5. SIRT4 è sovraespressa (K), così come MPC1 (I) e PDHA (G).
Norvegia vs Pisa
Ho parlato più volte dell studio sulla espressione genica nei mitocondri dei pazienti ME/CFS della Reumatologia di Pisa (vedi qui). Ricordo che nella prima parte di quello studio si verificò l’espressione genica di 194 proteine mitocondriali in due gemelli, uno sano e uno con ME/CFS. In questo caso le cellule usate sono le piastrine, non le PBMC. Confrontando questi dati con quelli dello studio norvegese si osserva quanto segue.
La subunità PDHA (PDH E1-alpha) è sovra espressa sia nello studio norvegese (figura 5, G) che nel gemello con ME/CFS, che ne esprime 1.6 volte di più del fratello sano (Tabella S1); il gemello malato sovraesprime anche la subunità E1 beta.
la PDK1 e la PDK3 sono sovraespresse nel gemello malato (1.54 volte e 1.73 volte, rispettivamente) (Tabella S1). Nello studio norvegese ad essere sovraespresse sono la PDK1, la PDK2 e la PDK4 (figura 3), quindi esiste una discrepanza fra i due risultati.
Nel complesso si riscontrano dunque sia somiglianze che incongruenze fra questi due studi.
Cosa emerge?
Riassumo i punti salienti di quanto qui visto:
i mitocondri consumano amminoacidi al posto dello zucchero per produrre energia, e questo fa pensare a un possibile blocco dell’enzima piruvato deidrogensi (PDH);
una serie di enzimi che inibiscono il piruvato deidrogenasi (PDK1,2,4) sono effettivamente sovra espressi nei pazienti ME/CFS, rispetto ai controlli sani;
una proteina (la PPARD) che aumenta l’espressione del PDK4 è sovra espressa nei pazienti ME/CFS;
la proteina SIRT4, che inibische la subunità PDHA del PDH, è sovra espressa nei pazienti ME/CFS e, come forma di compenso, la subunità PDHA è sovra espressa, probabilemnte grazie alla sovra espressione di MPC1, che ha funzione regolatoria.
Detto in termini più semplici:
qualcosa induce l’espressione di una serie di meccanismi che riducono la funzione del piruvato deidrogenasi, costringendo i pazienti ME/CFS a bruciare amminoacidi al posto dello zucchero. Ma dei sistemi di compenso intervengono per cercare di riportare il metabolismo energetico alla normalità, senza riuscirci. I sistemi di compenso sono diversi fra maschi e femmine, ma il difetto metabolico a monte è lo stesso nei due sessi.
In un prossimo articolo riporterò la terza parte dello studio, che ha esaminato l’effetto del sangue dei pazienti su cellule umane in coltura. In quella sede torneremo sulle conclusioni da trarre da questi dati, e da quelli di altri studi.
we obtain
paolo maccallini 