Archive for the ‘mathematical biology’ Category

SMB meeting in Montreal

July 27, 2019

This week I have been attending the SMB meeting in Montreal. There was a minisymposium on reaction networks and I gave a talk there on my work with Elisenda Feliu and Carsten Wiuf on multistability in the multiple futile cycle. There were also other talks related to that system. A direction which was new to me and was discussed in a talk by Elizabeth Gross was using a sophisticated technique from algebraic geometry (the mixed volume) to obtain an upper bound on the number of complex solutions of the equations for steady states for a reaction network (which is then of course also an upper bound for the number of positive real solutions). There were two talks about the dynamics of complex balanced reaction networks with diffusion. I have the impression that there remains a lot to be understood in that area.

At this conference the lecture rooms were usually big enough. An exception was the first session ‘mathematical oncology from bench to bedside’ which was completely overfilled and had to move to a different room. In that session there was a tremendous amount of enthusiasm. There is now a subgroup of the SMB for cancer modelling which seems to be very active with its own web page and blog. I should join that subgroup. Some of the speakers were so full of energy and so extrovert that it was a bit much for me. Nevertheless, it is clear that this is an exciting area and I would like to be part of it. There was also a session of cancer immunotherapy led by Vincent Lemaire from Genentech. He and two others described the mathematical modelling being done in cancer immunotherapy in three major pharmaceutical companies (Genentech, Pfizer and Glaxo-Smith-Kline). These are very big models. Lemaire said that at the moment that there are 2500 clinical trials going on for therapies related to PD-1. A recurring theme in these talks was the difference between mice and men.

This morning there was a talk by Hassan Jamaleddine concerning nanoparticles used to present antigen. These apparently primarily stimulate Tregs more than effector T cells and can thus be used as a therapy for autoimmune diseases. He showed some impressive pictures illustrating clearance of EAE using this technique. A central theme was interference between attempts to use the technique in animals with two autoimmune diseases in different organs, e.g. brain and liver. I was interested by the fact that for what he was doing steady state analysis was insufficient for understanding the biology.

This afternoon, the conference being over, I took to opportunity to visit Paul Francois at McGill, a visit which was well worthwhile.

The fold-Hopf and Hopf-Hopf bifurcations

June 30, 2019

Bifurcations of a dynamical system \dot x=f(x,\alpha) can be classified according to their codimension. The intuitive idea is that the set of vector fields exhibiting a given type of bifurcation form a submanifold of the space of all vector fields of the given codimension. Well-known examples of bifurcations of codimension one are the cusp bifurcation, which can already occur when the variables x and \alpha are one-dimensional, and the Hopf bifurcation where the variable x must be two-dimensional but the variable \alpha can be one-dimensional. In fact the minimal dimension of the variable \alpha required corresponds to the codimension of the bifurcation. In this post I want to discuss two bifurcations of codimension 2, the fold-Hopf bifurcation, where x must be at least three-dimensional and the Hopf-Hopf bifurcation where x must be at least four-dimensional

The first typical feature of these bifurcations is the configuration of eigenvalues of D_x f at the bifurcation point. For a fold there is an eigenvalue zero. For a Hopf bifurcation there is a complex conjugate pair of non-zero purely imaginary eigenvalues. For a fold-Hopf bifurcation there is a zero eigenvalue and a complex conjugate pair of non-zero purely imaginary eigenvalues. For a Hopf-Hopf bifurcation there are two complex conjugate pairs of non-zero purely imaginary eigenvalues. The basic hope now is that if some genericity conditions are satisfied the system can be locally reduced to a normal form by a transformation of variables. This is true for the fold and Hopf cases but for fold-Hopf and Hopf-Hopf it is no longer true. A weaker goal which can be attained is to reduce the system to an approximate normal form so that the right hand side is the sum of a simple explicit expression and a higher order error term. The genericity assumptions are as follows. For the cusp the steady state should move with non-zero velocity when the parameter is changed and the steady state of the system for the bifurcation value of the parameter should be as non-degenerate as possible. This means that although f_x=0 there (which is the bifurcation condition) f_{xx}\ne 0. For the Hopf case the eigenvalues which are on the imaginary axis at the bifurcation value should move off the axis with non-zero velocity when the parameter is changed. At the same time the steady state at the bifurcation value should be as non-degenerate as possible. Solutions close to this steady state circle it and a corresponding Poincare mapping can be defined which describes how the distance from the steady state changes when the solution circles it once. Call this p(x). The fact that there is a bifurcation means that p'(0)=0 and p''(0) is automatically zero. The non-degeneracy condition is that p'''(0)\ne 0.

