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

I have just attended the annual meeting of the Society for Mathematical Biology in Salt Lake City. Before reporting on my experiences there I will start with an apparently unrelated subject. I studied and did my PhD at the University of Aberdeen and I am very satisfied with the quality of the mathematical education I received there. There was one direction in mathematics in which I did not get as much training as I later needed, namely analysis and partial differential equations. This was not the fault of the lecturers, who had enough expertise and enthusiasm for teaching these things. It was the fault of the students. Since it was a small department the students chose which advanced courses were to be offered from a list of suggestions. Most of the students (all of them except me?) did not like analysis and so most of the advanced courses with a significant analysis content were not chosen. By the time I got my first postdoc position I had become convinced that in the area I was working in the research based on differential geometry was a region which was overgrazed and the thing to do was to apply PDE. Fortunately the members of the group of Jürgen Ehlers which I joined were of the same opinion. The first paper I wrote after I got there was on a parabolic PDE, the Robinson-Trautman equation. I had to educate myself for this from books and one of the sources which was most helpful was a book by Avner Friedman. Here is the connection to the conference. Avner Friedman, now 84 but very lively, gave a talk. Mathematically the subject was free boundary value problems for reaction-diffusion-advection equations, Friedman’s classic area. More importantly, these PDE problems came from the modelling of combination therapies for cancer. The type of therapies being discussed included antibodies to CTLA-4 and PD-1 and Raf inhibitors, subjects I have discussed at various places in this blog. I was impressed by how much at home Friedman seemed to be with these biological and medical themes. This is maybe not so surprising in view of the fact that he did found the Institute for Mathematical Biosciences in Ohio and was its director from 2002 to 2008. More generally I was positively impressed by the extent to which the talks I heard at this conference showed a real engagement with themes in biology and medicine and evidence of a lot of cooperations with biologists and clinicians. There were also quite a number of people there employed at hospitals and with medical training. As an example I mention Gary An from the University of Chicago who is trained as a surgeon and whose thoughtful comments about the relations between mathematics, biology and medicine I found very enlightening. There was a considerable thematic overlap with the conference on cancer immunotherapy I went to recently.

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.

A talk I found very interesting was by Sebastian Walcher. I already wrote briefly about a talk of his in Copenhagen in a previous post but this time I understood a lot more. The question he was concerned with is how to find interesting small parameters in dynamical systems which allow the application of geometric singular perturbation theory. In GSPT the system written in the slow time (with the smallness parameter included as a variable) contains a whole manifold of steady states, the critical manifold. The most straightforward theory is obtained when the eigenvalues of the linearization of the system transverse to the critical manifold lie away from the imaginary axis. This corresponds to the situation of a transversely hyperbolic manifold of steady states. The first idea of Walcher’s talk is that whenever we have a transversely hyperbolic manifold of steady states in a dynamical system this is an opportunity for identifying a small parameter. This may not sound very useful at first sight because it would seem reasonable that generic dynamical systems would never contain manifolds of steady states of dimension greater than zero. There is a reason why this observation is misleading for systems arising from reaction networks. In these systems the state space is defined by positivity conditions on the concentrations and there are also certain parameters (such as reaction constants and total amounts) which are required to be positive. To have a name let us call the region defined by these positivity conditions the conventional region of the spaces of states and parameters. In the conventional region manifolds of steady states are not to be expected. On the other hand it frequently happens that they arise when we go to the boundary of that region. A familiar example is the passage to the Michaelis-Menten limit in the system describing a reaction catalysed by an enzyme. This takes us from the extended mass action kinetics for substrate, free enzyme and substrate-enzyme complex to Michaelis-Menten kinetics for the substrate alone. Roughly speaking it is the limit where the amount of the enzyme is very small compared to the amount of the substrate. I often wondered whether there could not be a kind of ‘anti-Michaelis-Menten’ limit where the amount of enzyme is very large compared to the amount of the substrate. I asked Walcher whether he knew how to do this and how it fitted into his general scheme. He gave me a positive answer to this question and some references and I must look into this in detail when I get time. The reason for being interested in this is that if we can obtain suitable information about a limiting case on the boundary it may be possible obtain information on the part of the conventional region where a certain parameter is small but non-zero.

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.

