Archive for the ‘mathematical biology’ Category

The Higgins-Selkov oscillator, part 2

March 29, 2018

Some time ago I wrote about a mathematical model for glycolysis, the Higgins-Selkov oscillator. In fact I now prefer to call it the Selkov oscillator since it is a system of equations written down by Selkov, although it was obtained by modifying a previous model of Higgins. At that time I mentioned some of the difficulties involved in understanding the global properties of solutions of this system. I mentioned that I had checked that this system exhibits a supercritical Hopf bifurcation, thus proving that it has stable periodic solutions for certain values of the parameters. In the hope of getting more insights into this system I decided to make it the theme of a master’s thesis. Pia Brechmann, the student who got this subject, obtained a number of interesting results. After she submitted her thesis she and I decided to carry this a bit further so as to get a picture of the subject which was as comprehensive as possible. We have now written a paper on this subject. I cannot say that we were able to answer all the questions we would have liked to but at least we answered some of them. Here I will mention some of the things which are known and some which are not.

When written in dimensionless form the system contains two parameters \alpha (a positive real number) and \gamma (an integer which is at least two). I already mentioned a result about a Hopf bifurcation and this was originally only proved for \gamma=2. In the paper we showed that there is a supercritical Hopf bifurcation for any value of \gamma. I also previously mentioned doing a Poincaré compactification of the system and blowing up new stationary solutions which appear in this process and are too degenerate to analyse directly. The discussion of using blow-ups in polar coordinates previously given actually concerned the case \gamma=2 and does not seem practical for higher values of \gamma. It turns out that the technique of quasihomogeneous directional blow-ups, explained in the book ‘Qualitative theory of planar differential systems’ by Dumortier et al. can be used to treat the general case. This type of blow-up has the advantage that the transformations are given by monomials rather than trigonometric functions and that there is a systematic method for choosing good values of the exponents in the monomials.

We discovered a paper by d’Onofrio (J. Math. Chem. 48, 339) on the Selkov oscillator where he obtains interesting results on the uniqueness and stability of periodic solutions. It was suggested by Selkov on the basis of numerical calculations that for large \alpha there exist solutions with oscillations which grow arbitrarily large at late times (let us call these unbounded oscillations). We were not able to decide on a rigorous level whether such solutions exist or not. They cannot exist when a periodic solution does. When they exist they have to do with a heteroclinic cycle in the compactification. We showed that when a cycle of this kind exists it is asymptotically stable and in that case solutions with unbounded oscillations exist. However we were not able to decide whether a heteroclinic cycle at infinity ever exists for this system. What we did prove is that for all values of the parameters there exist unbounded solutions which are eventually monotone.

We also proved that when the steady state is stable each bounded solution converges to it and that when there exists a periodic solution it is unique and each bounded solution except the steady state converges to that. I find it remarkable that such an apparently harmless two-dimensional dynamical system is so resistent to a complete rigorous analysis.

Advertisements

Lotka’s system

March 11, 2018

The system of ODE \dot x=a-bxy, \dot y=bxy-cy was considered in 1910 by Lotka as a model for oscillatory chemical reactions (J. Phys. Chem. 14, 271). It exhibits damped oscillations but no sustained oscillations, i.e. no periodic solutions. It should not be confused with the famous Lotka-Volterra system for predator-prey interactions which was first written down by Lotka in 1920 (PNAS 6, 410) and which does have periodic solutions for all positive initial data. That the Lotka system has no periodic solutions follows from the fact that x^{-1}y^{-1} is a Dulac function. In other words, if we multiply the vector field defined by the right hand sides of the equations by the positive function x^{-1}y^{-1} the result is a vector field with negative divergence. This change of vector field preserves periodic orbits and it follows from the divergence theorem that the rescaled vector field has no periodic orbits. My attention was drawn to this system by the paper of Selkov on his model for glycolysis (Eur. J. Biochem. 4, 79). In his model there is a parameter \gamma which is assumed greater than one. He remarks that if this parameter is set equal to one the system of Lotka is obtained. Selkov obtains his system as a limit of a two-dimensional system with more complicated non-linearities. If the parameter \gamma is set to one in that system equations are obtained which are related to Higgins’ model of glycolysis. Selkov remarks that this last system admits a Dulac function and hence no sustained oscillations and this is his argument for discarding Higgins’ model and replacing it by his own. The equations are \dot x=a-b\frac{xy}{1+y+xy} and \dot y=b\frac{xy}{1+y+xy}-cy. (To obtain the limit already mentioned it is first necessary to do a suitable rescaling of the variables.) In this case the Dulac function is \frac{1+y+xy}{xy}.

