## Archive for the ‘dynamical systems’ Category

### Reaction networks in Copenhagen

July 9, 2015

Last 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 $\delta$ 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 $\delta'$. 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 $\delta\ge\delta'$, 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, 2015

In 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 $F_i$ 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 $s_i$. Thus if the time derivatives of the $s_i$ can be expressed in terms of the $F_i$ we are done. It is shown in the appendix to the paper how a linear combination of the time derivatives of the $s_i$ with coefficients only depending on the $s_i$ is equal to $F_i$. Moreover it is stated that the time derivatives of the $s_i$ 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, 2015

Photosynthesis 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, 2015

In 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 $i$ to vertex $j$ 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 $i$ to reaction $j$ if species $i$ is present on the LHS of reaction $j$. There is a directed edge from reaction $i$ to species $j$ if the net production of species $j$ in reaction $i$ is non-zero. The first type of edge is labelled by the partial derivative of flux $j$ with respect to species $i$. 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.

### Proof of the global attractor conjecture

May 14, 2015

In a previous post I discussed the global attractor conjecture which concerns the asymptotic behaviour of solutions of the mass action equations describing weakly reversible chemical reaction networks of deficiency zero (or more generally complex balanced systems). Systems of the latter class are sometimes called toric dynamical systems because of relations to the theory of toric varieties in algebraic geometry. I just saw a paper by Gheorghe Craciun which he put on ArXiv last January (arxiv.org/pdf/1501.0286) where he proves this conjecture, thus solving a problem which has been open for more than 40 years. The result says that for reaction networks of this class the long-time behaviour is very simple. For given values of the conserved quantities there is exactly one positive stationary solution and all other solutions converge to it. What needs to be done beyond the classical results on this problem is to show that a positive solution can have no $\omega$-limit points where some concentration vanishes. This property is sometimes known as persistence.

A central concept used in the proof of the result is that of a toric differential inclusion. This says that the time rate of change of the concentrations is contained in a special type of subset. The paper contains a lot of intricate geometric constructions. These are explained consecutively in dimensions one, two, three and four in the paper. This should hopefully provide a good ladder to understanding.

### The existence proof for Hopf bifurcations

September 22, 2014

In a Hopf bifurcation a pair of complex conjugate eigenvalues of the linearization of a dynamical system $\dot x=f(x,\alpha)$ at a stationary point pass through the imaginary axis. This has been discussed in a previous post. Often textbook results (e.g. Theorem 3.3 in Kuznetsov’s book) concentrate on the generic case where two additional conditions are satisfied. One of these is that the first Lyapunov coefficient is non-vanishing. The other is that the eigenvalues pass through the imaginary axis with non-zero velocity. The existence of periodic solutions can be obtained if only the second of these conditions are satisfied. This was already included in the original paper of Hopf in 1942. Hopf states his results only in the case of analytic systems but this should perhaps be seen as a historical accident. A similar result holds with mucher weaker regularity assumptions. It is proved under the assumption of $C^2$ dependence on $x$ and $C^1$ dependence on $\alpha$ in Hale’s book on ordinary differential equations. This has consequences for the case where the second genericity assumption is not satisfied. Let $\lambda$ be an eigenvalue which passes through the imaginary axis for $\alpha=0$ and suppose that the derivatives of ${{\rm Re}\lambda}$ with respect to $\alpha$ vanish up to order $2k$ for an integer $k$ but that the derivative of order $2k+1$ does not vanish. Then it is possible to replace $\alpha$ by $\alpha^{2k+1}$ as parameter and after this change the second genericity assumption is satisfied. Even if the original right hand side was analytic in $\alpha$ the transformed right hand side is in general not $C^2$. It is, however, $C^1$ and so the version of the theorem in Hale’s book applies to give the existence of periodic solutions. This theorem applies to a two-dimensional system but it then also evidently applies in general by a centre manifold reduction.

The theorem is proved as follows. The problem is transformed to polar coordinates $(\rho,\theta)$ and then $\rho$ is written as a function of $\theta$. In this way a non-autonomous scalar equation with $2\pi$-periodic coefficients is obtained and the aim is to find a $2\pi$-periodic solution. The first step is to reformulate the task as a fixed-point problem with the property that if a fixed point is periodic it will be a solution of the original problem. Then it is shown using the Banach fixed point theorem(in a minor variant of the local existence theorem for ODE using Picard iteration) that there always exists a fixed point depending on a certain new parameter. This fixed point is only periodic if the result of substituting it into the right hand side of the original equation has mean value zero. This condition can be written as $G(\alpha,a)=0$. Applying the implicit function theorem to $G$ shows the existence of a solution of $G(\alpha(a),a)=0$ for $a$ small. This completes the proof.

