Modelling the Calvin cycle

Some years ago the Max Planck Institute for Molecular Plant Physiology organized a conference on metabolic networks. I decided to see what was going on in the institute next to the one where I work and I went to some of the talks. The one which I found most interesting was by Zoran Nikoloski. His subject was certain models for the Calvin cycle, which is part of photosynthesis. A motivating question was whether photosynthesis can work in two different stable steady states. If that were the case it might be possible to influence the plant to move from one state to another and, in the best case, to increase its production of biomass. This is of interest for biotechnology. Mathematically the question is that of multistationarity, i.e. whether a system of evolution equations admits more than one stationary solution. Beyond this it is of interest whether there can be more than one stable stationary solution. In fact in this context the issue is not that of absolute uniqueness of stationary solutions but of uniqueness within a given stoichiometric compatibility class. This means that the solution is unique when certain conserved quantities are fixed. One thing I found attractive about the presentation was that the speaker was talking about rigorous mathematical results on the dynamics and not just about numerically calculating a few solutions.

If the system is modelled deterministically and diffusion is neglected there results a system of ordinary differential equations for the concentrations of the substances involved as functions of time. It is necessary to choose which substances should be included in the description. In a basic model of the Calvin cycle there are five substances. In the work discussed in the talk of Nikoloski and in a paper he wrote with Sergio Grimbs and others (Biosystems 303, 212) various ODE systems based on this starting point are considered. They differ by the type of kinetics used. They consider mass action kinetics (MA), extended Michaelis-Menten kinetics where the enzymes catalysing the reactions are included explicitly (MM-MA) and effective Michaelis-Menten (MM) obtained from the system MM-MA by a singular limit. The systems MA and MM consist of five equations while the system MM-MA consists of nineteen equations. In the paper of Grimbs et. al. they show among other things that the system MM never admits a stable stationary solution, whatever the reaction constants, while the system MM-MA can exhibit two different stationary solutions.

After the talk I started reading about this subject and I also talked to Nikoloski about it. Later I began doing some research on these systems myself. Some technical difficulties which arose (which I wrote about in a previous post) led me to consult Juan Velázquez and he joined me in this project. Now we have written a paper on models for the Calvin cycle. In a case where there is only one stationary solution and it is unstable it is of interest to consider the final fate of general solutions of the system. For some initial conditions the concentrations of all substances tend to zero at late times. For other data (a whole open set) we were able to show that all concentrations tend to infinity as . We called the latter class runaway solutions. These do not seem to be of direct biological relevance but they might be helpful in choosing between alternative models which are more or less appropriate. The proof of the existence of runaway solutions for the MA system is somewhat complicated since this turns out to be a system with two different timescales. The system MM-MA also admits runaway solutions. Although the system is larger than MA the existence proof is simpler and in fact can be carried out in the context of a larger class of systems. Runaway solutions are also found for the system MM.

In the paper of Grimbs et. al. one system is considered which includes the effect of diffusion. Restricting to homogeneous solutions of this system gives a system of ODE called MAdh which is different from the system MA. The difference is that while the concentration of ATP is a dynamical variable in MAdh it is taken to be constant in MA. We showed that the system MAdh has zero, one or two solutions depending on the values of the parameters and that all solutions are bounded. Thus runaway solutions are ruled out. Intuitively this is due to the fact that the supply of energy is bounded but this heuristic argument is far from providing a proof. There are many other models of the Calvin cycle in the literature. In general they consider the reactions between a larger class of substances. It is an interesting task for the future to extend the results obtained up to now to these more general models. This post has been very much concerned with the mathematics of the problem and has not said much about the biology. The reactions making up the Calvin cycle were determined experimentally by Melvin Calvin and I can highly recommend his Nobel lecture as a description of how this was achieved

Absolute concentration robustness

In the past years I have been on the committees for many PhD examinations. A few days ago, for the first time, I was was on the committee for a thesis on a subject belonging to the area of mathematical biology. This was the thesis of Jost Neigenfind and it was concerned with a concept called absolute concentration robustness (ACR).

