mathematical biology

Archive for the ‘mathematical biology’ Category

One-dimensional centre manifolds, part 4

April 30, 2024

In this post I consider the case of the transcritical bifurcation, starting with the case of a one-dimensional dynamical system depending on one parameter $\dot x=f(x,\alpha)$ . I consider the situation where for structural reasons $x=0$ is always a steady state. Hence $f(x,\alpha)=xg(x,\alpha)$ for a smooth function $g$ . The equation $g(x,\alpha)=0$ could potentially have a bifurcation at the origin but here I specifically consider the case where it does not. The assumptions are $g(0,0)=0$ and $g_x(0,0)\ne 0$ . Moreover it is assumed that $g_\alpha(0,0)\ne 0$ . By the implicit function theorem the zero set of $g$ is a function of $\alpha$ . By the last assumption made the derivative of this function at zero is non-zero. The graph crosses both axes transversally at the origin. It will be supposed, motivated by intended applications, that $g_x(0,0)<0$ . If $g_\alpha<0$ there exist positive steady states for $\alpha>0$ and if $g_\alpha>0$ there exist positive steady states with $\alpha<0$ . The conditions which have been required of $g$ can be translated to conditions on $f$ . Of course $f(0,0)=0$ and $f_\alpha (0,0)=0$ . Since $g(0,0)=0$ we also have $f_x(0,0)=0$ and the origin is a bifurcation point for $f$ . The other two conditions on $g$ translate to $f_{xx}(0,0)\ne 0$ and $f_{\alpha x}(0,0)\ne 0$ . It follows from the implicit function theorem, applied to $g$ that the bifurcation can be reduced by a coordinate transformation to the standard form $\dot x=x(\alpha\pm x)$ . Thus we see how many steady states there are for each value of $\alpha$ and what their stability properties are.

My primary motivation for discussing this here is to throw light on the concept of backward bifurcation. At a bifurcation point of this kind there is a one-dimensional centre manifold. I want to explain the relation of general centre manifold theory to Theorem 4 of in the paper of van den Driessche and Watmough. Background for this discussion can be found here. It is explicitly assumed in the discussion leading up to the theorem that the centre manifold at the bifurcation point is one-dimensional. All but one of the non-zero eigenvalues of the linearization have negative real parts close to the bifurcation point. The remaining one changes sign as the bifurcation point is passed. The main idea is to reduce the bifurcation to the two-dimensional extended centre manifold at the bifurcation point. In equations (23) and (24) the key diagnostic quantities $a$ and $b$ for the bifurcation are defined. They are expressed in invariant form both in terms of the extrinsic dynamical system and in a form intrinsic to the centre manifold. In the coordinate form used above $a=\frac12 f_{xx}(0,0)$ and $b=f_{\alpha x}(0,0)$ . Theorem 4 of the paper makes statements about the existence and stability of positive steady states which are essentially equivalent to those made above when it is taken into account that in the situation of the theorem the centre manifold is asymptotically stable with respect to perturbations in the transverse directions. It does not say more than that. The fact, often seen in pictures of backward bifurcation that the unstable branch undergoes a fold bifurcation is not included. In particular the fact that for some parameter values there exists more than one positive steady state is not included.

Posted in dynamical systems, mathematical biology | Leave a Comment »

Nobel lecture of Harvey Alter

March 1, 2024

In 2020 the Nobel prize for medicine was awarded to Harvey Alter, Michael Houghton and Charles Rice for their role in the discovery of the hepatitis C virus. I now watched the videos of the corresponding Nobel lectures. For my taste the lecture of Alter was by far the most interesting of the three. I think that he was also the one who played the most fundamental role in this discovery. At the beginning of his lecture he emphasizes the point that the most important discoveries in science often come as a complete surprise and not as a result of planned research programmes. Alter was 85 when he got the prize and so he had to wait a long time for it. The papers documenting his fundamental contributions were published in 1989. A central part of this work was the collection and preservation of blood samples from patients undergoing open heart surgery. Why was this group chosen? One of the most important modes of infection with hepatitis B and C used to be blood transfusions. This continued to be the case until tests were available to screen donors for these diseases. This kind of surgery involves extensive blood transfusions and so the chances of infection were relatively high in these patients. Also these patients suffered from relatively few other diseases which could have been confounding factors. These blood samples were an invaluable resource in the search for the virus. They were the basis of painstaking analysis over many years.

