Archive for the ‘mathematical biology’ Category

Codimension three bifurcations

July 13, 2021

In our recent work with Lisa Kreusser on Lck we observed that it can happen that to reach a certain goal it is easier to work with a bifurcation of higher codimension although this means working with an object which is intrinsically more complicated. In that specific case we used a Bogdanov-Takens bifurcation (codimension two) instead of a Hopf bifurcation (codimension one). The goal that we reached was to prove the existence of a periodic solution, which turned out to be unstable. We did not prove the existence of a stable periodic solution of the system being considered. I wondered if this strategy could be pushed further, going to a bifurcation of codimension three.

The general idea is as follows. We have a dynamical system \dot x=f(x,\alpha) depending on parameters. We want to compare it with a normal form which is a standard model dynamical system depending on parameters. Here everything is local in a neighbourhood of (0,0), which is a steady state. The number of parameters in the normal form is the codimension. In the end we want to conclude that dynamical features which are present in the normal form are also present in the original system. The normal form for the Hopf bifurcation contains periodic solutions. In fact there are two versions depending on a parameter which can take the value plus or minus one. These are known as supercritical and subcritical. In the supercritical case the periodic solutions are stable, in the subcritical case unstable. The normal form for the Bogdanov-Takens bifurcation also has two cases which may, by analogy with the Hopf case, be called super- and subcritical. The fact we used is that when a Bogdanov-Takens bifurcation exists there is always a Hopf bifurcation close to it. Moreover the super- or subcritical nature of the bifurcation is inherited. What I wanted to investigate is whether the presence of a suitable bifurcation of codimension three could imply the existence of a Bogdanov-Takens bifurcation close to it and whether it might even be the case that both super- and subcritical Bogdanov-Takens bifurcations are obtained. If this were the case then it would indicate an avenue to a statement of the kind that the existence of a generic bifurcation of codimension three of a certain type implies the existence of stable periodic solutions.

The successful comparison of a given system with the normal form relies on certain genericity conditions being satisfied. These are of two types. The first type is a condition on the system defined by setting \alpha=0. It says that the system does not belong to a certain degenerate set defined by the vanishing of particular functions of the derivatives of f(x,0) at x=0. In an abstract way this degenerate set may be thought of a subset of the set of k-jets of functions at the origin. This subset is a real algebraic variety which is a union of smooth submanifolds (strata) of different dimensions which fit together in a nice way. The set defining the bifurcation itself is also a variety of this type and the degenerate set consists of all strata of the bifurcation set except that of highest dimension. The family f(x,\alpha) induces a mapping from the parameter space into the set of k-jets. The second type of genericity condition says that this mapping should be transverse to the bifurcation set. Because of the genericity condition of the first type the image of x=0 lies in the stratum of highest dimension in the bifurcation set, which is a smooth submanifold. The transversality condition says that the sum of the tangent space to this manifold and the image of the linearization at x=0 of the mapping defined by f is the whole space. Writing about these things gives me a strange feeling since they involve concepts which I used in a very different context in my PhD, which I submitted 34 years ago, and not very often since then.

The following refers to a two-dimensional system. In the Bogdanov-Takens case the bifurcation set is defined by the condition that the linearization about (0,0) has a double zero eigenvalue. The genericity condition of the first type involves the 2-jet. The first part BT.0 (cf. this post for the terminology) says that the linearization should not be identically zero. The second part consists of conditions BT.1 and BT.2 involving the second derivatives. The remaining condition BT.3 is the transversality condition (genericity condition of the second type). Now I want to study bifurcations which satisfy the eigenvalue condition for a Bogdanov-Takens bifurcation and the condition BT.0 but which may fail to satisfy BT.1 or BT.2. At one point the literature on the subject appeared to me like a high wall which I did not have the resolve to try to climb. Fortunately Sebastiaan Janssens pointed me to a paper of Kuznetsov (Int. J. Bif. Chaos, 15, 3535) which was a gate through the wall and gave me access to things on the other side which I use in the following discussion. Following Dumortier et al. (Erg. Th. Dyn. Sys. 7, 375) we look at cases where the 2-jet of the right hand side of the equation for \dot y is of the form \alpha x^2+\beta xy. It turns out that we can divide the problem into two cases, in each of which one of the coefficients \alpha and \beta is zero and the other non-zero. We first concentrate on the case \beta=0 which is that studied in the paper just quoted. There it is referred to as the cusp case. It is the case where BT.2 is satisfied and BT.1 fails. As explained in the book of Dumortier et al. (Bifurcations of planar vector fields. Nilpotent singularites and Abelian integrals.) it is useful to further subdivide the case \alpha=0 into three subcases, known as the saddle, focus and elliptic cases. These cases give rise to different normal forms (partially conjectural). A bifurcation diagram for the cusp case can be found in the paper of Dumortier et al. Since we are dealing with a codimension 3 bifurcation the bifurcation diagram should be three-dimensional. It turns out, however, that the bifurcation set has a conical structure near the origin, so that its structure is determined by its intersection with any small sphere near the origin. This gives rise to a two-dimensional object which can be well represented in the plane. Note, however, that passing from the sphere to the plane necessarily involves discarding a ‘point at infinity’. It is this intersection which is represented in Fig. 3 of the paper of Dumortier et al. That this is an accurate representation of the bifurcation diagram in this case is the main content of the paper of Dumortier et al.