The concentration of a given substance in cells of a given type varies widely between the individual cells. (Cf. also this previous post). It is of interest to identify mechanisms which can ensure that the steady state concentration of a particular substance is independent of initial data. (This is a way in which the output of a system can be independent of background variation.) In saying this I am assuming implicitly that more general solutions converge to steady states. A more satisfactory formulation can be obtained as follows. In a chemical reaction network there are usually a number of conserved quanitities, say . These define affine subspaces of the state space, the stoichiometric compatibility classes. For many systems there is exactly one stationary solution in each stoichiometric compatibility class. The condition of interest here is that the value of one of the concentrations, call it , in the steady state solution is independent of the parameters . (The other concentrations will in general depend on the .) This property is ACR. I first heard of this in a talk by Uri Alon at the SMB conference in Krakow in the summer of 2011. The basic idea is explained clearly in a paper of Shinar and Feinberg (Science 327, 1389). They present a general theoretical approach but also describe some biological systems where ACR (in a suitable approximate sense) has been observed experimentally. In the terminology of Chemical Reaction Network Theory (CRNT) the examples they discuss have deficiency one. They mention that ACR is impossible in systems of deficiency zero. There is no reason why it should not occur in systems of deficiency greater than one but in those more complicated dynamics make it more difficult to decide whether the property holds or not.

The result of Shinar and Feinberg only covers a class of reaction networks which is probably very restricted. What Neigenfind does in his thesis is to develop more general criteria for ACR and computer algorithms which can check these criteria for given systems. The phenomenon of ACR is interesting since it is a feature which may be more common in reaction systems coming from biology than in generic systems. At least there is a good potential reason why this might be the case.

Talk on mathematical modelling in Karlstad

Yesterday I was in Karlstad in Sweden to give a talk on the uses of mathematical modelling in the natural sciences. I was invited to do this by Claes Uggla and I was very happy to have the opportunity to present some of my ideas on this subject. The talk was structured as a series of examples involving applications of different mathematical techniques. Many of these examples have been discussed in some form in this blog during the past few years and indeed a lot of my ideas on the subject were developed in conjunction with the blog posts. The subjects were William Harvey and the circulation of the blood, multidrug therapy for HIV-AIDS, the lizard Uta stansburiana, oscillations near the big bang, Liesegang rings, modelling oscillations in vole populations using a reaction-diffusion system, signal transduction in T cells.

As well as presenting a variety of applications of different types of mathematics I also wanted to explain some mathematical connections between these subjects. One central idea is that structural stability is an issue of key importance in modelling natural phenomena. Most phenomenological models involve parameters or other elements which are not known exactly. Thus to be of interest for applications features of the dynamics of the model should be invariant under arbitrary small perturbations of the system. More precisely, if a model does not possess an invariance of this type but is nevertheless useful this requires some explanation. One possible source of an explanation is the presence of what I call ‘absolute elements’ in the model. For instance, in population dynamics if a population is zero at some time then it will definitely remain zero. This fact is independent of the details of how the population grows when it is non-zero. Similarly a spacetime singularity can define an absolute element in cosmology. When the spacetime metric breaks down this ends the dynamics in a way which is independent of the details of the dynamics of the matter away from the singularity. Thus structural stability can be weakened to the condition of invariance under small perturbations which leave certain submanifolds fixed. This can lead to the appearance of relevant heteroclinic cycles although these are not structurally stable in the absolute sense. It explains the appearance of heteroclinic cycles in the models for lizards and for the big bang in a unified way. In a similar way, restricting the perturbations of a system of chemical reactions to those which leave a particular reaction irreversible can furnish the homoclinic orbit needed to model Liesegang rings.

I have now put a slightly extended version of this talk with references on my web page. On the same day there was a talk by Bernt Wennberg on models for the collective motion of birds and fish, concentrating mainly on kinetic models related to the Boltzmann equation. At the start of his talk he showed some of the well-known pictures of flocks of Starlings over Rome. In the evening I had my own pleasant experience with a flock of birds. A large number of Jackdaws (a couple of hundred) were flying around the central square in Karlstad and calling. For some reason I have become increasingly attached to the Jackdaw over the years. At this point, and without a good excuse, I want to tell a story about Jackdaws from the book ‘King Solomon’s Ring’ by Konrad Lorenz. It is a long time since I read the book and so I hope I do not distort the story too much. At one time Lorenz was living in a small village in Austria where he was regarded by the locals as a bit crazy. One of his interests was the social life of Jackdaws. There were Jackdaws living on the roofs of the houses and he climbed up to get close to them. In order to fit in better with his black subjects he decided to dress in black. The only ‘suitable’ black clothing he could find was a devil’s costume left over from a fancy dress party. No doubt the spectacle of him climbing over the roofs dressed as the devil perfected his reputation with the local inhabitants.

