Archive for the ‘dynamical systems’ Category

The method of averaging

June 24, 2018

Techniques of averaging in the theory of differential equations have interested me for a long time. It happens that when determining the asymptotics of certain solutions it is important to show that certain integrals are finite although these integrals are not absolutely convergent. For this there must be a suitable cancellation of positive and negative contributions. Back in 2007 I published a paper (Class. Quantum Grav. 24, 667) where I studied this kind of phenomenon in certain inflationary cosmological models. There I did not use any general techniques but instead I just derived estimates by hand. More recently I have spent some time learning what general techniques there are. One famous method is the Krylov-Bogoliubov averaging method. Next semester I will organize a seminar on the subject of the method of averaging.

An iconic example which is a good starting point for discussing this subject is the Kapitza pendulum. Suppose that a rigid rod is attached to a support about which it can rotate freely in a vertical plane. If the support is fixed we get an ordinary pendulum. The steady state where the rod is vertically above the support is obviously unstable. The Kapitza pendulum is obtained by supposing that instead the support undergoes oscillations in the vertical direction with small amplitude and large frequency. For suitable choices of the parameters this stabilizes the unstable steady state of the ordinary pendulum.

How can the situation just described be understood mathematically? I follow here the discussion in Hale’s book on ordinary differential equations. The basic equation of motion is second order. If friction is ignored the equations can be put in Hamiltonian form which means reducing them to a system of two first order equations in a certain way. Introducing a rescaled time coordinate leads to a system of the form x'=\epsilon f(t,x)+\epsilon h(\epsilon t,x) which is a standard form for the method of averaging. The system contains two time scales t and \epsilon t. We now average over the fast time t, defining f_0(x)=\frac{1}{T}\int_0^T f(t,x) dt. Then the original equation is replaced by the averaged equation x'=\epsilon f_0(x)+\epsilon h(\epsilon t,x). Here T is the period of the oscillation. The definition of f_0 given by Hale is more complicated since he wants to allow more general motions of the support which might be only almost periodic. The question is now to what extent solutions of the original equation can be approximated by solutions of the averaged equation for \epsilon small. Write the right hand side of the averaged equation in the form \epsilon G(\epsilon t,x). Suppose that the averaged equation has a periodic solution. We linearize the averaged equation about that solution. This gives rise to characteristic exponents, which are the exponential growth rates of linearized perturbations. If none of these characteristic exponents is purely imaginary then we obtain an existence theorem for solutions of the original equation which are small perturbations of the periodic solutions of the averaged equation. If in addition all characteristic exponents have negative real parts then the perturbed solution is asymptotically stable and if there is a characteristic exponent with positive real part it is unstable. It is a separate problem to prove the existence of a periodic solution of the averaged equation.


The Higgins-Selkov oscillator, part 2

March 29, 2018

Some time ago I wrote about a mathematical model for glycolysis, the Higgins-Selkov oscillator. In fact I now prefer to call it the Selkov oscillator since it is a system of equations written down by Selkov, although it was obtained by modifying a previous model of Higgins. At that time I mentioned some of the difficulties involved in understanding the global properties of solutions of this system. I mentioned that I had checked that this system exhibits a supercritical Hopf bifurcation, thus proving that it has stable periodic solutions for certain values of the parameters. In the hope of getting more insights into this system I decided to make it the theme of a master’s thesis. Pia Brechmann, the student who got this subject, obtained a number of interesting results. After she submitted her thesis she and I decided to carry this a bit further so as to get a picture of the subject which was as comprehensive as possible. We have now written a paper on this subject. I cannot say that we were able to answer all the questions we would have liked to but at least we answered some of them. Here I will mention some of the things which are known and some which are not.

When written in dimensionless form the system contains two parameters \alpha (a positive real number) and \gamma (an integer which is at least two). I already mentioned a result about a Hopf bifurcation and this was originally only proved for \gamma=2. In the paper we showed that there is a supercritical Hopf bifurcation for any value of \gamma. I also previously mentioned doing a Poincaré compactification of the system and blowing up new stationary solutions which appear in this process and are too degenerate to analyse directly. The discussion of using blow-ups in polar coordinates previously given actually concerned the case \gamma=2 and does not seem practical for higher values of \gamma. It turns out that the technique of quasihomogeneous directional blow-ups, explained in the book ‘Qualitative theory of planar differential systems’ by Dumortier et al. can be used to treat the general case. This type of blow-up has the advantage that the transformations are given by monomials rather than trigonometric functions and that there is a systematic method for choosing good values of the exponents in the monomials.

