In a previous post I mentioned Calvin’s Nobel lecture. Now I read it again and since I had learned a lot of things in the meantime I could profit from it in new ways. The subject of the lecture is the way in which Calvin and his collaborators discovered the mechanisms of the dark reactions of photosythesis. This involved years of experiments which I am not qualified to discuss. What I will do here is to describe some of the major conceptual components of this work. The first step was to discover which chemical substances are involved in the process. To make this a well-defined question it is necessary to fix a boundary between those substances to be considered and others. As their name suggests the dark reactions can take place in the dark and to start with the process was studied in the dark. It seems, however, that this did not lead to very satisfactory results and this led to a change of strategy. The dark reactions also take place in the light and the idea was to look at a steady state situation where photosynthesis is taking place in the presence of light. The dark reactions incorporate carbon dioxide into carbohydrates and the aim was to find the mechanism by which this occurs. At the end of the Second World War, when this work was done, carbon 14 had just become much more easily available due to the existence of nuclear reactors. Calvin also mentions that when doing difficult separations of compounds in his work on photosynthesis he used things he had learned when separating plutonium during the war. Given a steady state situation with ordinary carbon dioxide the radioactive form of the gas containing carbon 14 could be introduced. The radioactive carbon atoms became incorporated into some of the organic compounds in the plants used. (The principal subject of the experiment was the green alga Chlorella.) In fact the radioactive carbon atoms turned up in too many compounds – the boundary had been fixed too widely. This was improved on by looking what happened on sufficiently short time scales after the radioactive gas had been added, of the order of a few seconds. After this time the process was stopped, leading to a snapshot of the chemical concentrations. This meant that the labelled carbon had not had time to propagate too far through the system and was only found in a relatively small number of compounds. The compounds were separated by two-dimensional chromatography and those which were radioactive were located by the black spots they caused on photographic film. Calvin remarks ironically that the apparatus they were using did not label the spots with the names of the compounds giving rise to them. It was thus necessary to extract those compounds and analyse them by all sorts of techniques which I know very little about. It took about ten years. In any case, the endpoint of this process was the first major conceptual step: a set of relevant compounds had been identified. These are the carbon compounds which are involved in the reactions leading from the point where carbon dioxide enters the system and before too much of the carbon has been transferred to other systems connected to the Calvin cycle. While reading the text of the lecture I also had a modern picture of the reaction network in front of me and this was useful for understanding the significance of the elements of the story being told. From the point of view of the mathematician this step corresponds to determining the nodes of the reaction network. It remains to find out which compounds react with which others, with which stoichiometry. In looking for the reactions one useful source of information is the following. The carbon atoms in a given substance involved in the cycle are not equivalent to each other. By suitable experiments it can be decided which are the first carbon atoms to become radioactive. For instance, a compound produced in relatively large amounts right at the beginning of the process is phosphoglyceric acid (PGA) and it is found that the carbon in the carboxyl group is the one which becomes radioactive first. The other two carbons become radioactive at a common later time. This type of information provides suggestions for possible reaction mechanisms. Another type of input is obtained by simply counting carbon atoms in potential reactions. For instance, if the three-carbon compound PGA is to be produced from a precursor by the addition of carbon dioxide then the simple arthmetic relation indicates that there might be a precursor molecule with two carbons. However this molecule was never found and it turns out that the relevant arithmetic is . The reaction produces two molecules of PGA from a precursor with five carbon atoms, ribulose bisphosphate (RuBP). Combining the information about the order in which the carbon atoms were incorporated with the arithmetic considerations allowed a large part of the network to be reconstructed. Nevertheless the nature of one key step, that in which carbon dioxide is incorporated into PGA remained unclear. Further progress required a different type of experiment. The measurements used up to now are essentially measurements of concentrations at one time point (or very few time points). The last major step was taken using measurements of the dynamics. Here the concentrations of selected substances are determined at sufficiently many time points so as to get a picture of the time evolution of concentrations is certain circumstances. The idea is to first take measurements of PGA and RuBP in conditions of constant light. These concentrations are essentially time-independent. Then the light is switched off. It is seen that the concentration of PGA increases rapidly (it more than doubles within a minute) while that of RuBP rapidly decreases on the same time scale. This gives evidence that at steady state RuBP is being converted to PGA. This completes the picture of the reaction network. Further confirmation that the picture is correct is obtained by experiments where the amount of carbon dioxide available is suddenly reduced and the resulting transients in various concentrations monitored.

## Archive for the ‘mathematical biology’ Category

### Calvin on the Calvin cycle

July 31, 2015### Reaction networks in Copenhagen