One important feature of hepatitis C is that it becomes chronic in 70 per cent of cases. This looks like a failure of the immune system to handle this disease. What are the reasons for this failure? One concerns quasispecies. The hepatitis C virus has an RNA genome and the copying of RNA is very error-prone. This leads to a huge variety in the genomes of virions in a single patient. This in turn results in rapid mutations of the virus. If an antibody has developed to combat the virus then selective pressure will quickly cause a new form to become dominant which is not vulnerable to that antibody. It seems to me that if this type of effect is to be captured using mathematical model it will require a stochastic model. Deterministic models of the type I have studied in the past are probably not helpful for that. In the lecture it is also mentioned that the number of T cells (CD4+ and CD8+) declines very much in chronically infected hepatitis C patients. No explanation is offerred as to why that is the case. Deterministic mathematical models might be able to contribute some understanding in that case.

The lecture contains the following interesting story. There was a time at which liver cancer was much more common in Japan than in the West. The reason for this was that that cancer was in many cases a late stage effect of hepatitis C. During wars in the early part of the 20th century many Japanese soldiers injected drugs with shared needles and this was what spread the disease. It was observed that there were many cases of jaundice (the most striking symptom of hepatitis) on the battlefield. Decades later many of these men developed serious liver disease, including cancer. Japanese doctors predicted that a similar phenomenon would be seen in the West when the effects of recreational drug use became manifest. They were right.

Posted in diseases, immunology, mathematical biology, public health | Leave a Comment »

Positivity for systems of reaction-diffusion equations

January 29, 2024

Here I consider a system of reaction-diffusion equations of the form $\frac{\partial u_i}{\partial t}=d_i\Delta u_i+f_i(u)$ , $1\le i\le n$ . The functions $u_i(x,t)$ are defined on $\Omega\times [0,\infty)$ where $\Omega$ is a bounded domain in $R^n$ with smooth boundary. Let the $u_i$ be denoted collectively by $u$ . I assume that the diffusion coefficients $d_i$ are non-negative. If some of them are zero then the system is degenerate. In particular there is an ODE special case where all $d_i$ are zero. If this system really describes chemical reactions and the $u_i$ are concentrations then it is natural to assume that $u_i(0)\ge 0$ for all $i$ . It should then follow that in the presence of suitable boundary conditions $u(t)\ge 0$ for all $t\ge 0$ . I assume that $u$ is a classical solution and that it extends to the boundary with enough smoothness that the boundary conditions are defined pointwise. It is necessary to implement the idea that the system is defined by chemical reactions. This can be done by requiring that whenever $u\ge 0$ and $u_i=0$ it follows that $f_i(u)\ge 0$ . (This means that if a chemical species is not present it cannot be consumed in any reaction.) It turns out that this condition is enough to ensure positivity.

I will now explain a proof of positivity. The central ideas are taken from a paper of Maya Mincheva and David Siegel (J. Math. Chem. 42, 1135). Thanks to Maya for helpful comments on this subject. The argument is already interesting for ODE since important conceptual elements can already be seen. I will first discuss that case. Consider a solution of the equation $\dot u=-au+b$ on the interval $[0,\infty)$ where $a$ and $b$ are continuous functions with $b>0$ and suppose that $u(0)>0$ . I claim that $u(t)>0$ for all $t\ge 0$ . Let $t^*=\sup\{t_1:u(t)\ge 0\ {\rm on}\ [0,t_1]\}$ . Since $u(0)>0$ it follows by continuity that $t_1>0$ . Assume that $t_1<\infty$ . By continuity $u(t_1)\ge 0$ . If $u(t_1)$ were greater than zero then by continuity it would also be positive for $t$ slightly larger than $t_1$ , contradicting the definition of $t_1$ . Thus $u(t_1)=0$ . The evolution equation then implies that $\dot u(t_1)>0$ . This implies that $u(t)<0$ for $t$ slightly less than $t_1$ , which also contradicts the definition of $t_1$ . Hence in reality $t_1=\infty$ and this completes the proof of the desired result.