We discovered a paper by d’Onofrio (J. Math. Chem. 48, 339) on the Selkov oscillator where he obtains interesting results on the uniqueness and stability of periodic solutions. It was suggested by Selkov on the basis of numerical calculations that for large \alpha there exist solutions with oscillations which grow arbitrarily large at late times (let us call these unbounded oscillations). We were not able to decide on a rigorous level whether such solutions exist or not. They cannot exist when a periodic solution does. When they exist they have to do with a heteroclinic cycle in the compactification. We showed that when a cycle of this kind exists it is asymptotically stable and in that case solutions with unbounded oscillations exist. However we were not able to decide whether a heteroclinic cycle at infinity ever exists for this system. What we did prove is that for all values of the parameters there exist unbounded solutions which are eventually monotone.

We also proved that when the steady state is stable each bounded solution converges to it and that when there exists a periodic solution it is unique and each bounded solution except the steady state converges to that. I find it remarkable that such an apparently harmless two-dimensional dynamical system is so resistent to a complete rigorous analysis.

Lotka’s system

March 11, 2018

The system of ODE \dot x=a-bxy, \dot y=bxy-cy was considered in 1910 by Lotka as a model for oscillatory chemical reactions (J. Phys. Chem. 14, 271). It exhibits damped oscillations but no sustained oscillations, i.e. no periodic solutions. It should not be confused with the famous Lotka-Volterra system for predator-prey interactions which was first written down by Lotka in 1920 (PNAS 6, 410) and which does have periodic solutions for all positive initial data. That the Lotka system has no periodic solutions follows from the fact that x^{-1}y^{-1} is a Dulac function. In other words, if we multiply the vector field defined by the right hand sides of the equations by the positive function x^{-1}y^{-1} the result is a vector field with negative divergence. This change of vector field preserves periodic orbits and it follows from the divergence theorem that the rescaled vector field has no periodic orbits. My attention was drawn to this system by the paper of Selkov on his model for glycolysis (Eur. J. Biochem. 4, 79). In his model there is a parameter \gamma which is assumed greater than one. He remarks that if this parameter is set equal to one the system of Lotka is obtained. Selkov obtains his system as a limit of a two-dimensional system with more complicated non-linearities. If the parameter \gamma is set to one in that system equations are obtained which are related to Higgins’ model of glycolysis. Selkov remarks that this last system admits a Dulac function and hence no sustained oscillations and this is his argument for discarding Higgins’ model and replacing it by his own. The equations are \dot x=a-b\frac{xy}{1+y+xy} and \dot y=b\frac{xy}{1+y+xy}-cy. (To obtain the limit already mentioned it is first necessary to do a suitable rescaling of the variables.) In this case the Dulac function is \frac{1+y+xy}{xy}.

The fact that the Lotka-Volterra system admits periodic solutions can be proved by exhibiting a conserved quantity. At this point I recall the well-known fact that while conserved quantities and their generalizations, the Lyapunov functions, are very useful when you have them there is no general procedure for finding them. This naturally brings up the question: if I did not know the conserved quantity for the Lotka-Volterra system how could I find it? One method is as follows. First divide the equation for \dot y by that for \dot x to get a non-autonomous equation for dy/dx, cheerfully ignoring points where \dot x=0. It then turns out that the resulting equation can be solved by the method of separation of variables and that this leads to the desired conserved quantity.

One undesirable feature of the Lotka-Volterra system is that it has a one-parameter family of periodic solutions and must therefore be suspected to be structurally unstable. In addition, if we consider a solution where predators are initially absent the prey population grows exponentially. The latter feature can be eliminated by replacing the linear growth term in the equation for the prey by a logistic one. A similar term corresponding to higher death rates at high population densities can be added in the equation for the predators but the latter modification has no essential effect. This is a Lotka-Volterra model with intraspecific competition. As discussed in the book ‘Evolutionary Games and Population Dynamics’ by Josef Hofbauer and Karl Sigmund, when this model has a positive steady state that state is globally asymptotically stable. The proof uses the fact that the expression which defines the conserved quantity in the usual Lotka-Volterra model defines a Lyapunov function in the case with intraspecific competition. This is an example of the method of obtaining conserved quantities or Lyapunov functions by perturbing those which are already known in special cases.

It follows from Poincaré-Bendixson theory that the steady states in the Lotka model and the Higgins model are globally asymptotically stable. This raises the question whether we could not find Lyapunov functions for those systems. I do not know how. The method used for Lotka-Volterra fails here because the equation for dy/dx is not separable.