A simple system with two different timescales

While struggling with the proof of a result for a model for photosynthesis (which I intend to report on in more detail at a later date) I decided to apply the following principle which I heard about as an undergraduate and is apparently due to Paul Halmos: if there is a mathematical problem you cannot solve then there is also a simpler problem you cannot solve. With this in mind I developed a model problem for what I wanted to do which I still could not solve. In the meantime, with the help of some key input from Juan Velázquez, I can solve the model problem and I will report on that here. It seems that the right conceptual framework for this is that of systems with two different timescales. Consider two functions and which satisfy and . The fact that the power here is exactly five is not so important. This was the power that came up in the problem I was originally interested in and I wanted to rule out any misleading simplifications which might have arisen by replacing it by the power two, for instance. The first equation can be solved explicitly in the form .This reduces the original system to the scalar, but non-autonomous, equation . The question of interest is whether this equation has solutions which tend to infinity as and if so how many of these are there. A guess at the asymptotics of a solution of this kind which is formally consistent is . It is possible to take this further by looking for formal power series solutions where is a linear combination of integer powers of . It turns out that there is a solution of this type in the sense of formal power series. In other words, substituting the expression in the equations and comparing coefficients gives a consistent answer. The coefficients are determined uniquely. This means that if there is more than one solution with this type of asymptotics these coefficients cannot distinguish between them and hence cannot be used to parametrize them. In fact there is a one-parameter family of solutions having this type of asymptotics and the difference between any two of these is of order for a constant . This means that while these solutions decay on a timescale the difference between them decays on a timescale which is exponentially faster. I am reminded of the term ‘asymptotics beyond all orders’ which I have heard occasionally but I do not know exactly what that means.

How can these results be proved? First introduce a quantity since this can be expected to decay very fast. The evolution equation for can be rewritten as an evolution equation for of the form , where the quantity will eventually be small. This equation can be solved by variation of constants to give an integral equation for . It contains a parameter which distinguishes the different solutions. The problem of solving the integral equation can be reformulated as a fixed point problem for a mapping . The key step is to show that, for sufficiently large and sufficiently small, maps a suitable set of functions with a certain decay property to itself. The fixed point can then be obtained using the Arzela-Ascoli theorem. The simple system considered here can presumably be treated by simpler methods. My motivation for discussing it here that the technique of proof is of much wider applicability and the conclusions obtained are a model for other problems.

The multiple futile cycle

The multiple futile cycle is a simple type of network of chemical reactions which is often found in biological systems. In a previous post I mentioned it as a component of a slightly more complicated network found in many cells, the MAP kinase cascade. One concrete realization of the multiple futile cycle is a protein which can be phosphorylated at up to sites. All the phosphorylation steps are carried out by one kinase while all dephosphorylation steps are carried out by one phosphatase. Each step is modelled in the Michaelis-Menten way, including an enzyme-substrate complex as one of the species and using mass action kinetics. There results a system of ordinary differential equations with three conservation laws. These represent the conservation of the total amount of the two enzymes and of the substrate protein. In the case , which might be called the simple futile cycle, using the conservation laws to eliminate some of the variables leads to a three-dimensional dynamical system. A basic question is what can be said about the dynamics of solutions of this system.

It has been shown by Angeli and Sontag (Nonlinear Analysis RWA, 9, 128) that in the case every solution converges to a stationary solution and that this stationary solution is unique for given values of the conserved quantities. The proof uses the theory of monotone dynamical systems. The original dynamical system is not monotone and so the first step in their proof is to replace it by another system which is monotone and show that convergence properties of solutions of the second imply convergence properties of solutions of the first. The second step is to prove the convergence of solutions of monotone systems under the additional condition of the existence of a translational symmetry. The paper mentions that this result is dual to a previously known result due to Mierczyński about monotone systems with a conserved quantity. Up to now I thought that the only benefit of knowing that a dynamical system is monotone is the possibility of reducing it to an effective system of one dimension less. This is only interesting if the initial system is of dimension no more than three. What this work has shown me is that knowing that a system is monotone can sometimes be the key to concluding much more. One aspect of the paper of Angeli and Sontag which was a source of confusion for me was a difference in conventions to what I am familiar with from chemical reaction network theory. This seems to be essential for the monotonicity argument and not just a matter of taste. The stoichiometric matrix (or stoichiometry matrix) is defined differently because a reversible reaction is treated as a single reaction rather than as a pair.. I feel a spontaneous preference for the CRNT convention but here it seems that a different one can be a real advantage. In the case of the simple futile cycle an important effect is that the dimension of the kernel of the stoichiometric matrix is three with the CRNT convention and one with the Angeli-Sontag convention.