Suppose now that we weaken the assumptions to $b\ge 0$ and $u(0)\ge 0$ . We would like to conclude that $u(t)\ge 0$ for all $t$ . To do this we define a new quantity $v=u+\epsilon e^{\sigma t}$ for positive constants $\epsilon$ and $\sigma$ . Then $\dot v=\dot u+\epsilon\sigma e^{\sigma t}$ . Hence $\dot v=-au+[b+\epsilon\sigma e^{\sigma t}]$ and $v(0)=u(0)+\epsilon>0$ . Now $u=v-\epsilon e^{\sigma t}$ and so $\dot v=-av+[b+\epsilon(\sigma+a) e^{\sigma t}]$ . It follows that if $\sigma$ is large enough $v$ satisfies the conditions satisfied by $u$ in the previous argument and it can be concluded that $v(t)>0$ for all $t$ . Letting $\epsilon$ tend to zero shows that $u(t)\ge 0$ for all $t$ , the desired result.

This is different, and perhaps a bit more complicated than, the proof I know for this type of result. That proof involves considering the derivative of $\log u$ on $[0,t_1)$ . It also involves approximating non-negative data by positive data. A difference is that the proof just given does not use the continuous dependence of solutions of an ODE on initial data and in that sense it is more elementary. In Theorem 3 of the paper of Mincheva and Siegel the former proof is extended to a case involving a system of PDE.

Now I come to that PDE proof. The system of PDE concerned is the one introduced above. Actually the paper requires the $d_i$ be positive but that stronger condition is not necessary. This equation is solved with an initial datum $u_0(x)=u(x,0)$ and a boundary condition $\alpha u+\frac{\partial u}{\partial\nu}=g$ . Here $\alpha$ is a diagonal matrix with non-negative entries and the function $g$ is non-negative. The derivative $\frac{\partial}{\partial\nu}$ is that in the direction of the outward unit normal to the boundary. We assume that $u$ is a classical solution, i.e. all derivatives of $u$ appearing in the equation exist and are continuous. Moreover $u$ has a continuous extension to $t=0$ and a $C^1$ extension to $\bar\Omega\times (0,\infty)$ . We now replace the differential equation by the differential inequality $\frac{\partial u}{\partial t}-D\Delta u-f(u)\ge 0$ . We assume that the initial data are non-negative, $u_0(x)\ge 0$ . The assumption that $g$ is non-negative, together with the boundary condition, gives rise to the inequality $\alpha u+\frac{\partial u}{\partial\nu}\ge 0$ . The aim is to show that all solutions of the resulting system of inequalities are non-negative. We assume the condition for a system of chemical reactions already mentioned.

The proof is a generalization of that already given in the ODE case. The first step is to treat the case where each of the inequalities is replaced by the corresponding strict inequality. In contrast to the proof in the paper we assume that that $u_0$ is strictly positive on $\bar\Omega$ . We define $t_1$ as in the ODE case so that $[0,t_1]$ is the longest interval where the solution is non-negative. We suppose that $t_1$ is finite and obtain a contradiction. Note first that, as in the ODE case, $t_1>0$ by continuity. Now $u(t_1)\ge 0$ . If $u(t_1)$ were strictly positive on $\bar\Omega$ then by continuity $u$ would be strictly positive for $t$ slightly greater than $t_1$ , contradicting the definition of $t_1$ . Hence there is an index $i$ and a point $x_0\in\bar\Omega$ with $u_i(x_0,t_1)=0$ and $u_i(x_0,t)\ge 0$ for all $t<t_1$ . Suppose first that $x_0\in\Omega$ . Then $\frac{\partial u_i}{\partial t}(x_0,t_1)\le 0$ and $\Delta u_i(x_0,t_1)\ge 0$ . This contradicts the strict inequality related to the evolution equation for $u_i$ and so cannot happen. Suppose next that $x_0$ is on the boundary of $\Omega$ . Then $\alpha_i u_i(x_0,t_1)=0$ and $\frac{\partial u_i}{\partial\nu}\le 0$ . This contradicts the strict inequality related to the boundary condition for $u_i$ . Thus in fact $t_1=\infty$ and $u$ is strictly positive for all time.

Now we do a perturbation argument by considering $v=u+\bar\epsilon e^{\sigma t}H(x)$ . Here $\bar\epsilon$ is the vector all of whose components are $\epsilon$ and $H$ is a positive function. It is obvious that $v(0)>0$ . We now choose $H(x)=e^{h(x)}$ which ensures its positivity and require that the outward normal derivative of $h$ on the boundary is equal to one. (Here I will take for granted that a function $h$ of this kind exists. A source for this statement is cited in the paper. For me the positivity statement is already very interesting in the case that $\Omega$ is a ball and there the existence of the function $h$ is obvious.) Then the fact that $u$ satisfies the non-strict boundary inequality implies that $v$ satisfies the strict boundary inequality. It remains to derive an evolution equation for $v$ . A straightfoward calculation gives $\frac{\partial v_i}{\partial t}-d_i\Delta v_i-f_i(v)\ge e^{\sigma t}\epsilon[\sigma H-d_i\Delta H-L\sqrt{n}H]$ where $L$ is a Lipschitz constant for $f$ on the image of the compact region being considered. Choosing $\sigma$ large enough ensures that the right hand side and hence the left hand side of this inequality is strictly positive. We conclude that $v>0$ and, letting $\epsilon\to 0$ , that $u\ge 0$ .