The matrix-tree theorem

October 29, 2017

When I was an undergraduate graph theory was a subject which seemed to me completely uninteresting and I very consciously chose not to attend a lecture course on the subject. In the meantime I have realized that having become interested in reaction networks I could profit from knowing more about graph theory. One result which plays a significant role is the matrix-tree theorem. It asserts the equality of the number of subgraphs of a certain type of a given graph and a certain determinant. This could be used to calculate numbers of subgraphs. On the other hand, it could be used in the other direction calculate determinants and it is the second point of view which is relevant for reaction networks.

For the first part of the discussion here I follow these notes. If G is an (unoriented) graph then a spanning tree is a connected subgraph which includes all vertices of G and contains no cycles. The graph Laplacian L of G is the matrix for which L_{ii} is equal to the number of edges containing the vertex v_i, L_{ij}=-1 if i\ne j and there is an edge joining vertex v_i to vertex v_j and L_{ij}=0 otherwise. The first version of the matrix tree theorem, due to Kirchhoff, says that the number of spanning trees of G can be calculated in the following way. Choose any vertex v_j of G and remove the jth row and jth column from L to get a matrix L_j. Then compute the determinant of L_j. Surprisingly the value of the determinant is independent of j. The first version of the theorem can be obtained as a consequence of a second more sophisticated version proved a hundred years later by Tutte. This concerns directed graphs. A vertex b is said to be accessible from a if there is a directed path from a to b. A vertex in a directed graph is called a root if every other vertex is accessible from it. A directed tree is a directed subgraph which contains a root and which, when the orientation is forgotten, is a tree (i.e. it is connected and contains no unoriented cycles). There is then an obvious way to define a spanning rooted tree. To formulate the second version of the theorem we introduce the matrix Laplacian L of the directed graph. If i=j the entry L_{ij} is the number of edges with final vertex i. L_{ij}=-1 if i\ne j and there is a directed edge from v_i to v_j. L_{ij}=0 if i\ne j and there is no edge connecting v_i and v_j. The statement of the theorem is that the number of rooted trees with root v_j is equal to the determinant of L_j, a matrix derived from L as before. To see that the first version of the theorem follows the second first associate an oriented graph to an unoriented one by putting in oriented edges joining a pair of vertices in both directions whenever they are joined by an edge in the original graph. Then there is a bijection between trees rooted at v in the unoriented graph and oriented trees rootes at v in the oriented graph. On the other hand the two graphs have the same Laplacian since the number of edges ending at a vertex in the oriented graph is the same as the number of edges having that endpoint in the unoriented graph.

What is the connection of all this with reaction networks? Consider a chemical reaction network with only monomolecular reactions and reaction coefficients all equal to one. Then under mass action kinetics the evolution equations for the concentrations are \dot x=Lx, where L is the Laplacian matrix of the network. There is a conserved quantity, which is the sum of the concentrations of all species. A steady state is an element of the kernel of L. The next part of the discussion follows a paper of Gunawardena (PLoS ONE, 7(5), e36321). He allows general reaction constants. The notion of the Laplacian is extended correspondingly. If the network is strongly connected (in the terminology of graph theory) or weakly reversible (in the terminology of chemical reaction network theory) the kernel of the Laplacian matrix is one-dimensional. Thus there is precisely one steady state in each stoichiometric compatibility class. It is moreover possible to compute a vector spanning the kernel of N by graph theoretical methods. This also goes back to Tutte. To get the ith component first take the product of the reaction constants over a rooted tree and then sum these quantities over the rooted trees with root v_i. More generally the dimension of the kernel of N is equal to the number of terminal strong linkage classes. It is also revealed there that the Laplacian corresponds to the matrix A considered by Horn and Jackson. These ideas are also related to the King-Altman theorem of enzyme kinetics. I have the impression that I have as yet only scratched the surface of this subject. I hope I will be able to come back and understand more about it at a later date.

Sternberg’s theorem and smooth conjugacy

August 23, 2017

Suppose we have a dynamical system \dot x=f(x) and a steady state x_0. Let \phi be a homeomorphism defined on a neighbourhood of x_0 with \phi (x_0)=0. If we transform the system to coordinates y=\phi(x) under what circumstances can \phi be chosen so that the transformed system is equal to the linearization of the original system about x_0? The answer depends on the regularity required for f and \phi. If f is C^1 and x_0 is hyperbolic (i.e. all eigenvalues of the linearization at x_0 have non-vanishing real part) then there always exists a continuous mapping \phi with this property. (In that case we cannot necessarily transform the vector field but we can transform its integral curves.) This is the Grobman-Hartman theorem. Even if f is C^\infty it is not in general possible to choose \phi to be C^2.