Now we come to the fold-Hopf bifurcation. One non-degeneracy condition combines conditions from the two simpler bifurcations in a simple way. It says that the position of the steady state and the real part of the eigenvalue move independently as the two parameters are changed. This is condition ZH0.3 in Theorem 8.7 in the book of Kuznetsov. I was confused by the fact that this condition involves a quantity \gamma (\alpha) which is apparently nowhere defined in the book. It does occur on one other page. On that page there is also a \Gamma (\alpha) which is defined and I think that the solution to the problem is that these two are equal. Assuming that that is correct then \gamma (\alpha) is the projection of the position of the steady state onto the kernel of the linearization at the bifurcation point. The remaining non-degeneracy conditions are conditions on the system at the bifurcation value of the parameter. At the moment I do not have an intuition for the meaning of those conditions.

In the case of the Hopf-Hopf bifurcation a genericity assumption which is qualitatively different from those we have seen up to now is a non-resonance condition, condition HH.0 of Kuznetsov. It says that the two imaginary parts of the eigenvalues at the bifurcation point should not exhibit linear relations with integer coefficients. The next condition is that the real parts of the eigenvalues move independently as the two parameters are changed (HH.5 of Kuznetsov). As in the fold-Hopf case the remaining non-degeneracy conditions are conditions on the system at the bifurcation value of the parameter which I do not understand intuitively.

When analysing the Hopf bifurcation it turns out that after doing a suitable transformation of variables and discarding some terms which can be considered small the resulting system is rotationally invariant. In polar coordinates the angular component is constant and the radial component is cubic in the radius. For the fold-Hopf bifurcation it is natural to proceed as follows. We do a linear transformation so that the new axes belong to the eigenspaces of the linearization. Moreover this transformation is chosen so that the restriction of the linearization to the plane correponding to the complex eigenvalues is in standard form. Then the normal form is rotationally invariant and can be expressed in cylindrical polar coordinates. The component in the angular direction depends only on the coordinate \xi along the axis while the other two components depend only on \xi and the radial coordinate \rho. Thus the analysis of the phase portrait of the system near the bifurcation point can be reduced to the analysis of a two-dimensional dynamical system on the half-plane \rho\ge 0 called the amplitude system. A steady state of the amplitude system with \rho=0 corresponds to a steady state of the full system while a steady state of the amplitude system with \rho>0 corresponds to a periodic solution of the full system.

The amplitude system contains a parameter s=\pm 1 and a parameter \theta. The signs of these two quantities are of crucial importance. If s=1 and \theta>0 then the system can be reduced to normal form. If s=-1 and \theta<0 then the situation is still relatively simple but adding a small perturbation typically causes a heteroclinic orbit to break. The most difficult case is that where s\theta<0 and it is in that case that chaos may occur. A steady state of the amplitude system away from the axis can undergo a Hopf bifurcation and this corresponds to the occurrence of an invariant torus in the full system and is a Neimark-Sacker bifurcation of the limit cycle. This torus can break up for some value of the parameters and this is what leads to chaos.

In the Hopf-Hopf case the normal form involves two angular coordinates where the dynamics are trivial and two radial coordinates r_1 and r_2. Thus we again obtain a two-dimensional amplitude system, this time defined on a quadrant. Kuznetsov distinguishes between a ‘simple’ and a ‘difficult’ case according to the parameters but at my present level of understanding they both look very difficult. For Hopf-Hopf the truncated normal form is generically never topologically equivalent to the full system. A Neimark-Sacker bifurcation is always present.

As mentioned in a previous post, both bifurcations discussed here have been observed numerically in an ecological model and Hopf-Hopf bifurcations (but not fold-Hopf bifurcations, this was stated incorrectly in the previous post) in a model for the MAPK cascade.

Stability in the multiple futile cycle

April 8, 2019

In a previous post I described the multiple futile cycle, where a protein can be phosphorylated up to n times. About ten years ago Wang and Sontag proved that with a suitable choice of parameters this system has 2k+1 steady states. Here k denotes the integral part of n/2. The question of the stability of these steady states was left open. On an intuitive level it is easy to have the picture that stable steady states and saddle points should alternate. This suggests that there should be k+1 stable states and k saddles. On the other hand it is not clear where this intuition comes from and it is very doubtful whether it is reliable in a high-dimensional dynamical system. I have thought about this issue for several years now. I had some ideas but was not able to implement them in practise. Now, together with Elisenda Feliu and Carsten Wiuf, we have written a paper where we prove that indeed there are parameters for which there exist steady states with the stability properties suggested by the intuitive picture.

