Conference on quantitative principles in biology at EMBL

November 5, 2017

I just returned from a conference with the title ‘Quantitative Principles in Biology’ at EMBL in Heidelberg. There was quite a lot of ‘philosophical’ discussion there about the meaning attached by different people to the word ‘model’ and the roles of experiment and theory in making scientific progress, in particular in biology. The main emphasis of the conference was on specific scientific results and I will mention a few of them below. In general I noticed a rather strong influence of physics at this conference. In particular, ideas from statistical physics came up repeatedly. I also met several people who had moved their research area to biology from physics, for instance string theory. I was happy that the research field of immunology was very well represented at the conference and that Jeremy Gunawardena said to me that he believes that immunology is one area of biology where mathematics has a lot to contribute to biology.

In a talk of Susanna Manrubia I learned about a class of biological systems I had never heard of before. These are the so-called multipartite viruses. In this type of system a genome is distributed between two or more virus particles. To allow the production of new viruses in a cell it must be infected with all the components. The talk described experiments in which this type of system was made to evolve in the laboratory. They start with foot and mouth virus in a situation where infection takes place with very high frequency. This is followed through many generations. It was found that in this context a bipartite virus could evolve from a normal virus. This did not involve an intermediate situation where (as in influenza) a virus has several segments of genetic material in a single virus particles. When the bipartite virus was propagated in isolation under circumstances where the frequency of infection was much lower it evolved back to the ordinary virus state. There exists a large variety of multipartite viruses in nature. They seem to be most common in plants (where apparently multiple virus infections are very common) but are also found in other organisms, including animals.

I had a poster at the conference on my work on T cell activation with Eduardo Sontag. To my surprise Paul Francois, on whose results this work built, was at the conference and so I had the chance to discuss these things with him. In our work there is a strong focus on the T cell receptor and at the conference there were several other contributions related to the modelling of other receptors. Eugenia Lyaschenko talked about how receptors can sense relative levels of ligand concentration and how endocytosis plays a role in this. Nikolas Schnellbächer had a poster on the theme of how dimerization of receptors can lead to major qualitative changes in their response functions. There are also important differences between homo- and heterodimers. I learned something about the mechanisms which play a role there. Yaron Antebi talked about the dynamical significance of a situation where several related ligands can bind to several related receptors.

Turing instabilities came up repeatedly at the meeting and were getting positive comments. One ‘take-home message’ for me was that the usual story that a Turing instability requires different diffusion constants should be weakened. It is based on the analysis of a system with two components and as soon as there are more than two components no such clear statement can be made. In addition, taking into account cell growth can help to trigger Turing instabilities.

A talk by Pieter Rein ten Wolde deepened my understanding of circadian clocks in cyanobacteria. They have a clock on a post-translational level involving phosphorylations of the KaiABC proteins and also a clock which involves translation. In the talk it was explained how the bacterium needs both of these for its time-keeping. A key point is that the period of the oscillator defining the clock (around 24 hours) can be longer than the period of the cell cycle. Thus this is like a clock which can continue to tick while it is being disassembled into its constituent parts and put together again.

The matrix-tree theorem

October 29, 2017

When I was an undergraduate graph theory was a subject which seemed to me completely uninteresting and I very consciously chose not to attend a lecture course on the subject. In the meantime I have realized that having become interested in reaction networks I could profit from knowing more about graph theory. One result which plays a significant role is the matrix-tree theorem. It asserts the equality of the number of subgraphs of a certain type of a given graph and a certain determinant. This could be used to calculate numbers of subgraphs. On the other hand, it could be used in the other direction calculate determinants and it is the second point of view which is relevant for reaction networks.