The fact that the Lotka-Volterra system admits periodic solutions can be proved by exhibiting a conserved quantity. At this point I recall the well-known fact that while conserved quantities and their generalizations, the Lyapunov functions, are very useful when you have them there is no general procedure for finding them. This naturally brings up the question: if I did not know the conserved quantity for the Lotka-Volterra system how could I find it? One method is as follows. First divide the equation for \dot y by that for \dot x to get a non-autonomous equation for dy/dx, cheerfully ignoring points where \dot x=0. It then turns out that the resulting equation can be solved by the method of separation of variables and that this leads to the desired conserved quantity.

One undesirable feature of the Lotka-Volterra system is that it has a one-parameter family of periodic solutions and must therefore be suspected to be structurally unstable. In addition, if we consider a solution where predators are initially absent the prey population grows exponentially. The latter feature can be eliminated by replacing the linear growth term in the equation for the prey by a logistic one. A similar term corresponding to higher death rates at high population densities can be added in the equation for the predators but the latter modification has no essential effect. This is a Lotka-Volterra model with intraspecific competition. As discussed in the book ‘Evolutionary Games and Population Dynamics’ by Josef Hofbauer and Karl Sigmund, when this model has a positive steady state that state is globally asymptotically stable. The proof uses the fact that the expression which defines the conserved quantity in the usual Lotka-Volterra model defines a Lyapunov function in the case with intraspecific competition. This is an example of the method of obtaining conserved quantities or Lyapunov functions by perturbing those which are already known in special cases.

It follows from Poincaré-Bendixson theory that the steady states in the Lotka model and the Higgins model are globally asymptotically stable. This raises the question whether we could not find Lyapunov functions for those systems. I do not know how. The method used for Lotka-Volterra fails here because the equation for dy/dx is not separable.

Conference on mechanobiology and cell signalling in Oberwolfach

March 2, 2018

I just attended a conference in Oberwolfach in an area rather far away from my usual interests, although it was about mathematical biology. I did meet some interesting (known and unkown) people and encountered some new ideas. Here I will just discuss two talks which particularly caught my attention.

In a talk of Takashi Hiiragi I learned a number of interesting things about embryology. The specific subject was the mouse embryo but similar things should apply to the human case. On the other hand it is very far away from what happens in Drosophila, for instance. There is a stage where the first eight cells are essentially identical. More precisely there is a lot of random variation in these cells, but no systematic differences. The subsequent divisions of these cells are not temporally correlated. By the time the number of cells has reached thirty-two an important differentiation step has taken place. By that time there are some of the cells which belong to the embryo while the others will be part of the placenta. If the first eight cells are separated then each one is capable of giving rise to a complete mouse. (The speaker did seem to indicate some restriction but did not go into details.) In order to understand the development process better one of these cells is studied in isolation. The cell contains a clock and so it ‘knows’ that it is in the eight cell stage. It then develops into a group of four cells in the same way that the eight cells would normally develop into 32. Differentiation takes place. The key symmetry-breaking step takes place when one end of the cell (at the eight-cell stage) develops an area at one end where actin is absent. This polarization then influences the further motion of the cells. It is interesting that the interaction between the cells in these processes seems to have more to do with mechanical signals then with chemical ones.

There was a talk of Fredric Cohen about cholesterol. His claim was that the concentration of cholesterol as usually measured is not a useful quantity and that the quantity which should be measured is the chemical potential of cholesterol. This has to do with the fact that cholesterol is hardly soluble in water or in hydrophobic liquids. I must say that the term ‘chemical potential’ was something which was always very opaque for me. As a result of this talk I think I am beginning to see the light. The cholesterol in a cell is mainly contained in the cell membrane. However it is not simply dissolved there as single molecules. Instead most of the molecules are interacting with proteins or with other cholesterol molecules. The chemical potential has to do with how many molecules get transferred when the system is connected to a reservoir. Only those molecules which are free are available to be transferred. So the issues seem to be the relationship between the amounts of free and bound molecules and what the real significance of the concentration of free molecules is for understanding a system.

Energy budget models

February 26, 2018

In a previous post I mentioned a talk I gave on dinosaurs and promised more information at a later date. Now a paper related to this has appeared and I will keep the promise. A basic issue is the way in which dinosaurs regulated their body temperature. The traditional idea was that they were cold-blooded (exotherm) like crocodiles. Later it was suggested that they might have been warm-blooded (endotherm) like birds or mammals. Then it was claimed that this was unrealistic and that they were mesotherm. This means something in between exotherm and endotherm, with a limited control of body temperature. There do exist organisms like this, a notable example being the tuna. I was involved in a project with Jan Werner, Eva Maria Griebeler and Nikos Sfakianakis on this subject. The first published result coming from this effort has now appeared (J. Theor. Biol. 444, 83).