How can information about the relative stability of steady states be obtained by analytical calculations? For this it is good if the steady states are close together so that their stability can be investigated by local calculations. One way they can be guaranteed to be close together is if they all originate in a single bifurcation as a parameter is varied. This is the first thing we arranged in our proof. The next observation is that the intuition I have been talking about is based on thinking in one dimension. In a one-dimensional dynamical system alternating stability of steady states does happen, provided degenerate situations are avoided. Thus it is helpful if the centre manifold at the bifurcation is one-dimensional. This is the second thing we arranged in our proof. To get the particular kind of alternating stability mentioned above we also need the condition that the flow is contracting towards the centre manifold. I had previously solved the case n=2 of this problem with Juliette Hell but we had no success in extending it to larger values of n. The calculations became unmanageable. One advantage of the case n=2 is that the bifurcation there was a cusp and certain calculations are done in great generality in textbooks. These are based on the presence of a one-dimensional centre manifold but it turns out to be more efficient for our specific problem to make this explicit.

The general structure of the proof is that we first reduce the multiple futile cycle, which has mass action kinetics, to a Michaelis-Menten system which is much smaller. This reduction is well-behaved in the sense of geometric singular perturbation theory (GSPT), since the eigenvalues of a certain matrix are negative. With this in place steady states can be lifted from the Michaelis-Menten system to the full system while preserving their stability properties. The bifurcation arguments mentioned above are then applied to the Michaelis-Menten system.

The end result of the ideas discussed so far is that the original analytical problem is reduced to three algebraic problems. The first is the statement about the eigenvalues required for the application of GSPT. This was obtained for the case n=2 in my work with Juliette but we had no idea how to extend it to higher values of n. The second is to analyse the eigenvalues of the linearization of the system about the bifurcation point. What we want is that two eigenvalues are zero and that all others have negative real parts. (One zero eigenvalue arises because there is a conservation law while the second corresponds to the one-dimensional centre manifold.) There are many parameters which can be varied when choosing the bifurcation point and a key observation is that this choice can be made in such a way that the linearization is a symmetric matrix, which is very convenient for studying eigenvalues. The third problem is to determine the leading order coefficient which determines the stability of the bifurcation point within the centre manifold.

I started to do parts of the algebra and I would describe it as being like entering a jungle with a machete. I was able to find a direction to proceed and show that some progress could be made but I very soon got stuck. Fortunately my coauthors came and built a reliable road to the final goal.

Herd immunity

February 14, 2019

I have a long term interest in examples where mathematics has contributed to medicine. Last week I heard a talk at a meeting of the Mainzer Medizinische Gesellschaft about vaccination. The speaker was Markus Knuf, director of the pediatric section of the Helios Clinic in Wiesbaden. In the course of his talk he mentioned the concept of ‘herd immunity’ several times. I was familiar with this concept and I have even mentioned it briefly in some of my lectures and seminars. It never occurred to me that in fact this is an excellent example of a case where medical understanding has benefited from mathematical considerations. Suppose we have a population of individuals who are susceptible to a particular disease. Suppose further that there is an endemic state, i.e. that the disease persists in the population at a constant non-zero level. It is immediately plausible that if we vaccinate a certain proportion \alpha of the population against the disease then the proportion of the population suffering from the disease will be lower than it would have been without vaccination. What is less obvious is that if \alpha exceeds a certain threshold \alpha_* then the constant level will be zero. This is the phenomenon of herd immunity. The value of \alpha_* depends on how infectious the disease is. A well-known example with a relatively high value is measles, where \alpha is about 0.95. In other words, if you want to get rid of measles from a population then it is necessary to vaccinate at least 95% of the population. It occurs to me that this idea is very schematic since measles does not occur as a constant rate. Instead it occurs in large waves. This idea is nevertheless one which is useful when making public health decisions. Perhaps a better way of looking at it is to think of the endemic state as a steady state of a dynamical model. The important thing is that this state is asymptotically stable in the dynamic set-up so that it recovers from any perturbation (infected individuals coming in from somewhere else). It just so happens that in the usual mathematical models for this type of phenomenon whenever a positive steady state (i.e. one where all populations are positive) exists it is asymptotically stable. Thus the distinction between the steady state and dynamical pictures is not so important. After I started writing this post I came across another post on the same subject by Mark Chu-Carroll. I am not sad that he beat me to it. The two posts give different approaches to the same subject and I think it is good if this topic gets as much publicity as possible.