For the first part of the discussion here I follow these notes. If $G$ is an (unoriented) graph then a spanning tree is a connected subgraph which includes all vertices of $G$ and contains no cycles. The graph Laplacian $L$ of $G$ is the matrix for which $L_{ii}$ is equal to the number of edges containing the vertex $v_i$, $L_{ij}=-1$ if $i\ne j$ and there is an edge joining vertex $v_i$ to vertex $v_j$ and $L_{ij}=0$ otherwise. The first version of the matrix tree theorem, due to Kirchhoff, says that the number of spanning trees of $G$ can be calculated in the following way. Choose any vertex $v_j$ of $G$ and remove the $j$th row and $j$th column from $L$ to get a matrix $L_j$. Then compute the determinant of $L_j$. Surprisingly the value of the determinant is independent of $j$. The first version of the theorem can be obtained as a consequence of a second more sophisticated version proved a hundred years later by Tutte. This concerns directed graphs. A vertex $b$ is said to be accessible from $a$ if there is a directed path from $a$ to $b$. A vertex in a directed graph is called a root if every other vertex is accessible from it. A directed tree is a directed subgraph which contains a root and which, when the orientation is forgotten, is a tree (i.e. it is connected and contains no unoriented cycles). There is then an obvious way to define a spanning rooted tree. To formulate the second version of the theorem we introduce the matrix Laplacian $L$ of the directed graph. If $i=j$ the entry $L_{ij}$ is the number of edges with final vertex $i$. $L_{ij}=-1$ if $i\ne j$ and there is a directed edge from $v_i$ to $v_j$. $L_{ij}=0$ if $i\ne j$ and there is no edge connecting $v_i$ and $v_j$. The statement of the theorem is that the number of rooted trees with root $v_j$ is equal to the determinant of $L_j$, a matrix derived from $L$ as before. To see that the first version of the theorem follows the second first associate an oriented graph to an unoriented one by putting in oriented edges joining a pair of vertices in both directions whenever they are joined by an edge in the original graph. Then there is a bijection between trees rooted at $v$ in the unoriented graph and oriented trees rootes at $v$ in the oriented graph. On the other hand the two graphs have the same Laplacian since the number of edges ending at a vertex in the oriented graph is the same as the number of edges having that endpoint in the unoriented graph.

What is the connection of all this with reaction networks? Consider a chemical reaction network with only monomolecular reactions and reaction coefficients all equal to one. Then under mass action kinetics the evolution equations for the concentrations are $\dot x=Lx$, where $L$ is the Laplacian matrix of the network. There is a conserved quantity, which is the sum of the concentrations of all species. A steady state is an element of the kernel of $L$. The next part of the discussion follows a paper of Gunawardena (PLoS ONE, 7(5), e36321). He allows general reaction constants. The notion of the Laplacian is extended correspondingly. If the network is strongly connected (in the terminology of graph theory) or weakly reversible (in the terminology of chemical reaction network theory) the kernel of the Laplacian matrix is one-dimensional. Thus there is precisely one steady state in each stoichiometric compatibility class. It is moreover possible to compute a vector spanning the kernel of $N$ by graph theoretical methods. This also goes back to Tutte. To get the $i$th component first take the product of the reaction constants over a rooted tree and then sum these quantities over the rooted trees with root $v_i$. More generally the dimension of the kernel of $N$ is equal to the number of terminal strong linkage classes. It is also revealed there that the Laplacian corresponds to the matrix $A$ considered by Horn and Jackson. These ideas are also related to the King-Altman theorem of enzyme kinetics. I have the impression that I have as yet only scratched the surface of this subject. I hope I will be able to come back and understand more about it at a later date.

Sternberg’s theorem and smooth conjugacy

August 23, 2017

Suppose we have a dynamical system $\dot x=f(x)$ and a steady state $x_0$. Let $\phi$ be a homeomorphism defined on a neighbourhood of $x_0$ with $\phi (x_0)=0$. If we transform the system to coordinates $y=\phi(x)$ under what circumstances can $\phi$ be chosen so that the transformed system is equal to the linearization of the original system about $x_0$? The answer depends on the regularity required for $f$ and $\phi$. If $f$ is $C^1$ and $x_0$ is hyperbolic (i.e. all eigenvalues of the linearization at $x_0$ have non-vanishing real part) then there always exists a continuous mapping $\phi$ with this property. (In that case we cannot necessarily transform the vector field but we can transform its integral curves.) This is the Grobman-Hartman theorem. Even if $f$ is $C^\infty$ it is not in general possible to choose $\phi$ to be $C^2$.