Sternberg’s theorem shows that if f is C^\infty and there are no resonances then \phi can be chosen to be C^\infty. The definition of a resonance in this context is as follows. Let \lambda_i be the eigenvalues of the linearization. A resonance is a relation of the form \lambda_j=\sum_i m_i\lambda_i with positive integers m_i satisfying \sum_i m_i\ge 2. For example, if the eigenvalues are \lambda_1=1 and \lambda_2=2 we have the resonance \lambda_2=2\lambda_1. Consider the two-dimensional case. If there the eigenvalues are real with opposite signs then there can be no resonance and Sternberg’s theorem applies to all saddles in two dimensions. If both eigenvalues are positive then it is known that \phi can be chosen to be C^1 but resonances can occur and there are cases where although f is analytic \phi cannot be chosen to be C^2. Even when there are no resonances it is in general not the case that if f is analytic the mapping \phi can also be chosen analytic. It is necessary to make a further assumption (no small divisors). This says roughly speaking that the expressions \lambda_j-\sum_i m_i\lambda_i should not only be non-zero but that in addition they should not be able to approach zero too fast as the number of summands increases. This is the content of a theorem of Siegel.

All the results discussed up to now concern the hyperbolic case. What happens in the presence of purely imaginary eigenvalues? The Grobman-Hartman theorem has a generalization, the theorem of Shoshitaishvili, also known as the reduction theorem. In words it says that there is a continuous mapping \phi which reduces the system to the product of the flow on any centre manifold and a standard saddle. A standard saddle is linear and so if the centre manifold is trivial this reduces to Grobman-Hartman. The next question is what can be achieved with mappings \phi of higher regularity. A fact which is implicit in the reduction theorem is that linear systems where the origin is a hyperbolic fixed point whose stable and unstable manifolds have the same dimension are topologically equivalent. Evidently they are in general not smoothly equivalent since any smooth mapping preserving the origin will not change the eigenvalues of the linearization. However once this has been taken into account the theorem does generalise under the assumptions that there are no resonances among the eigenvalues with non-zero real parts.More precisly the eigenvalues \lambda_i in the definition of a resonance must have non-zero real parts. By contrast quantity \lambda_j should either be an eigenvalue with non-zero real part or zero. The result being referred to here is Takens’ theorem.

Suppose that for a steady state we introduce coordinates x, y and z corresponding to the stable, unstable and centre manifolds. Then the form achieved by the reduction in Takens theorem is that we get the equations \dot x=A(z)x, \dot y=A(z)y and \dot z=g(z) for matrices A and B and a function g depending on the central variables z. One technical point is that while for a smooth system the reduction can be done by a mapping which is C^k for any finite k it cannot in general be done by a mapping which is itself infinitely differentiable. As an example, consider a two-dimensional system and a fixed point where the linearization has one zero and one non-zero eigenvalue. In this case there can be no resonances and the theorem says that there is a C^k coordinate transformation which reduces the system to the form \dot x=a(y)x, \dot y=b(y) with a(0)\ne 0, b(0)=0 and b'(0)=0.

Conference on reaction networks and population dynamics in Oberwolfach

July 9, 2017

This is a belated report on a conference in Oberwolfach I attended a couple of weeks ago. The title includes two elements. The first held no suprises for me but the second was rather different from what I had expected. My expectation was that it would be about the evolution of populations of organisms. In fact it was rather focussed on models related to genetics, in other words with the question of how certain genetic traits spread through a population.