Coming back to the talk I mentioned, a valuable aspect of it was that the speaker could report on things from his everyday experience in the clinic. This makes things much more immediate than if someone is talking about the subject on an abstract level. Let me give an example. He showed a video of a small boy with an extremely persistent cough. (He had permission from the child’s parents to use this video for the talk.) The birth was a bit premature but the boy left the hospital two weeks later in good health. A few weeks after that he returned with the cough. It turned out that he had whooping cough which he had caught from an adult (non-vaccinated) relative. The man had had a bad cough but the cause was not realised and it was attributed to side effects of a drug he was taking for a heart problem. The doctors did everything to save the boy’s life but the infection soon proved fatal. It is important to realize that this is not an absolutely exceptional case but a scenario which happens regularly. It brings home what getting vaccinated (or failing to do so) really means. Of course an example like this has no statistical significance but it can nevertheless help to make people think.

Let me say some more about the mathematics of this situation. A common type of model is the SIR model. The dependent variables are S, the number of individuals who are susceptible to infection by the disease, I, the number of individuals who are infected (or infectious, this model ignores the incubation time) and R, the number of individuals who have recovered from the disease and are therefore immune. These three quantities depend on time and satisfy a system of ODE containing a number of parameters. There is a certain combination of these parameters, usually called the basic reproductive rate (or ratio) and denoted by R_0 whose value determines the outcome of the dynamics. If R_0\le 1 the infection dies out – the solution converges to a steady state on the boundary of the state space where I=0. If, on the other hand, R_0>1 there exists a positive steady state, an endemic equilibrium. The stability this of this steady state can be examined by linearizing about it. In fact it is always stable. Interestingly, more is true. When the endemic steady state exists it is globally asymptotically stable. In other words all solutions with positive initial data converge to that steady state at late time. For a proof of this see a paper by Korobeinikov and Wake (Appl. Math. Lett. 15, 955). They use a Lyapunov function to do so. At this point it is appropriate to mention that my understanding of these things has been improved by the hard work of Irena Vogel, who recently wrote her MSc thesis on the subject of Lyapunov functions in population models under my supervision.

 

Mathematics enters the canon of molecular biology

December 18, 2018

On the shelf in my office I have a heavy red tome. It is the fifth edition of ‘Molecular Biology of the Cell’ by Alberts et al., a standard work. While reading an editorial by Arup Chakraborty I learned that the sixth edition of the book contains a section on mathematical models. Its title is ‘Mathematical Analysis of Cell Functions’ and it is almost 20 pages long. It is perhaps not very deep mathematically but I am nevertheless delighted that it is there. I find it very important that a textbook of central importance in cell biology takes time to discuss mathematics in that context and presents arguments why mathematics is valuable for biology. I wonder how much of what is written there would be understood by a mathematician with no background in molecular biology. Of course that is not the intended audience of the book and it is just my idle curiosity that makes me ask the question.

Let me list some of the main themes treated in the section. I give them as informal statements. (1) Negative feedback can stabilize a steady state. (2) Negative feedback plus delay can give rise to sustained oscillations. (3) Positive feedback plus cooperativity can give rise to multistability. (4) An incoherent feed-forward loop can give rise to a transient response to a signal. (5) A coherent feed-forward loop can give rise to a delay. I should emphasize that in all cases mentioned here the models involved use ODE. The delay in (2) does not have to do with a delay equation but just with a sufficiently large number of steps which take place one after another.

I would like to make some connections between the points (1)-(5) above and mathematical theorems. I will start with (2) since that seems the easiest case. I will replace (2) with another statement which is even easier: (2a) negative feedback is a necessary condition for the existence of an attracting periodic solution. This is related to ideas of Rene Thomas which I disussed in a previous post. Suppose we have a system of ODE \dot x=f(x) where the signs of the elements of the linearization Df(x) are independent of x. Here a ‘sign’ may be positive, negative or zero. As explained elsewhere a system of this kind defines a graph called the species graph or influence graph. A feedback loop is an oriented cycle in this graph and its sign is the product of the signs of the edges making up the cycle. The system exhibits negative feedback if its influence graph contains a negative feedback loop. The claim is then that there exist no attracting periodic solutions. Before I go further I should acknowledge my debt to Frederic Beck, who wrote his MSc under my supervision on the subject of feedback. The discussion which follows has benefitted a lot from his exposition of this subject. A system without negative feedback loops as defined above is called coherent. It was proved by Angeli, Hirsch and Sontag (J. Diff. Eq. 246, 3058) that by reversing the signs of some of the variables the system can be made quasicooperative. This means that along any feedback loop all signs are equal. Note that this does not mean that signs need be equal along an unoriented cycle. This can be illustrated by the case of an incoherent feedforward loop. It follows from the results of Angeli et al. that in a coherent system there exist no attracting periodic solutions. This result can be strengthened by replacing ‘attracting’ by ‘stable’. This has been proved by Richard and Comet (J. Math. Biol 63, 593). They also claim that if the condition that the signs of the entries in the linearization are spatially constant is dropped then the analogous statement is false. Their proof of this latter statement contained an error, as noticed by Frederic Beck. They repaired it in an erratum (J. Math. Biol. 70, 957).