Sternberg’s theorem shows that if $f$ is $C^\infty$ and there are no resonances then $\phi$ can be chosen to be $C^\infty$. The definition of a resonance in this context is as follows. Let $\lambda_i$ be the eigenvalues of the linearization. A resonance is a relation of the form $\lambda_j=\sum_i m_i\lambda_i$ with positive integers $m_i$ satisfying $\sum_i m_i\ge 2$. For example, if the eigenvalues are $\lambda_1=1$ and $\lambda_2=2$ we have the resonance $\lambda_2=2\lambda_1$. Consider the two-dimensional case. If there the eigenvalues are real with opposite signs then there can be no resonance and Sternberg’s theorem applies to all saddles in two dimensions. If both eigenvalues are positive then it is known that $\phi$ can be chosen to be $C^1$ but resonances can occur and there are cases where although $f$ is analytic $\phi$ cannot be chosen to be $C^2$. Even when there are no resonances it is in general not the case that if $f$ is analytic the mapping $\phi$ can also be chosen analytic. It is necessary to make a further assumption (no small divisors). This says roughly speaking that the expressions $\lambda_j-\sum_i m_i\lambda_i$ should not only be non-zero but that in addition they should not be able to approach zero too fast as the number of summands increases. This is the content of a theorem of Siegel.

All the results discussed up to now concern the hyperbolic case. What happens in the presence of purely imaginary eigenvalues? The Grobman-Hartman theorem has a generalization, the theorem of Shoshitaishvili, also known as the reduction theorem. In words it says that there is a continuous mapping $\phi$ which reduces the system to the product of the flow on any centre manifold and a standard saddle. A standard saddle is linear and so if the centre manifold is trivial this reduces to Grobman-Hartman. The next question is what can be achieved with mappings $\phi$ of higher regularity. A fact which is implicit in the reduction theorem is that linear systems where the origin is a hyperbolic fixed point whose stable and unstable manifolds have the same dimension are topologically equivalent. Evidently they are in general not smoothly equivalent since any smooth mapping preserving the origin will not change the eigenvalues of the linearization. However once this has been taken into account the theorem does generalise under the assumptions that there are no resonances among the eigenvalues with non-zero real parts.More precisly the eigenvalues $\lambda_i$ in the definition of a resonance must have non-zero real parts. By contrast quantity $\lambda_j$ should either be an eigenvalue with non-zero real part or zero. The result being referred to here is Takens’ theorem.

Suppose that for a steady state we introduce coordinates $x$, $y$ and $z$ corresponding to the stable, unstable and centre manifolds. Then the form achieved by the reduction in Takens theorem is that we get the equations $\dot x=A(z)x$, $\dot y=A(z)y$ and $\dot z=g(z)$ for matrices $A$ and $B$ and a function $g$ depending on the central variables $z$. One technical point is that while for a smooth system the reduction can be done by a mapping which is $C^k$ for any finite $k$ it cannot in general be done by a mapping which is itself infinitely differentiable. As an example, consider a two-dimensional system and a fixed point where the linearization has one zero and one non-zero eigenvalue. In this case there can be no resonances and the theorem says that there is a $C^k$ coordinate transformation which reduces the system to the form $\dot x=a(y)x$, $\dot y=b(y)$ with $a(0)\ne 0$, $b(0)=0$ and $b'(0)=0$.

SMB conference in Utah

July 21, 2017

Now Subgroups are being set up within the Society to concentrate on particular subjects. One of these, the Immunobiology and Infection Subgroup had its inaugural meeting this week and of course I went. There I and a number of other people learned a basic immunological fact which we found very surprising. It is well known that the thymus decreases in size with age so that presumably our capacity to produce new T cells is constantly decreasing. The obvious assumption, which I had made, is that this is a fairly passive process related to the fact that many systems in our bodies run down with age. We learned from Johnna Barnaby that the situation may be very different. It may be that the decrease in the size of the thymus is due to active repression by sexual hormones. She is involved in work on therapy for prostate cancer and said that it has been found that in men with prostate cancer who are getting drugs to reduce their testosterone levels it is seen that their thymus increases in size.