In another paper (J. Math. Biol. 61, 581) Angeli, De Leenheer and Sontag present a more general theory related to this. Here the hypotheses needed to obtain a suitable monotone system involve the properties of certain graphs constructed from the reaction network. In this theory the notion of persistence of the dynamical system plays an important role. This is the property that a positive solution can never have any limit points on the boundary of the positive region. The case (dual futile cycle) has been considered in a paper of Wang and Sontag (J. Nonlin. Sci. 18, 527). There they are able to show that for certain ranges of the parameters generic solutions converge to stationary solutions. To emphasize the power of the techniques developed in these papers it should be pointed out that they can be applied to systems with arbitarily large numbers of unkowns and parameters and that when they apply they give strong conclusions.

D. Angeli and E. D. Sontag (2008). Translation-invariant monotone systems, and a global convergence result for enzymatic futile cycles Nonlinear Analysis: Real World Applications, 9 DOI: 10.1016/j.nonrwa.2006.09.006

D. Angeli, P. De Leenheer and E. D. Sontag (2010). Graph-theoretic characterizations of monotonicity of chemical networks in reaction coordinates. Journal of mathematical biology, 61 (4) PMID: 19949950
 

The Goodwin oscillator

During a talk by Jae Kyoung Kim at last week’s SMB meeting the speaker showed a system of equations and called it ‘a system you all know’. This revealed to me a gap in my knowledge of mathematical biology. The system is the Goodwin oscillator. It is described in Murray’s book on mathematical biology and I am sure I have read the relevant section on some level. This just shows that there a big difference between reading something and understanding its significance and being able to situate in a wider context. Now I have done my homework on this and I will write something about it here. The system, in the form it is given in Murray’s book, is an example of systems of the form . Here the labels on the are supposed to be interpreted modulo . In other words there are equations and is interpreted as . In the Goodwin model itself and the functions are linear for . The function is equal to for constants , and . This function is positive and its derivative is negative. Thus it can be interpreted as representing a negative feedback on the production of . In the context in which it was introduced by Goodwin the quantities , and represent concentrations of mRNA, the enzyme it codes for and the product of a reaction it catalyzes. The substrate of the enzyme is assumed present at a constant level and is not modelled explicitly.

It is known that the system admits a periodic solution if the Hill coefficient is greater than eight and not otherwise. Since this number is considered unrealistically large for the application which inspired the model modifications of this have been considered where periodic solutions can be obtained for lower values of . It is proved by Hastings, Tyson and Webster (J. Diff. Eq. 25, 39) that for the Goodwin system and a larger class of similar models the following is true. The system has a unique steady state and if the linearization of the system at that point has no repeated eigenvalues and at least one eigenvalue with positive real part there exist periodic solutions. This reduces the existence question to the analysis of the linearization. The existence proof relies on the Brouwer fixed point theorem and is similar to a proof I described in a previous post. Although the Goodwin system is three dimensional the method is not restricted to that case. The proof does not give information about the stability of the periodic solutions. In the paper of Hastings et. al. they indicate that an alternative analysis using a Hopf bifurcation can give stability in some cases. However no details of the stability argument are given in that paper.

The Goodwin model was inspired by the fundamental work of Monod and Jacob on gene regulation.  Various things have given me an appetite for learning more about gene regulatory networks and this was increased by some of the talks I heard last week.

Stability of heteroclinic cycles

A heteroclinic orbit is a solution of a dynamical system which converges to one stationary solution in the past and to another stationary solution in the future. A heteroclinic chain is a sequence of heteroclinic orbits where the past limit of each orbit is the future limit of the preceding one. If this sequence is periodic we get what is called a heteroclinic cycle. Given such an object it is of interest to ask about its stability. For an initial datum sufficiently close to the cycle, when does the corresponding solution converge to the cycle at late times? In particular, when is the -limit set of the solution of interest equal to the entire cycle? To obtain information about this question it is useful to consider the linearization of the system about the vertices of the cycle. For a solution of the kind we are looking for, if it exists, will spend most of its time near the stationary points which are the vertices. If during time periods near a vertex it tends to approach the cycle then this is a good sign that the whole solution will approach the cycle. The behaviour of the solution near the stationary solution is determined by the linearization, at least if the stationary solution is hyperbolic.