In the bifurcation diagram in the paper of Dumortier et al. we see that there are two Bogdanov-Takens points b_1 and b_2 which are joined by a curve of Hopf bifurcations. The points b_1 and b_2 are stated to be codimension 2 and this indicates that they are non-degenerate Bogdanov-Takens points. There is one exceptional point c_2 on the curve of Hopf bifurcations which is stated to be of codimension 2. This indicates that it is a non-degenerte Bautin bifurcation point. This in turn implies that one of the Bogdanov-Takens points is supercritical and the other subcritical. Thus in this case there exist both stable and unstable periodic solutions in each neighbourhood of the bifurcation point. This indicates that the strategy outlined at the beginning of this post can in principle work. Whether it can work in practise in a given system depends on how difficult it is to check the relevant non-degeneracy conditions. I end with a brief comment on the case \alpha=0. The book of Dumortier et al. presents conjectured bifurcation diagrams for all subcases but does not claim complete proofs that they are all correct. In all cases we have two Bogdanov-Takens points joined by a curve of Hopf bifurcations on thich there is precisely one Bautin point. Thus in this sense we have the same configuration as in the case \beta=0. In a discussion of these matters in the paper of Kuznetsov mentioned above (which is from 2005) it is stated that a complete proof is only available in the saddle case. I do not know if there has been any further progress on this since then.

A model for the Calvin cycle with diffusion

June 29, 2021

In the past I have spent a lot of time thinking about mathematical models for the Calvin cycle of photosynthesis. These were almost exclusively systems of ordinary differential equations. One thing which was in the back of my mind was a system involving diffusion written down in a paper by Grimbs et al. and which I will describe in more detail in what follows. Recently Burcu Gürbüz and I have been looking at the dynamics of solutions of this system in detail and we have just produced a preprint on this. The starting point is a system of five ODE with mass action kinetics written down by Grimbs et al. and which I call the MA system (for mass action). That system has only one positive steady state and on the basis of experimental data the authors of that paper were expecting more than one. To get around this discrepancy they suggested a modified system where ATP is introduced as an additional variable and the diffusion of ATP is included in the model. I call this the MAd system (for mass action with diffusion). Diffusion of the other five chemical species is not included in the model. The resulting system can be thought of as a degenerate system of six reaction reaction-diffusion equations or as a system where five ODE are coupled to one (non-degenerate) diffusion equation. The idea was that the MAd system might have more than one steady state, with an inhomogenous steady state in addition to the known homogeneous one. Experiments which measure only spatial averages would not detect the inhomogeneity. To be relevant for experiments the steady states should be stable.

For simplicity we consider the model in one space dimension. The spatial domain is then a finite closed interval and Neumann boundary conditions are imposed for the concentration of ATP. We prove that for suitable values of the parameters in the model there exist infinitely many smooth inhomogeneous steady states. It turns out that all of these are (nonlinearly) unstable. This is not a special feature of this system but in fact, as pointed out in a paper of Marciniak-Czochra et al. (J. Math. Biol. 74, 583), it is a frequent feature of systems where ODE are coupled to a diffusion equation. This can be proved using a method discussed in a previous post which allows nonlinear instability to be concluded from spectral instability. We prove the spectral instability in our example. There may also exist non-smooth inhomogeneous steady states but we did not enter into that theme in our paper. If stable inhomogeneous steady states cannot be used to explain the experimental observations what alternatives are available which still keep the same model? If experiments only measure time averages then an alternative would be limit sets other than steady states. In this context it would be interesting to know whether the system has spatially homogeneous solutions which are periodic or converge to a strange attractor. A preliminary investigation of this question in the paper did not yield definitive results. With the help of computer calculations we were able to show that it can happen that the linearization of the (spatially homogeneous) system about a steady state has non-real eigenvalues, suggesting the presence of at least damped oscillations. We proved global existence for the full system but we were not able to show whether general solutions are bounded in time, not even in L^1. It can be concluded that there are still many issues to be investigated concerning the long-time behaviour of solutions of this system.

Dynamics of the activation of Lck

January 26, 2021