July 9, 2015Last week I attended a workshop on reaction network theory organized by Elisenda Feliu and Carsten Wiuf. It took place in Copenhagen from 1st to 3rd July. I flew in late on the Tuesday evening and on arrival I had a pleasant feeling of being in the north just due to the amount and quality of the light. Looking at the weather information for Mainz I was glad I had got a reduction in temperature of several degrees by making this trip. A lot of comments and extra information on the talks at this conference can be found on the blog of John Baez and that of Matteo Polettini. Now, on my own slower time scale, I will write a bit about things I heard at the conference which I found particularly interesting. The topic of different time scales is very relevant to the theme of the meeting and the first talk, by Sebastian Walcher, was concerned with it. Often a dynamical system of interest can be thought of as containing a small parameter and letting this parameter tend to zero leads to a smaller system which may be easier to analyse. Information obtained in this way may be transported back to the original system. If the parameter is a ratio of time scales then the limit will be singular. The issue discussed in the talk is that of finding a suitable small parameter in a system when one is suspected. It is probably unreasonable to expect to find a completely general method but the talk presented algorithms which can contribute to solving this type of problem.

In the second talk Gheorghe Craciun presented his proof of the global attractor conjecture, which I have mentioned in a previous post. I was intrigued by one comment he made relating the concept of deficiency zero to systems in general position. Later he explained this to me and I will say something about the direction in which this goes. The concept of deficiency is central in chemical reaction network theory but I never found it very intuitive and I feel safe in claiming that I am in good company as far as that is concerned. Gheorghe’s idea is intended to improve this state of affairs by giving the deficiency a geometric interpretation. In this context it is worth mentioning that there are two definitions of deficiency on the market. I had heard this before but never looked at the details. I was reminded of it by the talk of Jeanne Marie Onana Eloundou-Mbebi in Copenhagen, where it played an important role. She was talking about absolute concentration robustness. The latter concept was also the subject of the talk of Dave Anderson, who was looking at the issue of whether the known results on ACR for deterministic reaction networks hold in some reasonable sense in the stochastic case. The answer seems to be that they do not. But now I return to the question of how the deficiency is defined. Here I use the notation for the deficiency as originally defined by Feinberg. The alternative, which can be found in Jeremy Gunawardena’s text with the title ‘Chemical reaction network theory for in silico biologists’ will be denoted by . Gunawardena, who seems to find the second definition more natural, proves that the two quantities are equal provided a certain condition holds (each linkage class contains precisely one terminal strong linkage class). This condition is, in particular, satisfied for weakly reversible networks and this is perhaps the reason that the difference in definitions is not often mentioned in the literature. In general , so that deficiency zero in the sense of the common definition implies deficiency zero in the sense of the other definition.

For a long time I knew very little about control theory. The desire to change this motivated me to give a course on the subject in the last winter semester, using the excellent textbook of Eduardo Sontag as my main source. Since that I had never taken the time to look back on what I learned in the course of doing this and this became clearer to me only now. In Copenhagen Nicolette Meshkat gave a talk on identifiability in reaction networks. I had heard her give a talk on a similar subject at the SIAM life science conference last summer and not understood much. I am sure that this was not her fault but mine. This time around things were suddenly clear. The reason is that this subject involves ideas coming from control theory and through giving the course I had learned to think in some new directions. The basic idea of identifiability is to extract information on the parameters of a dynamical system from the input-output relation.

There was another talk with a lot of control theory content by Mustafa Khammash. He had brought some machines with him to illustrate some of the ideas. These were made of Lego, driven by computers and communicating with each other via bluetooth devices. One of these was a physical realization of one of the favourite simple examples in control theory, stabilization of the inverted pendulum. Another was a robot programmed to come to rest 30 cm in front of a barrier facing it. Next he talked about an experiment coupling living cells to a computer to form a control system. The output from a population of cells was read by a combination of GFP labeling and a FACS machine. After processing the signal the resulting input was done by stimulating the cells using light. This got a lot of media attention unter the name ‘cyborg yeast’. After that he talked about a project in which programmes can be incorporated into the cells themselves using plasmids. In one of the last remarks in his talk he mentioned how cows use integral feedback to control the calcium concentration in their bodies. I think it would be nice to incorporate this into popular talks or calculus lectures in the form ‘cows can do integrals’ or ‘cows can solve differential equations’. The idea would be to have a striking example of what the abstract things done in calculus courses have to do with the real world.

My talk at the conference was on phosphorylation systems and interestingly there was another talk there, by Andreas Weber, which had a possibly very significant relation to this. I only became aware of the existence of the corresponding paper (Errami et. al., J. Comp. Phys. 291, 279) a few weeks ago and since it involves a lot of techniques I am not too familiar with and has a strong computer science component I have only had a limited opportunity to understand it. I hope to get deeper into it soon. It concerns a method of finding Hopf bifurcations.

This conference was a great chance to maintain and extend my contacts in the community working on reaction networks and get various types of inside information on the field

### Models for photosynthesis, part 2