There were some recurrent themes at the conference. One was oncolytic viruses. These are genetically engineered viruses intended to destroy cancer cells. In modelling these it is common to use extensions of the fundamental model of virus dynamics which is very familiar to me. For instance Dominik Wodarz talked about some ODE models for oncolytic viruses in vitro where the inclusion of interferon production in the model leads to bistability. (In reponse to a question from me he said that it is a theorem that without the interferon bistability is impossible.) I was pleased to see how, more generally, a lot of people were using small ODE models making real contact to applications. Another recurrent theme was that there are two broad classes of macrophages which may be favourable or unfavourable to tumour growth. I should find out more about that. Naveen Vaidya talked about the idea that macrophages in the brain may be a refuge for HIV. Actually, even after talking to him I am not sure if it should not rather be microglia than macrophages. James Moore talked about the question of how T cells are eliminated in the thymus or become Tregs. His talk was more mathematical than biological but it has underlined once again that I want to understand more about positive and negative selection in the thymus and the related production of Tregs.

On a quite different subject there were two plenary talks related to coral reefs. A theme which is common in the media is that of the damage to coral due to climate change. Of course this is dominated by politics and usually not accompanied by any scientific information on what is going on. The talk of Marissa Blaskett was an excellent antidote to this kind of thing and now I have really understood something about the subject. The other talk, by Mimi Koehl, was less about the reefs themselves but about the way in which the larvae of snails which graze on the coral colonize the reef. I found the presentation very impressive because it started with a subject which seemed impossibly complicated and showed how scientific investigation, in particular mathematical modelling, can lead to understanding. The subject was the interaction of microscopic swimming organisms with the highly turbulent flow of sea water around the reefs. Investigating this involved among other things the following. Measuring the turbulent flow around the reef using Doppler velocimetry. Reconstructing this flow in a wave tunnel containing an artificial reef in order to study the small-scale structure of the transport of chemical substances by the flow. Going out and checking the results by following dye put into the actual reef. And many other things. Last but not least there was the mathematical modelling. The speaker is a biologist and she framed her talk by slides showing how many (most?) biologists hate mathematical modelling and how she loves it.

Conference on reaction networks and population dynamics in Oberwolfach

July 9, 2017

This is a belated report on a conference in Oberwolfach I attended a couple of weeks ago. The title includes two elements. The first held no suprises for me but the second was rather different from what I had expected. My expectation was that it would be about the evolution of populations of organisms. In fact it was rather focussed on models related to genetics, in other words with the question of how certain genetic traits spread through a population.

There was one talk which did have a connection to population biology in way closer to what I had expected. It happens all the time that ecosystems are damaged by exotic species imported, deliberately or by accident, from other parts of the world. There are also well-known stories of the type that to try to control exotic species number one exotic species number two is introduced and is itself very harmful. It is nice to hear an example where this kind of introduction of an exotic species was very successful. It is the case of the cassava plant which was introduced from South America to Africa and became a staple food there. Then an insect from South America (species number one) called the mealy bug was introduced accidentally and caused enormous damage. Finally an ecologist called Hans Herren introduced a parasitic wasp (species number two) from South America, restoring the food supply and saving numerous lives (often the number 20 million is quoted). More details of this story can be found here.

I want to mention one statement made in the talk of Gheorghe Craciun in Oberwolfach which I found intriguing. I might have heard this before but it did not stick in my mind properly. The statement is that the set of dynamical systems which possess a complex balanced steady state is a variety of codimension $\delta$, where $\delta$ is the deficiency. There seemed to be some belief in the audience that this variety is actually a smooth manifold. On one afternoon we had something similar to the breakout sessions in Banff. I suggested the topic for one of these, which was Lyapunov functions. The idea was to compare classes of Lyapunov functions which people working on different classes of dynamical systems knew. This certainly did not lead to any breakthrough but I think it did lead to a useful exchange of information. I documented the discussion for my own use and I think I could profit by following some of the leads there.

To finish I want to mention a claim made by Ankit Gupta in his talk. It did not sound very plausible to me but I expect that it at least contains a grain of truth. He said that these days more papers are published on $NF\kappa B$ than on all of mathematics.

Conference on mathematical analysis of biological interaction networks at BIRS

June 9, 2017

I have previously written a post concerning a meeting at the Banff International Research Station (BIRS). This week I am at BIRS again. Among the topics of the talks were stochastic chemical reaction networks, using reaction networks in cells as computers and the area of most direct relevance to me, multiple steady states and their stability in deterministic CRN. Among the most popular examples occurring in the talks in the latter area were the multiple futile cycle, the MAPK cascade and the EnvZ/OmpR system. In addition to the talks there was a type of event which I had never experienced before called breakout sessions. There the participants split into groups to discuss different topics. The group I joined was concerned with oscillations in phosphorylation cycles.