A talk I found very interesting was by Sebastian Walcher. I already wrote briefly about a talk of his in Copenhagen in a previous post but this time I understood a lot more. The question he was concerned with is how to find interesting small parameters in dynamical systems which allow the application of geometric singular perturbation theory. In GSPT the system written in the slow time (with the smallness parameter included as a variable) contains a whole manifold of steady states, the critical manifold. The most straightforward theory is obtained when the eigenvalues of the linearization of the system transverse to the critical manifold lie away from the imaginary axis. This corresponds to the situation of a transversely hyperbolic manifold of steady states. The first idea of Walcher’s talk is that whenever we have a transversely hyperbolic manifold of steady states in a dynamical system this is an opportunity for identifying a small parameter. This may not sound very useful at first sight because it would seem reasonable that generic dynamical systems would never contain manifolds of steady states of dimension greater than zero. There is a reason why this observation is misleading for systems arising from reaction networks. In these systems the state space is defined by positivity conditions on the concentrations and there are also certain parameters (such as reaction constants and total amounts) which are required to be positive. To have a name let us call the region defined by these positivity conditions the conventional region of the spaces of states and parameters. In the conventional region manifolds of steady states are not to be expected. On the other hand it frequently happens that they arise when we go to the boundary of that region. A familiar example is the passage to the Michaelis-Menten limit in the system describing a reaction catalysed by an enzyme. This takes us from the extended mass action kinetics for substrate, free enzyme and substrate-enzyme complex to Michaelis-Menten kinetics for the substrate alone. Roughly speaking it is the limit where the amount of the enzyme is very small compared to the amount of the substrate. I often wondered whether there could not be a kind of ‘anti-Michaelis-Menten’ limit where the amount of enzyme is very large compared to the amount of the substrate. I asked Walcher whether he knew how to do this and how it fitted into his general scheme. He gave me a positive answer to this question and some references and I must look into this in detail when I get time. The reason for being interested in this is that if we can obtain suitable information about a limiting case on the boundary it may be possible obtain information on the part of the conventional region where a certain parameter is small but non-zero.

There was one talk which did have a connection to population biology in way closer to what I had expected. It happens all the time that ecosystems are damaged by exotic species imported, deliberately or by accident, from other parts of the world. There are also well-known stories of the type that to try to control exotic species number one exotic species number two is introduced and is itself very harmful. It is nice to hear an example where this kind of introduction of an exotic species was very successful. It is the case of the cassava plant which was introduced from South America to Africa and became a staple food there. Then an insect from South America (species number one) called the mealy bug was introduced accidentally and caused enormous damage. Finally an ecologist called Hans Herren introduced a parasitic wasp (species number two) from South America, restoring the food supply and saving numerous lives (often the number 20 million is quoted). More details of this story can be found here.

I want to mention one statement made in the talk of Gheorghe Craciun in Oberwolfach which I found intriguing. I might have heard this before but it did not stick in my mind properly. The statement is that the set of dynamical systems which possess a complex balanced steady state is a variety of codimension \delta, where \delta is the deficiency. There seemed to be some belief in the audience that this variety is actually a smooth manifold. On one afternoon we had something similar to the breakout sessions in Banff. I suggested the topic for one of these, which was Lyapunov functions. The idea was to compare classes of Lyapunov functions which people working on different classes of dynamical systems knew. This certainly did not lead to any breakthrough but I think it did lead to a useful exchange of information. I documented the discussion for my own use and I think I could profit by following some of the leads there.

To finish I want to mention a claim made by Ankit Gupta in his talk. It did not sound very plausible to me but I expect that it at least contains a grain of truth. He said that these days more papers are published on NF\kappa B than on all of mathematics.

Conference on mathematical analysis of biological interaction networks at BIRS

June 9, 2017

I have previously written a post concerning a meeting at the Banff International Research Station (BIRS). This week I am at BIRS again. Among the topics of the talks were stochastic chemical reaction networks, using reaction networks in cells as computers and the area of most direct relevance to me, multiple steady states and their stability in deterministic CRN. Among the most popular examples occurring in the talks in the latter area were the multiple futile cycle, the MAPK cascade and the EnvZ/OmpR system. In addition to the talks there was a type of event which I had never experienced before called breakout sessions. There the participants split into groups to discuss different topics. The group I joined was concerned with oscillations in phosphorylation cycles.

In the standard dual futile cycle we have a substrate which can be phosphorylated up to two times by a kinase and dephosphorylated again by a phosphatase. It is assumed that the (de-)phosphorylation is distributive (the number of phosphate groups changes by one each time a substrate binds to an enzyme) and sequential (the phosphate groups are added in one order and removed in the reverse order). A well-known alternative to this is processive (de-)phosphorylation where the number of phosphate groups changes by two in one encounter between a substrate and an enzyme. It is known that the double phosphorylation system with distributive and sequential phosphorylation admits reaction constants for which there are three steady states, two of which are stable. (From now on I only consider sequential phosphorylation here.) By contrast the corresponding system with processive phophorylation always has a unique steady state. Through the talk of Anne Shiu here I became aware of the following facts. In a paper by Suwanmajo and Krishnan (J. R. Soc. Interface 12:20141405) it is stated that in a mixed model with distributive phosphorylation and processive dephosphorylation periodic solutions occur as a result of a Hopf bifurcation. The paper does not present an analytical proof of this assertion.