May 22, 2015In my previous post on this subject I discussed the question of the status of the variables in the Poolman model of photosynthesis and in the end I was convinced that I had understood which concentrations are to be considered as dynamical unknowns and which as constants. The Poolman model is a modified version of the Pettersson model and the corresponding questions about the nature of the variables have the same answers in both cases. What I am calling the Pettersson model was introduced in a paper of Pettersson and Ryde-Pettersson (Eur. J. Biochem 175, 661) and there the description of the variables and the equations is rather complete and comprehensible. Now I will go on to consider the second question raised in the previous post, namely what the evolution equations are. The evolution equations in the Poolman model are modifications of those in the Pettersson model and are described relative to those in the original paper on the former model. For this reason I will start by describing the equations for the Pettersson model. As a preparation for that I will treat a side issue. In a reaction network a reaction whose rate depends only on the concentrations of the substances consumed in the reaction is sometimes called NAC (for non-autocatalytic). For instance this terminology is used in the paper of Kaltenbach quoted in the previous post. The opposite of NAC is the case where the reaction rate is modulated by the concentrations of other substances, such as activators or inhibitors.

The unknowns in the Pettersson model are concentrations of substances in the stroma of the chloroplast. The substances involved are 15 carbohydrates bearing one or more phosphate groups, inorganic phosphate and ATP, thus 17 variables in total. In addition to ordinary reactions between these substances there are transport processes in which sugar phosphates are moved from the stroma to the cytosol in exchange for inorganic phosphate. For brevity I will also refer to these as reactions. The total amount of phosphate in the stroma is conserved and this leads to a conservation law for the system of equations, a fact explicitly mentioned in the paper. On the basis of experimental data some of the reactions are classified as fast and it is assumed that they are already at equilibrium. They are also assumed to be NAC and to have mass-action kinetics. This defines a set of algebraic equations. These are to be used to reduce the 17 evolution equations which are in principle there to five equations for certain linear combinations of the variables. The details of how this is done are described in the paper. I will now summarize how this works. The time derivatives of the 16 variables other than inorganic phosphate are given in terms of linear combinations of 17 reaction rates. Nine of these reaction rates, which are not NAC, are given explicitly. The others have to be treated using the 11 algebraic equations coming from the fast reactions. The right hand sides of the five evolution equations mentioned already are linear combinations of those reaction rates which are given explicitly. These must be expressed in terms of the quantities whose time derivatives are on the left hand side of these equations, using the algebraic equations coming from the fast reactions and the conservation equation for the total amount of phosphate. In fact all unknowns can be expressed in terms of the concentrations of RuBP, DHAP, F6P, Ru5P and ATP. Call these quantities . Thus if the time derivatives of the can be expressed in terms of the we are done. It is shown in the appendix to the paper how a linear combination of the time derivatives of the with coefficients only depending on the is equal to . Moreover it is stated that the time derivatives of the can be expressed in terms of these linear combinations.

Consider now the Poolman model. One way in which it differs from the Pettersson model is that starch degradation is included. The other is that while the kinetics for the ‘slow reactions’ (i.e. those which are not classified as fast in the Pettersson model) are left unchanged, the equilibrium assumption for the fast reactions is dropped. Instead the fast reactions are treated as reversible with mass action kinetics. In the thesis of Sergio Grimbs (Towards structure and dynamics of metabolic networks, Potsdam 2009) there is some discussion of the models of Poolman and Pettersson. It is investigated whether information about multistability in these models can be obtained using ideas coming from chemical reaction network theory. Since the results from CRNT considered require mass action kinetics it is implicit in the thesis that the systems are being considered which are obtained by applying mass action to all reactions in the networks for the Poolman and Pettersson models. These systems are therefore strictly speaking different from those of Pettersson and Poolman. In any case it turned out that these tools were not useful in this example since the simplest results did not apply and for the more complicated computer-assisted ones the systems were too large.

In the Pettersson paper the results of computations of steady states are presented and a comparison with published experimental results looks good in a graph presented there. So whay can we not conclude that the problem of modelling the dynamics of the Calvin cycle was pretty well solved in 1988? The paper contains no details on how the simulations were done and so it is problematic to repeat them. Jablonsky et. al. set up simulations of this model on their own and found results very different from those reported in the original paper. In this context the advantage of the Poolman model is that it has been put into the BioModels database so that the basic data is available to anyone with the necessary experience in doing simulations for biochemical models. Forgetting the issue of the reliability of their simulations, what did Petterson and Ryde-Pettersson find? They saw that depending on the external concentration of inorganic phosphate there is either no positive stationary solution (for high values of this parameter) or two (for low values) with a bifurcation in between. When there are two stationary solutions one is stable and one unstable. It looks like there is a fold bifurcation. There is a trivial stationary solution with all sugar concentrations zero for all values of the parameter. When the external phosphate concentration tends to zero the two positive stationary solutions coalesce with the trivial one. The absence of positive stationary solutions for high phosphate concentrations is suggested to be related to the concept of ‘overload breakdown’. This means that sugars are being exported so fast that the production from the Calvin cycle cannot keep up and the whole system breaks down. It would be nice to have an analytical proof of the existence of a fold bifurcation for the Pettersson model but that is probably very difficult to get.

### Models for photosynthesis