In the standard dual futile cycle we have a substrate which can be phosphorylated up to two times by a kinase and dephosphorylated again by a phosphatase. It is assumed that the (de-)phosphorylation is distributive (the number of phosphate groups changes by one each time a substrate binds to an enzyme) and sequential (the phosphate groups are added in one order and removed in the reverse order). A well-known alternative to this is processive (de-)phosphorylation where the number of phosphate groups changes by two in one encounter between a substrate and an enzyme. It is known that the double phosphorylation system with distributive and sequential phosphorylation admits reaction constants for which there are three steady states, two of which are stable. (From now on I only consider sequential phosphorylation here.) By contrast the corresponding system with processive phophorylation always has a unique steady state. Through the talk of Anne Shiu here I became aware of the following facts. In a paper by Suwanmajo and Krishnan (J. R. Soc. Interface 12:20141405) it is stated that in a mixed model with distributive phosphorylation and processive dephosphorylation periodic solutions occur as a result of a Hopf bifurcation. The paper does not present an analytical proof of this assertion.

It is a well-known open question, whether there are periodic solutions in the case that the modificiations are all distributive. It has been claimed in a paper of Errami et. al. (J. Comp. Phys. 291, 279) that a Hopf bifurcation had been discovered in this system but the claim seems to be unjustified. In our breakout sessions we looked at whether oscillations might be exported from the mixed model to the purely distributive model. We did not get any definitive results yet. There were also discussions on effective ways of detecting Hopf bifurcations, for instance by using Hurwitz determinants. It is well-known that oscillations in the purely distributive model, if they exist, do not persist in the Michaelis-Menten limit. I learned from Anne Shiu that it is similarly the case that the oscillations in the mixed model are absent from the Michaelis-Menten system. This result came out of some undergraduate research she supervised. Apart from these specific things I learned a lot just from being in the environment of these CRN people.

Yesterday was a free afternoon and I went out to look for some birds. I saw a few things which were of interest to me, one of which was a singing Tennessee warbler. This species has a special significance for me for the following reason. Many years ago when I still lived in Orkney I got an early-morning phone call from Eric Meek, the RSPB representative. He regularly checked a walled garden at Graemeshall for migrants. On that day he believed he had found a rarity and wanted my help in identifying it and, if possible, catching it. We did catch it and it turned out to be a Tennessee warbler, the third ever recorded in Britain. That was big excitement for us. I had not seen Eric for many years and I was sad to learn now that he died a few months ago at a relatively young age. The name of this bird misled me into thinking that it was at home in the southern US. In fact the name just came from the fact that the first one to be described was found in Tennessee, a migrant. The breeding range is much further north, especially in Canada. Thus it is quite appropriate that I should meet it here.

Becoming a German citizen

May 24, 2017

I first moved to Germany in 1987 and I have spent most of the time since then here. The total time I have spent elsewhere since my first arrival in Germany does not add up to more than two years. There is every reason to expect I will spend the rest of my life here. I am married to a German, I have a job here I like and a house. I could have applied for German citizenship a long time ago but I never bothered. Being an EU citizen living in Germany I had almost almost all privileges of a native. The only exception is that I could not vote except in local elections but since I am not a very political person that was not a big issue for me. It was also the case that for a long time I might have moved to another country. For instance I applied for a job in Vienna a few years ago and I might well have taken it if it had been offered to me. Now the chances of my moving are very small and so there is no strong argument left against becoming a German citizen.

What is more important is that there are now arguments in favour of doing so. With the EU showing signs of a possible disintegration the chance that I could lose the privileges I have here as an EU citizen is not so small that it should be neglected. The referendum in which the Scots voted on the possibility of leaving the UK was the concrete motivation for my decision to start the application process. Scotland stayed in the UK but then the Brexit confirmed that I had made the right decision. At the moment there is no problem with keeping British citizenship when obtaining German citizenship and I am doing so. This may change sometime, meaning that I will have to give up my British citizenship to keep the German one, but I see this as of minor importance.