This discussion shows that a lot is understood about negative feedback as a necessary condition for periodic solutions although there may still be more to be discovered on that subject. I think that much less is understood about the possibility of negative feedback being a sufficient condition for periodic solutions. As a concrete example let us consider the Selkov oscillator. It does contain a negative feedback loop consisting of one positive and one negative edge. For some parameter values it exhibits periodic solutions but for others there are only damped oscillations. Can the assumption of a negative feedback loop be supplemented in some way by an assumption of delay to give a sufficient condition? What should that assumption of delay be?

Hahn’s minimal model for the Calvin cycle

December 4, 2018

As has been discussed elsewhere in this blog there are many mathematical models for the Calvin cycle of photosynthesis in the literature and it is of interest to try to understand the relations between these models. One extreme case which may be considered is to look at the simplest possible models of this system which still capture some of the essential biological features. In 1984 Brian Hahn published a model of the Calvin cycle with 19 chemical species. Like many of the models in the literature this does not include a description of photorespiration. In 1987 he extended this to a model which does include photorespiration and has 31 chemical species. Later (1991) Hahn followed the strategy indicated above, looking for simplified models. He found one two-dimensional model and one three-dimensional model. Like many models in the literature the three-dimensional model contains the fifth power of the concentration of GAP (glyceraldehyde phosphate). This arises because elementary reactions have been lumped together to give an effective reaction where five molecules of the three-carbon molecule GAP go in and three molecules of a five-carbon sugar come out. Assuming mass action kinetics for this effective reaction then gives the fifth power. In an attempt to simplify the analysis yet further Hahn replaced the fifth power by the second power in his two-dimensional model.

Sadly, Brian Hahn is no longer around to study models of photosynthesis. He was a mathematics professor at the University of Cape Town. As head of department he was asked to give an opinion on whether the position of one member of the department should be extended. He found he could not make a positive recommendation and the matter was referred to a committee. At the meeting of this committee the candidate became very aggressive and, at least partly as a consequence of this, it was decided that his position should not be extended. Some time later, in the mathematics building, the candidate beat Hahn to death. At his trial it was decided that he was mentally ill and therefore could not be convicted of murder. After less than two years he was released from confinement. From what I have read it seems that the public discussion of this case has been influenced by issues of racism: Hahn was white and the man who killed him was black.

Let me return to photosynthesis. I thought that the small models of Hahn should be subjected to detailed scrutiny. This became the project of my PhD student Hussein Obeid. Now we have written a paper on the two-dimensional model of Hahn (with quadratic nonlinearity) and the related model with the fifth power. These can be studied with and without photorespiration. Without photorespiration it turns out that there are solutions which tend to infinity at late times. These are similar to solutions previously discovered in another model of the Calvin cycle by Juan Velazquez and myself, which we called runaway solutions. The difference is that in the two-dimensional model we can see how the runaway solutions fit into the global dynamics of the system. In this case there is a unique positive steady state which is unstable. The one-dimensional stable manifold of this point divides the space into two regions. In one region all solutions tend to infinity and have the same late-time asymptotics. In the other region all solutions tend to zero. The existence of these solutions is related to the phenomenon of overload breakdown studied by Dorothea Möhring and myself in another paper. When photorespiration is added a stable steady state comes in from infinity. Convergence to the stable steady state replaces the runaway asymptotics and all solutions are bounded.