The enzyme Lck (lymphocyte-associated tyrosine kinase) is of central importance in the function of immune cells. I hope that mathematics can contribute to the understanding and improvement of immune checkpoint therapies for cancer. For this reason I have looked at the relevant mathematical models in the literature. In doing this I have realized the importance in this context of obtaining a better understanding of the activation of Lck. I was already familiar with the role of Lck in the activation of T cells. There are two tyrosines in Lck, Y394 and Y505, whose phosphorylation state influences its activity. Roughly speaking, phosphorylation of Y394 increases the activity of Lck while phosphorylation of Y505 decreases it. The fact that these are influences going in opposite directions already indicates complications. In fact the kinase Lck catalyses its own phosphorylation, especially on Y394. This is an example of autophosphorylation in trans, i.e. one molecule of Lck catalyses the phosphorylation of another molcule of Lck. It turns out that autophosphorylation tends to favour complicated dynamics. It is already the case that in a protein with a single phosphorylation site the occurrence of autophosphorylation can lead to bistability. Normally bistability in a chemical reaction network means the existence of more than one stable positive steady state and this is the definition I usually adopt. The definition may be weakened to the existence of more than one non-negative stable steady state. That autophosphorylation can produce bistability in this weaker sense was already observed by Lisman in 1985 (PNAS 82, 3055). He was interested in this as a mechanism of information storage in a biochemical system. In 2006 Fuss et al. (Bioinformatics 22, 158) found bistability in the strict sense in a model for the dynamics of Src kinases. Since Lck is a typical member of the family of Src kinases these results are also of relevance for Lck. In that work the phosphorylation processes are embedded in feedback loops. In fact the bistability is present without the feedback, as observed by Kaimachnikov and Kholodenko (FEBS J. 276, 4102). Finally, it was shown by Doherty et al. (J. Theor. Biol. 370, 27) that bistability (in the strict sense) can occur for a protein with only one phosphorylation site. This is in contrast to more commonly considered phosphorylation systems. These authors have also seen more complicated dynamical behaviour such as periodic solutions.

All these results on the properties of solutions of reaction networks are on the level of simulations. Recently Lisa Kreusser and I set out to investigate these phenomena on the level of rigorous mathematics and we have just put a paper on this subject on the archive. The model developed by Doherty et al. is one-dimensional and therefore relatively easy to analyse. The first thing we do is to give a rigorous proof of bistability for this system together with some information on the region of parameter space where this phenomenon occurs. We also show that it can be lifted to the system from which the one-dimensional system is obtained by timescale separation. The latter system has no periodic solutions. To obtain bistability the effect of the phosphorylation must be to activate the enzyme. It does not occur in the case of inhibition. We show that when an external kinase is included (in the case of Lck there is an external kinase Csk which may be relevant) and we do not restrict to the Michaelis-Menten regime bistability is restored.

We then go on to study the dynamics of the model of Kaimachnikov and Kholodenko, which is three-dimensional. These authors mention that it can be reduced to a two-dimensional model by timescale separation. Unfortunately we did not succeed in finding a rigorous framework for their reduction. Instead we used another related reduction which gives a setting which is well-behaved in the sense of geometric singular perturbation theory (normally hyperbolic) and can therefore be used to lift dynamical features from two to three dimensions in a rather straightforward way. It then remains to analyse the two-dimensonal system. It is easy to deduce bistability from the results already mentioned. We go further and show that there exist periodic and homoclinic solutions. This is done by showing the existence of a generic Bogdanov-Takens bifurcation, a procedure described more generally here and here. This system contains an abundance of parameters and the idea is to fix these so as to get the desired result. First we choose the coordinates of the steady state to be fixed simple rational numbers. Then we fix all but four of the parameters in the system. The four conditions for a BT bifurcation are then used to determine the values of the other four parameters. To get the desired positivity for the computed values the choices must be made carefully. This was done by trial and error. Establishing genericity required a calculation which was complicated but doable by hand. When a generic BT bifurcation has been obtained it follows that there are generic Hopf bifurcations nearby in parameter space and the nature of these (sub- or supercritical) can be determined. It turns out that in our case they are subcritical so that the associated periodic solutions are unstable. Having proved that periodic solutions exist we wanted to see what a solution of this type looks like by simulation. We had difficulties in finding parameter values for which this could be done. (We know the parameter values for the BT point and that those for the periodic solution are nearby but they must be chosen in a rather narrow region which we do not know explicitly.) Eventually we succeeded in doing this. In this context I used the program XPPAUT for the first time in my life and I learned to appreciate it. I see this paper as the beginning rather than the end of a path and I am very curious as to where it will lead.

Paper on hepatitis C

November 13, 2020