It is a well-known open question, whether there are periodic solutions in the case that the modificiations are all distributive. It has been claimed in a paper of Errami et. al. (J. Comp. Phys. 291, 279) that a Hopf bifurcation had been discovered in this system but the claim seems to be unjustified. In our breakout sessions we looked at whether oscillations might be exported from the mixed model to the purely distributive model. We did not get any definitive results yet. There were also discussions on effective ways of detecting Hopf bifurcations, for instance by using Hurwitz determinants. It is well-known that oscillations in the purely distributive model, if they exist, do not persist in the Michaelis-Menten limit. I learned from Anne Shiu that it is similarly the case that the oscillations in the mixed model are absent from the Michaelis-Menten system. This result came out of some undergraduate research she supervised. Apart from these specific things I learned a lot just from being in the environment of these CRN people.

Yesterday was a free afternoon and I went out to look for some birds. I saw a few things which were of interest to me, one of which was a singing Tennessee warbler. This species has a special significance for me for the following reason. Many years ago when I still lived in Orkney I got an early-morning phone call from Eric Meek, the RSPB representative. He regularly checked a walled garden at Graemeshall for migrants. On that day he believed he had found a rarity and wanted my help in identifying it and, if possible, catching it. We did catch it and it turned out to be a Tennessee warbler, the third ever recorded in Britain. That was big excitement for us. I had not seen Eric for many years and I was sad to learn now that he died a few months ago at a relatively young age. The name of this bird misled me into thinking that it was at home in the southern US. In fact the name just came from the fact that the first one to be described was found in Tennessee, a migrant. The breeding range is much further north, especially in Canada. Thus it is quite appropriate that I should meet it here.

The Routh-Hurwitz criterion

May 7, 2017

I have been aware of the Routh-Hurwitz criterion for stability for a long time and I have applied it in three dimensions in my research and tried to apply it in four. Unfortunately I never felt that I really understood it completely. Here I want to finally clear this up. A source which I found more helpful than other things I have seen is One problem I have had is that the Hurwitz matrices, which play a central role in this business, are often written in a form with lots of … and I was never sure that I completely understood the definition. I prefer to have a definite algorithm for constructing these matrices. The background is that we would like to understand the stability of steady states of a system of ODE. Suppose we have a system \dot x=f(x) and a steady state x_0, i.e. a solution of f(x_0)=0. It is well-known that this steady state is asymptotically stable if all eigenvalues \lambda of the linearization A=Df(x_0) have negative real parts. This property of the eigenvalues is of course a property of the roots of the characteristic equation \det(A-\lambda I)=a_0\lambda^n+\ldots+a_{n-1}\lambda+a_n=0. It is always the case here that a_0=1 but I prefer to deal with a general polynomial with real coefficients a_i, 0\le i\le n and a criterion for the situation where all its roots have negative real parts. It is tempting to number the coefficients in the opposite direction, so that, for instance, a_n becomes a_0 but I will stick to this convention. Note that it is permissible to replace a_k by a_{n-k} in any criterion of this type since if we multiply the polynomial by \lambda^{-n} we get a polynomial in \lambda^{-1} where the order of the coefficients has been reversed. Moreover, if the real part of \lambda is non-zero then it has the same sign as the real part of \lambda^{-1}. I find it important to point this out since different authors use different conventions for this. It is convenient to formally extend the definition of the a_i to the integers so that these coefficients are zero for i<0 and i>n.

For a fixed value of n the Hurwitz matrix is an n by n matrix defined as follows. The jth diagonal element is a_j, with 1\le j\le n. Starting from a diagonal element and proceeding to the left along a row the index increases by one in each step. Similarly, proceeding to the right along a row the index decreases by one. In the ranges where the index is negative or greater than n the element a_n can be replaced by zero. The leading principal minors of the Hurwitz matrix, in other words the determinants of the submatrices which are the upper left hand corner of the original matrix, are the Hurwitz determinants \Delta_k. The Hurwitz criterion says that the real parts of all roots of the polynomial are negative if and only if a_0>0 and \Delta_k>0 for all 1\le k\le n. Note that a necessary condition for all roots to have negative real parts is that all a_i are positive. Now \Delta_n=a_n\Delta_{n-1} and so the last condition can be replaced by a_n>0. Note that the form of the \Delta_k does not depend on n. For n=2 we get the conditions a_0>0, a_1>0 and a_2>0. For n=3 we get the conditions a_0>0, a_1>0, a_1a_2-a_0a_3>0 and a_3>0. Note that the third condition is invariant under the replacement of a_j by a_{n-j}. When a_0a_3-a_1a_2>0, a_0>0 and a_3>0 then the conditions a_1>0 and a_2>0 are equivalent to each other. In this way the invariance under reversal of the order of the coefficients becomes manifest. For n=4 we get the conditions a_0>0, a_1>0, a_1a_2-a_0a_3>0, a_1a_2a_3-a_1^2a_4-a_0a_3^2>0 and a_4>0.