In a previous post I described a result of Stefan Liebscher and collaborators which provides detailed information on the nature of the initial singularities of some spatially homogeneous spacetimes which are vacuum or where the matter content is described by a perfect fluid with linear equation of state . In that situation the Einstein equations can be reduced to a system of ordinary differential equations, the Wainwright-Hsu system, which treats all Bianchi class A models in a unified way. In particular it includes the type IX models. There is a heteroclinic cycle consisting of three Bianchi type I vacuum solutions. The main theorem of the paper is that there is a codimension one submanifold of initial data for which the -limit set of the corresponding solution is the heteroclinic cycle just described. The qualitative nature of this result is just as in the general discussion above except that the direction of time has been reversed. The system for vacuum solutions is four-dimensional. The vertices of the cycle are embedded in a one-dimensional manifold of stationary solutions and so the linearization must have at least one zero eigenvalue. As a consequence these vertices are not hyperbolic but the problem can be overcome. Of the remaining eigenvalues one is negative and the others are positive. The theorem makes use of the fact that for . In the presence of a fluid there is additional positive eigenvalue . The same idea of proof applies provided . This inequality is equivalent to an inequality for the parameter in the equation of state of the fluid.

An analogue of the vacuum solutions of type IX is given by solutions of type with a magnetic field. The dynamics of these solutions near the singularity was studied a long time ago by Marsha Weaver. In this situation there is a heteroclinic cycle essentially identical to that in the vacuum case. It is then natural to ask whether an analogue of the known theorem in the vacuum case in the paper by Stefan and collaborators holds. Together with Stefan and Blaise Tchapnda we have now written a paper on this subject. It turns out that there is a closely analogous result but that it is a lot harder to prove. The reason is that the eigenvalues of the linearization are in a less favourable configuration. Fortunately a weaker condition on the eigenvalues suffices. Suppose that denotes the eigenvalue at a vertex corresponding to the outgoing orbit in the cycle at that point. Then it suffices to assume that without imposing conditions of the other , provided that another condition on the existence of invariant manifolds is satisfied. The existence of these manifolds is a consequence of the geometric nature of the problem which gives rise to the dynamical system being considered. In this way we get a result on the stability of the heteroclinic cycle in the model with magnetic field. We are also able to remove the undesirable restriction on in the case with fluid. This work gives rise to a number of new questions on possible generalizations of this result. For more information on this I refer to the discussion section of our paper.

Dynamics of the MAP kinase cascade

The MAP kinase cascade is a group of enzymes which can iteratively add phosphate groups to each other. More specifically, when a suitable number of phosphate groups have been added to one enzyme in the cascade it becomes activated and can add a phosphate to the next enzyme in the row. I found this kind of idea of enzymes modifying each other with the main purpose of activating each other fascinating when I first came across it. (The first example I saw was actually the complement cascade which occurs in immunology.) This type of structure is just asking to be modelled mathematically and not surprisingly a lot of work has been done on it. Here I will survey some of what is known.

The MAP kinase cascade is a structure which occurs in many types of cells. It has three layers. The first layer consists of a protein which can be phosphorylated once. The second layer consists of a protein which can be phosphorylated twice by the same enzyme. This enzyme is the phosphorylated form of the protein in the first layer. The third layer also consists of a protein which can be phosphorylated twice by the same enzyme. This enzyme is the doubly phosphorylated form of the protein in the second layer. The protein in the third layer is the one which is called MAP kinase (mitogen activated protein kinase, MAPK). A kinase is an enzyme which phosphorylates something else and so it is not suprising that the protein in the second layer is called a MAP kinase kinase (MAPKK). The protein in the first layer is accordingly called a MAP kinase kinase kinase (MAPKKK). The roles of the players in this scheme can be taken by different enzymes. For concreteness I name those which occur in the case of human T cells. There the proteins in the first, second and third layers are called Raf, MEK and ERK, respectively. The protein which phophorylates Raf, and hence starts the whole cascade, is Ras. It, or rather the corresponding gene ras, is famous as an oncogene. This means that when the gene is not working properly cancer can result. In fact many drugs used in cancer treatment target proteins belonging to the MAP kinase cascade.