Together with Alexis Nangue and other colleagues from Cameroon we have just completed a paper on mathematical modelling of hepatitis C. What is being modelled here is the dynamics of the amount of the virus and infected cells within a host. The model we study is a modification of one proposed by Guedj and Neumann, J. Theor. Biol. 267, 330. One part of it is a three-dimensional model, which is related to the basic model of virus dynamics. It is coupled with a two-dimensional model which gives a simple description of the way in which the virus replicates inside a host cell. This intracellular process is related to the mode of action of the modern drugs used to treat hepatitis C. I gave some information about the disease itself and its treatment in a previous post. In the end the object of study is a system of five ODE with many parameters.

For this system we first proved global existence and boundedness of the solutions, as well as positivity (positive initial data lead to a positive solution). One twist here is a certain lack of regularity of the coefficients of the system. When some of the unknowns become zero the right hand side of the equations is discontinuous. This means that it is necessary to prove that this singular set is not approached during the evolution of a positive solution. The source of the irregularity is the use of something called the standard incidence function instead of simple mass action kinetics. The former type of kinetics has a long history in epidemiology and I do not want to try to explain the background here. In any case, there are arguments which say that mass action kinetics leads to unrealistic results in within-host models of hepatitis and that the standard incidence function is better.

We show that the model has up to two virus-free steady states and determine their stability. The study of positive steady states is more difficult and, at the moment, incomplete. We have proved that there cannot be more than three steady states but we do not know if there is ever more than one. Under increasingly restrictive assumptions on the parameters (restrictions which are unfortunately not all biologically motivated) we show that there is at least one positive steady state, that there is exactly one and that that one is asymptotically stable. Under certain other assumptions we can show that every solution converges to a steady state. This last proof uses the method of Li and Muldowney discussed in a previous post. Learning about this method was one of the (positive) side effects of working on this project. Another was an improvement of my understanding of the concept of the basic reproductive number as discussed here and here. During this project I have learned a lot of new things, mathematical and biological, and I feel that I am now in a stronger position to tackle other projects modelling hepatitis C and other infectious diseases.


Talk by David Ho on COVID-19

August 22, 2020

On Thursday David Ho gave a keynote lecture at the SMB conference. He talked about work to develop monoclonal antibodies against SARS-CoV-2. He started by apologising, in view of the given audience, that there would be no mathematics in his talk but he did make make clear his continuing belief in the importance of applying mathematics to biology. He has been leading an effort with a precise medical goal – to find effective neutralising antibodies against this new virus. Antibodies were obtained from five patients severely ill with COVID-19. Four of them survived while one later died of the disease. These antibodies were then analysed by biochemical and bioinformatic means to find those which bound best to the spike protein of the virus. In this context I learned some basic things about the virus. The spike, which is used by the virus to enter cells is considered the number one target for antibodies which could be effective in combating the disease. More precisely there are two different subdomains which are possible targets, one more at the tip of the spike (the receptor-binding domain) and another more on the sides (the N-terminal domain), which is a trimer. A number of antibodies were found which bind to the first subdomain or to one of the subunits of the second. Another was found whose binding site is somewhat less local. This whole process was carried out in just few weeks, a remarkable achievement.

The antibodies just mentioned are the therapeutic candidates. The idea is to either produce monoclonal antoibodies with these sequences or possibly versions which are improved so as to be longer-lived. Monoclonal antibodies are known to be extremely expensive when used to treat other diseases, such as cancer. They are also expensive in the present context, but the speaker said that the retail cost depends very much on the quantity produced. In other applications the number of patients is relatively small and the cost correspondingly high. If the antobodies were being used for a very large number of patients the cost would be lower. It would remain problematic for low and middle income countries. It has been discussed that the Gates foundation might make it possible to offer this treatment in poorer countries for fifty dollars a dose. The main advantage of this method compared with that of trying to use antibodies from the serum of patients directly is that it is much more practical to apply on a very large scale. The effectiveness of the antibodies against the disease has been tested in hamsters. Ho made the impression of someone tackling a major problem of humanity head on with some of the best tools available. He said that an article giving an account of the work had appeared in Nature (Potent neutralizing antibodies against multiple epitopes on SARS-CoV-2 spike, Nature 584, 450). In response to one question on one aspect of the treatment he said that the answer was not known but he would just be continuing to a Zoom meeting of researchers leading the attempts to develop therapies which was to discuss exactly that question.

My first virtual conference (SMB 2020)

August 20, 2020

At the moment I am attending the annual conference of the Society for Mathematical Biology, which is taking place online. This is my first experience of this kind of format. The conference has many more participants than in any previous year, more than 1700. It takes place in a virtual building which is generated by the program Sococo. I find this environment quite disorienting and a bit stressful. This reaction probably has to do with the facts that I am no longer so young and that I have always tried to avoid social media as much as possible. I am sure that younger generations (and members of older generations with an enthusiasm for new technical developments) have far fewer problems getting used to it. In advance I was a bit worried about setting up the necessary computer requirements to be able to give my talk or even to go to others. In the end it worked out and my talk, given via Zoom, went smoothly. I got some good feedback, I am already convinced that it was worth joining this meeting and I may be less sceptical about joining others of this type in the future. There have been technical hitches. For instance the start of one big talk was delayed by about 20 minutes for a reason of this kind. Nevertheless, many things have gone well. Of course it is much preferable to meet people personally but when that is not possible virtual meetings with old friends are also pleasant.