The long-term goal of this project is to understand the evolution of warm-blooded animals in connection with the evolution of birds from dinosaurs. This involves understanding the way in which animals allocate energy to different tasks. How much do they use for generating body heat and how much do they use for other tasks such as maintenance or reproduction? A first step is to find a parametrization of the possible energy allocation strategies. In other words we want to identify suitable variables which could be used to describe the evolutionary process we want to understand. This is the content of the paper which has just appeared. At this point it might be asked how we can find out about the energy consumption of dinosaurs at all. It turns out that there are general relations known between the energy consumption of an animal and its growth rate over its lifetime. Thus the growth curve of a dinosaur gives indirect information about its energy consumption. But how can we get information about the growth curve? This is something I learned in the course of this project. The large bones of dinosaurs exhibit annual growth rings like those known for trees. The rings are of different thicknesses and thus give information on the growth rate in different years.

The paper does not contain a detailed dynamic model of the energy use of an animal over its lifetime. Instead it introduces a set of possible time evolutions depending on a finite number of parameters and then tests (by numerical methods) whether these suffice to reproduce the experimental growth curves of a number of animals with sufficient accuracy. It is also checked that the results obtained are consistent with known facts about the age of sexual maturation of the different species. It turns out that the mathematical model is successful in fitting the experimental constraints. It is found that, as expected, the model predicts that exotherms continue to grow as long as they live while endotherms stop growing at a time comparable to the age of sexual maturity.

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 jth row and jth 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 ith 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.

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

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.

Kestrels and Dirichlet boundary conditions

February 28, 2017

The story I tell in this post is based on what I heard a long time ago in a talk by Jonathan Sherratt. References to the original work by Sherratt and his collaborators are Proc. R. Soc. Lond. B269, 327 (more biological) and SIAM J. Appl. Math. 63, 1520 (more mathematical). There are some things I say in the following which I did not find in these sources and so they are based on my memories of that talk and on things which I wrote down for my own reference at intermediate times. If this has introduced errors they will only concern details and not the basic story. The subject is a topic in population biology and how it relates to certain properties of reaction-diffusion equations.

In the north of England there is an area called the Kielder Forest with a lake in the middle and the region around the lake is inhabited by a population of the field vole Microtus\ agrestis. It is well known that populations of voles undergo large fluctuations in time. What is less known is what the spatial dependence is like. There are two alternative scenarios. In the first the population density of voles oscillates in a way which is uniform in space. In the second it is a travelling wave of the form U(x-ct). In that case the population at a fixed point of space oscillates in time but the phase of the oscillations is different at different spatial points. In general there is relatively little observational data on this type of thing. The voles in the Kielder forest are an exception to this since in that case a dedicated observer collected data which provides information on both the temporal and spatial variation of the population density. This data is the basis for the modelling which I will now describe.

The main predators of the voles are weasels Mustela\ nivalis. It is possible to set up a model where the unknowns are the populations of voles and weasels. Their interaction is modelled in a simple way common in predator-prey models. Their spatial motion is described by a diffusion term. In this way a system of reaction-diffusion equations is obtained. These are parabolic equations and to the time evolution is non-local in space. The unknowns are defined on a region with boundary which is the complement of a lake. Because of this we need not only initial values to determine a solution but also boundary conditions. How should they be chosen? In the area around the lake there live certain birds of prey, kestrels. They hunt voles from the air. In most of the area being considered there is very thick vegetation and the voles can easily hide from the kestrels. Thus the direct influence of the kestrels on the vole population is negligible and the kestrels to not need to be included in the reaction-diffusion system. They do, however, have a striking indirect effect. On the edge of the lake there is a narrow strip with little vegetation and any vole which ventures into that area is in great danger of being caught by a kestrel. This means that the kestrels essentially enforce the vanishing of the population density of voles at the edge of the lake. In other words they impose a homogeneous Dirichlet boundary condition on one of the unknowns at the boundary. Note that this is incompatible with spatially uniform oscillations. On the boundary oscillations are ruled out by the Dirichlet condition. When the PDE are solved numerically what is seen that the shore of the lake generates a train of travelling waves which propagate away from it. This can also be understood theoretically, as explained in the papers quoted above.