May 15, 2015Photosynthesis is a process of central importance in biology. There is a large literature on modelling this process. One step is to identify networks of chemical reactions involved. Another is to derive mathematical models (usually systems of ODE) from these networks. Here when I say ‘model’ I mean ‘mathematical model’ and not the underlying network. In a paper by Jablonsky et. al. (BMC Systems Biology 5: 185) existing models are surveyed and a number or errors and inconsistencies in the literature are pointed out. This reminded me of the fact that a widespread problem in the biological literature is that the huge amount of data being generated these days contains very many errors. Here I want to discuss some issues related to this, concentrating on models for the Calvin cycle of photosynthesis and, in particular, on what I will call the Poolman model.

A point which might seem obvious and trivial to the mathematician is that a description of a mathematical model (I consider here only ODE models) should contain a clear answer to the following two questions. 1) What are the unknowns? 2) What are the evolution equations? One source of ambiguity involved in the first question is the impossibility of modelling everything. It is usually unreasonable to model a whole organism although this has been tried for some simple ones. Even if it were possible, the organism is in interaction with other organisms and its environment and these things cannot also be included. In practise it is necessary to fix a boundary of the system we want to consider and cut there. One way of handling the substances outside the cut in a mathematical model is to set their concentrations to constant values, thus implicitly assuming that to a good approximation these are not affected by the dynamics within the system. Let us call these external species and the substances whose dynamics is included in the model internal species. Thus part of answering question 1) is to decide on which species are to be treated as internal. In this post I will confine myself to discussing question 1), leaving question 2) for a later date.

Suppose we want to answer question 1) for a model in the literature. What are potential difficulties? In biological papers the equations (and even the full list of unknowns) are often banished to the supplementary material. In addition to being less easy to access and often less easy to read (due to typographical inferiority) than the main text I have the feeling that this supplementary material is often subjected to less scrutiny by the referees and by the authors, so that errors or incompleteness can occur more easily. Sometimes this information is only contained in some files intended to be read by a computer rather than a human being and it may be necessary to have, or be able to use, special software in order to read them in any reasonable way. Most of these difficulties are not absolute in nature. It is just that the mathematician embarking on such a process should ideally be aware of some of the challenges awaiting him in advance.

How does this look in the case of the Poolman model? It was first published in a journal in a paper of Poolman, Fell and Thomas (J. Exp. Botany, 51, 319). The reaction network is specified by Fig. 1 of the paper. This makes most of the unknowns clear but leaves the following cases where something more needs to be said. Firstly, it is natural to take the concentration of ADP to be defined implicitly through the concentration of ATP and the conservation of the total amount of adenosine phosphates. Secondly, it is explictly stated that the concentrations of NADP and NADPH are taken to be constant so that these are clearly external species. Presumably the concentration of inorganic phosphate in the stroma is also taken to be constant, so that this is also an external variable, although I did not find an explicit statement to this effect in the paper. The one remaining possible ambiguity involves starch – is it an internal or an external species in this model? I was not able to find anything directly addressing this point in the paper. On the other hand the paper does refer to the thesis of Poolman and some internet resources for further information. In the main body of the thesis I found no explicit resolution of the question of external phosphate but there it does seem that this quantity is treated as an external parameter. The question of starch is particularly important since this is a major change in the Poolman model compared to the earlier Pettersson model on which it is based and since Jablonsky et. al. claim that there is an error in the equation describing this step. It is stated in the thesis that ‘a meaningful concentration cannot be assigned to’ … ‘the starch substrate’ which seems to support my impression that the concentration of starch is an external species. Finally a clear answer confirming my suppositions above can be found in Appendix A of the thesis which describes the computer implementation. There we find a list of variables and constants and the latter are distinguished by being preceded by a dollar sign. So is there an error in the equation for starch degradation used in the Poolman model? My impression is that there is not, in the sense that the desired assumptions have been implemented successfully. The fact that Jablonsky et. al. get the absurd result of negative starch concentrations is because they compute an evolution for starch which is an external variable in the Poolman model. What could be criticised in the Poolman model is that the amount of starch in the chloroplast varies a lot over the course of the day. Thus a model with starch as an external variable could only be expected to give a good approximation to reality on timescales much shorter than one day.

### The species-reaction graph

May 14, 2015In the study of chemical reaction networks important features of the networks are often summarised in certain graphs. Probably that most frequently used is the species graph (or interaction graph), which I discussed in a previous post. The vertices are the species taking part in the network and the edges are related to the non-zero entries of the Jacobian matrix of the vector field defining the dynamics. Since the Jacobian matrix depends in general on the concentrations at which it is evaluated there is a graph (the local graph) for each set of values of the concentrations. Sometimes a global graph is defined as the union of the local graphs. A sign can be attached to each edge of the local graph according to the sign of the corresponding partial derivative. In the case, which does occur quite frequently in practise, that the signs of the partial derivatives are independent of the concentrations the distinction between local and global graphs is not necessary. In the general case a variant of the species graph has been defined by Kaltenbach (arXiv:1210.0320). In that case there is a directed edge from vertex to vertex if there is any set of concentrations for which the corresponding partial derivative is non-zero and instead of being labelled with a sign the edge is labelled with a function, namely the partial derivative itself.