As prerequisites for my application I had to do a number of things. Of course it was necessary to submit a number of documents but I have the feeling that the amount of effort was less than when obtaining the documents needed to get married here. I had to take an examination concerning my knowledge of the German language, spoken and written. It was far below my actual level of German and so from that point of view it was a triviality. It was just a case of investing a bit of time and money. I also had to do a kind of general knowledge test on Germany and on the state where I live. This was also easy in the sense that the questions were not only rather simple for anyone who has lived in the country for some time but they are also taken from a list which can be seen in advance. Again it just meant an investment of time and money. At least I did learn a few facts about Germany which I did not know before. In my case these things were just formalities but I think it does make sense that they exist. It is important to ensure that other applicants with a background quite different from mine have at least a minimal knowledge of the language and the country before they are accepted.

After all these things had been completed and I had submitted everything it took about a year before I heard that the application had been successful. This time is typical here in Mainz – I do not know how it is elsewhere in Germany – and it results from the huge backlog of files. People are queueing up to become German citizens, attracted by the prospect of a strong economy and a stable political system. Yesterday I was invited to an event where the citizenship of the latest group of candidates was bestowed in a ceremony presided over by the mayor. There were about 60 new citizens there from a wide variety of countries. The most frequent nationality by a small margin was Turkish, followed by people from other middle eastern countries such as Iraq and Iran. There were also other people from the EU with the most frequent nationality in that case being British. My general feeling was one of being slightly uneasy that I was engaged in a futile game of changing horses. It is sad that the most civilised countries in the world are so much affected by divisive tendencies instead of uniting to meet the threats confronting them from outside.

The Routh-Hurwitz criterion

May 7, 2017

I have been aware of the Routh-Hurwitz criterion for stability for a long time and I have applied it in three dimensions in my research and tried to apply it in four. Unfortunately I never felt that I really understood it completely. Here I want to finally clear this up. A source which I found more helpful than other things I have seen is https://www.math24.net/routh-hurwitz-criterion/. One problem I have had is that the Hurwitz matrices, which play a central role in this business, are often written in a form with lots of … and I was never sure that I completely understood the definition. I prefer to have a definite algorithm for constructing these matrices. The background is that we would like to understand the stability of steady states of a system of ODE. Suppose we have a system $\dot x=f(x)$ and a steady state $x_0$, i.e. a solution of $f(x_0)=0$. It is well-known that this steady state is asymptotically stable if all eigenvalues $\lambda$ of the linearization $A=Df(x_0)$ have negative real parts. This property of the eigenvalues is of course a property of the roots of the characteristic equation $\det(A-\lambda I)=a_0\lambda^n+\ldots+a_{n-1}\lambda+a_n=0$. It is always the case here that $a_0=1$ but I prefer to deal with a general polynomial with real coefficients $a_i, 0\le i\le n$ and a criterion for the situation where all its roots have negative real parts. It is tempting to number the coefficients in the opposite direction, so that, for instance, $a_n$ becomes $a_0$ but I will stick to this convention. Note that it is permissible to replace $a_k$ by $a_{n-k}$ in any criterion of this type since if we multiply the polynomial by $\lambda^{-n}$ we get a polynomial in $\lambda^{-1}$ where the order of the coefficients has been reversed. Moreover, if the real part of $\lambda$ is non-zero then it has the same sign as the real part of $\lambda^{-1}$. I find it important to point this out since different authors use different conventions for this. It is convenient to formally extend the definition of the $a_i$ to the integers so that these coefficients are zero for $i<0$ and $i>n$.