Next we look at the issue of loss of stability. If H is the region in matrix space where the Routh-Hurwitz criteria are satisfied, what happens on the boundary of H? One possibility is that at least one eigenvalue becomes zero. This is equivalent to the condition a_n=0. Let us look at the situation where the boundary is approached while a_n remains positive, in other words the determinant of the matrix remains non-zero. Now a_0=1 and so one of the quantities \Delta_k with 1\le k\le n-1 must become zero. In terms of eigenvalues what happens is that a number of complex conjugate pairs reach the imaginary axis away from zero. The generic case is where it is just one pair. An interesting question is whether and how this kind of event can be detected using the \Delta_k alone. The condition for exactly one pair of roots to reach the imaginary axis is that \Delta_{n-1}=0 while the \Delta_k remain positive for k<n-1. In a paper of Liu (J. Math. Anal. Appl. 182, 250) it is shown that the condition for a Hopf bifurcation that the derivative of the real part of the eigenvalues with respect to a parameter is non-zero is equivalent to the condition that the derivative of \Delta_{n-1} with respect to the parameter is non-zero. In a paper with Juliette Hell (Math. Biosci. 282, 162), not knowing the paper of Liu, we proved a result of this kind in the case n=3.

Poincaré, chaos and the limits of predictability

March 5, 2017

In the past I was surprised that there seemed to be no biography of Henri Poincaré. I recently noticed that a biography of him had appeared in 2013. The title is ‘Henri Poincaré. A scientific biography’ and the author is Jeremy Gray. At the moment I have read 390 of the 590 pages. I have learned interesting things from the book but in general I found it rather disappointing. One of the reasons is hinted at by the subtitle ‘A scientific biography’. Compared to what I might have hoped for the book concentrates too much on the science and too little on the man. Perhaps Poincaré kept his private life very much to himself and thus it was not possible to discuss these aspects more but if this is so then I would have found it natural that the book should emphasize this point. I have not noticed anything like that. I also found the discussion of the scientific topics of Poincaré’s work too technical in many places. I would have preferred a presentation of the essential ideas and their significance on a higher level. There are other biographies of great mathematicians which made a better impression on me. I am thinking of the biography of Hilbert by Constance Reid and even of the slim volumes (100 pages each) on Gauss and Klein written in East Germany.

On important discovery of Poincaré was chaos. He discovered it in the context of his work on celestial mechanics and indeed that work was closely connected to his founding the subject of dynamical systems as a new way of approaching ordinary differential equations, emphasizing qualitative and geometric properties in contrast to the combination of complex analysis and algebra which had dominated the subject up to that point. The existence of chaos places limits on predictability and it is remarkable that these do not affect our ability to do science more than they do. For instance it is known that there are chaotic phenomena in the motion of objects belonging to the solar system. This nevertheless does not prevent us from computing the trajectories of the planets and those of space probes sent to the other end of the solar system with high accuracy. These space probes do have control systems which can make small corrections but I nevertheless find it remarkable how much can be computed a priori, although the system as a whole includes chaos.

This issue is part of a bigger question. When we try to obtain a scientific understanding of certain phenomena we are forced to neglect many effects. This is in particular true when setting up a mathematical model. If I model something using ODE then I am, in particular, neglecting spatial effects (which would require partial differential equations) and the fact that often the aim is not to model one particular object but a population of similar objects and I neglect the variation between these objects which I do not have under control and for whose description a stochastic model would be necessary. And of course quantum phenomena are very often neglected. Here I will not try to address these wider issues but I will concentrate on the following more specific question. Suppose I have a system of ODE which is a good description of the real-world situation I want to describe. The evolution of solutions of this system is uniquely determined by initial data. There remains the problem of sensitive dependence on initial data. To be able to make a prediction I would like to know that if I make a small change in the initial data the change in some predicted quantity should be small. What ‘small’ means in practice is fixed by the application. A concrete example is the weather forecast whose essential limits are illustrated mathematically by the Lorenz system, which is one of the icons of chaos. Here the effective limit is a quantitative one: we can get a reasonable weather forecast for a couple of days but not more. More importantly, this time limit is not set by our technology (amount of observational data collected, size of the computer used, sophistication of the numerical programs used) but by the system itself. This time limit will not be relaxed at any time in the future. Thus one way of getting around the effects of chaos is just to restrict the questions we ask by limits on the time scales involved.

