Is mathematics being driven out by computers?

In the past two weeks I attended two conferences. The first was the annual meeting of the Deutsche Mathematikervereinigung (DMV, the German mathematical society) in Berlin. The second was the joint annual meeting of the ESMTB (European Society for Mathematical and Theoretical Biology) and the SMB (Society for Mathematical Biology) in Heidelberg. I had the impression that the participation of the SMB was relatively small compared to previous years. (Was this mainly due to the pandemic or due to other problems in international travel?) There were about 500 participants in total who were present in person and about another 100 online. I was disappointed with the plenary talks at both conferences. The only one which I found reasonably good was that of Benoit Perthame. One reason I did not like them was the dominance of topics like machine learning and artificial intelligence. This brings me to the title of this post. I have the impression that mathematics (at least in applied areas) is becoming ever weaker and being replaced by the procedure of developing computer programmes which could be applied (and sometimes are) to the masses of data which our society produces these days. This was very noticeable in these two conferences. I would prefer if we human beings would continue to learn something and not just leave it to the machines. The idea that some day the work of mathematicians might be replaced by computers is an old one. Perhaps it is now happening, but in a different way from that which I would have expected. Computers are replacing humans but not because they are doing everything better. There is no doubt there are some things they can do better but I think there are many things which they cannot. The plenary talks at the DMV conference on topics of this kind were partly critical. There occurred examples of a type I had not encountered before. A computer is presented with a picture of a pig and recognizes it as a pig. Then the picture is changed in a very specific way. The change is quantitatively small and is hardly noticeable to the human eye. The computer identifies the modified picture as an aeroplane. In another similar example the starting picture is easily recognizable as a somewhat irregular seven and is recognized by the computer as such. After modification the computer recognizes it as an eight. This seems to provide a huge potential for mistakes and wonderful opportunities for criminals. I feel that the trend to machine learning and related topics in mathematics is driven by fashion. It reminds me a little of the ‘successes’ of string theory in physics some years ago. Another aspect of the plenary talks at these conferences I did not like was that the speakers seemed to be showing off with how much they had done instead of presenting something simple and fascinating. At the conference in Heidelberg there were three talks by young prizewinners which were shorter than the plenaries. I found that they were on average of better quality and I know that I was not the only one who was of that opinion.

In the end there were not many talks at these conferences I liked much but let me now mention some that I did. Amber Smith gave a talk on the behaviour of the immune system in situations where bacterial infections of the lung arise during influenza. In that talk I really enjoyed how connections were made all the way from simple mathematical models to insights for clinical practise. This is mathematical biology of the kind I love. In a similar vein Stanca Ciupe gave a talk about aspects of COVID-19 beyond those which are common knowledge. In particular she discussed experiments on hamsters which can be used to study the infectiousness of droplets in the air. A talk of Harsh Chhajer gave me a new perspective on the intracellular machinery for virus production used by hepatitis C, which is of relevance to my research. I saw this as something which is special for HCV and what I learned is that it is a feature of many positive strand RNA viruses. I obtained another useful insight on in-host models for virus dynamics from a talk of James Watmough.

Returning to the issue of mathematics and computers another aspect I want to mention is arXiv. For many years I have put copies of all my papers in preprint form on that online archive and I have monitored the parts of it which are relevant for my research interests for papers by other people. When I was working on gravitational physics it was gr-qc and since I have been working on mathematical biology it has been q-bio (quantitative biology) which I saw as the natural place for papers in that area. q-bio stands for ‘quantitative biology’ and I interpreted the word ‘quantitative’ as relating to mathematics. Now the nature of the papers on that archive has changed and it is also dominated by topics strongly related to computers such as machine learning. I no longer feel at home there. (To be fair I should say there are still quite a lot of papers there which are on stochastic topics which are mathematics in the classical sense, just in a part of mathematics which is not my speciality.) In the past I often cross-listed my papers to dynamical systems and maybe I should exchange the roles of these two in future – post to dynamical systems and cross-list to q-bio. If I succeed in moving further towards biology in my research, which I would like to I might consider sending things to bioRxiv instead of arXiv.

In this post I have written a lot which is negative. I feel the danger of falling into the role of a ‘grumpy old man’. Nevertheless I think it is good that I have done so. Talking openly about what you are unsatisfied with is a good starting point for going out and starting in new positive directions.

Advertisement

Lyapunov-Schmidt reduction and stability