We were able to show that the dynamics for the system with the fifth power is qualitatively similar to that with the second power. The only differences are that it not possible to obtain such detailed information about the parameter ranges for which certain qualitative behaviour occurs and that it is not ruled out that non-hyperbolic behaviour might happen in exceptional cases. Hahn derives his two-dimensional model by an informal reduction of the three-dimensional one. We showed how this can be made rigorous using geometric singular perturbation theory. This allows some limited information to be obtained about the dynamics of solutions of the three-dimensional system. We hope that this paper will serve as a good foundation for later investigations of the the three- and higher dimensional models of Hahn in the future.

Visit to Aachen

September 13, 2018

In the past few days I was in Aachen. I had only been there once before, which was for a job interview. It was just before Christmas. I only remember that it was dark, that there was a big Christmas market which I had trouble navigating through and that the building of the RWTH where the job was was large and imposing. This time the visit was more relaxed. I was there for a small meeting organized by Sebastian Walcher centred around the topic of singular perturbation theory and with an emphasis on models coming from biology.

I visited the cathedral of which a major part is older than all comparable buildings in northern Europe (eighth century). Even for someone with such rudimentary knowledge of architecture as myself it did appear very different from other cathedrals, with an oriental look. If I had been asked beforehand to imagine what the church of Charlemagne might look like I would never have guessed anything like the reality.

I gave a talk about my work with Juliette Hell about the existence of periodic solutions of the MAP kinase cascade. The emphasis was not only on the proof of this one result but also on explaining some general techniques with the help of this particular example. A theme running through this meeting was how, given a system of ODE depending on parameters, it is possible to find a small parameter \epsilon such that the limit \epsilon\to 0, while singular, nevertheless has some nice properties. In particular the aim is to obtain a limiting system with less equations and less parameters. Then we would like to lift some features of the dynamics of the limiting system to the full system. One talk, which I would like to mention, was by Satya Samal of the Joint Research Centre for Computational Biomedicine on tropical geometry. This is a term which I have often heard without understanding what it means intrinsically or how it can be applied to the study of chemical reaction networks. A lot of what was in the introduction of the talk can also be found in a paper by Radulescu et al. (Math. Modelling of Natural Phenomena, 10, 124). There was one thing in the talk which I found enlightening, which is not in the paper and which I have not seen elsewhere. I will try to explain it now. In this tropical business an algebra is introduced with ‘funny’ definitions of addition and multiplication for real numbers. It seems that often the formal properties of these operations are explained and then the calculations start. This is perhaps satisfactory for people with a strong affinity for algebra. I would like something more, an intuitive explanation of where these things come from. Let us consider multiplication first. If I start with two real numbers a and b, multiply them with \log\epsilon and exponentiate, giving the powers of \epsilon with these exponents, then multiply the results and finally take the logarithm and divide by \log\epsilon I get the sum a+b. This should explain why the tropical equivalent of the product is the sum. Now let us consider the sum. I start with numbers a and b, take powers of a quantity \epsilon with these exponents as before, then sum the results and take the logarithm. Thus we get \log (\epsilon^a+\epsilon^b). Now we need an extra step, namely to take the limit \epsilon\to 0. If b>a then the second term is negligible with respect to the first for \epsilon small and the limit is a\log\epsilon. If a>b then the limit is b\log\epsilon. Thus after dividing by \log\epsilon I get \min\{a,b\}. This quantity is the tropical equivalent of the sum. I expect that this discussion might need some improvement but it does show on some level where the ‘funny’ operations come from. It also has the following intuitive interpretation related to asymptotic expansions. If I have two functions whose leading order contributions in the limit \epsilon\to 0 are given by the powers a and b then the leading order of the product is the power given by the tropical product of a and b while the leading order of the sum is the power given by the tropical sum of a and b.

The technique of matched asymptotic expansions has always seemed to me like black magic or, to use a different metaphor, tight-rope walking without a safety net. At some point I guessed that there might be a connection between (at least some instances of) this technique and geometric singular perturbation theory. I asked Walcher about this and he made an interesting suggestion. He proposed that what is being calculated in leading order is the point at which a solution of a slow-fast system hits the critical manifold (asymptotically) thus providing an initial condition for the slow system. I must try to investigate this idea more carefully.

SIAM Life Science Conference in Minneapolis

August 11, 2018