A model for the MAP kinase cascade was written down by Huang and Ferrell (PNAS, 93, 10078). They used a description of Michaelis-Menten type where for each basic substance three species are included in the network. These are the substance itself (free substrate), the enzyme and the complex of the two. Of course since in the MAP kinase cascade certain proteins act both as substrate and enzyme in different reactions there is some overlap between these. For clarity this may be called the ‘extended Michaelis-Menten’ description to contrast it with the ‘effective Michaelis-Menten’ description arising from the extended version by a quasi-steady state limiting process. Note that for a given basic reaction network with species the extended MM description has more than species but still has mass-action kinetics whereas the effective MM description has species but kinetics more complicated than mass action. Phosphatases catalysing the reverse reactions are also included in the model. The phosphatase which removes both phosphate groups of ERK is called MKP3.

In the paper the steady states of the model are studied and an input-output relation is computed numerically. The activity of the MAPK is plotted as a function of the concentration of the first enzyme (Ras in the example). A sigmoidal curve is found which corresponds to what is called ultrasensitivity. The dynamical properties of the model are not discussed. In particular it is not discussed whether there might be multistability (more than one stable stationary solution for fixed values of the parameters) or periodic solutions. The authors also did experiments whose results agreed well with the theoretical predictions. The experiments were done with extracts from the oocytes (immature egg cells) of the frog Xenopus laevis.

The possible dynamic behaviour was investigated in later papers. In some of these the effect of adding an additional feedback was considered. This kind of feedback is probably important in real biological systems. It may, for instance, explain why the results of experiments on whole oocytes are different from those done with extracts. Here, for mathematical simplicity, I will restrict to the case without additional feedback, in other words to the original Huang-Ferrell model. Multistability in this type of model was found in a paper of Markevich, Hoek and Kholodenko (J. Cell Biol. 164, 353). They investigate both extended and effective MM dynamics numerically and find bistability for both. In the extended MM model, which is the one I am most interested in here, the phosphorylation is supposed to be distributive. In other the words the kinase is released between the two phosphorylation steps. The alternative to this is called processive phosphorylation. In a paper of Conradi et. al. this result is compared with chemical reaction network theory (CRNT). It is found that while techniques from CRNT yield results agreeing with those of Markevich et. al. for the case where both the kinase and the phosphatase act in a distributive way, if one of these is replaced by a processive mechanism it can be proved using the Deficiency One Algorithm of CRNT that there is no multistationarity. The case with distributive phosphorylation is the special case of what is called a multiple futile cycle with steps. Wang and Sontag (J. Math Biol. 57, 29) proved upper and lower bounds for the number of steady states in this type of system under certain assumptions on the parameters. In particular this confirms that there can be three steady states (without determining their stability). Going beyond the single layer to the full cascade opens up more possibilities. Numerical evidence has been presented by Qiao et. al. (PLOS Comp. Biol. 9, 2007) that there are periodic solutions. To understand why these should exist it might be best to think of them as relaxation oscillations.

Entrainment by oscillations

In the book of Goldbeter which I have mentioned in several recent posts a concept which occurs repeatedly is that of entrainment. While looking for some more information about this topic I found a paper of Russo, di Bernardo and Sontag (PloS Computational Biology 4, e1000739) which gives an insightful treatment of the subject. The basic idea is to consider two systems which are coupled in some way and to consider the influence of oscillations in one system on the behaviour of the other. It is easy to see how this might be translated into a problem expressed in terms of dynamical systems. A classical example related to this is contained in a story about Christiaan Huygens who was, among other things, the inventor of the pendulum clock in the mid 17th century. Apparently he did not construct clocks himself but had them made by others according to his plans. The well-known story is that he noticed that when two pendulum clocks were placed next to each other the phase of their oscillations became synchronized with, say, one always at the leftmost point of its swing when the other was at the rightmost. Another example is that of the circadian rhythm. There is a 24 hour rhythm in our body and it is interesting to know whether it comes from an intrinsic oscillator or not. Experiments with subjects isolated from the usual rhythm of day and night show that there is an intrinsic oscillator but that its period is closer to 25 hours. Under normal circumstances its period is brought to 24 hours due to the cycle of day and night by entrainment.