Another aspect of this question is that even when we are in a regime where a system of ODE is fully chaotic there will be some aspects of its behaviour which will be predictable. This is why is is possible to talk of ‘chaos theory’- I know too little about this subject to say more about it here. One thing I find intriguing is the question of model reduction. Often it is the case that starting from a system of ODE describing something we can reduce it to an effective model with less variables which still includes essential aspects of the behaviour. If the dimension of the reduced model is one or two then chaos is lost. If there was chaos in the original model how can this be? Has there been some kind of effective averaging? Or have we restricted to a regime (subset of phase space) where chaos is absent? Are the questions we tend to study somehow restricted to chaos-free regions? If the systems being modelled are biological is the prevalence of chaos influenced by the fact that biological systems have evolved? I have seen statements to the effect that biological systems are often ‘on the edge of chaos’, whatever that means.

This post contains many questions and few answers. I just felt the need to bring them up.

Models for photosynthesis, part 4

September 19, 2016

In previous posts in this series I introduced some models for the Calvin cycle of photosynthesis and mentioned some open questions concerning them. I have now written a paper where on the one hand I survey a number of mathematical models for the Calvin cycle and the relations between them and on the other hand I am able to provide answers to some of the open questions. One question was that of the definition of the Pettersson model. As I indicated previously this was not clear from the literature. My answer to the question is that this system should be treated as a system of DAE (differential-algebraic equations). In other words it can be written as a set of ODE \dot x=f(x,y) coupled to a set of algebraic equations g(x,y)=0. In general it is not clear that this type of system is locally well-posed. In other words, given a pair (x_0,y_0) with g(x_0,y_0)=0 it is not clear whether there is a solution (x(t),y(t)) of the system, local in time, with x(0)=x_0 and y(0)=y_0. Of course if the partial derivative of g with repect to y is invertible it follows by the implicit function theorem that g(x,y)=0 is locally equivalent to a relation y=h(x) and the original system is equivalent to \dot x=f(x,h(x)). Then local well-posedness is clear. The calculations in the 1988 paper of Pettersson and Ryde-Pettersson indicate that this should be true for the Pettersson model but there are details missing in the paper and I have not (yet) been able to supply these. The conservative strategy is then to stick to the DAE picture. Then we do not have a basis for studying the dynamics but at least we have a well-defined system of equations and it is meaningful to discuss its steady states.

I was able to prove that there are parameter values for which the Pettersson model has at least two distinct positive steady states. In doing this I was helped by an earlier (1987) paper of Pettersson and Ryde-Pettersson. The idea is to shut off the process of storage as starch so as to get a subnetwork. If two steady states can be obtained for this modified system we may be able to get steady states for the original system using the implicit function theorem. There are some more complications but the a key step in the proof is the one just described. So how do we get steady states for the modified system? The idea is to solve many of the equations explicitly so that the problem reduces to a single equation for one unknown, the concentration of DHAP. (When applying the implicit function theorem we have to use a system of two equations for two unknowns.) In the end we are left with a quadratic equation and we can arrange for the coefficients in that equation to have convenient properties by choosing the parameters in the dynamical system suitably. This approach can be put in a wider context using the concept of stoichiometric generators but the proof is not logically dependent on using the theory of those objects.

Having got some information about the Pettersson model we may ask what happens when we go over to the Poolman model. The Poolman model is a system of ODE from the start and so we do not have any conceptual problems in that case. The method of construction of steady states can be adapted rather easily so as to apply to the system of DAE related to the Poolman model (let us call it the reduced Poolman model since it can be expressed as a singular limit of the Poolman model). The result is that there are parameter values for which the reduced Poolman model has at least three steady states. Whether the Poolman model itself can have three steady states is not yet clear since it is not clear whether the transverse eigenvalues (in the sense of GSPT) are all non-zero.

By analogy with known facts the following intuitive picture can be developed. Note, however, that this intuition has not yet been confirmed by proofs. In the picture one of the positive steady states of the Pettersson model is stable and the other unstable. Steady states on the boundary where some concentrations are zero are stable. Under the perturbation from the Pettersson model to the reduced Poolman model an additional stable positive steady state bifurcates from the boundary and joins the other two. This picture may be an oversimplification but I hope that it contains some grain of truth.