This meeting has been organized around the subgroups of the society. I am a member of the subgroups for immunology and oncology. In fact my talk got scheduled in the group for mathematical modelling. My greatest allegiance is to the immunology subgroup and so I was happy to see that it was represented so strongly at this conference. It has seven sessions of lectures, made up of 28 talks. If I should choose my favourite from those talks in this section I have heard so far (it is not finished yet) then I pick the talk of Ruy Ribeiro on CD8 cells in HIV infection. I did feel I was missing some necessary background but I nevertheless found the talk useful in introducing me to important ideas which were new to me. The immunology subgroup have managed to get a very prominent speaker for their keynote talk, David Ho. I wrote about his work and its relation to mathematics in one of my first ever posts on this blog, way back in 2008. I am looking forward to hearing his talk later today (which will be on COVID-19). I will now mention some of my other personal highlights from the conference so far. I noticed on the program that Stas Shvartsman was giving a talk. His work has made a positive impression on me in the past and I did have a little e-mail contact with him. On the other hand I never met him personally and I had never heard a talk by him. This was my chance and I went in with high expectations. They were not disappointed. He talked about developmental defects arising from mutations in a single gene. In particular he concentrated on mutations in the Raf-MEK-ERK MAPK cascade. I was familiar with the role of mutations in this cascade in cancer but I had never heard about this other role. Stas described experiments in Drosophila where one base in MEK is mutated. This produces flies with a particular small change in the pattern of the veins in their wings. Interestingly, this does not occur in all flies with the mutation but only in 30% of them (I hope I am remembering the right number). Another talk yesterday which I appreciated was that by Robert Insall, who was talking about aspects of chemotaxis. He started with the phenomenon of the spread of melanoma. Melanoma cells do have a very strong tendency to spread in space and the question is what controls that. Is it chemotaxis along a chemical gradient?He showed pictures of melanoma cells moving fast up a chemical gradient. Then he showed a similar picture showing them moving just as fast without the chemical gradient. So what is going on? This kind of experiment starts with cells concentrated on one side of a region and a spatially homogeneous distribution of a relevant substance. The cells consume this substance, create their own gradient and then undergo chemotaxis along it. The speaker explained how there are many instances of chemotaxis in biology which can only be explained by a self-created gradient and not be a pre-existing one. He also made some interesting remarks about how mathematical modelling can lead to insights in biology which would not be possible with the usual verbal approaches of the biologists.

I also want to mention an interesting conversation I had in the poster session. The poster concerned was that of Daniel Koch. The theme of his work (which has already been published in a paper in J. Theor. Biol.) is that the formation of oligomers of proteins (or their posttranslationally modified variants) can lead to interesting dynamics. At first sight this may sound too simple to be interesting but in fact in mathematics it is often the careful consideration of apparently simple situations which leads to fundamental progress. I imagine that this principle also applies to other disciplines (such as biology) but it is perhaps strongest in mathematics. In any case, I am strongly motivated to study this work carefully. The only question in when it will be, given the many other directions I want to pursue.

The basic reproduction number for infectious diseases, part 2

July 9, 2020

Here I continue the discussion of the previous post. There I mentioned that when defining a basic reproduction number it is necessary, in addition to defining a set of ODE, to make a choice of infected and uninfected compartments. This means partitioning the unknowns into two subsets. In addition it is necessary to distinguish processes regarded as new infections from others. This means splitting the right hand side of the equations into a sum of two terms. There are also some other biologically motivated assumptions on the right hand side. As indicated in the previous post, a reproduction number is a feature of a disease-free steady state of the model (i.e. a steady state where the unknowns belonging to the infected compartment are zero). In fact it is a feature of the linearization of the model about that point. Intuitively, it has to do with a situation where a population contains a very small number of infected individuals. The matrix defining the right hand side of the linearization is a sum of two terms, each of which is partitioned in a certain way. We now consider a situation where the disease-free steady state is perturbed by introducing a small number of individuals into each of the infected compartments while supposing that they cannot cause a significant number of secondary infections. The numbers of individuals in the infected compartments then satisfy a linear system. Under the given assumptions all eigenvalues of the matrix defining this system have negative real parts. Thus the solution tends to zero exponentially. We are interested in the average number of individuals in these compartments over all positive times. To understand this consider first the simpler example of a single compartment with an exponential behaviour x(t)=x_0e^{-at}. The average value of x is \int_0^\infty x_0e^{-at}dt=\left[-a^{-1}x_0e^{-at}\right]_0^\infty=x_0a^{-1}. In fact what we are interested in is the number of infections expected from the individuals in the infectious compartment. We have the equation \psi'(t)=-V\psi (t), where all eigenvalues of V have positive real parts and we want to calculate F\int_0^\infty e^{-Vt}dt=FV^{-1}\psi (0), where F describes the process of new infections. The matrix FV^{-1} is the next generation matrix mentioned in the previous post. It is now a matter of some rather involved linear algebra to show that the modulus of the largest eigenvalue of this matrix is greater than one if and only if the greatest real part of an eigenvalue of the linearization of the system of ODE is greater than zero.