If we want to prove an inequality for solutions of a PDE it is common to proceed as follows. We deform the problem continuously as a function of a small parameter $\epsilon$ so as to get a simpler problem. When that has been solved we let $\epsilon$ tend to zero to get a solution of the original problem. Often it is the equations which are deformed. Then we need a theorem on existence and continuous dependence to get a continuous deformation of the solution. The above proof is different. We perturb the solution in a way whose continuity is obvious and then derive a family of equations of which that is a family of solutions. This is easy and comfortable. The hard thing is to guess a good deformation of the soluion.

Posted in analysis, mathematical biology, partial differential equations | Leave a Comment »

Conference on mathematics and immunology in Blagoevgrad

September 15, 2023

I recently attended a conference on mathematics and immunology in Blagoevgrad in Bulgaria. I must admit that I knew very little about Bulgaria. I knew the name, I knew where it was and I knew that it had a Slavic language. That was about all. Thus it was educational for me to be in the country and also have the chance to talk to Bulgarians. It was in particular interesting to learn that there was a large Bulgarian empire centuries ago. I flew to Sofia and took a train south to the site of the conference. The train was quite old and slow but it passed through interesting countryside, some of which was extremely green and lush. I saw several Bee-eaters from the train. On the evening I arrived there were a lot of swallows flying around. They were flying quite high and fast and I was not able to identify them definitely. They looked and sounded different to the species I was familiar with. While they were eating mosquitos the mosquitos were eating me and that made it difficult to concentrate. A couple of days later I was able to convince myself that one suspicion I had had at the beginning was correct – they were Red-rumped Swallows.

The conference itself was small and I much prefer small and thematically focussed meetings such as this to large conferences. It was very pleasant to be among people who are interested in and knowledgeable about both mathematics and immunology. The first speaker was Becca Asquith who talked about KIRs, the inhibitory receptors found on natural killer cells. What I did not know before was that KIRs are also presented on T cells. It is well known that the susceptibility of individuals to infectious or auto-immune diseases is dependent on which MHC molecules they have. I was aware that, for instance, certain HLA-DR types are associated to an increased risk of multiple sclerosis. In fact HLA molecules may have positive or negative effects on disease risk. It turns out that the types of KIRs a given individual has can also have analogous effects. A central theme of the talk was how these statistically observed effects can be explained mechanistically and how mathematical models can help to distinguish between different mechanisms. A number of talks were related to epidemiology, considering in-host models, population models and the coupling of the two. One concerned a virus I did not know previously, the Usutu virus. Stanca Ciupe described modelling for experiments where house sparrows were exposed to mosquitos carrying this virus. The virus (and its name) originated in Africa but it has been spread to parts of Europe by migratory birds. Doing some more reading on the subject I discovered that it is reponsible for some cases where a lot of blackbirds have been seen to die. Humans can become infected with the virus but it does not seem to be directly harmful to us. The virus is nevertheless interesting as a candidate for one which might mutate and cause a serious disease of humans in the future. Jonathan Forde gave an interesting talk on the question of the best way to invest limited (financial or human) resources in fighting COVID-19 or other infectious diseases. Sometimes there is an optimal balance between vaccination and testing.

There were also talks on subjects other than infectious diseases. Vladimira Suvandjieva talked about some work she has done on modelling the role of NETs in lupus. Here mathematics is applied to a biological subject which I wrote about in this blog a long time ago. In this phenomenon the immune cells involved are neutrophils. Doron Levy talked about using mathematical models to better understand cancer immunotherapy. There were talks on the ways in which the time of administration can influence the effectiveness of chemotherapy against cancer and the strength of side effects. What is notable is a strong dependence on the sex of the patient. Apart from the science it was interesting to talk to the other participants from various parts of the world about political themes, thus being exposed to different facts and opinions. I have not had so many opportunities for this kind of conversation since the pandemic.

Posted in diseases, mathematical biology, public health | Leave a Comment »