Another more complicated graph is the species-reaction graph or directed species-reaction graph (DSR graph). As explained in detail by Kaltenbach the definition (and the name of the object) are not consistent in the literature. The species-reaction graph was introduced in a paper of Craciun and Feinberg (SIAM J. Appl. Math. 66, 1321). In a parallel development which started earlier Mincheva and Roussel (J. Math. Biol. 55, 61) developed results using this type of graph based on ideas of Ivanova which were little known in the West and for which published proofs were not available. In the sense used by Kaltenbach the DSR graph is an object closely related to his version of the interaction graph. It is a bipartite graph (i.e. there are two different types of vertices and each edge connects a vertex of one type with a vertex of the other). In the DSR graph the species define vertices of one type and the reactions the vertices of the other type. There is a directed edge from species to reaction if species is present on the LHS of reaction . There is a directed edge from reaction to species if the net production of species in reaction is non-zero. The first type of edge is labelled by the partial derivative of flux with respect to species . The second type is labelled by the corresponding stoichiometric coefficient. The DSR graph determines the interaction graph. The paper of Soliman I mentioned in a recent post uses the DSR graph in the sense of Kaltenbach.

A type of species-reaction graph has been used in quite a different way by Angeli, de Leenheer and Sontag to obtain conditions for the montonicity of the equations for a network written in reaction coordinates.

### Systems and synthetic biology in Strasbourg

March 26, 2015This past week I was at a conference entitled ‘Advances in Systems and Synthetic Biology’ in Strasbourg. The first talk was by Albert Goldbeter and anyone who has read many of the posts on my blog will realize that I went there with high expectations. I was not disappointed. It was a broad talk with the main examples being glycolysis, circadian rhythms and the cell cycle. A lot of the things he talked about were already rather familiar. It was nevertheless rewarding to hear the inside story. There were also enough themes in the talk which were quite new to me. For instance he mentioned that oscillations have been observed in red blood cells, where transcription is ruled out. I enjoyed listening to him, perhaps even more than I did reading his book. Another talk on Monday was by Jacques Demongeot. I am sure that he is a brilliant, versatile and highly knowledgeable person. Unfortunately he made no concessions in his talk for the benefit of non-experts. He jumped into the talk without saying where it was going and I did not have the background knowledge to be able to supply that information on my own. I felt as if I was flattened by a blast wave of information and unfortunately I understood essentially nothing.

The first talk on Tuesday was by Nicolas Le Novère and it had to do with engineering-type approaches to molecular biology. This is very far from my world but gave me some fascinating glimpses as to how this kind of thing works. Incidentally, I found out that Le Novère has a blog with a number of contributions which I found enlightening. The next talk was by François Fages and was focussed on computer science issues. It nevertheless contained more than one aspect which I found very interesting. At this point I should give some background. There are influential ideas on the relation of feedback loops to the qualitative behaviour of the solutions of dynamical systems due to René Thomas. They have been developed over many years by several people and a number of them are at this conference (El Houssine Snoussi, Marcelle Kaufman, Jacques Demongeot) and today I attended a tutorial held by Gilles Bernot on related themes. The basic idea is ‘positive feedback loops are necessary for bistability, negative feedback loops for periodic solutions’. I will not get into this more deeply here but I will just mention that some of the conjectures of Thomas have been made into theorems over the years. For instance a rather definitive version of the result on multistability was proved by Christophe Soulé. In his talk Fages mentioned a recent generalization of this result due to Sylvain Soliman. In the past I had the suspicion that the interest of the conjectures of Thomas was severely limited by the fact that the hypotheses rule out certain configurations which are very widespread in reaction networks of practical importance. It seemed to me that Fages made exactly this point and was saying that the improved results of Soliman overcome this difficulty. I must go into the matter in more detail as soon as I have time. Another point mentioned in the talk was an automated way to find siphons. This is a concept in reaction networks which I should know more about and I have the impression that in a couple of cases I have discovered these objects in examples without realizing that they were instances of this general concept.

On Wednesday there was an extended presentation by Oliver Ebenhöh. One speaker had cancelled and Oliver extended his presentation to fill the resulting extra time. I felt that listening to the presentation was time well spent and I did not feel my attention waning. He explained many things related to plant science and, in particular, the use of starch by plants. One key topic was the way in which a particular enzyme acts on chains of glucose monomers (generalization of maltose, which is the case of two units). It creates a kind of maximal entropy distribution of different lengths. The talk presented both a theoretical analysis and precise experimental results. The theoretical part involved an application of elementary thermodynamic ideas. I liked this and it brought me to a realization about my relation to physics. In the past I have been exposed to too much theoretical physics of a very pure kind, remote from applications to phenomena close to everyday life. It was refreshing to see basic physical ideas being applied in a down to earth way to the analysis of real experiments, in this case in biology.