In the current reporting of epidemiological data the quantity R is time-dependent. I am far away from being able to link this with the above discussion. I looked at the web page of the Robert Koch Institute to see what they say about how they produce their curves. This led me to a paper of Cintron-Arias et al. (Math. Biosci. Eng. 6, 261), which might be able to provide a gateway into this but it is clear that it involves many things which I do not understand at present. Perhaps I should first look at the classic book ‘Infectious Diseases of Humans’ by Anderson and May. Robert May died in April 2020 so that he did not have the opportunity to observe the present pandemic.

The basic reproduction number for infectious diseases

June 12, 2020

These days reproduction numbers for epidemiology are prominent in the popular media. Many people are familiar with the idea that stopping a resurgence of COVID-19 infections in a region has to do with making and keeping something called R less than one. They may also be familiar with the informal definition of R that it is the number of new infections caused by an infected individual. But how is R (or R_0 as it is more commonly called by scientists) defined? A mathematician might not expect to find an answer in the media but it might be reasonable to expect one in the scientific literature on epidemiology. In the past I have been frustrated by the extent to which this fails to be the case. What is typically given is a description in words which I never found possible to convert into a precise mathematical account, despite considerable effort. Now, in the context of a project on hepatis C which I have been working on with colleagues from Cameroon, my attention was drawn to a paper of van den Driessche and Watmough (Math. Biosci. 180, 29) which contains some answers of the type I was looking for. I was vaguely aware of this paper before but I had never seriously tried to read it because I did not realise its nature.

The context in which I would have liked to find answers is that of models given by systems of ordinary differential equations where the unknowns are the numbers of individuals in different categories (susceptible, infected, recovered etc.) as functions of time. How the numbers reported in the media are calculated (on the basis of discrete data) is something I have not yet tried to find out. At the moment I would like an answer in the context which bothered me in the past and this is the context treated in the paper mentioned above. A typical situation is that found in the basic model of virus dynamics, a system of three ODE describing the dynamics of a virus within a host, with the unknowns being uninfected cells, infected cells and virions. There is a quantity R_0 which can be expressed in terms of the coefficients of the system. If R_0\le 1 then the only non-negative steady state is virus-free. This is the uninfected state and it is globally asymptotically stable. If R_0>1 there is an uninfected state which is unstable and an infected state which is positive and globally asymptotically stable. This kind of situation is not unique to this example and similar things are seen in many models of infection. There is a reproductive number (or perhaps more than one) which defines a threshold between different types of late-time behaviour.

It is not obvious that the analysis of van den Driessche and Watmough applies to models of in-host dynamics of a pathogen since it is necessary to make a choice of infected and uninfected compartments which is related to the biological interpretation of the variables and not just to the mathematical structure of the model. Their analysis does apply to the basic model of virus dynamics if the infected compartments are chosen to be the infected cells and virions and the reaction fluxes are partitioned in a suitable way. The simple picture of the significance of the reproductive number given above does not always hold. There is also another scenario which can occur and does so in many practical examples and involves the notion of a backward bifurcation. It goes as follows. For R_0 sufficiently small the disease-free steady state is globally asymptotically stable but as R_0 is increased this property breaks down before R_0=1 is reached. A fold bifurcation occurs which creates a stable and an unstable positive steady state. The unstable steady state moves so as to meet the disease-free state when R_0=1. For R_0>1 there are exactly two steady states and the positive one is globally asymptotically stable. There is bifurcation for R_0=1 but it has a different structure from that in the classical scenario (which is a transcritical bifurcation). It bears some resemblance to a sub-critical Hopf bifurcation.