Fomites and backward bifurcations

July 25, 2023

I first met the curious word ‘fomite’ a few months ago. It means a carrier of infection which is not a human being or an animal. It evolved in an interesting way, as recounted in the Wikipedia article on this concept. It started out as the Latin word ‘fomes’. The plural of that is ‘fomites’, which was taken over into English. There it was interpreted as an English plural and gave rise to the corresponding singular ‘fomite’. Let me now get away from the word and consider the concept it designates. During the COVID-19 pandemic there was a lot of discussion about the modes of transmission of infection. In particular people wanted to understand the importance of the disinfection of surfaces (and of hands). These inanimate surfaces, where the virus might be present and act as a source of infection, are then (the surfaces of) fomites. As far as I can judge the consensus developed that this mode of transmission was not the key factor in the case of COVID-19 with attention being concentrated on aerosols. However important they may be in that particular case there is no doubt that there are diseases for which transmission by fomites is very important. The theme of hospital infections is a very important one these days and there fomites are likely to play a central role.

I came across the word fomite in the course of some research I did recently with Aytül Gökce and Burcu Gürbüz on some population models for infectious diseases. We just produced a preprint on that. One aspect of infection which was included in our model and which is not very common in the literature is that of infections coming from virus in the environment and here we can think of fomites. In modelling this phenomenon a choice has to be made of a function which describes the force of infection. We can consider an expression for the rate of infection of the from $Sf(C)$ , where $S$ is the number of susceptibles and $C$ is the concentration of virus in the environment which may cause infection. What should we choose for the function $f$ ? This type of issue comes up in other modelling settings, in particular in those related to biology, where a choice of response function has to be made. In a biochemical reaction $f$ could describe the reaction rate as a function of the concentration of the substrate. Typical choices are the Michaelis-Menten function $f(C)=\frac{V_{\max}C}{K+C}$ and the Hill function where $C$ is replaced in this expression by a power $C^n$ . The Michaelis-Menten function has a mechanistic basis, introduced by the people who it is now named after. Hill also presented a mechanistic basis for his function. It incorporates the phenomenon of cooperative binding, with the original example being the binding of oxygen to haemoglobin. Another example of the choice of response functions occurs in predator-prey models in ecology. Here $f$ describes the dependence of the rate of uptake of prey as a function of the density of prey. The functions which are mathematically identical to the Michaelis-Menten and Hill functions are called Holling type II and Holling type III. Holling gave mechanistic interpretations of these. Type II is associated with the fact that the predator needs some time to process prey and that that time is lost for the search for prey. Type III is associated with learning by the predator. In this list type I is missing. Type I refers to a linear function, except that to account for the fact that the rate of uptake of prey is in reality limited the linear function is usually cut off at some level, after which it is constant.

It seems that in epidemiology these issues have not been studied so much. There is an extended discussion of the subject in Chapter 10 of the book ‘Mathematical Epidemiology of Infectious Diseases’ by Diekmann and Heesterbeek. I mentioned this book in a previous post but unfortunately forgot to mention the names of the authors. At the beginning of the chapter it is stated clearly that it will not give a definitive answer to the questions it raises. I appreciate this kind of modesty. In that chapter there is a mechanistic derivation of a response function in a particular epidemiological setting. It is more complicated that the conventional ones and is not a rational function. The function corresponding to Holling type II has been been discussed in relation to epidemiology, apparently for the first time by Dietz. His discussion is purely phenomenological as opposed to mechanistic. In other words he chooses the simplest function $f$ satisfying some basic desired properties. In our paper we consider phenomenological models containing a power $n$ as in Holling type II ( $n=1$ ) and Holling type III ( $n\ge 2$ ). We found that the case $n=2$ gives rise to backward bifurcations while the case $n=1$ does not. Thus here the details of the response function can be of considerable mathematical or even medical importance. (The presence of backward bifurcations has implications for the strategy of therapies.) Another topic of the paper is the different ways in which imperfect vaccination can be described mathematically but I will not go any deeper into this theme here.

Posted in diseases, dynamical systems, mathematical biology, public health | Leave a Comment »

The formalism of van den Driessche and Watmough

May 10, 2023