In his talk on Thursday Joel Bader talked about his work on engineering yeast in such a way as to find out which combinations of genes are essential for survival. The aim is to look for a kind of minimal genome within yeast. One of the techniques of gathering information is to do a random recombination using a Cre-Lox system and looking to see which of the mutants produced are viable. The analysis of these experiments leads to consideration of self-avoiding or non-self-avoiding random walks and at this point I had a strange feeling of deja vu. A few weeks ago I gave a talk in the Hausdorff colloquium in Bonn. In this event two speakers are invited on one afternoon and their themes are not necessarily correlated. On the day I was there the other speaker was Hugo Duminil-Copin and he was talking about self-avoiding random walks, a topic which I knew very little about. Now I was faced with (at least superficially) similar ideas in the context of DNA recombination. At the end of his talk Bader spent a few minutes on a quite different topic, namely bistability in the state of M. tuberculosis. I would have liked to have heard more about that. He is collaborating on this together with Denise Kirschner whose work on modelling tuberculosis I have discussed in a previous post.

This meeting had the advantages of relatively small conferences (in this case of the order of 50 participants) and has served the purpose of opening new perspectives for me.

### SIAM Conference on the Life Sciences in Charlotte

August 7, 2014This week I have been attending the SIAM Conference on the Life Sciences in Charlotte. Here I want to mention some highlights from my personal point of view. First I will mention some of the plenary talks. John Rinzel talked about mathematical modelling of certain perceptual phenomena. We are all familiar with the face-vase picture which switches repeatedly between two forms. I had never considered the question of trying to predict how often the picture switches. Rinzel presented models for this and for other related auditory phenomena which he demonstrated in the lecture. I find it remarkable that such apparently subjective phenomena can be brought into such close connection with precise mathematical models. Kristin Swanson talked about her work on modelling the brain cancer known as glioma and its various deadly forms. I had heard her talk on the same theme at the meeting of the Society for Mathematical Biology in Dundee in 2003. Of course there has been a lot of progress since then. This was long before I started this blog but if the blog had existed I would certainly have written about the topic. I will not try to resurrect the old stories from that distant epoch. Instead I will just say that Kristin is heavily involved in using computer simulations to optimize the treatment (surgery, radiotherapy, chemotherapy) of individual patients. One of the main points in her talk this week is that it seems to be possible to divide patients into two broad categories (with nodular or diffuse growth of the tumour) and that this alone may have important implications for therapeutic decisions. Oliver Jensen talked about a multiscale model for predicting plant growth, for instance the way in which a root manages to sense gravity and move downwards. This involves some very sophisticated continuum mechanics which the speaker illustrated by everyday examples in a very effective and sometimes humorous way. The talk was both impressive and entertaining. Norman Mazer talked about the different kinds of cholesterol (LDL, HDL etc.). According to what he said lowering LDL levels is an effective means for avoiding risks of cardiovascular illness but the alternative strategy of raising HDL levels has not been successful. He explained how mathematical modelling can throw light on this phenomenon. My understanding is that the link between high HDL level and lower cardiovascular risks is a correlation and not a sign of a causal influence of HDL level on risk factors. The last talk was by James Collins, a pioneer of synthetic biology. The talk was full of good material, both mathematical and non-mathematical. Maybe I should invest some time into learning about that field.

There was one very interesting subject which was not the subject of a talk at the conference (at least not of one I heard – it was briefly referred to in the talk of Collins mentioned above) but was a subject of conversation. It is a paper called ‘Paradoxical Results in Perturbation-Based Signaling Network Reconstruction’ by Sudhakaran Prabakaran, Jeremy Gunawardena and Eduardo Sontag which appeared in Biophys. J. 106, 2720. It suggests that the ways in which biologists deduce the influence of substances on each other on the basis of experiments are quite problematic. The mathematical content of the paper is rather elementary but its consequences for the way in which theoretical ideas are applied in biology may be considerable. The system studied in the paper is an in vitro reconstruction of part of the MAP kinase cascade and so not so far from some of my research.

Among the parallel sessions those which were most relevant for me were one entitled ‘Algebra in the Life Sciences’ and organized by Elisenda Feliu, Nicolette Meshkat and Carsten Wiuf and one called ‘Developments in the Mathematics of Biochemical Reaction Networks’ organized by Casian Pantea and Maya Mincheva. My talk was in the second of these. These sessions were very valuable for me since they allowed me to meet a considerable number of people working in areas close to my own research interests, including several whose papers were well known to me but whom I had never met. I think that this will bring me to a new level in my work in mathematical biology due to the various interactions which took place. I will not discuss the contents of individual talks here. It is rather the case that what I learned form them will flow into my research effort and hence indirectly influence future posts in this blog. I feel that this conference has gained me entrance into a (for me) new research community which could be the natural habitat for my future research. I am very happy about that. The whole conference was an enjoyable and stimulating experience for me. I noticed no jet lag at all but I must be suffering from a lack of sleep due to the fact that the many things going on here just did not leave me the eight hours of sleep per night I am used to.

### Breaking waves in Madrid

