Cytokine dynamics, part 2

In a previous post I discussed a paper of Lev Bar-Or about the dynamics of the interactions between T cells and macrophages via the cytokines they produce. The model of that paper is a four-dimensional dynamical system with eighteen parameters. In that post I also discussed a paper of mine where I proved some mathematical results about this system. The features I captured were that for certain values of the parameters there is one stationary solution which is a global attractor while for other special values of the parameters there are three stationary solutions, two of which are stable. Both these cases were already seen in phase portraits included in the original paper. Some numerical simulations I did at that time did not reveal any more complicated behaviour.

More recently experiments I did using Mathematica led me to discover some new phenomena. For certain values of the parameters it is possible to extract a two-dimensional model system with just two parameters and and this already displays interesting dynamical phenomena. The parameter encodes the effect of antigen presentation. A two-dimensional system has the advantage that its qualitative properties can be displayed effectively in the plane. I found the command ‘Manipulate’ in Mathematica to be very helpful in searching through parameter space. First I used it on a phase portrait. This did not get me very far for the following reasons. I found it quite hard to recognize the important dynamical features in the phase portrait and I found that searching in parameter space quickly became tiring. Next I tried displaying nullclines as a method to look for stationary solutions. There it was easier to see what was going on since it is only necessary to see where two curves intersect. My first attempt with this approach also gave nothing new but the fact that it was less tiring meant that I was encouraged to try again. The second attempt revealed parameter values where there are seven stationary solutions. Once I had found these parameter values I could go back and compute the corresponding phase portrait. It shows that four of the seven steady states are sinks while the other two are saddle points.

Having seen these pictures I had the ambition to prove the existence of seven stationary solutions, four of which are stable, under conditions on the parameters which are as general as possible. I succeeded in doing this. In addition I was able to obtain generalizations of these results to the original four-dimensional system. These theorems and their proofs are contained in a new paper. I will now describe a key technique which was used. I first thought that this method was quite different from anything I had used before but I later saw that it does bear some resemblance to the method of Fuchsian equations. I leave it to the interested reader to think about what the relations are. The idea for proving the existence of stationary solutions is as follows. Often there is an obvious way of defining a mapping whose fixed points are in one-to-one correspondence with stationary solutions of the dynamical system. Suppose now that for a given system it is possible to guess that there is a stationary solution near some point with coordinates . Consider a small box of size containing this point. If it is possible to show that this box is invariant under then the Brouwer fixed point theorem implies the existence of a stationary solution inside that box. If it can further be shown that the restriction of to the box is a contraction the stationary solution is the only one in the box. In the example the necessary estimates are obtained with the help of a parameter which has the property that the stationary solution tends to as it gets large. The procedure is then to choose the parameter large and the size of the box to be neither too small nor too large. With the obvious mapping this is enough to capture the four stable stationary solutions. The other three are more difficult and it is necessary to define a different mapping which has the same fixed points as but has other good mapping properties in addition. This more refined technique can be used to get the other three stationary solutions, at least in certain cases.

On the biological side these results suggest a scenario where the alternative states of the immune system are not just distinguished by the dominance of Th1 or Th2 cytokines but also different discrete possibilities for the strength of the dominance. Changing from one of these steady states to the other leads to a 50% increase in the total cytokine concentration and so it is by no means a small effect. I am not aware that this kind of phenomenon has been observed experimentally. One thing I will take away from this project is that in the future I plan to try to combine analytical and numerical approaches in an increasingly sophisticated way when trying to investigate the properties of dynamical systems.

Advertisement

Fuchsian equations

Fuchsian equations are a class of differential equations which have played a big role in my research and which have now found a niche in the mathematical relativity community. In fact they have a very wide range of applications. The story begins with Lazarus Fuchs in the mid nineteenth century. He was interested in the case of a higher order linear scalar ODE for a complex-valued function whose coefficients are analytic except for isolated singularities. The cases which I have mainly been concerned with are first order nonlinear systems of PDE for real-valued functions whose coefficients may only be smooth. It sounds like almost everything is different in the two cases and so what is the similarity? The basic idea is that a Fuchsian problem concerns a system of PDE with a singularity in the coefficients where the desire is to describe the behaviour of solutions near the singularity. If certain structural conditions are satisfied then for any formal power series solution of the problem there is a corresponding actual solution. In the original situation of Fuchs the series converged and the solution was its limit. More generally there need be no convergence and the series is only asymptotic. A related task is to find singular solutions of a regular system. By suitable changes of variables this can sometimes be transformed to the problem of finding a regular solution of a singular system.