For a fixed value of $n$ the Hurwitz matrix is an $n$ by $n$ matrix defined as follows. The $j$th diagonal element is $a_j$, with $1\le j\le n$. Starting from a diagonal element and proceeding to the left along a row the index increases by one in each step. Similarly, proceeding to the right along a row the index decreases by one. In the ranges where the index is negative or greater than $n$ the element $a_n$ can be replaced by zero. The leading principal minors of the Hurwitz matrix, in other words the determinants of the submatrices which are the upper left hand corner of the original matrix, are the Hurwitz determinants $\Delta_k$. The Hurwitz criterion says that the real parts of all roots of the polynomial are negative if and only if $a_0>0$ and $\Delta_k>0$ for all $1\le k\le n$. Note that a necessary condition for all roots to have negative real parts is that all $a_i$ are positive. Now $\Delta_n=a_n\Delta_{n-1}$ and so the last condition can be replaced by $a_n>0$. Note that the form of the $\Delta_k$ does not depend on $n$. For $n=2$ we get the conditions $a_0>0$, $a_1>0$ and $a_2>0$. For $n=3$ we get the conditions $a_0>0$, $a_1>0$, $a_1a_2-a_0a_3>0$ and $a_3>0$. Note that the third condition is invariant under the replacement of $a_j$ by $a_{n-j}$. When $a_0a_3-a_1a_2>0$, $a_0>0$ and $a_3>0$ then the conditions $a_1>0$ and $a_2>0$ are equivalent to each other. In this way the invariance under reversal of the order of the coefficients becomes manifest. For $n=4$ we get the conditions $a_0>0$, $a_1>0$, $a_1a_2-a_0a_3>0$, $a_1a_2a_3-a_1^2a_4-a_0a_3^2>0$ and $a_4>0$.

Next we look at the issue of loss of stability. If $H$ is the region in matrix space where the Routh-Hurwitz criteria are satisfied, what happens on the boundary of $H$? One possibility is that at least one eigenvalue becomes zero. This is equivalent to the condition $a_n=0$. Let us look at the situation where the boundary is approached while $a_n$ remains positive, in other words the determinant of the matrix remains non-zero. Now $a_0=1$ and so one of the quantities $\Delta_k$ with $1\le k\le n-1$ must become zero. In terms of eigenvalues what happens is that a number of complex conjugate pairs reach the imaginary axis away from zero. The generic case is where it is just one pair. An interesting question is whether and how this kind of event can be detected using the $\Delta_k$ alone. The condition for exactly one pair of roots to reach the imaginary axis is that $\Delta_{n-1}=0$ while the $\Delta_k$ remain positive for $k. In a paper of Liu (J. Math. Anal. Appl. 182, 250) it is shown that the condition for a Hopf bifurcation that the derivative of the real part of the eigenvalues with respect to a parameter is non-zero is equivalent to the condition that the derivative of $\Delta_{n-1}$ with respect to the parameter is non-zero. In a paper with Juliette Hell (Math. Biosci. 282, 162), not knowing the paper of Liu, we proved a result of this kind in the case $n=3$.

Mathematical models for T cell activation

May 2, 2017

The proper functioning of our immune system is heavily dependent on the ability of T cells to identify foreign substances and take appropriate action. For this they need to be able to distinguish the foreign substances (non-self) from those coming from substances belonging to the host (self). In the first case the T cell should be activated, in the second not. The process of activation is very complicated and takes days. On the other hand it seems that an important part of the distinction between self and non-self only takes a few seconds. A T cell must scan the surface of huge numbers of dendritic cells for the presence of the antigen it is specific for and it can only spare very little time for each one. Within that time the cell must register that there is something relevant there and be induced to stay longer, instead of continuing with its search.

A mathematical model for the initial stages of T cell activation (the first few minutes) was formulated and studied by Altan-Bonnet and Germain (PloS Biol. 3(11), e356). They were able to use it successfully to make experimental predictions, which they could then confirm. The predictions were made with the help of numerical simulations. From the point of view of the mathematician a disadvantage of this model is its great complexity. It is a system of more than 250 ordinary differential equations with numerous parameters. It is difficult to even write the definition of the model on paper or to describe it completely in words. It is clear that such a system is difficult to study analytically. Later Francois et. el. (PNAS 110, E888) introduced a radically simplified model for the same biological situation which seemed to show a comparable degree of effectiveness to the original model in fitting the experimental data. In fact the simplicity of the model even led to some new successful experimental predictions. (Altan-Bonnet was among the authors of the second paper.) This is the kind of situation I enjoy, where a relatively simple mathematical model suffices for interesting biological applications.