July 19, 2014Last week I was at the 10th AIMS Conference on Dynamical Systems, Differential Equations and Applications in Madrid. It was very large, with more than 2700 participants and countless parallel sessions. This kind of situation necessarily generates a somewhat hectic atmosphere and I do not really like going to that type of conference. I have heard the same thing from many other paople. There is nevertheless an advantage, namely the possibility of meeting many people. To do this effectively it is necessary to proceed systematically since it is easy to go for days without seeing a particular person of interest. This aspect was of particular importance for me since I am still at a relatively early stage in the process of entering the field of mathematical biology and I have few contacts there in comparison to my old field of mathematical relativity. In any case, the conference allowed me to meet a lot of interesting people and learn a lot of interesting things.

The first plenary talk, by Charles Fefferman, was on a subject related to a topic I was interested in many years ago. I learned that a lot has happened since I last thought about this. The attempt to model a body of fluid with a free surface leads to considerable mathematical difficulties. When I started working on dynamical models for this kind of situation few people seemed to be interested in proving theorems on the subject. The source of my interest in the subject was the influence of Jürgen Ehlers, who always had a clear vision of what were the important problems. In this way I found myself in the position of a pioneer in a certain research area. Being in that situation has the advantage of not being troubled by strong competition. On the other hand it can also mean that whatever you achieve can be largely ignored and it is not the best way to get wide recognition. Often finishing mathematical research directions gets more credit than starting them. This could no doubt be compensated by suitable advertizing but that was never my strong point. This is a configuration which I have often found myself in and in fact, comparing advantages and disadvantages, I do not feel I need to change it. Coming back to the fluids with free surface, this is now a hot topic and played a prominent role at the conference. When I was working on this the issue of local existence in the case of inviscid fluids was still open. A key step was the work of Sijue Wu on water waves. I learned from the talk of Fefferman that this has been extended in the meantime to global existence for small data. The question which is now the focus of interest is formation of singularities (i.e. breakdown of classical solutions) for large data. Instead of considering the breaking of one wave the idea is to consider two waves which are approaching each other while turning over until they meet. There are already analytical results on parts on this process by Fefferman and collaborators and they plan to extend this to a more global picture by using a computer-assisted proof. Another plenary was by Ingrid Daubechies, who talked about applications of image processing to art history. I must admit that beforehand the theme did not appear very attractive to me but in fact the talk was very entertaining and I am glad I went to hear it.

I gave a talk on my recent work with Juliette Hell on the MAPK cascade in a session organized by Bernold Fiedler and Atsushi Mochizuki. I found the session very interesting and the highlight for me was Mochizuki’s talk on his work with Fiedler. The subject is how much information can be obtained about a network of chemical reactions by observing a few nodes, i.e. by observing a few concentrations. What I find particularly interesting are the direct connections to biological problems. Applied to the gene regulatory network of an ascidian (sea squirt) this theoretical approach suggests that the network known from experimental observations is incomplete and motivates searching for the missing links experimentally. Among the many other talks I heard at the conference, one which I found particularly impressive concerned the analysis of successive MRT pictures of patients with metastases in the lung. The speaker was using numerical simulations with these pictures as input to provide the surgeon with indications which of the many lesions present was likely to develop in a dangerous way and should therefore be removed. One point raised in the talk is that it is not really clear what information about the tissue is really contained in an MRT picture and that this could be an interesting mathematical problem in itself. In fact there was an encouragingly (from my point of view) large number of sessions and other individual talks at the conference on subjects related to mathematical biology.

The conference took place on the campus of the Universidad Autonoma somewhat outside the city. A bonus for me was hearing and seeing my first bee-eater for many years. It was quite far away (flying high) but it gave me real pleasure. I was grateful that the temperatures during the week were very moderate, so that I could enjoy walking through the streets of Madrid in the evening without feeling disturbed by heat or excessive sun.

### The Higgins-Selkov oscillator

May 14, 2014In a previous post I wrote about glycolytic oscillations and mentioned a mathematical model for them, the Higgins-Selkov oscillator. Higgins introduced this as a chemical model while Selkov also wrote about some mathematical aspects of modelling this system. When I was preparing my course on dynamical systems I wanted to present an example where the existence of periodic solutions can be concluded by using the existence of a confined region in a two-dimensional system and Poincare-Bendixson theory. An example which is frequently treated in textbooks is the Brusselator and I wanted to do something different. I decided to try the Higgins-Selkov oscillator. Unfortunately I came up against difficulties since that model has unbounded solutions and it is hard to show that there are any bounded solutions except a stationary solution which can be calculated explicitly. For the purposes of the course I went over to considering the Schnakenberg model, a modification of the Higgins-Selkov oscillator where it is not hard to see that all solutions are bounded.

More recently I decided to try to finally find out what happens with the Higgins-Selkov oscillator itself. Reading Selkov’s paper I originally had the impression that he had proved the essential properties of the solutions. This turned out to be mistaken. One obstacle for me was that Selkov cited a theorem from a famous Russian textbook of Andronov et. al. and I did not know what the theorem was. An English translation of the book exists in the university library here but since Selkov only cites a page number I did not know how to find the theorem. I was able to get further when Jan Fuhrmann got hold of a copy of the page in question from the Russian original. This page has an easily identifiable picture on it and this allowed me to identify the corresponding page of the English translation and hence the theorem. I found that, as far as it is applicable to the oscillator problem this was something I could prove myself by a simple centre manifold argument. Thus I realized that the results quoted by Selkov only resolve some of the simpler issues in this problem.