Here are some impressions from the SIAM Life Science Conference. On Monay I heard a talk by Nathan Kutz about the structure and function of neural networks in C. elegans and a moth whose name I have forgotten. Coincidentally, when I was in Boston I heard a presentation by someone in the Gunawardena group about genes and transcription factors related to different types of neurons in C. elegans. On that occasion I asked how many of the few hundred neurons of C. elegans are regarded as being of different types and I was suprised to hear that it of the order of one hundred. Returning to the talk of Kutz, the key thing is that all the neurons in C. elegans and all their connections are known. Thus it is possible to simulate the whole network and reproduce central aspects of the behaviour of the worm. Another simplifying circumstance is that the motion of the worm can be described by only four parameters. A central message of the talk was that the behaviour of the nervous system of the worm itself can be reduced to a low dimensional dynamical system. This reminded me, with how much justification I do not know, of what I heard about the beaks of Darwin’s finches in Lisbon recently. As for the moth, the question described in the talk was how it learns to identify odours. A lot is known about the architecture of the neurons in the moth’s olfactory system. The speaker has compared the learning ability of this biological neural network with that of artifical neural networks used in machine learning. In reaching quite a high accuracy (70%) on the basis a small amount of training data the moth network beat all the artifical networks clearly. When given more data the result of the moth network hardly improved while those of the artificial network got better and better and eventually overtook the moth. The speaker suggested that the moth system is very effectively optimized for the task it is supposed to perform. I also heard a nice talk by Richard Bertram about canards in a three-dimensional ODE system describing the electrical activity of heart cells. I found what he explained about canards enlightening and this has convinced me that it is time for me to finally understand this subject properly. In particular he explained that the oscillations associated with a canard have to do with twisting of the stable and unstable manifolds. I also found what he said about the relations between the canard and concrete biological observations (the form of electrical signals from neurons) very interesting.

On Wednesday I heard a talk by Benjamin Ribba, who after an academic past now works with Roche. One of the main themes of his talk was cancer immunotherapy but I did not feel I learned too much there. What I found more interesting was what he said in the first part of his talk about chemotherapy for low-grade glioma. This is a disease which only develops very slowly, so that even after it has been discovered the tumour may grow very slowly over a period of several years. He showed patient data for this. He explained why it was reasonable to leave patients so long without therapy. The reason is that the standard course of therapy at that time was something which could only be given to a patient once in a lifetime (presumably due to side effects). Thus it made sense to wait with the therapy until the time which was likely to increase the patients lifetime as much as possible. During the slow growth phase the disease is typically asymptomatic, so that it is not a question of years of physical suffering for the patient. The speaker develeloped a mathematical model for the evolution of the disease (a model with a few ODE) and this could be made to fit the data well. The time course of the tumour growth is as follows. Before treatment it grows steadily. After the treatment is started the tumour starts to shrink. After the treatment is stopped it continues to shrink for a long time but usually eventually takes off again after a few years. In the model it was possible to try out alternatives to the standard treatment protocol and one was found which performed a lot better. The talk did not include any indication of whether this information had been used in paractise. I asked that question at the end of the talk. The answer was that by the time the analysis had been done the standard treatment has been replaced by a very different one and so the theory was too late to have clinical consequences.

On Thursday I heard a talk by Dominik Wodarz. One of his themes was oncolytic viruses but here I want to concentrate on another one. This had to do with chronic lymphocytic leukemia. This is a B cell malignancy and is treated using ibrutinib, an inhibitor of BTK (Bruton’s tyrosine kinase). Most of the malignant B cells form a tumour in the bone marrow while a small percentage of them circulate in the blood. Measuring the latter is the easiest way of monitoring the course of the disease. When the drug is given the tumour shrinks but it was not clear whether this is because the cells leave the bone marrow for the blood, where they die, or whether a large proportion die in the marrow. Using a simple mathematical model it could be shown that the second of these is the correct explanation. There are several things I like about this work. First, very simple models are used. Second, after comparison with data, these give convincing explanations for what is happening. Third, these conclusions may actually influence clinical practise. This is a good example of what mathematical biology should be like.

Visit to Boston

August 5, 2018

At the moment I am visiting Boston for the first time in my life. On Thursday I gave a talk at Harvard Medical School, at the invitation of Jeremy Gunawardena. The fact was not lost on me that this address was a quite exceptional one among the long list of places where I have given talks in my life. The subject of my talk was T cell signalling and, in particular, my work on this with Eduardo Sontag. I had a look around in the surrounding area and was impressed to see the variety of prestigious medical institutions which I knew by reputation, such as the Boston Children’s Hospital, The Dana-Farber Cancer Center and the Brigham and Women’s Hospital. I was happy how many pedestrians there were, indicating that at least in Boston the pedestrian is not an endangered species.