I discussed the process of Lyapunov-Schmidt reduction in a previous post. Here I give an extension of that to treat the question of stability. I again follow the book of Golubitsky and Schaeffer. My interest in this question has the following source. Suppose we have a system of ODE depending on a parameter . The equation for steady state solutions is then . Sometimes we can eliminate all but one of the variables to obtain an equation for a real variable whose solutions are in one to one correspondence with those of the original equation for steady states. Clearly this situation is closely related to Lyapunov-Schmidt reduction in the case where the linearization has corank one. Often the reduced equation is much easier to treat that the original one and this can be used to obtain information on the number of steady states of the ODE system. This can be used to study multistationarity in systems of ODE arising as models in molecular biology. In that context we would like more refined information related to multistability. In other words, we would like to know something about the stability of the steady states produced by the reduction process. Stablity is a dynamical property and so it is not a priori clear that it can be investigated by looking at the equation for steady states on its own. Different ODE systems can have the same set of steady states. Note, however, that in the case of hyperbolic steady states the stability of a steady state is determined by the eigenvalues of the linearization of the function at that point. Golubitsky and Schaeffer prove the following remarkable result. (It seems clear that results of this kind were previously known in some form but I did not yet find an earlier source with a clear statement of this result free from many auxiliary complications.) Suppose that we have a bifurcation point where the linearization of has a zero eigenvalue of multiplicity one and all other eigenvalues have negative real part. Let be the corresponding zero of . The result is that if and then for close to the linearization of about the steady state has a unique eigenvalue close to zero and its sign is the same as that of . Thus the stability of steady states arising at the bifurcation point is determined by the function .

I found the proof of this theorem hard to follow. I can understand the individual steps but I feel that I am still missing a global intuition for the strategy used. In this post I describe the proof and present the partial intuition I have for it. Close to the bifurcation point the unique eigenvalue close to zero, call it , is a smooth function of and because it is of multiplicity one. The derivative is also a smooth function. The aim is to show that they have the same sign. This would be enough to prove the desired stability statement. Suppose that the gradient of at is non-zero. Then the zero set of is a submanifold in a neighbourhood of . It turns out that vanishes on that manifold. If we could show that the gradient of is non-zero there then it would follow that the sign of off the manifold is determined by that of . With suitable sign conventions they are equal and this is the desired conclusion. The statement about the vanishing of is relatively easy to prove. Differentiating the basic equations arising in the Lyapunov-Schmidt reduction shows that the derivative of applied to the gradient of a function arising in the reduction process is zero. Thus the derivative of has a zero eigenvalue and it can only be equal to . For by the continuous dependence of eigenvalues no other eigenvalue can come close to zero in a neighbourhood of the bifurcation point.

After this the argument becomes more complicated since in general the gradients of and could be zero. This is got around by introducing a deformation of the original problem depending on an additional parameter and letting tend to zero at the end of the day to recover the original problem. The deformed problem is defined by the function . Lyapunov-Schmidt reduction is applied to to get a function . Let be the eigenvalue of which is analogous to the eigenvalue of . From what was said above it follows that, in a notation which is hopefully clear, implies . We now want to show that the gradients of these two functions are non-zero. Lyapunov-Schmidt theory includes a formula expressing in terms of . This formula allows us to prove that . Next we turn to the gradient of , more specifically to the derivative of with respect to . First it is proved that for all . I omit the proof, which is not hard. Differentiating and evaluating at the point shows that is an eigenvalue of there with eigenvalue . Hence for all . Putting these facts together shows that and the derivative of with respect to at the bifurcation point is equal to one.

We now use the following general fact. If and are two smooth functions, vanishes whenever does and the gradients of both functions are non-zero then extends smoothly to the zero set of and the value of the extension there is given by the ratio of the gradients (which are necessarily proportional to each other). In our example we get with . Setting in the first equation then gives the desired conclusion.

Another paper on hepatitis C: absence of backward bifurcations

In a previous post I wrote about a paper by Alexis Nangue, myself and others on an in-host model for hepatitis C. In that context we were able to prove various things about the solutions of that model but there were many issues we were not able to investigate at that time. Recently Alexis visited Mainz for a month, funded by an IMU-Simons Foundation Africa Fellowship. In fact he had obtained the fellowship a long time ago but his visit was delayed repeatedly due to the pandemic. Now at last he was able to come. This visit gave us the opportunity to investigate the model from the first paper further and we have now written a second paper on the subject. In the first paper we showed that when the parameters satisfy a certain inequality every solution converges to a steady state as . It was left open, whether this is true for all choices of parameters. In the second paper we show that it is not: there are parameters for which periodic solutions exist. This is proved by demonstrating the presence of Hopf bifurcations. These are obtained by a perturbation argument starting from a simpler model. Unfortunately we could not decide analytically whether the periodic solutions are stable or unstable. Simulations indicate that they are stable at least in some cases.

Another question concerns the number of positive steady states. In the first paper we showed under a restriction on the parameters that there are at most three steady states. This has now been extended to all positive parameters. We also show that the number of steady states is even or odd according to the sign of , where is a basic reproductive ratio. It was left open, whether the number of steady states is ever greater than the minimum compatible with this parity condition. If there existed backward bifurcations (see here for the definition) it might be expected that there are cases with and two positive solutions. We proved that in fact this model does not admit backward bifurcations. It is known that a related model for HIV with therapy (Nonlin. Anal. RWA 17, 147) does admit backward bifurcations and it would be interesting to have an intuitive explanation for this difference.

In the first paper we made certain assumptions about the parameters in order to be able to make progress with proving things. In the second paper we drop these extra restrictions. It turns out that many of the statements proved in the first paper remain true. However there are also new phenomena. There is a new type of steady state on the boundary of the positive orthant and it is asymptotically stable. What might it mean biologically? In that case there are no uninfected cells and the state is maintained by infected cells dividing to produce new infected cells. This might represent an approximate description of a biological situation where almost all hepatocytes are infected.