At this stage I decided to follow the direction pointed out by Selkov on my own. The first stage, which can be used to obtain information about solutions which tend to infinity, is to do a Poincare compactification. This leads to a dynamical system on a compact subset of Euclidean space. In general it leads to the creation of new stationary points on the boundary which are not always hyperbolic. In this particular example two new stationary points are created. One of these has a one-dimensional centre manifold and it is relatively easy to determine its qualitative nature. This is what Selkov was quoting the result of Andronov et. al. for. The other new stationary solution is more problematic since the linearization of the system at that point is identically zero. More information can be obtained by transforming to polar coordinates about that point. This creates two new stationary points. One is hyperbolic and hence unproblematic. The linearization about the other is identically zero. Passing to polar coordinates about that point creates three new stationary points. One of them is hyperbolic while the other two have one-dimensional centre manifolds. The process comes to an end. When trying this kind of thing in the past I was haunted by the nightmare that the process would never stop. Is there a theorem which forbids that? In any case, in this example it is possible to proceed in this way and determine the qualitative behaviour near all points of the boundary. The problem is that this does not seem to help with the original issue. I see no way in which, even using all this information, it is possible to rule out that every solution except the stationary solution tends to infinity as tends to infinity.

Given that this appeared to be a dead end I decided to try an alternative strategy in order to at least prove that there are some parameter values for which there exists a stable periodic solution. It is possible to do this by showing that a generic supercritical Hopf bifurcation occurs and I went to the trouble of computing the Lyapunov coefficient needed to prove this. I am not sure how much future there is for the Higgins-Selkov oscillator since there are more modern and more complicated models for glycolysis on the market which have been studied more intensively from a mathematical point of view. More information about this can be found in a paper of Kosiuk and Szmolyan, SIAM J. Appl. Dyn. Sys. 10, 1307.

Finally I want to say something about the concept of feedback, something I find very confusing. Often it is said in the literature that oscillations are related to negative feedback. On the other hand the oscillations in glycolysis are often said to result from positive feedback. How can this be consistent? The most transparent definition of feedback I have seen is the one from a paper of Sontag which I discussed in the context of monotone systems. In that sense the feedback in the Higgins-Selkov oscillator is definitely negative. An increase in the concentration of the substrate leads to an increase in the rate of production of the product. An increase in the concentration of the product leads to an increase of the rate of consumption of the substrate. The combination of a positive and a negative sign gives a negative loop. The other way of talking about this seems to be related to the fact that an increase in the concentration of the product leads to an increase in the reaction rate between substrate and product. This is consistent with what was said before. The difference is what aspects of the system are being regarded as cause and effect, which can lead to a different assignment of the labels positive and negative. The problem as I see it is that feedback is frequently invoked but rarely defined, with the implicit suggestion that the definition should be obvious to anyone with an ounce of understanding. I seem to be lacking that ounce.

### Proofs of dynamical properties of the MAPK cascade

April 3, 2014The MAP kinase cascade, which I mentioned in a previous post, is a biochemical network which has been subject to a lot of theoretical and experimental study. Although a number of results about mathematical models for this network have been proved, many widely accepted results are based on numerical and/or heuristic approaches. Together with Juliette Hell we set out to extend the coverage of rigorous results in this area. Our first results on this can be found in a paper we just posted on q-bio.

The system of equations which is fundamental for this work is that of Huang and Ferrell discussed in my previous post on the subject. I call it the MM-MA system (for Michaelis-Menten via mass action). When this system can be reduced to a smaller system by means of a quasistationary approximation the result will be called the MM system (for Michaelis-Menten) (cf. this post). With a suitable formulation the MM system is a singular limit of the MM-MA system. The MAPK cascade consists of three coupled layers. The first main result of our paper concerns the dual futile cycle, which can be thought of as the second layer of the cascade in isolation (cf. this post). We proved that the MM system for the dual futile cycle exhibits a generic cusp bifurcation and hence that for suitable values of the parameters there exist two different stable stationary solutions (bistability). Using the fact that this system is a singular limit of the system arising from the MM-MA description of the same biological system we then used geometric singular perturbation theory (cf. this post) to conclude that the MM-MA system also shows bistability.

The second main result concerns the system obtained by truncating that of Huang-Ferrell by keeping only the first two layers. It is subtle to find a useful quasistationary approximation for this system and we were put on the right track by a paper of Ventura et. al. (PLoS Comp. Biol. 4(3):e1000041). This allowed us to obtained an MM system which is a limit of the MM-MA system in a way which allows geometric singular perturbation theory to be applied. This leads to the following relative statement: if the MM system for the truncated MAPK cascade has a hyperbolic periodic solution then the same is true for the MM-MA system. To get an absolute statement it remains to prove the existence of periodic solutions of the MM system, which in this case is of dimension three. That there are solutions of this kind is indicated by numerical work of Ventura et. al.