With hindsight my first contact with a problem of Fuchsian type was in some work I did with Bernd Schmidt (Class. Quantum Grav. 8, 985) where we proved the existence of certain models for spherically symmetry stars in general relativity. Writing the equations in polar coordinates leads to a system of ODE with a singularity corresponding to the centre of symmetry. We proved the existence near the centre of solutions satisfying the appropriate boundary conditions by doing a direct iteration. In later years I was interested in the problem of mathematical models of the big bang in general relativity. I spent the academic year 1994-95 at IHES (I have written about this in a previous post) and in that period I often went to PDE seminars in Orsay. On one occasion the speaker was Satyanad Kichenassamy and his subject was Fuchsian equations. This is a long term research theme of his and he has written about it extensively in his books ‘Nonlinear Wave Equations’ and ‘Fuchsian Reduction: Applications to Geometry, Cosmology and Mathematical Physics’. I found the talk very stimulating and after some time I realized that this technique might be useful for studying spacetime singularities. Kichenassamy and I cooperated on working out an application to Gowdy spacetimes and this resulted in a joint paper. A general theorem on Fuchsian systems proved there has since (together with some small later modifications) been a workhorse for investigating spacetime singularities in various classes of spacetimes.

One of the highlights of this development was the proof by Lars Andersson and myself (Commun. Math. Phys. 218, 479) that there are very general classes of solutions of the Einstein equations coupled to a massless scalar field whose initial singularities can be described in great detail. In later work with Thibault Damour, Marc Henneaux and Marsha Weaver (Ann. H. Poincare 3, 1049) we were able to generalize this considerably, in particular to the case of solutions of the vacuum Einstein equations in sufficiently high dimensions.  For these results no symmetry assumptions were necessary. More recently these results were generalized in another direction by Mark Heinzle and Patrik Sandin (arXiv:1105.1643). Fuchsian systems have also been applied to the study of the late-time behaviour of cosmological models with positive cosmological constant. In the meantime there are more satisfactory results on this question useing other methods (see this post) but this example does show that the Fuchsian method can be applied to problems in general relativity which have nothing to do with the big bang. In general this method is a kind of machine for turning heuristic calculations into theorems.

The Fuchsian method works as follows. Suppose that a system of PDE is given and write it schematically as . I consider the case that the equation itself is regular and the aim is to find singular solutions. Let be an explicit function which satisfies the equation up to a certain order in an expansion parameter and which is singular on a hypersurface defined by . Look for a solution of the form where is less singular than . The original equation can be rewritten as , where is singular. Now the aim is to show that there is a unique solution of the transformed equation which vanishes as . A theorem of this kind was proved in the analytic case in the paper mentioned above which I wrote with Kichenassamy. Results on the smooth case are harder to prove and there are less of them known. A variant of the procedure is to define as a solution (not necessarily explicit) of an equation which is a simplified version of the equation . In the cases of spacetime singularities which have been successfully handled the latter system is the so-called velocity-dominated system.

When is a dynamical system of mass action type?

At the moment I am at the European Conference on Mathematical and Theoretical Biology in Krakow. This is a joint conference with the Society for Mathematical Biology and it is very large, with more than 900 participants. This leads to a huge number of parallel sessions and the need to choose very carefully in order to get the most profit from the conference. On one day, for instance, there were two cases with two sessions on immunology occurring simultaneously.

On Tuesday I went to a session on biochemical reaction networks. This included a talk by Gheorghe Craciun with a large expository component which I found enlightening. He raised the question of when a system of ODE with polynomial coefficients can be interpreted as coming from a system of chemical reactions with mass action kinetics. He mentioned a theorem about this and after asking him for details I was able to find a corresponding paper by Hars and Toth. This is in the Colloquia Mathematica Societatis Janos Bolyai, which is a priori not easily accessible. The paper is, however, available as a PDF file on the web page of Janos Toth. A chemical reaction network gives rise to a system of equations of the form where the and are polynomials with positive coefficients. They represent the contributions from reactions where the species with concentration is on the right and left side respectively. The result of Hars and Toth is that any system of this algebraic form can be obtained from a reaction network. It was pointed out by Craciun in his talk that this means that arbitrarily complicated dynamics can be incorporated into systems coming from reaction networks. If we have a system of the form we can replace it by . This changes the system but does not change the orbits of solutions. If, for instance, we start with the Lorenz system with unknowns , and we can simply translate the coordinates so as to move the interesting dynamics into the region where all coordinates are positive and then multiply the result by . This preserves the strange attractor structure. This result may be compared with the Perelson-Wallwork theorem discussed in a previous post.

The construction of the reaction network reproducing the given equations is not very complicated. The main problem is keeping track of the notation. Suppose we start with a system of equations in variables which is polynomial and satisfies the necessary condition already mentioned. The reaction network can be constructed in the following way. (It is not at all unique.) Introduce one species for each . The right hand side of each equation is a sum of terms of the form and one reaction is introduced for each of these terms. To explain what it is suppose without generality that it belongs to the first equation. If hen the reaction transforms the complex to the complex with rate constant . The only species where there is a net production is and so this reaction only contributes to the first equation. Moreover it does so with the desired term. On the other hand if the reaction transforms the complex to the complex with rate constant . The assumption on the system assures that .