The particular mathematical set-up considered in the paper of Russo et. al. is the following. Consider an autonomous dynamical system containing some parameters. Now replace one or more of those parameters by functions of time with period . If solutions of the original system have a suitable tendency to converge to a stationary solution for a given choice of the parameters then solutions of the resulting non-autonomous system converge to periodic solutions of period . In the papers there are nice plots of numerical simulations which give a striking picture of this behaviour. The central result of the paper is a theorem which guarantees this type of behaviour under certain hypotheses. As pointed out in the paper verifying these hypotheses has some similarity to finding a Lyapunov function for an autonomous system. The positive side is that if it can be done it is possible to get strong conclusions. The negative side is that verifying the hypotheses is generally a matter of trial and error. There is no algorithm available for doing that.

The criterion is dependent on the choice of a matrix norm. This is used to define a quantity called the matrix measure of a matrix . The criterion is that the Jacobian of the function defining the dynamical system should have a matrix measure which is bounded above by a negative constant. In that case the system is said to be infinitesimally contracting. The matrix measure is defined by a limiting procedure, , but for particular choices of the matrix norm it is possible to calculate in a purely algebraic way. I have no intuitive feeling for what this definition means.

La vie oscillatoire

I have continued reading the book ‘La Vie Oscillatoire’ and I have learned many interesting things. Some of them were things I was aware of on some level already but they have now become clearer. Others were quite new to me. Some of them have to do with biology, some with mathematics. Chapter 5 is concerned with the secretion of hormones. I had the naive view that the effect of a hormone was due to its overall level. In reality frequency encoding is important for many hormone signals. For instance the triggering of ovulation is dependent on having the right kind of oscillatory signals and there is a therapy for infertility based on delivering a hormone in an appropriate oscillatory way. The menstrual cycle makes it natural to think about oscillations in that context but the oscillations just mentioned are on timescales of one hour rather than one month. Similar control mechanisms seem to apply to many other hormones. In the book biological systems are often compared across very different animals or other living organisms. One interesting conclusion is that the fact that the human menstrual cyle takes about one month is an accident. This deals a blow to more or less mystical ideas relating this cycle with the moon. A recurring theme in the book is the way in which frequency modulated signals are decoded in biological systems. Here I feel that I am at the edge of a part of the theory of dynamical systems which I should learn a lot more about.

The theme of Chapter 6 is rhythms in the brain. Elsewhere in this blog I have written about the Hodgkin-Huxley model several times. This can be used to describe the propagation of an action potential along an axon. However this is not its only application. It can be also be used to describe the oscillatory behaviour of individual neurons. The basic phenomena involved are the flow of sodium and potassium ions across the cell membrane. Calcium ions also play a role in some cases. A mathematical phenomenon which comes up in this discussion is that of bursting oscillations. Since I was not previously familiar with that I now read up on it a bit. My main source was the book ‘Mathematical Physiology’ by Keener and Sneyd. A variable in a dynamical system is said to display bursting oscillations when it has the following type of behaviour. There is a period where it changes very little followed by a period where it oscillates with high frequency and fairly constant amplitude. Then it returns to the quiescent phase it started in. It goes through this cycle repeatedly. The minimum ingredients required for a dynamical system to show this type of behaviour are dimension at least three and equations with different timescales. One example occurs in a model for the production of insulin by the -cells of the pancreas. In this case there are one slow and two fast variables. Heuristically the slow variable is first thought of as a parameter. As it is varied the dynamics of the fast system changes. For certain parameter values the fast system has a stable steady state. Starting from this point the slow variable changes in such a way that this steady vanishes in a fold bifurcation. The solution then moves quickly to being close to a limit cycle of the fast system and stays there for a while. The slow variable then changes in such a way that the limit cycle is destroyed in a homoclinic bifurcation. By that time the stable steady state has reappeared and the solution can jump back to it. Burst oscillations are classified into three types and the scenario just sketched is Type I. Coming back to the brain (or at least to neurons) a popular experimental system is the sea slug Aplysia. An isolated neuron of this organism can exhibit bursting oscillations. In this case (Type II) there are two slow variables. The oscillatory phase starts and ends with a homoclinic bifurcation. In Type III a period of oscillations is bounded by two subcritical Hopf bifurcations.

Another class of oscillations arises from the collective behaviour of a small number (up to about thirty) neurons. This has been studied for instance in the context of the motion of the snail Lymnaea. In the human brain there is a great variety of oscillations, producing characteristic traces in the EEG. They include a number of type of waves named after letters of the Greek alphabet, starting with , some of which have entered popular culture. I would like to learn more about them sometime, but the day for that has not yet arrived.