The most useful insights I got from reading and thinking about the paper of van den Driessche and Watmough are as follows. The primary significance of R_0 concerns the disease-free steady state and its stability. The fact that it can sometimes characterise the stability of an infected steady state is a kind of bonus which does not apply to all models. What it can do is to provide information about the stability of a positive steady state in a regime where it is close to the bifurcation point where it separates from the disease-free steady state. This circumstance is analysed in the paper using centre manifold theory. The significance for the stability of a steady state which is far away is a weak one. Continuity arguments can be used to propagate information about stability through parameter space but only as long as no bifurcations happen. When this is the case depends on the details of the particular example being considered. What is the definition of R_0 given in the paper? It is the largest modulus of an eigenvalue of a certain matrix (the next generation matrix) constructed from the linearization of the system about the disease-free steady state, whereby the construction of this matrix incorporates information about the biological meaning of the variables. Consider the example of the basic model of virus dynamics with the choices of infected and uninfected compartments as above. There is more than way of partitioning the reaction fluxes. I first tried to put both the production of infected cells and the production of virions into the category of fluxes called \cal F in the paper. Applying the definition of the reproduction number given there leads to \sqrt{R_0}, where R_0 is the reproduction number usually quoted for this model. If instead only the production of infected cells is put into the category \cal F then the general definition gives the conventional answer R_0. The two quantities defining the threshold are different but the definition of being above or below the threshold are the same. (\sqrt{x}<1\Leftrightarrow x>1). That this kind of phenomenon can occur is shown by example in the paper.

The Higgins-Selkov oscillator, part 3

March 17, 2020

Here I report on some further progress in understanding the dynamics of mathematical models for glycolysis. In a previous post I described some results on the Selkov model obtained in a paper by Pia Brechmann and myself. The main question left open there is whether this model admits unbounded oscillations. We have now been able to resolve this issue in a new paper, showing that there is precisely one value of the parameter \alpha for which solutions of this type exist. This progress was stimulated by a paper of Merkin, Needham and Scott (SIAM J. Appl. Math. 47, 1040). It was a piece of luck that we found this paper since its title gives no indication that it has any relation to the problem we were interested in. It concerns a model for certain chemical reactors which turns out to be mathematically identical to the Selkov model with the parameter \gamma equal to two. The authors claim that there are unbounded solutions for exactly one value of \alpha and present a derivation of this which is an intricate argument involving matched asymptotic expansions. I find this argument impressive but I have no idea how it could be transformed into a rigorous proof. Despite this it helped us indirectly. We had previously tried, without success, to find some way to transform the equations so as to obtain a well-behaved limit \alpha\to\infty. Combining some transformations in the paper of Merkin et al. with some we had used in our previous paper allowed this goal to be achieved. Once this has been done the dynamics is under control for \alpha sufficiently large. The next step is to do a shooting argument to show that there is a parameter value \alpha_1 for which a heteroclinic cycle at infinity exists. It had already been shown in our first paper that this is enough to conclude the occurrence of unbounded oscillations. The paper of Merkin et al. also helped us to understand what was the right set-up for the shooting argument. The uniqueness of \alpha_1 was obtained using a monotonicity property of the parameter dependence which appears to be new. At this point I think it is justified to say that the analysis of the main qualitative features of solutions of the Selkov oscillator is essentially complete. There is just one more thing I would like to know although I see it just as a curiosity. We showed that for \alpha a bit less than \alpha_1 there exists a stable periodic solution and that its diameter tends to infinity as \alpha\to\alpha_1. Merkin et al. give an expression for the leading order contribution to the diameter in this limit for \gamma=2. How could this be proved rigorously? How could an analogous expression be derived for other values of \gamma?

In the original paper of Selkov the system we have been discussing up to now is derived from a system with kinetics of Michaelis-Menten type by letting a parameter \nu tend to zero. He suggests that in the system with \nu>0 the unbounded oscillations are eliminated so that they could be thought of as an artefact of setting \nu=0. In our new paper we investigated this question although we were unfortunately only able to get partial results. One feature of the Selkov model not mentioned by Selkov in the original paper is that there are solutions which are unbounded and tend to infinity in a monotone manner. We showed that solutions of this type also exist for the Michaelis-Menten model and that they have exactly the same leading order asymptotics as in the case \nu=0. The periodic solutions of the Selkov oscillator are stable and arise in a supercritical Hopf bifurcation. That system admits no unstable periodic solutions whatsoever. By contrast, we found that the Michaelis-Menten system also exhibits subcritical Hopf bifurcations and correspondingly unstable periodic solutions.

I see a number of interesting directions in which this work could be extended but we will not attempt to do so in the foreseeable future since we have too many other projects to work on. I would only return to this if there turned out to be intriguing connections to other research projects or to the themes of master theses I was supervising.

Monotone systems revisited

December 4, 2019

There are some topics in mathematics and physics which are a lasting source of dissatisfaction for me since I feel that I have not properly understood them despite having made considerable efforts to do so. In the case of physics the reason is often that the physicists who understand the subject are not able to explain it in a way which provides what a mathematician sees as a comprehensible account. In mathematics the problem is a different one. Mathematicians frequently have a tendency (often justified) to discuss things on a level which is as general as possible. This leads to theorems which are loaded down with detail and where the many technical conditions make it difficult to see the wood for the trees. When confronted with such things I sometimes feel exhausted and give up. I prefer an account which builds up ideas step by step from simple beginnings. Here I return to a subject which I have written about more than once in this blog before but where the sense of dissatisfaction remains. I hope to reduce it here.