In the paper of Francois et. al. they not only do simulations but also carry out interesting analytical calculations for their model. On the other hand they do not follow the direction of attempting to use these calculations to formulate and prove mathematical theorems about the solutions of the model. Together with Eduardo Sontag we have now written a paper where we obtain some rigorous results about the solutions of this system. In the original paper the only situation considered is that where the system has a unique steady state and any other solution converges to that steady state at late times. We have proved that there are parameters for which there exist three steady states. A numerical study of these indicates that two of them are stable. A parameter in the system is the number $N$ of phosphorylation sites on the T cell receptor complex which are included in the model. The results just mentioned on steady states were obtained for $N=3$.

An object of key importance is the response function. The variable which measures the degree of activation of the T cell in this model is the concentration $C_N$ of the maximally phosphorylated state of the T cell receptor. The response function describes how $C_N$ depends on the important input variables of the system. These are the concentration $L$ of the ligand and the constant $\nu$ describing the rate at which the ligand unbinds from the T cell receptor. A widespread idea (the lifetime dogma) is that the quantity $\nu^{-1}$, the dissociation time, determines how strongly an antigen signals to a T cell. It might naively be thought that the response should be an increasing function of $L$ (the more antigen present the stronger the stimulation) and a decreasing function of $\nu$ (the longer the binding the stronger the stimulation). However both theoretical and experimental results lead to the conclusion that this is not always the case.

We proved analytically that for certain values of the parameters $C_N$ is a decreasing function of $L$ and an increasing function of $\nu$. Since these rigorous results give rather poor information on the concrete values of the parameters leading to this behaviour and on the global form of the function we complemented this analytical work by simulations. These show how $C_N$ can have a maximum as a function of $\nu$ within this model and that as a function of $L$ it can have the following form in a log-log plot. For $L$ small the graph is a straight line of slope one. As $L$ increases it switches to being a straight line of slope $1-N/2$ and for still larger values it once again becomes a line of slope one, shifted with respect to the original one. Finally the curve levels out as it must do, since the function is bounded. The proofs do not make heavy use of general theorems and are in general based on doing certain estimates by hand.

All of these results were of the general form ‘there exist parameter values for the system such that $X$ happens’. Of course this is just a first step. In the future we would like to understand better to what extent biologically motivated restrictions on the parameters lead to restrictions on the dynamical behaviour.

New hope for primary progressive multiple sclerosis?

April 12, 2017

Multiple sclerosis is generally classified into three forms. The relapsing-remitting form is the most common initial form. It is characterized by periods when the symptoms get much worse separated by periods where they get better. The second form is the primary progressive form where the symptoms slowly and steadily get worse. It is generally thought to have a worse prognosis than the relapsing-remitting form. In many cases the relapsing-remitting form converts to a progressive form at some time. This is then the secondary progressive form. In the meantime there is a big variety of drugs on the market which are approved for the treatment of the RR form of MS. They cannot stop the disease but they can slow its progression. Until very recently there was no drug approved for the treatment of progressive MS. This has now changed with the approval of ocrelizumab, an antibody against the molecule CD20 which is found on the surface of B cells. It has been approved for both the RR form and some cases of the progressive form of MS.

Ocrelizumab acts by causing B cells to be killed. It has been seen to have strong positive effects in combatting MS in some cases. This emphasizes the fact that T cells, usually regarded as the main culprit causing damage during MS, are not alone. B cells also seem to play an important role although what role that is is not so clear. There previously existed an antibody against CD20, rituximab, which was used in the therapy of diseases other than MS. Ocrelizumab has had problemtic side effects, with a high frequency of infections and a slightly increased cancer risk. For this reason it has been abandoned as a therapy for rheumatoid arthritis. On the other hand the trial for MS has less problems with side effects.

One reason not to be too euphoric about this first treatment for progressive MS is the following. It has been shown to be effective against patients in the first few years of illness and those where there are clear signs of inflammatory activity in MRT scans. This suggests to me a certain suspicion. The different types of MS are not clearly demarcated. Strong activity in the MRT is typical of the RR form. So I wonder if the patients where this drug is effective are perhaps individuals with an atypical RR form where the disease activity just does not cross the threshold to becoming manifest on the symptomatic level for a certain time. This says nothing against the usefuleness of the drug in this class of patients but it might be a sign that its applicability will not extend to a wider class of patients with the progressive form in the future. It also suggests caution in hoping that the role of B cells in this therapy might help to understand the mechanism of progressive MS.