Another conference on biological oscillators at EMBL in Heidelberg

I recently attended a conference at EMBL in Heidelberg and I very much enjoyed experiencing a live conference for the first time in a couple of years. I heard similar sentiments expressed by many of the other participants at the meeting. This conference at EMBL was a sequel to one which I previously wrote about here. The present event was in hybrid form with many of the speakers remote. There were nevertheless more than a hundred people attending on site. The conference started with a presymposium. This was intended to teach some mathematics to biologists. I attended it since I saw it as an opportunity to learn more about what kind of mathematics is really of interest to biologists. Among the main themes discussed were the relationships between positive feedback and multistability and between negative feedback and oscillations. First there was a one-hour talk by Hanspeter Herzel. Then there was a practical part where we were supposed to play with a computer programme. I had downloaded the necessary programmes (R and RStudio) as recommended but this part of the event was a failure for me. I am simply lacking in basic computer competence. It was not explained to us how to begin using the programme and I was not able to supply this missing information on my own. The first part, the lecture, was more interesting for me. The speaker mentioned a paper which he wrote with others about circadian oscillations in the number of lymphocytes in different tissues (D. Druzd et al., Immunity 46, 120). I had previously wondered about the possible roles of oscillations in immunology but I never thought of that direction. I spoke to Herzel about this in a coffee break. This demonstrates a huge advantage of live versus online conferences. I am sure that the information he and I exchanged over coffee would never have been communicated if the conference had been only online. There is a standard picture in immunology in which antigen is being continuously transported to lymph nodes, where it can activate lymphocytes. A key point of the paper is that this does not happen at a constant rate. Instead the process is highly oscillatory. Lymphocytes reach their highest level in the lymph nodes mainly at the beginning of the active phase  (i.e. the beginning of the dark phase in the mice in which these observations were carried out). This means that the effectiveness of a vaccination or another chemical intervention may depend strongly on the time at which it is administered. Herzel told me about an example where this has been seen in practise in cancer immunotherapy. I decided that I wanted to investigate this more closely. Before I could do that I heard the talk of Francis Levi, which was exactly on this topic. Returning to the paper quoted above, according to Herzel the mathematical content was very elementary, using a linear model. I am happy that simple mathematical models and ideas can lead to useful biological insights. What I do not find so good is that the information on the mathematics presented in the paper is so minimal, even in the supplementary material. There is one aspect of this story which is unclear to me. It is important for the functioning of the immune system that a given T cell visits many lymph nodes in a day. Thus the delays to the entrance or exit from lymph nodes which are supposed to implement the rhythm must act in some kind of averaged sense.

I also had a chance to talk to Levi over coffee and get some additional insights about some aspects of his lecture. He has been working on chronotherapy in oncology for many years. This means the idea that the effectiveness of a cancer therapy can be very dependent on the time of day it is administered. He has applied these ideas in practise but the ideas have not gained wide acceptance in the community of oncologists. There is a chance that this may change soon due to the appearance of two papers on this subject in the prestigious journal ‘The Lancet Oncology’ in November 2021. One of the papers (22, 1648) is by Levi, the other (22, 1777) by Qian et al.

Now let me mention a couple of the other contributions I liked best. On Monday there was a (remote) talk by Albert Goldbeter on the coupling between the cell cycle and the circadian clock. Here, as elsewhere in this conference, entrainment was a central theme. There was a discussion of the role of multiple limit cycles in these models. There was also a (remote) talk by Jim Ferrell. His subject was cataloguing certain aspects of an organism called the mouse lemur. The idea was to have a list of cell types and hormones and to know which cell types produce and are affected by which hormones. There is a preprint on this subject on BioRxiv. One feature of these primates which I found striking is the following. They are much fatter in winter than in summer and this is related to a huge difference in thyroid hormones. If I remember correctly it is a factor of ten. For comparison, in humans thyroid hormones also vary with the time of year but only on the scale of a couple of per cent. In a talk by Susan Golden (live) on the Kai system in cyanobacteria I was able to experience one of the pioneers in that field.

Codimension three bifurcations

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 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 , 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 . It says that the system does not belong to a certain degenerate set defined by the vanishing of particular functions of the derivatives of at . In an abstract way this degenerate set may be thought of a subset of the set of -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 induces a mapping from the parameter space into the set of -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 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 of the mapping defined by 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 has a double zero eigenvalue. The genericity condition of the first type involves the -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 -jet of the right hand side of the equation for is of the form . It turns out that we can divide the problem into two cases, in each of which one of the coefficients and is zero and the other non-zero. We first concentrate on the case 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 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 and which are joined by a curve of Hopf bifurcations. The points and are stated to be codimension 2 and this indicates that they are non-degenerate Bogdanov-Takens points. There is one exceptional point 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 . 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 . 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

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

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

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.

This

Talk by David Ho on COVID-19

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)

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.