On Friday I gave a talk on T cell activation at Northeastern University at the invitation of Eduardo Sontag. This is an institution which has recently jumped up in the rankings by cutting its student numbers drastically and increasing its fees correspondingly. The basic subject of the talk was the same as on Thursday but it was modified in order to try to cater for a different audience. The talk on Friday included more mathematics (since I expected that on average the Friday audience would have a stronger mathematical background than that on Thursday) and more biology (since on Thursday I left out more of the biology which I assumed would be known to most of the audience). There were two other talks by Michael Margaliot and Yoram Zarai. They were on the subject of certain models for the way that ribosomes are allocated to mRNA. A key idea is that many ribosomes moving along an RNA molecule could get into traffic jams. Indeed there is a close relation between these and models for traffic jams in the literal sense. The question is then how the machinery of the cell should be organized to avoid such traffic jams. Typical mathematical techniques which were applied in the work explained in the talk were theorems about monotone systems of ODE (or corresponding control systems) satisfying some extra conditions leading to a simplification of their dynamics.

Today, Saturday, I spent most of the day walking around in Boston and Cambridge. In the morning it rained quite a bit but with an umbrella that was no problem. I was blissfully unaware of the fact that there had been a tornado warning and I also only heard later that a tornado did today hit a place not too far from Boston. In Cambridge I observed a turkey sitting on the pavement and looking into a shop window as if it was considering what it might buy. After a while it crossed the road and seemed to be well-educated since it was very careful to stay on the zebra crossing. By chance and without expecting anything special I happened to enter the palatial Boston Public Library. Later in the afternoon I was on the waterfront near the aquarium and the pleasant experience of sitting by the sea was heightened by the fact that the gull on top of the lamppost was a Ring-Billed Gull, a treat for my European eyes. I also spontaneously decided to take a boat trip around the harbour. In coming to Boston I had no special expectations about the city itself. Despite my short and superficial acquaintance with the place I can say that it has gained a secure place on my personal list of the most attractive cities I know.

ECMTB 2018 in Lisbon, part 2

July 27, 2018

In the mini-workshops at the conference related to chemical reaction network theory the most striking new result to be announced was by Balazs Boros. His preprint on this is arXiv:1710.04732. In fact it is necessary to say in what sense this is new but I will postpone that point and first discuss the mathematics. This result is very easy to formulate and I will try to make the discussion here as self-contained as possible. We start with a chemical reaction network consisting of reactions and complexes (the expressions on the left and right hand sides of the reactions like X+Y). This network can be represented by a directed graph where the nodes are the complexes and the edges are the reactions, oriented from the complex on the LHS to that on the RHS. The network is called weakly reversible if whenever we can get from node X to node Y by following directed edges we can get from Y to X. If we assume mass action kinetics and choose a positive reaction rate for each reaction we get a system of ODE describing the evolution of the concentrations of the substances belonging to the network in a standard way. Because of the interpretation we are interested in positive solutions, i.e. solutions for which all concentrations are positive. The theorem proved by Boros says: if the network is weakly reversible then the corresponding ODE system with mass action kinetics has at least one positive steady state. Actually he proves that the stronger (and more natural) statement holds that there is a solution in each positive stoichiometric compatibility class. Evidently the hypotheses only involve the graph of the network and require no details of the form of the complexes or the values of the reaction constants. Thus it is a remarkably strong result. In contrast to the statement of the theorem the proof is not at all easy. It involves reducing the desired statement to an application of the Brouwer fixed point theorem. Returning to the question of the novelty of the result, it was announced in a preprint of Deng et al. in 2011 (arXiv:1111.2386). It has never been published and it seems that the proof proposed by the authors is far from complete. Furthermore, the authors do not seem to be willing and able to repair the argument. Thus this result has been blocked for seven years. For anyone else it is an ungrateful task to provide a proof since a positive reaction from the authors of the original paper is doubtful. Furthermore other people not familiar with the background may also fail to give due credit to the author of the new paper. I think that with this work Balazs has done a great service to the reaction network community and we who belong to this community should take every opportunity to express our gratitude for this.

There was a nice talk by Ilona Kosiuk on her work with Peter Szmolyan on NF\kappaB. She expressed doubts about the derivation of the three-dimensional system mentioned in a previous post from the four-dimensional system. The work she explained in some detail concerned the four-dimensional system and uses GSPT to investigate the existence of periodic solutions of that system.

I feel that I got a lot more out of this conference than that I did out of that in Nottingham two years ago. I found more points of contact with my own research. This fact perhaps has less to do with the conference itself than it does with me. It is simply that I have penetrated a lot more deeply into mathematical biology during the last two years.