Here I continue the discussion of the paper of van den Driessche and Watmough from the point I left it in the last post. We have a system of ODE with unknowns $x_i$ , $1\le i\le n$ . The first $m$ of these have the status of infected variables. The right hand side of the evolution equation for $x_i$ is the sum of three terms ${\cal F}_i$ , ${\cal V}_i^+$ and $-{\cal V}_i^-$ where ${\cal F}_i$ , ${\cal V}_i^+$ and ${\cal V}_i^-$ are all non-negative. This is condition (A1) of the paper. On the right hand side of this type of system there are positive and negative terms. If we take ${\cal V}_i^-$ to be the total contribution of the negative terms then it is uniquely defined by the system of ODE. On the other hand splitting the positive terms involves a choice. The notation ${\cal V}_i={\cal V}_i^--{\cal V}_i^+$ is also used. As is usual in reaction networks ${\cal V}_i^-=0$ when $x_i=0$ . This is condition (A2) of the paper. The ${\cal F}_i$ represent new infections and hence must be zero for $i>m$ . This is condition (A3) of the paper. It restricts the ambiguity in the splitting. In condition (A4) the set of disease-free steady states is supposed to be invariant and this leads to the vanishing of some non-negative terms at those points. In general this condition restricts the choice of the infected compartments. In the fundamental model of virus dynamics the infected variables are a subset of $\{y,v\}$ . However choosing it to be a proper subset would violate condition (A4). Thus for this model there is no choice at this point. There does remain a choice in splitting the positive terms, as already discussed in a previous post. It remains to consider the last condition (A5). Its interpretation in words is that if the infection is stopped the disease-free steady state is asymptotically stable. Mathematically, stopping the infection means setting some parameters in the model to zero.

Consider now the linearization of the system at a disease-free steady state. The choice of infectious compartments induces a block structure. Splitting the RHS of the system into ${\cal F}_i$ and ${\cal V}_i$ also splits the linearization into a sum. It is an elementary consequence of the assumptions made that some of the entries of the matrices occurring in these decompositions are zero. It can also be shown that some entries are non-negative. To describe a new infection we can take an initial datum close to the disease-free steady state. Then we may reasonably suppose (leaving aside questions of proof) that the actual dynamics is well-approximated by the linearized dynamics during an initial period. Under the given assumptions the linearization is block lower triangular. The linearized evolution is partially decoupled and the uninfected variables decay exponentially. So the dynamics essentially reduces to that of the infected variables and is generated by the upper left block, called $-V$ in the paper. We get an exponential phase of decay with populations proportional to $e^{-tV}$ . During this phase the number of susceptibles is approximately constant, equal to its value at the disease-free steady state. Hence we get an explicit expression for the rate of infection in this phase. Integrating this with respect to time gives an expression for the ratio of total number of individuals infected in that phase as a function of the original number of individuals in the different compartments. It is given by the matrix $FV^{-1}$ , which is the NGM previously mentioned.

I have not succeeded in making a connection between the formalism of van den Driessche and Watmough and the discussion on p. 19 of the book of May and Nowak. Towards the bottom of that page it is easy to see the connection with the unique positive eigenvalue of the linearization at the boundary steady state if $R_0$ is simply interpreted as the known combination of parameters. After I had written this I noticed a source cited by van den Driessche and Watmough, the book ‘Mathematical Epidemiology of Infectious Diseases’ by Diekmann and Heesterbeek and, in particular, Theorem 6.13 of that book. (To avoid any confusion, I am talking about the original edition of this book, not the later extended version.) I had never seen the book before and I borrowed it from the library. I found a text which pays attention to the relations between mathematical rigour and applications in an exemplary way. I feel that if I read the whole book carefully and did all the exercises (which is the recommendation of the authors) I might be able to get rid of all the problems which are the subject of this and the previous related posts. Since this ideal engagement with the book is something I will not achieve soon, if ever, I restrict myself to some limited comments. On p. 105 the authors write ‘Many modellers … even distrust whether the threshold value one of $R_0$ … corresponds exactly to their stability condition. … The remainder of this section is intended for modellers who recognise themselves in the above description.’ I recognise myself in this sense. The section referred to contains Theorem 6.13.

Posted in books, diseases, dynamical systems, mathematical biology | 1 Comment »

The basic reproduction number, yet again

May 4, 2023