Poincaré, chaos and the limits of predictability

March 5, 2017

In the past I was surprised that there seemed to be no biography of Henri Poincaré. I recently noticed that a biography of him had appeared in 2013. The title is ‘Henri Poincaré. A scientific biography’ and the author is Jeremy Gray. At the moment I have read 390 of the 590 pages. I have learned interesting things from the book but in general I found it rather disappointing. One of the reasons is hinted at by the subtitle ‘A scientific biography’. Compared to what I might have hoped for the book concentrates too much on the science and too little on the man. Perhaps Poincaré kept his private life very much to himself and thus it was not possible to discuss these aspects more but if this is so then I would have found it natural that the book should emphasize this point. I have not noticed anything like that. I also found the discussion of the scientific topics of Poincaré’s work too technical in many places. I would have preferred a presentation of the essential ideas and their significance on a higher level. There are other biographies of great mathematicians which made a better impression on me. I am thinking of the biography of Hilbert by Constance Reid and even of the slim volumes (100 pages each) on Gauss and Klein written in East Germany.

On important discovery of Poincaré was chaos. He discovered it in the context of his work on celestial mechanics and indeed that work was closely connected to his founding the subject of dynamical systems as a new way of approaching ordinary differential equations, emphasizing qualitative and geometric properties in contrast to the combination of complex analysis and algebra which had dominated the subject up to that point. The existence of chaos places limits on predictability and it is remarkable that these do not affect our ability to do science more than they do. For instance it is known that there are chaotic phenomena in the motion of objects belonging to the solar system. This nevertheless does not prevent us from computing the trajectories of the planets and those of space probes sent to the other end of the solar system with high accuracy. These space probes do have control systems which can make small corrections but I nevertheless find it remarkable how much can be computed a priori, although the system as a whole includes chaos.

This issue is part of a bigger question. When we try to obtain a scientific understanding of certain phenomena we are forced to neglect many effects. This is in particular true when setting up a mathematical model. If I model something using ODE then I am, in particular, neglecting spatial effects (which would require partial differential equations) and the fact that often the aim is not to model one particular object but a population of similar objects and I neglect the variation between these objects which I do not have under control and for whose description a stochastic model would be necessary. And of course quantum phenomena are very often neglected. Here I will not try to address these wider issues but I will concentrate on the following more specific question. Suppose I have a system of ODE which is a good description of the real-world situation I want to describe. The evolution of solutions of this system is uniquely determined by initial data. There remains the problem of sensitive dependence on initial data. To be able to make a prediction I would like to know that if I make a small change in the initial data the change in some predicted quantity should be small. What ‘small’ means in practice is fixed by the application. A concrete example is the weather forecast whose essential limits are illustrated mathematically by the Lorenz system, which is one of the icons of chaos. Here the effective limit is a quantitative one: we can get a reasonable weather forecast for a couple of days but not more. More importantly, this time limit is not set by our technology (amount of observational data collected, size of the computer used, sophistication of the numerical programs used) but by the system itself. This time limit will not be relaxed at any time in the future. Thus one way of getting around the effects of chaos is just to restrict the questions we ask by limits on the time scales involved.

Another aspect of this question is that even when we are in a regime where a system of ODE is fully chaotic there will be some aspects of its behaviour which will be predictable. This is why is is possible to talk of ‘chaos theory’- I know too little about this subject to say more about it here. One thing I find intriguing is the question of model reduction. Often it is the case that starting from a system of ODE describing something we can reduce it to an effective model with less variables which still includes essential aspects of the behaviour. If the dimension of the reduced model is one or two then chaos is lost. If there was chaos in the original model how can this be? Has there been some kind of effective averaging? Or have we restricted to a regime (subset of phase space) where chaos is absent? Are the questions we tend to study somehow restricted to chaos-free regions? If the systems being modelled are biological is the prevalence of chaos influenced by the fact that biological systems have evolved? I have seen statements to the effect that biological systems are often ‘on the edge of chaos’, whatever that means.

This post contains many questions and few answers. I just felt the need to bring them up.