Brain normalization is the rigid rotation and/or deformation of a 3D brain scan so that it can match a template brain. This is a necessary step for the analysis of brain magnetic resonance data (whether it is morphological or functional) as well as of brain PET data. It allows introducing a set of spatial coordinates such that each triplet (x,y,z) identifies the same anatomical region both in the brain we are studying and in the template (Mandal P.K. et al. 2012). So, for instance, once the normalization has been performed on an fMRI of the brain of a patient and on a set of fMRIs from a suited control group, we can compare the BOLD activation of each anatomical region of the patient’s brain with the activation of the corresponding anatomical region of the control group.
This part can be skipped, it is not necessary to read these passages for the understanding of the following paragraphs. If we assume that P is a point of the brain before normalization and we call S(P) its position after the normalization, we can consider the vectorial function:
which gives the new position of each point of the brain, after normalization. If P_0 is a point of the brain before the normalization, then we can write:
and remembering the expression of the differential of a vectorial function, we have
With a few passages we can write:
From the above formula, we realize that in order to define the configuration of the brain after normalization, we have to define, for each point P, a set of 12 parameters. Of these parameters, 6 describe the rigid movement and can be considered the same for each point. The other 6 (the coefficients of the matrix) are those that describe the change of shape and size and, in general, they are different for each point. The job of SPM is to calculate these parameters.
There are several templates (also called brain atlases) that have been developed across the years. The very first one was published in 1957 (Talairach J. et al. 1957). The same researcher then participated in building one of the most used brain atlas ever, based on a single female brain, published in 1988 (Talairach J. et Tournoux P. 1988). Another widely used template is the so-called MNI-152. It was built adapting the 3D MRI brain scans of 152 healthy individuals to the Talairach and Tournoux template. The adaptation was achieved using both a rigid roto-translation and a deformation (Maintz J.B. et Viergever M.A. 1988). The second step is required for overcoming the problem of differences in brain shape and dimension that we encounter within the human population.
Available brain atlases have some limitations. One of them being the fact that despite diseased brains are the most widely studied, they are also the most difficult to register to a template built from healthy individuals, because of usually marked differences in shape and/or size. This is true for instance for brains of patients with Alzheimer’s disease (Mandal P.K. et al. 2012). Another important limitation is that registration algorithms perform poorly for the brainstem (particularly for pons and medulla) (Napadow V. et al. 2006). This might have represented a problem for the study of diseases where a possible involvement of the brainstem is suspected, like for instance ME/CFS (VanElzakker M. et al. 2019).
The SPM software package is one of the most widely used instruments for the analysis of functional brain imaging data (web page). It is freely available for download, but it requires that you have a MatLab copy in your computer. Those who don’t have a MatLab license can request and install a standalone version of SPM12 by following the instructions of this page.
Importing a DICOM file
Once you have installed SPM12 in your computer, the first step in order to register a brain is to convert the format the series of images are written in, to a format that SPM12 can read. MRI images are usually in .dcm format (DICOM) while SPM12 reads .nii files. In order to do that, click DICOM import (figure below, on the left, red ellipse), then click on DICOM files (on the right, red ellipse), then select your .dcm file and click DONE (below, red ellipse). If you then click DISPLAY (blue ellipse, left) you will see your MRI scan in another window (see next paragraph). A video tutorial on these operations is available here.
Setting the origin
To start the normalization process, it is highly recommended to set manually the origin of the coordinates. If this is done properly, the registration will not only take less time but, even more importantly, the chances of a successful normalization will increase. The origin is set at the level of the anterior commissure (figure below). To find this anatomical structure, you can follow this video. Once you have put the cross on the right place in the sagittal, coronal and axial windows, just click SET ORIGIN (red ellipse) and then save your work clicking REORIENT (blue ellipse).
In this step, SPM12 calculates the set of distortions that have to be applied to the source brain to adapt it to the template in MNI space. On the main menu select NORMALIZE (ESTIMATE) (figure, on the left, red). This will open the batch editor where you are asked to load the subject you want to apply normalization to (figure, right, red). You have also a set of estimation options (blue), that we leave as they are. Then you click the RUN button, the arrow on the top of the batch editor.
At this point, your PC will perform a set of calculations that will require a few minutes. At the end of this process, a new .nii file will be saved in the spm12 folder. This is the set of distortions that will allow your subject’s brain to be registered to the template.
Now click on NORMALIZE (WRITE) on the main menu. The batch editor will then ask you for the deformation field, which is the file generated in the previous step, and for the images to write, which is the scan of your subject (figure below). Select them, then press the RUN button on the batch editor. A new .nii file will be written in the spm12 folder. This is the normalized brain!
In the next figure, you have the normalized brain on the left and the initial scan of the same subject on the right. As you can see, there is an overall change of shape.
Now that we have normalized our brain in the MNI space, we can easily find anatomical regions within its sections. We can, for instance, load the normalized brain with MRIcron and overlay a template with Brodmann’s areas highlighted in various colours (figure below).
La versione in inglese di questo articolo è disponibile qui.
Molti dei miei lettori sono probabilmente a conoscenza dei tentativi attualmente fatti per simulare matematicamente il metabolismo energetico dei pazienti ME/CFS, integrando i dati metabolici con i dati genetici.In particolare, il dr.Robert Phair ha sviluppato un modello matematico delle principali vie metaboliche coinvolte nella conversione dell’energia, dall’energia immagazzinata nei legami chimici di grandi molecole come glucosio, acidi grassi e amminoacidi, all’energia immagazzinata nell’adenosina trifosfato (ATP), pronta per l’uso. Phair, che è un ingegnere, ha determinato le equazioni differenziali che regolano questa enorme quantità di reazioni chimiche e le ha adattate al profilo genetico trovato nei pazienti ME/CFS. Ma già alcuni anni fa due fisici pubblicarono un interessante modello matematico del metabolismo energetico durante e dopo l’esercizio, nei pazienti ME/CFS (Lengert N. et Drossel B. 2015). In quanto segue descriverò questo modello e le sue previsioni e vedremo da vicino queste equazioni differenziali.
Le vie metaboliche che sono state analizzate
Il modello di Lengert e Drossel estende due sistemi di equazioni differenziali precedentemente pubblicati che descrivono il comportamento della glicolisi, del ciclo di Krebs (enormemente semplificato come una singola reazione!), della catena di trasporto degli elettroni mitocondriale (descritta in dettaglio), del sistema della creatina chinasi e della conversione di adenosina difosfato (ADP) in ATP, nei muscoli scheletrici (Korzeniewski B. et Zoladz JA. 2001), (Korzeniewski B. et Liguzinski P. 2004). Gli autori hanno aggiunto equazioni per l’accumulo di lattato e il suo efflusso fuori dalla cellula, per la sintesi de novo di inosina monofosfato (IMP) durante il recupero, per la degradazione dell’adenosina monofosfato (AMP) in IMP, per la degradazione di IMP in inosina e ipoxantina. Tutte le vie coinvolte sono raccolte nella figura 1. Queste reazioni sono descritte da 15 equazioni differenziali e la soluzione è un insieme di 15 funzioni del tempo che rappresentano la concentrazione dei principali metaboliti coinvolti (come il lattato, il piruvato, l’ATP, ecc.). Diamo ora uno sguardo più da vicino a una di queste equazioni e alla struttura generale dell’intero sistema di equazioni.
Figura 1. Questa è una rappresentazione schematica dei percorsi metabolici descritti dal modello matematico sviluppato da Lengert e Drossel. In dettaglio: sintesi citosolica e degradazione di ADP, AMP e IMP (a sinistra), via della protein chinasi e glicolisi (centro), catena di trasporto degli elettroni e ciclo TCA (a destra). Da Lengert N. et Drossel B. 2015.
Equazioni differenziali per reazioni chimiche
Consideriamo l’equazione utilizzata dagli autori per la reazione catalizzata dalla lattato deidrogenasi (la trasformazione del piruvato in lattato, figura 2) dove si è anche tenuto conto dell’efflusso di lattato dal citosol. L’equazione differenziale è la seguente:
dove i tre parametri sono determinati sperimentalmente e i loro valori sono
Il primo descrive l’attività dell’enzima lattato deidrogenasi: più questo parametro è elevato, più l’enzima è attivo. Il secondo descrive la reazione inversa (dal lattato al piruvato). Il terzo è una misura di quanto lattato la cellula è in grado di trasportare al di fuori della sua membrana. Forse il lettore si è reso conto che l’equazione del lattato è una equazione differenziale ordinaria del primo ordine. Si dice “primo ordine” perché nell’equazione compare solo la derivata prima della funzione che dobbiamo determinare (lattato, in questo caso); “ordinario” si riferisce al fatto che il lattato è funzione di una sola variabile (il tempo, in questo caso). Si vede immediatamente che un’equazione come questa può essere scritta come segue:
Supponiamo ora di avere altre due equazioni differenziali di questo tipo, una per il piruvato e una per i protoni (le altre due funzioni del tempo che sono presenti nell’equazione):
Allora avremmo un sistema di tre equazioni differenziali ordinarie come questo:
I valori iniziali delle funzioni che dobbiamo determinare sono raccolti nell’ultima riga: questi sono i valori che le funzioni incognite assumono all’inizio della simulazione (t = 0). In questo caso, questi valori sono le concentrazioni di lattato, piruvato e protoni nel citosol, a riposo. Le tre funzioni del tempo sono chiamate la soluzione del sistema. Questo tipo di sistema di equazioni è un esempio di problema di Cauchy, e sappiamo dalla teoria matematica che non solo ha una soluzione, ma che questa soluzione è unica. Inoltre, mentre questa soluzione può non essere sempre facilmente trovata con metodi rigorosi, è abbastanza facile risolvere il problema con metodi approssimati, come il metodo di Runge-Kutta o il metodo di Heun. Detto questo, il sistema di equazioni differenziali ordinarie proposto da Lengert e Drossel per il metabolismo energetico è proprio come quello qui sopra, con l’eccezione che comprende 15 equazioni anziché tre. Quindi, la principale difficoltà in questo tipo di simulazione non è l’aspetto computazionale, ma la determinazione dei parametri (come quelli enzimatici) e dei valori iniziali, che devono essere raccolti dalla letteratura medica o devono essere determinati sperimentalmente, se non sono già disponibili. L’altro problema è come progettare le equazioni: esistono spesso diversi modi per costruire un modello matematico di una reazione chimica o di qualsiasi altro processo biologico.
Il modello matematico della ME/CFS
Come adattiamo ai pazienti ME/CFS un modello del metabolismo energetico che è stato impostato con parametri presi da esperimenti condotti su soggetti sani? Questa è un’ottima domanda, e abbiamo visto che Robert Phair ha dovuto usare i dati genetici dei pazienti ME/CFS relativi agli enzimi chiave del metabolismo energetico, al fine di impostare il suo modello. Ma questi dati non erano disponibili quando Lengert e Drossel hanno progettato le loro equazioni. E allora? I due fisici hanno cercato studi sulla fosforilazione ossidativa nei pazienti ME/CFS e hanno scoperto che qusto processo cellulare era stato misurato con diverse impostazioni sperimentali e da diversi gruppi e che il denominatore comune di tuti gli studi era una riduzione di funzione che andava da circa il 35% (Myhill S et al. 2009), (Booth, N et al 2012), (Argov Z. et al. 1997), (Lane RJ. et al. 1998) a circa il 20% (McCully KK. et al. 1996), (McCully KK. et al. 1999). Quindi l’idea degli autori è stata di moltiplicare i parametro enzimatici di ciascuna reazione appartenente alla fosforilazione ossidativa per un numero compreso tra 0,6 (grave ME / CFS) a 1,0 (persona sana). In particolare, i due fisici hanno scelto un valore di 0,7 per la ME/CFS, nei loro esperimenti in silico (cioè esperimenti virtuali condotti nel processore di un computer).
Previsioni del modello matematico
Il modello matematico è stato utilizzato per eseguire prove di esercizio in silico con varie lunghezze e intensità. Quello che Lengert e Drossel hanno trovato è stato che il tempo di recupero nel paziente ME/CFS medio era sempre maggiore se confrontato con quelli di una persona sana. Il tempo di recupero è definito come il tempo necessario affinché una cellula ripristini il suo contenuto di ATP (97% del livello in stato di riposo) dopo lo sforzo. Nella figura 3 si vedono i risultati della simulazione per un esercizio molto breve (30 secondi) e molto intenso. Come potete vedere, nel caso di una cellula sana (a sinistra) il tempo di recupero è di circa 600 minuti (10 ore) mentre una cellula di una persona con ME/CFS (a destra) richiede più di 1500 minuti ( 25 ore) per recuperare.
Un altro risultato interessante della simulazione è un aumento di AMP nei pazienti rispetto al controllo (figura 3, linea arancione). Ciò è dovuto all’uso compensativo delle due vie metaboliche in figura 4: la reazione catalizzata dall’adenilato chinasi, in cui due molecole di ADP sono utilizzate per produrre una molecola di ATP e una molecola di AMP; e la reazione catalizzata dalla deaminasi AMP, che degrada AMP in IMP (che viene quindi convertito in inosina e ipoxantina). Queste due reazioni sono utilizzate dai pazienti ME/CFS più che dal controllo sano, al fine di aumentare la produzione di ATP al di fuori dei mitocondri.
Se diamo un’occhiata più da vicino alle concentrazioni di AMP e IMP nelle 4 ore successive allo sforzo (figura 5), vediamo effettivamente una maggiore produzione di IMP (linea verde) e AMP (linea arancione) nei muscoli scheletrici dei pazienti (destra) rispetto ai controlli (sinistra).
Un’ulteriore via di compensazione utilizzata dai pazienti (secondo questo modello) è la produzione di ATP da ADP da parte dell’enzima creatina chinasi (figura 6). Questo è un altro modo che abbiamo per produrre ATP nel citosol senza l’aiuto dei mitocondri. In questo modello di ME/CFS, vi è un aumento nell’uso di questo percorso, che porta a una diminuzione della concentrazione cellulare di fosfocreatina e un aumento della concentrazione cellulare di creatina (figura 7).
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.
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.
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.
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).
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).
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.
I have performed a set of analysis on experimental data already published about autoimmunity to muscarinic receptors in ME/CFS. My predictions are that extracellular loop 2 and 3, and also transmembrane helix 5 of both muscarinic cholinergic receptors 4 and 3, are main autoantigens in a subset of ME/CFS patients. Moreover, I have found that autoimmunity to M4 and M3 ChR is independent of autoimmunity to beta 2 adrenergic receptor, also reported in ME/CFS patients.
Myalgic encephalomyelitis/chronic fatigue syndrome (ME/CFS) is a debilitating disease characterized by cognitive deficits, fatigue, orthostatic intolerance with symptoms exacerbated after exertion (post-exertional malaise, PEM) (IOM, 2015). This disease has no known cause but several abnormalities have been observed in energy metabolism (Tomas C. and Newton J. 2018), immune system, gut flora (Blomberg J. et al. 2018), brain (Zeineh MM. et al. 2014). A possible role for autoantibodies in the pathogenesis of the disease has been suggested by the finding of reactivity of patient sera to two nuclear antigens (Nishikai, et al., 1997), (Nishikai, et al., 2001), to cardiolipin (Hokama, et al., 2009), to HSP60 (Elfaitouri A. et al. 2013), and to muscarinic cholinergic (M ChR) and beta adrenergic receptors (ß AdR) (Tanaka S et al. 2003), (Loebel M et al. 2016); reactivity that was significantly elevated when compared to healthy contols. Reactivity to adrenergic and muscarinic Ch receptors has been confirmed by two independent groups, but these results have not been published yet (R). A role for autoantibodies in at least a subgroup of patients has also been suggested by a response to rituximab, a CD20 B cells depleting agent (Fluge Ø. et al. 2011), (Fluge Ø. et al. 20115), and to immunoadsorption (Scheibenbogen C. et al 2018). Sera response to muscarinic cholinergic receptors is confirmed in two studies but both of them used an immune assay with proteins coated on a plate. This kind of test does not allow to identify the exact autoepitope on the receptor and – even more importantly – it is subjected to false positive results because it exposes to sera surfaces of receptors that are hidden when they are in their physiological position (Ramanathan S et al 2016). Nevertheless, the amount of data provided in the study by Loebel et al. where reactivity of sera to 5 subtypes of muscarinic cholinergic receptors have been measured simultaneously, has – in our opinion – the potential to unveil the exact autoepitope(s). Thus, I performed a bioinformatical analysis on experimental data from this study in order to extract hidden information. I used a software for the in silico study of B cell epitope cross-reactivity (Maccallini P. et al. 2018) and a software for amino acid protrusion index calculation (Ponomarenko J. et al., 2008). Our prediction is that patients sera mainly react to three epitopes that belong to the second and third extracellular loop of M3 and M4 ChR, but also to a hidden epitope of the same two receptors, leading to possible false positive results of this test. I have also found that the reactivity to beta 2 adrenergic receptor (ß2 AdR) found in the study by Loebel et al. is not due to the same antibody that reacts to muscarinic cholinergic receptors.
Search for cross-reactive epitopes. Cross-reactivity between muscarinic cholinergic receptors M4 and M3, and between M4 and M1 has been studied in silico using EPITOPE, a software already described (Maccallini P. et al. 2018). Briefly, EPITOPE searches for cross-reactive epitopes shared between two proteins (let’s say protein A and protein B) by comparing each possible 7-mer peptide of A with each possible 7-mer peptide of B. The comparison is made using the algorithm by Needleman and Wunsch (Needleman SB. and Wunsch CD. 1970) with a gap model a + b·x, where a is the opening gap penalty, b is the extending one, and x is the extension of the gap. A penalty for gaps at the end of the alignment was also assumed. The choice for gap penalties and substitution matrix were done according to the theory already developed for peptide alignments (Altschul SF. 1991), (Karlin S. and Altschul SF. 1990). Available experimental data on cross-reactivity between γ enolase and α enolase (McAleese SM. et al. 1988) have been used for EPITOPE calibration: a score >60 was considered the cut-off for cross-reactivity, a score below 50 indicates non-cross-reactive epitopes; a score between 50 and 60 defines a borderline result. A simpler version of EPITOPE has been used for single local alignments. The main program used for M4-M3 comparison, its subroutine NeWalign and the substitution matrix employed are available for download. Primary structures used in this work have been downloaded from UniProt and are the following ones: M1 ChR (P11229), M3 ChR (P20309), M4 ChR (P08173), B2 AdR (P07550).
Surface exposure. In order to select only those 7-mer peptides that are on the surface of proteins, I have considered their mean protrusion indexes. A protrusion index of at least 0.5 has been considered the cut-off for surface exposure. Protrusion indexes of single amino acids have been calculated with ElliPro. A protrusion index of 0.5 means that the amino acid is outside the ellipsoid of inertia which includes 50% of the centers of mass of all the amino acids of the protein (Ponomarenko J. et al., 2008). For M4 ChR I have used the crystal structure 5DSG (Thal DM. et al. 2016). The 3D structure of human M3 ChR has not been experimentally determined yet, so I have used a theoretical model built using murine M3 ChR (PDB ID: 4DAJA) as a template, provided by ModBase.
Selection criteria. Our purpose is to predict to what epitopes of M3 and M4 ChRs sera from ME/CFS patients react. So I search for M4 ChRs 7-mer peptides that are cross-reactive to M3 ChR, but non-cross-reactive to M1 ChR. Moreover, they have to present surface exposure both on M4 and on M3 ChR (otherwise antibodies can’t reach them). So, selection criteria for M4 ChR epitopes are as follows:
they have to be cross-reactive to M3 ChR;
they have to be non-cross reactive to M1 ChR or borderline;
they have to present a mean protrusion index ≥0.5;
M3 ChR peptides to which thy cross-react 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.
The search for 7-mer peptides of M4 ChR that are cross-reactive to M3 ChR found 108 sequences. We then studied cross-reactivity to M1 ChR for each of these peptides and we found that 11 of them are non-cross-reactive and that other 9 peptides have borderline reactivity. None of these 20 peptides presented a cross-reactivity to B2 AdR (Table 1S, column 1). Scores between peptides of M4 ChR and the other two muscarinic cholinergic receptors are plotted in Figure 1. The distribution of scores from the comparison of M3 ChR with M1 ChR and with M3 ChR are reported in Figure 2. For the M4 ChR 20 epitopes mentioned above, we calculated the mean protrusion indexes and we did the same calculation for their cross-reactive peptides on M3 ChR. We also indicated their position with respect to the plasma membrane. All these data are collected in Table 2S. Once we apply selection criteria on these 20 peptides, we obtain 9 epitopes (Table 1). Of these selected epitopes, one belongs to a transmembrane helix: peptide 186-192 of M4 ChR, which cross-reacts to peptide 231-237 of M3 ChR. Peptide 418-431 of M4 ChR is partially immersed in the plasmatic membrane, even though its cross-reactive peptide of M3 ChR is entirely exposed to the extracellular space, and the same applies to the other two epitopes found (figure 1). Peptide 175-181 of M4 ChR cross-reacts to peptide 211-217 of M3 ChR; peptide 186-192 of M4 Chr cross-reacts to peptide 222-228 of M3 ChR; peptide 418-431 of M4 Chr cross-reacts to peptide 513-522 of M3 ChR. Sequences that fulfill selection criteria and their respective inverted sequences are collected in Table 2.
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 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.
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 have been described in Ss (Segal B. et al. 2008). A pathogenic role of these antibodies in fatigue for both ME/CFS and Sjögren syndrome could perhaps be due to their vasodilatory effect.
Our analysis unveiled reactivity to a hidden autoepitope, which belongs to transmembrane helix 5 of M3 and of M4 ChR. This epitope is buried inside the plasma membrane when these two receptors are in their physiological position, so this reactivity can’t contribute to the pathogenesis of ME/CFS.
None of the 7-mer peptides of M4 ChR that cross-react to M3 ChR and at the same time don’t cross-react to M1 ChR presents in silico reactivity to B2 AdR. This means that in those patients whose sera present reactivity to both M4-M3 ChR and B2 AdR, there are two distinct autoantibodies. This prediction of our model is consistent with the low correlation found by Loebel and colleagues between anti-M4 ChR and anti-B2 AdR antibodies (Loebel M et al. 2016).
Most B cells epitopes on non-denaturated proteins (i.e. proteins that conserve their tertiary structure) are believed to be conformational (Morris, 2007), so a significant limitation of this study is due to the fact that our analysis considers only linear epitopes. Nevertheless, the main limitation of this study remains by far my encephalopathy.
This analysis of previously published data suggests a role for the second and the third extracellular loop of M4 and M3 ChR as autoantigens in ME/CFS. It also predicts the presence of a hidden autoantigen and thus a risk of false-positive results with standard ELISA. The eight peptides found by this analysis and their inverse sequences (Table 2) should be employed as query sequences for the search for possible triggering pathogens and for other autoantigens. These predictions should be tested using both cell-based assays and ELISA tests with these 8 peptides coated on the plate.
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.
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.
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).
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 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).
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.
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 questo articolo dimostro che un test LTT per malattia di Lyme che utilizzi come uno degli antigeni la OspC (proteina integra) di B. burgdorferi sensu stricto può teoricamente risultare positivo (falso positivo) in soggetti con aumentata permeabilità intestinale.
Some lymphocyte transformation tests (LTT) popular in Europe for the diagnosis of Lyme disease, use full-length OspC of B. burgdorferi as one of their antigens and request a positive stimulation index against only one or two antigens, in order to be considered positive. In what follows, we demonstrate that, in the case of patients with gut bacteria translocation, such a test has a theoretical risk of false positive results.
Lymphocyte transformation test
Lymphocyte transformation test (LTT) is an assay which allows measuring the activity of peripheral blood Th cells against specific antigens. T cell activation starts shortly after infection, with T cells proliferation and the production of cytokines (such as INF-γ) which regulate the adaptive immune response (Sompayrac, 2012). As T cell response vanishes after the resolution of the infection (Kaech, et al., 2007)↑, LTT may be useful in providing a proof of active infection. When an LTT assay is performed, Th cells from peripheral blood of a patient are exposed to proteins from a particular pathogen. If a significant reaction is noted, which could be either Th cells proliferation or INF-γ expression, the assay is considered positive and suggestive of an active infection by that particular pathogen. The response is expressed through a number, often referred to as stimulatory index (SI). In Lyme disease, several attempts have been made in order to obtain such a tool, either by T cells proliferation assays or by INF-γ measures (Dressler, et al., 1991)↑, (Chen, et al., 1999)↑, (Valentine-Thon, et al., 2007)↑, (von Baehr, et al., 2012)↑, (Callister, et al., 2016 May)↑. Nevertheless, this procedure has not been fully recognized as useful at present and neither the European guidelines (Stanek, et al., 2011)↑ nor the CDC (Centers for disease control and prevention, 2015)↑ recommend the use of this kind of test.
Th cells activation and cross-reactive T cell epitopes
Th cells are activated when their T cell receptors (TCR) recognize a complementary antigen presented by MHC II molecules (see Figure 1) (Sompayrac, 2012). Peptides presented by MHC II to T helper cells are exclusively linear epitopes, and they have a length between 13 and 17 amino acids (Rudensky, et al., 1991)↑. Various experiments have demonstrated that peptides with 5 identical amino acids in a sequence of 10 have good chances to represent cross-reactive T cell epitopes (Root-Bernstein, 2014)↑. That said, the algorithm described above for the LTT test is not free from the risk of false positive results, as each protein used as antigen could present one or more linear epitopes of 10 amino acids which share at least 5 amino acids with some epitope of 10 amino acids from another pathogen. This risk is particularly high when the assay uses complete proteins as antigens, and when a high SI for only a few antigens is required in order to have a positive result of the test.
OspC and Pseudomonas aeruginosa
We have used BLAST from NCBI (National Library of Medicine), with OspC from Borrelia burgdorferi (strain ATCC 35210 / B31 / CIP 102532 / DSM 4680) identified by the swiss-prot ID Q07337 (↑) as query sequence, settings being as follows: expected threshold of 10, BLOSUM62 as substitution matrix, and a word of 3 amino acids. We have built a custom database with the main Phyla of the human gut microbiome observed in a healthy population, which are Bacteroides, Firmicutes, Proteobacteria, Verrucomicrobia, Actinobacteria, Tenericutes, and Euryarchaeota (Giloteaux, et al., 2016)↑. One of the possible matches that BLAST gives back is the following alignment between the query sequence and the outer membrane protein G (OprG) of Pseudomonas aeruginosa (PDB ID: 2X27):
As you can see, we have 6 identical amino acids in a peptide 10 amino acids long. This means that this peptide from Borrelia burgdorferi could theoretically bind a Th cell previously activated by P. aeruginosa. Peptide 111-120 from OspC is reported in Figure 2. Peptide 51-60 of OrpG is in Figure 3. The 3D structure of OspC from B. burgdorferi strain B31 used for that picture has been experimentally determined with X rays and a resolution of 2,51 Å in 2001 (Kumaran, et al., 2001)↑ and its MMDB ID is 15958 (↑). The conclusion from this data is that Th cells from a patient with an active infection by P. aeruginosa could proliferate and produce INF-γ when exposed to OspC from B. burgdorferi. In other words, a patient with an active P. aeruginosa infection would come out to have a positive LTT test for OspC.
Gut bacteria translocation
A disrupted mucosal barrier of the bowel, with consequent translocation of bacteria from the gut to the peripheral blood, has been described in patients with liver diseases (Zhu, et al., 2013)↑, chronic HIV infection (Openshaw, 2009)↑, Crohn’s disease (Wyatt, et al., 1993)↑, and in myalgic encephalomyelitis/chronic fatigue syndrome (ME/CFS) (Giloteaux, et al., 2016)↑. In ME/CFS it has been possible, in particular, to demonstrate the translocation of Pseudomonas aeruginosa, among other gram-negative enterobacteria. In fact serum concentration of IgA against lipopolysaccharides from P. aeruginosa and other enterobacteria has been found to be significantly greater in ME/CFS patients than in normal volunteers (Maes, et al., 2007)↑. Thus in ME/CFS patients the adaptive immune system usually reacts against pathogens which exit from the gut, and in particular, we know that it reacts against P. aeruginosa.
ME/CFS patients are among the main users of this kind of tests, because of the similarities between Lyme disease and the clinical picture of ME/CFS (Gaudino, et al., 1997)↑. ME/CFS patients have a high prevalence of increased gut permeability and gut microbiome translocation (Giloteaux, et al., 2016)↑, and their immune system reacts against P. aeruginosa in many cases (Maes, et al., 2007)↑. Thus, each LTT for Lyme disease which uses full-length OspC from B. burgdorferi ss as the antigen could theoretically lead to a high rate of false positive results in this population of patients. The Lyme disease LTT discussed above, which is currently popular in Europe, is one of such tests. More researches are warranted in order to confirm or exclude the theoretical danger of cross-reaction of Lyme disease LTT with gut microbiome. Moreover, on the basis of what here presented, it might be possible to develop an LTT specific for the diagnosis of gut bacteria translocation.