In this post I come back to the basic reproduction number, a concept which seems to me like a thorn in my flesh. Some of what I write here might seem rather like a repetition of things I have written in previous posts but I feel a strong need to consolidate my knowledge on this issue so as to possibly finally extract the thorn. The immediate reason for doing so is that I have been reading the book ‘Infectious diseases of humans’ by Anderson and May. I find this book very good and pleasant to read but my progress was halted when I came to Section 2.2 on p. 17 whose title is ‘The basic reproductive rate of a parasite’. (The question of whether it is better to use the word ‘rate’ or the word ‘number’ in this context is not the one that interests me.) What is written there is ‘The basic reproductive rate, $R_0$ is essentially the average number of successful offspring a parasite is capable of producing’. It is clear that this is not a precise mathematical definition, as shown in particular by the use of the word ‘essentially’. In that section there is further discussion of the concept without a definition. As a mathematician I am made uneasy by a discussion using what is apparently a mathematical concept without giving a definition. A possible reaction to this unease would be to say, ‘Be patient, maybe a definition will be given later in the book’. Unfortunately I have not been able to find a definition elsewhere in the book. I also checked whether it might not be included in some supplementary material or appendices. Later in the book this quantity is calculated in some models. This might provide an opportunity for reverse engineering the definition. This book is not aimed at mathematicians, which is an excuse for the lack of a definition. Next I turned to the book ‘Virus dynamics’ by Nowak and May which is related but a bit more mathematical. Here the place I get stuck is on p. 19. Again there is a definition of $R_0$ in words, ‘the number of newly infected cells that arise from any one infected cell when almost all cells are uninfected’. Just before that a specific model, the ‘basic model of virus dynamics’ has been introduced.

The basic model is clearly defined in the book. It is a system of ordinary differential equations depending on some parameters. Depending on the values of these parameters there is one positive steady state or none. These two cases can be characterized by a certain ratio of the parameters. If this ratio is greater than one there is a positive steady state. If it is less than or equal to one there is none. We can introduce the notation $R_0$ for this quantity and call it the basic reproduction number. This is a clear and permissible mathematical definition and it is an interesting statement about the model that such a quantity exists and can be written down explicitly in terms of the parameters of the system. This definition does, however, have two disadvantages. There are many models similar to this one in the literature. Thus the question comes up: if I have a specific model does there exist a quantity analogous to $R_0$ and if so how can I calculate it? The other is the idea that $R_0$ is not just a mathematical quantity. It potentially contains useful biological insights and that is something I would like to understand.

Is there more to be learned about this issue from the book of Nowak and May? They do have a kind of derivation in words of the expression for $R_0$ . It is obtained as the product of three quantities, each of which comes with a certain interpretation. Thus a strategy would be to try and understand these three quantities individually and then understand why it is relevant to consider their product. At this point I will introduce some insights which I mentioned in a previous post. For many models the existence of a positive steady state is coupled to the stability of a steady state on the boundary of the non-negative orthant (disease-free steady state). I believe that $R_0$ is really a quantity related to the disease-free steady state and that in a sense it is just a coincidence that it is related to the question of the existence of a positive steady state. This is illustrated by the fact that there are models exhibiting a so-called backward bifurcation where there do exist positive steady states in some cases with $R_0<1$ . Beyond this, I believe that $R_0$ depends only on the linearization of the system about the disease-free steady state.

As I have discussed elsewhere, the best source I know for understanding $R_0$ mathematically is a paper of van den Driessche and Watmough. There it is explained how, subject to making some choices, the linearization of the system at the disease-free steady state gives rise to something called the ‘next generation matrix’, let me abbreviate it by NGM. In that paper $R_0$ is defined to be the spectral radius of the NGM. Thus we obtain a rigorous definition of $R_0$ provided we have a rigorous definition of the NGM. In addition it is proved in that paper that the linearization at the disease-free steady state has an eigenvalue with positive real part precisely when $R_0>1$ , a fact with obvious consequences for stability of that steady state. It is clearly stated in the paper that the NGM does not only depend on the system of ODE defining the model but also on the choices already mentioned. The first choice is that of the $n$ variables $m$ are supposed to correspond to infected individuals. A disease-free steady state is then by definition one at which these $m$ variables vanish. This means that if we have a specific model and a specific choice of boundary steady state where some unknowns vanish and we want to apply this approach the infected variables must be a subset of the unknowns which vanish there. It is not ruled out that they might be a proper subset. This discussion has made it clear to me what strategy I should pursue to make progress in this area. I should look in great detail at the part of the paper of van den Driessche and Watmough where the method is set up. Since I think this post is already long enough I will do that on another occasion.

Posted in books, diseases, dynamical systems, mathematical biology | Leave a Comment »

Is mathematics being driven out by computers?

September 28, 2022

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.

Posted in diseases, dynamical systems, immunology, life, mathematical biology | 1 Comment »