Summing up, there are two types of theorem about Hopf bifurcation, a ‘coarse’ theorem of the type just sketched with weak hypotheses and a weak but still very interesting conclusion and a ‘fine’ theorem which gives stronger conclusions but needs a stronger hypothesis (non-vanishing of the Lyapunov coefficient and its sign). In his original paper Hopf proved both types. Are there also ‘rough’ versions of theorems about other bifurcations?

### SIAM Conference on the Life Sciences in Charlotte

August 7, 2014

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

### The Higgins-Selkov oscillator

May 14, 2014

In 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 $t$ 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, 2014

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

### Lyapunov exponents

March 25, 2014

In the literature on chaos and strange attractors one concept which plays a prominent role is that of Lyapunov exponents. I came across this repeatedly but never understood the definition. I think I have now understood the reason that I had problems. The fact that I have spent so much time with differential geometry in the past years sometimes makes me see the mathematical world through ‘differential geometric spectacles’. I felt that dynamical systems were objects which naturally live on smooth manifolds and that the definition of important concepts should not be dependent on coordinates or the presence of a preferred metric. The usual definitions of Lyapunov exponents appeared to me strongly tied to Euclidean space although I had seen a couple of comments in the literature (without further justification) that the definition was coordinate independent. In what follows I will give a manifestly coordinate independent definition. Incidentally, the definition of Lyapunov exponents is not invariant under changes of the time coordinate and at one time this led to confusion among people studying chaos in cosmological models. The source of the problem was the lack of a preferred time coordinate in general relativity.

Consider a dynamical system defined by a smooth vector field on a manifold $M$. Let $A$ be a compact subset of the manifold which is invariant under the flow $\phi$ generated by the vector field. The aim here is to define the maximum Lyapunov exponent of a point $x_0\in A$. The derivative of the flow, $J_t=D_x\phi (t,x_0)$ is a linear mapping from $T_{x_0}M$ to $T_{\phi (t,x_0)}M$. In the Euclidean space picture $J_t$ is treated as a matrix and this matrix is multiplied by its transpose. What is this transpose in an invariant setting? It could be taken to be the mapping from $T^*_{\phi (t,x_0)}M$ to $T^*_{x_0}M$ naturally associated to $J_t$ by duality. The product of the matrices could be associated with the composition of the linear mappings but unfortunately the domains and ranges do not match. To overcome this I introduce a Riemannian metric $g$ on a neighbourhood of $A$. It is then necessary to show at the end of the day that the result does not depend on the metric. The key input for this is that since $A$ is compact the restrictions of any two metrics $g_1$ and $g_2$ to $A$ are uniformly equivalent. In other words, there exists a positive constant $C$ such that $C^{-1}g_1(v,v)\le g_2(v,v)\le Cg_1(v,v)$ for all tangent vectors $v$ at points of $A$. Once the metric $g$ has been chosen it can be used to identify the tangent and cotangent spaces with each other at the points $x_0$ and $\phi (t,x_0)$ and thus to compose $J_t$ and its ‘transpose’ to get a linear mapping $B(t)$ on the vector space $T_{x_0}M$. This vector space does not depend on $t$. The eigenvalues $\lambda_i (t)$ of the mapping $B(t)$ are easily shown to be positive. The maximum Lyapunov exponent is the maximum over $i$ of the limes superior for $t\to\infty$ of $\frac{1}{t}$ times the logarithm of $\sqrt{\lambda_i (t)}$. Note that the ambiguity of a multiplicative constant in the definition of $B(t)$ becomes an ambiguity of an additive constant in the definition of the logarithms and because of the factor $\frac{1}{t}$ this has no effect on the end result.

In general if the maximum Lyapunov exponent at a point $x_0$ is positive this is regarded as a sign of instability of the solution starting at that point (sensitive dependence on initial conditions) and if the exponent is negative this is regarded as a sign of stability. Unfortunately in general these criteria are not reliable, a fact which is known as the Perron effect. This is connected with the question of reducing the study of the asymptotic behaviour of a non-autonomous linear system of ODE to that of the autonomous systems obtained by freezing the coefficients at fixed times.