I start with a system of ordinary differential equations \dot x_i=f_i(x). It should be defined on the n-dimensional Euclidean space or on one of its orthants. (An orthant is the subset of Euclidean space defined by making a choice of the signs of its components. It generalises a quadrant in the two-dimensional case.) The system is said to be cooperative if \frac{\partial f_i}{\partial x_j}>0 for all i\ne j. The name comes from the fact that the equations for the population dynamics of a set of species has this property if each species benefits the others. Suppose we now have two solutions x and \bar x of the system and that x_i(t_0)\le\bar x_i(t_0) for all i at some time time t_0. We may abbreviate this relation by x(t_0)\le\bar x(t_0). Here we see a partial order on Euclidean space defined by the ordering of the components. A theorem of Müller and Kamke says that if the initial data for two solutions of a cooperative system at time t_0 satisfies this relation then x(t)\le\bar x(t) for all t\ge t_0. Another way of saying this is that the time-t flow of the system is preserves the partial order. A system of ODE with this property is called monotone. Thus the Müller-Kamke theorem says that a cooperative system is monotone.

The differential condition for monotonicity can be integrated. If x and \bar x are two points in Euclidean space with x_i=\bar x_i for a certain i and x_j\le\bar x_j for j\ne i then f_i(x)\le f_i(\bar x). To see this we join x to \bar x by a piecewise linear curve where the coordinates other than the ith are increased successively from x_j to \bar x_j. On each segment of this curve the value of f_i does not decrease, as a consequence of the fundamental theorem of calculus. Hence its value at the end of the entire path is at least as big as its value at the beginning. We now want to prove that a certain inequality holds at all times t\ge t_0. In order to do this we would like to consider the first time t_*>t_0 where the inequality fails and get a contradiction. Unfortunately there might be no such time – in principle the condition might fail immediately. To get around this we deform the system for the solution \bar x to \frac{d\bar x_i}{dt}=f_i(\bar x)+\epsilon. If we can prove the result for the deformed system the result for the initial system follows by continuous dependence of the solution on \epsilon. For the deformed system let t_* be the supremum of the times where the desired inequality holds. If the inequality does not hold globally then the system is still defined at t=t_*. For t=t_* we have x_i=\bar x_i for some i and we can assume w.l.o.g. that x_j<\bar x_j for some j since otherwise the two solutions would be equal and the result trivial. The integrated form of the cooperativity condition implies that at t_* the right hand side of the evolution equation for \bar x_i-x_i is positive. On the other hand the fact that it just reached zero coming from positive values implies that the right hand side of the evolution equation is non-positive and we get a contradiction.

A key source of information about monotone dynamical systems is the book of Hal Smith with this title. I have repeatedly looked at this book but always got bogged down quite quickly. Now I realise that for my purposes it would have been much better if I had started with chapter 3. The Müller-Kamke theorem is discussed in section 3.1. The range of application of this theorem can be extended considerably by the following trick, discussed in section 3.5. Suppose that we define y_i=(-1)^{m_i}x_i where each of the m_i are zero or one. This transforms the signs of Df in a certain way and so cooperativity of the system for y corresponds to a certain sign pattern for the entries of Df. A first important condition is that each off-diagonal element of Df(x) should be either non-negative or non-positive. Next, the sign of \frac{\partial f_i}{\partial x_j}\frac{\partial f_j}{\partial x_i} is not changed be the transformation and must thus be non-negative. In the context of population models this can be interpreted as saying that there is no pair of species which are in a predator-prey relationship. Given that these two conditions are satisfied we consider a labelled graph where the nodes are the numbers from 1 to n and there is an edge between two nodes if at least one of the corresponding partial derivatives is non-zero at some point. The edge is then labelled with the sign of this non-zero value. A loop in the graph can be assigned the sign which is the product of those of its edges. It turns out that a system can be transformed to a cooperative system in the way indicated if and only if the graph contains no negative loops. I will call a system of this type ‘cooperative up to sign reversal’. The system can be transformed by a permutation of the variables into one where Df has diagonal blocks with non-negative entries and off-diagonal elements with non-positive entries.

If all elements of Df are required to be non-positive we get the class of competitive systems. It should be noted that being competitive leads to less restrictions on the dynamics of a system (towards the future) than being cooperative. We can define a class of systems which are competitive up to sign reversal. An example of such a system is the basic model of virus dynamics. In that system the unknowns are the populations of uninfected cells x, infected cells y and virus particles v. The transformation y\mapsto -y makes it into a competitive system. In various models of virus dynamics including the immune response the target cells of the virus and the immune cells are in a predator-prey relationship and so these systems can be neither cooperative up to sign or competitive up to sign.