Lyapunov-Schmidt reduction and stability

August 30, 2022

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 $\dot y+F(y,\alpha)=0$ depending on a parameter $\alpha$ . The equation for steady state solutions is then $F(y,\alpha)=0$ . Sometimes we can eliminate all but one of the variables to obtain an equation $g(x,\alpha)=0$ for a real variable $x$ 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 $F$ 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 $(y_0,\lambda_0)$ where the linearization of $F$ has a zero eigenvalue of multiplicity one and all other eigenvalues have negative real part. Let $x_0$ be the corresponding zero of $g$ . The result is that if $g(x)=0$ and $g'(x)\ne 0$ then for $x$ close to $x_0$ the linearization of $F$ about the steady state has a unique eigenvalue close to zero and its sign is the same as that of $g'(x)$ . Thus the stability of steady states arising at the bifurcation point is determined by the function $g$ .

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 $\mu$ , is a smooth function of $x$ and $\alpha$ because it is of multiplicity one. The derivative $g'$ 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 $g'$ at $x_0$ is non-zero. Then the zero set of $g$ is a submanifold in a neighbourhood of $x_0$ . It turns out that $\mu$ vanishes on that manifold. If we could show that the gradient of $\mu$ is non-zero there then it would follow that the sign of $\mu$ off the manifold is determined by that of $g'$ . With suitable sign conventions they are equal and this is the desired conclusion. The statement about the vanishing of $\mu$ is relatively easy to prove. Differentiating the basic equations arising in the Lyapunov-Schmidt reduction shows that the derivative of $F$ applied to the gradient of a function $\Omega$ arising in the reduction process is zero. Thus the derivative of $F$ has a zero eigenvalue and it can only be equal to $\mu$ . 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 $g'$ and $\mu$ could be zero. This is got around by introducing a deformation of the original problem depending on an additional parameter $\beta$ and letting $\beta$ tend to zero at the end of the day to recover the original problem. The deformed problem is defined by the function $\tilde F(y,\alpha,\beta)=F(y,\alpha)+\beta y$ . Lyapunov-Schmidt reduction is applied to $\tilde F$ to get a function $\tilde g$ . Let $\tilde\mu$ be the eigenvalue of $D\tilde F$ which is analogous to the eigenvalue $\mu$ of $DF$ . From what was said above it follows that, in a notation which is hopefully clear, $\tilde g_x=0$ implies $\tilde\mu=0$ . We now want to show that the gradients of these two functions are non-zero. Lyapunov-Schmidt theory includes a formula expressing $\tilde g_{x\beta}$ in terms of $F$ . This formula allows us to prove that $\tilde g_{x\beta}(0,0,0)=\langle v_0^*,v_0\rangle>0$ . Next we turn to the gradient of $\tilde \mu$ , more specifically to the derivative of $\tilde \mu$ with respect to $\beta$ . First it is proved that $\tilde\Omega (0,0,\beta)=0$ for all $\beta$ . I omit the proof, which is not hard. Differentiating $\tilde F$ and evaluating at the point $(0,0,\beta)$ shows that $v_0$ is an eigenvalue of $D\tilde F$ there with eigenvalue $\beta$ . Hence $\tilde\mu (0,0,\beta)=0$ for all $\beta$ . Putting these facts together shows that $\tilde\mu (\tilde \Omega,0,\beta)=\beta$ and the derivative of $\tilde\mu (\tilde \Omega,0,\beta)$ with respect to $\beta$ at the bifurcation point is equal to one.

We now use the following general fact. If $f_1$ and $f_2$ are two smooth functions, $f_2$ vanishes whenever $f_1$ does and the gradients of both functions are non-zero then $f_2/f_1$ extends smoothly to the zero set of $f_1$ 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 $\tilde\mu(\tilde\Omega,\alpha,\beta)=\tilde a(x,\alpha,\beta)\tilde g_x(x,\alpha,\beta)$ with $\tilde a(0,0,0)=[\frac{\partial}{\partial\beta}(\tilde\Omega,0,0)]/\tilde g_{x\beta}(0,0,0)]>0$ . Setting $\beta=0$ in the first equation then gives the desired conclusion.

Posted in analysis, dynamical systems, mathematical biology | Leave a Comment »

Another paper on hepatitis C: absence of backward bifurcations

June 13, 2022

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 $t\to\infty$ . 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 $R_0-1$ , where $R_0$ 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 $R_0<1$ 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.